All experiments were performed in accordance with the Norwegian Animal Welfare Act and the European Convention for the Protection of Vertebrate Animals used for Experimental and Other Scientific Purposes, Permit numbers 18011 and 29893.
Subjects
Male C57/Bl6 mice were housed in social groups of 2–6 individuals per cage (calcium imaging experiments) or individually (electrophysiology experiments, after implantation). The mice had access to nesting material and a planar running wheel and were kept on a 12 h light/12 h darkness schedule in a temperature and humiditycontrolled vivarium. Food and water were provided ad libitum. Twophoton calcium imaging data were collected from a cohort of 12 mice (5 implanted in MEC, 4 in PaS, and 3 in VIS). Electrophysiological data from the MEC were collected from 2 mice.
Surgeries
For all surgeries, anaesthesia was induced by placing the subjects in a plexiglass chamber filled with isoflurane vapour (5% isoflurane in medical air, flow of 1 l min^{−1}). Surgery was performed on a heated surgery table (38 °C). Air flow was kept at 1 l min^{−1} with 1–3% isoflurane as determined from physiological monitoring of breathing and heartbeat. The mice were allowed to recover from surgery in a heated chamber (33 °C) until they regained complete mobility and alertness. Postoperative analgesia was given in the form of subcutaneous injections of Metacam (5 mg kg^{−1}) 24 and 48 h after the first Metacam injection as long as was deemed necessary. Additionally, the mice were given subcutaneous injections or oral administration of Temgesic (0.05–0.1 mg kg^{−1}) with 6 to 8h (injections) or 12h (oral) intervals for the first 36 h after the first Temgesic injection.
Surgeries for calcium imaging
Surgeries were performed according to a twostep protocol. During the first procedure, newborn pups or adult mice were injected in MEC or PaS, or adult mice were injected in VIS with a virus carrying a construct for the expression of the calcium indicator GCaMP6m. The virus (for all injections: AAV1SynGcaMP6m; titre 3.43 × 10^{13} genome copies per ml, AV1PV2823, UPenn Vector Core, University of Pennsylvania, USA) was diluted 1:1 in sterile DPBS (1× Dulbecco’s Phosphate Buffered Saline, Gibco, ThermoFisher). During the second procedure, two weeks later, a microprism was implanted to gain optical access to infected neurons located in MEC and PaS, or a glass window was inserted to obtain similar access in VIS.
Virus injection and microprism implantation in MEC and PaS
In the first surgical procedure, newborn pups received injections of AAV1SynGCaMP6m one day after birth^{51}. An analgesic was provided immediately before the surgery (Rymadil, Pfizer, 5 mg kg^{−1}). Preheated ultrasound gel (39 °C, Aquasonic 100, Parker) was generously applied on the pup’s head in order to create a large medium for the transmission of ultrasound waves. Realtime ultrasound imaging (Vevo 1100 System, Fujifilm Visualsonics) allowed for targeted delivery of the viral mixture to specific areas of the brain. During ultrasound imaging, the pup was immobilized through a custommade mouth adapter. The ultrasound probe (MS550S) was lowered to be in close contact with the gel and thus the pup’s head to allow visualization of the targeted structures. The probe was kept in place for the whole duration of the procedure via the VEVO injection mount (VEVO Imaging Station. Imaging in BMode, frequency: 40 MHz; power: 100%; gain: 29 dB; dynamic range: 60 dB). Target regions were identified by structural landmarks: the MEC or PaS were identified in the antero–posterior and medio–lateral axis by the appearance of the aqueduct of Sylvius and the lateral sinus. The target area for injection was comparable to a coronal section at ∼−4.7 mm from bregma in the adult mouse. The solution containing the virus (250 ± 50 nl per injection) was injected in the target regions via beveled glass micropipettes (Origio, custom made; outer tip opening: 200 μm; inner tip opening: 50 μm) using a pressurepulse system (Visualsonics, 5 pulses, 50 nl per pulse). The pipette tip was pushed through the brain without any incision on the skin, or a craniotomy, and, to reduce the duration of the procedure, retracted immediately after depositing the virus in the target area. The anatomical specificity of the infection was verified by imaging serial sections of the infected hemispheres after experiment completion (see ‘Histology of calcium imaging mice and reconstruction of fieldofview location’).
Two weeks after the viral injection, we performed a second procedure, in which a microprism was implanted in the left hemisphere to gain optical access to the superficial layers of MEC and PaS^{52}. The implanted microprism was a rightangle prism with 2 mm side length and reflective enhanced aluminium coating on the hypotenuse (Tower Optical). The prism was glued to a 4mmdiameter (CS4R, thickness no. 1) round coverslip with UVcurable adhesive (Norland). On the day of surgery, mice were anaesthetized with isoflurane (IsoFlo, Zoetis, 5% isoflurane vapourised in medical air delivered at 0.8–1 l min^{−1}) after which two analgesics were provided through intraperitoneal injection (Metacam, Boehringer Ingelheim, 5 mg kg^{−1} or Rimadyl, Pfizer, 5 mg kg^{−1}, and Temgesic, Indivior, 0.05–0.1 mg kg^{−1}) and one local analgesic was applied underneath the skin covering the skull (Marcain, Aspen, 1–3 mg kg^{−1}). Their scalp was removed with surgical scissors and the surface of the bone was dried before being generously covered with optibond (Kerr). To increase the thickness and stability of the skull and overall preparation, a thin layer of dental cement (Charisma, Kulzer) was applied on the exposed skull, except in the location above the implant, where a 4mmwide circular craniotomy was made. The craniotomy was positioned over the dorsal surface of the cortex and cerebellum, with the centre positioned ∼ 4 mm lateral from the centre of the medial sinus, and above the transverse sinus just above the MEC and PaS. After the dura was removed above the cerebellum, the lower edge of the prism was slowly pushed in the empty space between the forebrain and the cerebellum, just posterior to the transverse sinus. The edges of the coverslip were secured to the surrounding skull with UVcurable dental cement (Venus Diamond Flow, Kulzer). A customdesigned steel headbar was attached to the dorsal surface of the skull, centred upon and positioned parallel to the top face of the microprism. All exposed areas of the skull, including the headbar, were finally covered with dental cement (Paladur, Kulzer) and made opaque by adding carbon powder (Sigma Aldrich) until the dental cement powder became dark grey.
Virus injection and glass window implantation in VIS
In a different cohort of mice than those used for MEC/PaS imaging, we induced the expression of GCaMP6m in neurons of the adult VIS for subsequent imaging. We targeted the injection of the same AAV1SynGCaMP6m viral solution used in the developing MEC and PaS to the primary visual cortex. On the day of surgery, 3 to 5monthold mice were anaesthetized with isoflurane (IsoFlo, Zoetis, 5 % isoflurane vapourized in medical air delivered at 0.8–1 l min^{−1}) after which two analgesics were provided through intraperitoneal injection (Metacam, Boehringer Ingelheim, 5 mg kg^{−1} or Rimadyl, Pfizer, 5 mg kg^{−1}, and Temgesic, Indivior, 0.05–0.1 mg kg^{−1}) and one local anaesthetic was applied underneath the skin covering the skull (Marcain, Aspen, 1–3 mg kg^{−1}). The virus was injected at three locations in VIS, all of which were within the following anatomical ranges in the right hemisphere: 2.3–2.5 mm lateral from the midline, 0.9–1.3 mm anterior from lambda^{53}. At each injection site, 50 nl of the virus was injected 0.5 mm below the dura and the pipette was left in place for 3–4 min to enable the virus to diffuse. The pipette was then brought to 0.3 mm below the dura and another 50 nl was injected. The pipette was then left in place for 5–10 min before retracting it completely. The speed of the injections was 5 nl s^{−1}.
Two weeks after the viral injection, a surgery to chronically implant a glass window over VIS was performed. The mice were handled as previously described for the prism surgery in MEC/PaS, including anaesthesia, delivery of analgesics, and scalp removal. Optibond was applied to the exposed skull except in the location of the craniotomy. A 4mmwide craniotomy was made, centred on the virus injection coordinates, and a 4mm glass window was placed underneath the skull edges of the craniotomy. The glass was slightly larger than the craniotomy, so after it was manoeuvred in place, the upward pressure exerted by the brain secured it in place against the skull, thereby minimizing the presence of empty gaps that might favour tissue and bone regrowth. The edges of the window were secured with UVcurable dental cement and superglue before the positioning of the headbar as described for the MEC–PaS implantation. All exposed areas of the skull, including the headbar, were finally covered with dental cement (Paladur, Kulzer) that was made opaque by adding carbon powder (Sigma Aldrich) until the dental cement powder became dark grey.
Neuropixels probe implants
Two adult mice (4 to 5 months old) were implanted with fourshank Neuropixels 2.0 silicon probes^{54} targeting the superficial layers of MEC in the left hemisphere. Prior to the surgery, the mice were given general analgesics (Metacam, Boehringer Ingelheim, 5 mg kg^{−1} and Temgesic, Indivior, 0.05–0.1 mg kg^{−1}) subcutaneously and one local anaesthetic was applied underneath the skin covering the skull (Marcain, Aspen, 1–3 mg kg^{−1}). After incision, a hole was drilled over the cerebellum for an anchor screw connected to a ground wire. Craniotomies were then drilled. Probes targeting the MEC were lowered from the surface to depths between 2.5 mm and 2.7 mm relative to the dura mater. They were implanted with the most medial shank placed on the brain surface 3.2 mm lateral to the midline and 0.4 mm anterior to the transverse sinus edge. The four shanks were oriented with the electrode sites on the posterior side. In one of the two mice (no. 104638), the probe was first rotated 7° in the horizontal plane (angle with reference to the coronal plane), with the most lateral shank in the most posterior position such that the shanks were parallel to the transverse sinus. The four shanks were then lowered vertically from this position.
The Neuropixels probe of the second mouse (no. 102335) was not rotated in the horizontal plane—that is, all shanks had the same anterior–posterior coordinates. The electrode shanks of this mouse were lowered from the surface with a 2° angle relative to the coronal plane, such that the shank tips were the most posterior. The shanks remained within the same sagittal plane as they were lowered. This second mouse was also implanted with a probe targeting the CA1 region in the right hemisphere, 1.225–1.975 mm relative to the midline, at a depth of 3 mm relative to dura mater, with all shanks 2.1 mm posterior to bregma. The hippocampal data were not used in the present study. The probes were secured to the skull using an adhesive (OptiBond, Kerr), UVcurable dental cement (Venus Diamond Flow, Kulzer), and dental cement (Meliodent, Kulzer). A headbar was attached as described above for the calcium imaging studies.
Selfpaced running behaviour under sensoryminimized conditions
Training of mice began 2 days after the prism implantation in MEC and PaS, 12 days after the implantation of a cranial window in VIS, and 5–7 days after Neuropixels probe implantation. All mice used for calcium imaging recordings and one Neuropixelsimplanted mouse (no. 104638) were headrestrained by a headbar with their limbs resting on a freely rotating styrofoam wheel with a metal shaft fixed through the centre. The radius of the wheel was ∼85 mm and the width 70 mm. Low friction ball bearings (HK 0608, Kulelager) were affixed to the ends of the metal shaft and held in place on the optical table using a custom mount. This arrangement allowed the mice to selfregulate their movement. The position of the mouse on the rotating wheel was measured using a rotary encoder (E6B2CWZ3E, YUMO) attached to its centre axis. Step values of the encoder (4,096 per full revolution, ∼130 µm resolution) were digitized by a microcontroller (Teensy 3.5, PJRC) and recorded using custom Python scripts at 40–50 Hz. Wheel tracking was triggered at the start of imaging and synchronized to the ongoing image acquisition through a digital input from the 2photon microscope. In a subset of mice recorded with calcium imaging (3 out of 12; 2 implanted in MEC, 1 implanted in PaS), the precise synchronization was not available to us and these data were hence not used for comparison of movement and imaging data. A Tslot photo interrupter (EESX672, Omron) served as a lap (full revolution) counter. Design and code of the wheel are publicly available under https://github.com/kavlintnu/wheel_tracker.
The other Neuropixels probeimplanted mouse (no. 102335) was headrestrained by a headbar while resting on a circular disc coated with rubber spray. The radius of this wheel was ∼85 mm. The mouse was allowed selfpaced movement on the wheel. Threedimensional motion capture (OptiTrack Flex 6 cameras and Motive recording software) was used to track the rotation of the wheel by tracking retroreflective markers placed on the wheel edge. Digital pulses were generated using an Arduino microcontroller which were used to align the Neuropixels acquisition system and the OptiTrack system via direct TTL input and infrared LEDs.
In all mice, the selfpaced task was performed under conditions of minimal sensory stimulation, in darkness, and with no rewards to signal elapsed time or distance run^{16,17}. Prior to the imaging sessions, the calcium imaging mice were accustomed to the setup through daily exposures over the course of between 5 and 15 sessions, one session per day. Neuropixelsimplanted mice were habituated to the setup by gradually increasing the time spent on the wheel over four days. In each session, after the mice were positioned on the wheel, they were gently headrestrained and free to run or rest^{55,56} for 30, 45 or 60 min.
Recording sessions of Neuropixelsimplanted mice also consisted of trials where the mice were freely foraging in a 80 cm × 80 cm open field arena for 30 min. These open field trials preceded the selfpaced wheel trials and were not used in the present study.
Twophoton imaging in headfixed mice
A custombuilt 2photon benchtop microscope (Femtonics, Hungary) was used for 2photon imaging of the target areas (that is, superficial layers of MEC, PaS and VIS). A Ti:Sapphire laser (MaiTai Deepsee eHP DS, SpectraPhysics) tuned to a wavelength of 920 nm was used as the excitation source. Average laser power at the sample (after the objective) was 50–120 mW. Emitted GCaMP6m fluorescence was routed to a GaAsP detector through a 600 nm dichroic beamsplitter plate and 490–550 nm bandpass filter. Light was transmitted through a 16×/0.8 NA waterimmersion objective (MRP07220, Nikon) carefully lowered in close contact to the coverslip glued to the microprism (for MEC–PaS imaging) or above the coverslip in contact with the brain surface (for VIS imaging). For the microprismimplanted mice, the objective lens was aligned to the ventro–lateral corner of the prism, to consistently identify the position of MEC and PaS across mice. Ultrasound gel (Aquasonic 100, Parker) or water was used to fill the gap between the objective lens and the glass coverslips. The software MESc (v 3.3 and 3.5, Femtonics, Hungary) was used for microscope control and data acquisition. Imaging time series of either ∼30 min or ∼60 min were acquired at 512 × 512 pixels (sampling frequency: 30.95 Hz, frame duration: ∼32 ms; pixel size: either 1.78 × 1.78 µm^{2} or 1.18 × 1.18 µm^{2}). Time series acquisition was initiated arbitrarily after the mouse was headrestrained on the setup.
Neuropixels recordings in headfixed mice
Signals were recorded using a Neuropixels acquisition system as described previously^{25,57}. In short, the electrophysiological signal was amplified with a gain of 80, lowpassfiltered at 0.5 Hz, highpassfiltered at 10 kHz, and digitized at 30 kHz on the probe circuit board. The digitized signal was then multiplexed by the ‘headstage’ circuit board and transmitted along a 5 m tether cable using twisted pair wiring to a Neuropixels PXIe acquisition module. The data was visualized and recorded using SpikeGLX version 20201103 software (https://billkarsh.github.io/SpikeGLX).
Histology
Histology of calcium imaging mice and reconstruction of fieldofview location
On the last day of imaging, after the imaging session, the mice were anaesthetized with isoflurane (IsoFlo, Zoetis) and then received an overdose of sodium pentobarbital before transcardial perfusion with freshly prepared PFA (4% in PBS). After perfusion, the brain was extracted from the skull and kept in 4% PFA overnight for postfixation. The PFA was then exchanged with 30% sucrose to cryoprotect the tissue.
To verify the anatomical location of the imaged FOVs in the microprismimplanted mice, we used small, custommade pins, derived from a thin piano wire coated with a solution of 1,1′dioctadecyl3,3,3′,3′tetramethylindocarbocyanine perchlorate (DiI; DiIC18(3)) (ThermoFischer), to mark the location of the imaged tissue in relation to the prism footprint. A DiIcoated pin was inserted into the brain tissue at the location left empty by the prism footprint, and specifically targeted to the ventro–lateral corner of the footprint (see ‘Surgeries’). The pin was left in place to favour transfer of DiI from the metal pin to the brain tissue, and to leave a fluorescent mark on the location of the imaged FOV. After 30 to 60 s, the pin was removed and the brain was sliced on a cryostat in 30–50 µm thick sagittal sections. All slices were collected sequentially in a 24well plate filled with PBS, before being mounted in their appropriate anatomical order on a glass slide in custommade mounting medium. For confocal imaging, a Zeiss LSM 880 microscope (Carl Zeiss) was used to scan through the whole series of slices and locate the position of the DiI fluorescent mark. Images were then acquired using an EC PlanNeofluar 20×/0.8 NA air immersion, 40×/1.3 oil immersion, or 63×/1.4 oil immersion objective (Zeiss, laser power: 2–15%; optical slice: 1.28–1.35 airy units, step size: 2 µm). Before acquisition, gain and digital offset were established to optimize the dynamic range of acquisition to the dynamic range of the GCaMP6m and DiI signals. Settings were kept constant during acquisition across brains. Based on the location of the red fluorescent mark, we could infer where, on the medio–lateral and dorso–ventral extent of the brain, the ventro–lateral corner of the microprism (and hence the 2photon FOV aligned to it) was located.
We used the Paxinos mouse brain atlas^{53} to produce a reference flat map representing the medio–lateral and dorso–ventral extent of the MEC and PaS. Flat maps helped delineate the extent of the FOV that fell within the anatomical boundaries of either the MEC and adjacent PaS, and allowed for a standardized comparison across mice. For each imaged mouse, we mapped the dorso–ventral and medio–lateral location of the DiI mark on the refence flat map (Extended Data Fig. 1c). Mice were assigned to ‘MEC imaging’ or ‘PaS imaging’ groups depending on the location of the FOV: a mouse would be further analysed as being part of the MEC imaging group if more than 50% of the area of the FOV occupied by GCaMP6mexpressing cells could be located in the MEC.
To verify the anatomical location of the FOVs in VIS in the glass window implanted mice, we sliced the brain until we reached the anatomical coordinates at which the virus was infused (see ‘Surgeries’). Coronally cut slices of 50 µm thickness were collected sequentially in a 24well plate, and immediately mounted in their appropriate anatomical order on a glass slide in custommade mounting medium. For confocal imaging, a Zeiss LSM 880 microscope (Carl Zeiss) was used according to the same specification as described above for MEC/PaS.
Histology and reconstruction of Neuropixels probe placement
After the end of experiments, the mice were anaesthetized and received an overdose of isoflurane (IsoFlo, Zoetis) before transcardial perfusion with saline followed by 4% formaldehyde. The brain was either extracted after perfusion or kept overnight in 4% formaldehyde for postfixation before extraction. The brains were then stored in 4% formaldehyde. Frozen 30 µm thick sagittal sections were cut on a cryostat, mounted on glass, and stained with Cresyl violet (Nissl). To estimate the shank locations, we used an Axio Scan.Z1 (Carl Zeiss) slide scanner microscope for brightfield detection at 20x magnification. We used Paxinos mouse brain atlas^{53} and the Allen Mouse Brain Common Coordinate Framework^{58} version 3 through the siibraexplorer (Forschungszentrum Juelich, https://atlases.ebrains.eu/viewer/) to estimate anatomical location of recording sites. A map of the probe shank was aligned to the histology assuming that the cutting plane was nearparallel to the sagittal plane. When possible, the anatomical locations were calculated using the tip of the probe shanks and the intersection of the shank with the brain surface as reference frames. When this was not possible, the profile of a nearby brain region (for example, the hippocampus) was used to estimate the MEC implant site. We observed thetarhythmicity of neural activity on all recorded shanks, as expected for recording locations in the MEC.
Analysis of imaging time series
Imaging time series data were analysed using the Suite2p^{59} Python library (https://github.com/MouseLand/suite2p). We used its builtin routines for motion correction, region of interests (ROI) extraction, neuropil signal estimation, and spike deconvolution. Nonrigid motion correction was chosen to align each frame iteratively to a template. Quality was assessed by visual inspection of the corrected stacks and builtin motion correction metrics. The Suite2p GUI was used to manually subselect putative neurons based on anatomical and signal characteristics and to discard obvious artefacts that accumulated during the analysis—for example, ROIs with footprints spanning large areas of the FOV, ROIs that did not have clearly delineated circumferences in the generated maximum intensity projection, or ROIs that were extracted automatically but showed no visible calcium transients.
Raw fluorescence calcium traces of each ROI were neuropilcorrected to create a fluorescence calcium signal F_{corr} by subtracting 0.7 times the neuropil signal from the raw fluorescence traces. We used the Suite2p integrated version of nonnegative deconvolution^{60} with tau = 1 s to deconvolve F_{corr}, yielding the basis for the binarized sequences that we refer to as the calcium activity (see ‘Binary deconvolved calcium activity and matrix of calcium activity’). Due to the absence of ground truth data for our combination of indicator, region, and imaging conditions, we used a decay tau that was at the lower end of biologically plausible values (tau = 1 s), which allowed even short and low amplitude spiking responses to be picked up by the analysis and therefore did not bias our analysis towards largeamplitude calcium transients (presumed bursting responses). To estimate the signaltonoise ratio (SNR) of each cell individually, we further thresholded the calcium activity (without binarization) at 1 s.d. over the mean, yielding filtered calcium activity, and classified the remaining activity as noise. We additionally ensured that noise was temporally well segregated from filtered calcium activity by requiring data points classified as noise to be separated by at least one second before and ten seconds after filtered calcium activity. The SNR of the cell was then estimated as the ratio of the mean amplitude of F_{corr} during episodes of filtered calcium activity over the s.d. of F_{corr} during episodes of noise. If no data points remained after the filtering of calcium activity, the cell was assigned a SNR of zero.
Binary deconvolved calcium activity and matrix of calcium activity
In order to denoise the recorded fluorescence calcium signals and have good temporal resolution, all analyses in the study were performed using the deconvolved calcium activity of the recorded cells. For each cell whose SNR was larger than 4, the deconvolved calcium activity (see ‘Analysis of imaging time series’) was downsampled by a factor of 4 by calculating the mean over time windows of ∼129 ms (original sampling frequency = 30.95 Hz, sampling frequency used in the analyses = 7.73 Hz). Because the ultraslow oscillations and periodic sequences unfolded at the time scales of seconds to minutes, this downsampling step gave a good temporal resolution for all quantifications while allowing us to work with smaller arrays (ultraslow oscillations and the oscillatory sequences were also detectable when using the original sampling frequency), which in some of the analyses reduced the computing time. Next, the downsampled deconvolved calcium activity was averaged over time and its s.d. was calculated. A threshold equal to this average plus 1.5 times the s.d. was used to convert the deconvolved calcium activity into a binary deconvolved calcium activity, such that all values above the threshold were set to 1 (calcium events), and all values below or equal to that threshold were set to 0. Unless stated otherwise, for all analyses throughout the study we used the deconvolved and binary calcium activity, to which for simplicity we refer to as ‘deconvolved calcium activity’ or simply ‘calcium activity’. The calcium activity of all cells in a session with SNR > 4 was stacked to construct a binary matrix of calcium activity which had as many rows as neurons, and as many columns as time bins sampled at 7.73 Hz. The population vectors are the columns of the matrix of calcium activity.
Note that the recorded calcium signals likely reflect a combination of groups of single spikes and higherfrequency bursts, although it was not possible to distinguish between the two types of firing. The sensitivity of the calcium indicator was likely not high enough to detect subthreshold potentials.
Spike Sorting and singleunit selection
Spike sorting of Neuropixels data was performed using a version of KiloSort 2.5 (ref. ^{54}) with some customizations to improve performance on recordings from the MEC region as described previously^{25}. All trials in a session were clustered together. Single units were discarded from analysis based on a < 20% estimated contamination rate with spikes from other neurons. These units were automatically labelled by the KiloSort 2.5 algorithm as ‘good’ units. In the example session from mouse no. 104638 only good units were considered. In the example session of mouse no. 102335, because the number of good units was lower (<250), we also used multiunit activity (MUA).
Autocorrelations and spectral analysis of singlecell calcium activity
To determine if the calcium activity of single cells displays ultraslow oscillations, for each neuron the PSD was calculated on the autocorrelation of its calcium activity. The PSD was computed using Welch’s method (pwelch, builtin Matlab function), with Hamming windows of 17.6 min (8,192 bins of 129 ms in each window) and 50% of overlap between consecutive windows. Note that when calculating the PSD a large window was needed to identify oscillation frequencies ≪0.1 Hz.
To visualize whether specific oscillatory patterns at fixed frequencies were present in the neural population, all autocorrelations from one session were sorted and stacked into a matrix, where rows are cells and columns are time lags. The sorting of autocorrelations was performed according to the maximum power of each PSD in a descending manner. The frequency at which the PSD peaked was used as an estimate of the oscillatory frequency of the cell’s calcium activity.
In order to determine significance for the peak of the PSD, we considered two extreme and opposite shuffling procedures: On the one hand, given that circularly shuffling the data preserves all inter calcium events (Extended Data Fig. 3c,d), taking this approach would preserve the shape and the position of the peak in the PSD calculated on experimental data. On the other hand, destroying the inter calcium event intervals by assigning a random position to each calcium event in the time series would lead to a flat PSD (Extended Data Fig. 3c,d). In the latter approach, all cells would be classified as oscillatory. To bridge these two approaches we developed a new shuffling procedure. For each cell we divided its calcium activity vector into n epochs of length W, with \(n=T/(W\,\bullet \,{\rm{SF}})\), where T is the total number of time bins sampled at a frequency SF = 7.73 Hz (that is, bin size = 129 ms). We next shuffled those epochs (and preserved the ordering of the time bins within each epoch). This method preserved the inter calcium event interval, but at the same time disrupted the periodicity. In the limit where W = 129 ms, this method coincides with shuffling all calcium events without preserving the inter calcium event intervals; in the limit where W = T/SF, this method is equivalent to circularly shuffling the data. For each of the 200 shuffled realizations we calculated the PSD and the fraction of cells for which the peak of the PSD in experimental data was above the 95th percentile of a shuffled distribution built with the values of the PSDs calculated on shuffled data (and at the frequency at which the PSD computed on experimental data peaked). Here we present the results for 5 different epoch lengths:
W = 1 s: 6226 oscillatory cells out of 6231 (99%)
W = 10 s: 6153 oscillatory cells out of 6231 (99%)
W = 20 s: 5695 oscillatory cells out of 6231 (91%)
W = 50 s: 4642 oscillatory cells out of 6231 (74%)
W = 100 s: 3521 oscillatory cells out of 6231 (56%)
When W is below the typical duration of the sequences (W < 50 s), the great majority of cells are classified as having a peak in the PSD. As expected, when W is similar to the duration of the sequences (W ≥ 50 s), the fraction of oscillatory cells quickly drops. This fraction is no longer significantly above a chance level of 5%.
This approach was used for determining the fraction of oscillatory cells both in calcium imaging and in Neuropixels data. In the main text we present the results corresponding to W = 20 s.
Finally, we note that there was some variability in the frequency at which the PSD peaked across cells within a session. For example, in the example session shown in Fig. 1b–d and Fig. 2a, some singlecell PSDs peaked at a frequency of 0.0066 Hz, while others did so at a frequency of 0.0075 Hz. However, in many cases the PSDs were wide enough to exhibit high power in neighbouring frequencies too, providing support to the frequencies being rather clustered among a subset of values, with some slight variability around those values. When all cells were analysed (n = 6,231 cells pooled across 15 oscillatory sessions, 5 mice), in approximately half of the MEC data the oscillatory frequency at the singlecell level was very similar to the frequency at the population level (Extended Data Fig. 7e). This finding points to a small variability in the frequency of singlecell activity in MEC, as expected in the presence of recurring sequences.
Correlation and PCA sorting methods
To determine whether neural population activity exhibits temporal structure we visualized the population activity by means of raster plots in which we sorted all cells according to different methods.
Correlation method
This method sorts cells such that those that are nearby in the sorting are more synchronized than those that are further away. First, each calcium activity was downsampled by a factor 4 by calculating the mean over counts of calcium events in bins of 0.52 s. The obtained calcium activity was then smoothed by convolving it with a gaussian kernel of width equal to four times the oscillation bin size, a bin size that was representative of the temporal scale of the population dynamics (see ‘Oscillation bin size’). The cross correlations between all pairs of cells were calculated using time bins as data points, and a maximum time lag of 10 time points, equivalent to ∼ 5 s. This small time lag allowed us to identify near instantaneous correlation while keeping information about the temporal order of activity between cell pairs. The maximum value of the crosscorrelation between cell i and cell j was stored in the entry (i,j) of the correlation matrix C, which was a square matrix of N rows and N columns, where N was the total number of recorded neurons in the session with SNR > 4. If the crosscorrelation peaked at a negative time lag the value in the entry (i,j) was multiplied by −1. The entry with the highest crosscorrelation value was identified and its row, denoted by i_{max,} was used as the ‘seed’ cell for the sorting procedure and chosen to be the first cell in the sorting. Cells were then sorted according to the values in the entries \({(i}_{\max },j)\), \(j=\mathrm{1,2},\ldots ,N\), j ≠ i_{max}, that is, their correlations with the seed cell, in a descending manner.
PCA method
Computing correlations from the calcium activity or the calcium signals can be noisy due to fine tuning of hyperparameters (for example, the size of the kernel used to smooth the calcium activity of all cells). To avoid this, we leveraged the fact that the periodic sequences of neural activity constitute lowdimensional dynamics with intrinsic dimensionality equal to 1, and sorted the cells based on an unsupervised dimensionality reduction^{61} approach (a similar approach was used in ref. ^{62}). For each recording session, PCA was applied to the matrix of calcium activity (bin size = 129 ms; using Matlab’s builtin pca function), including all epochs of movement and immobility and using the rows (neurons) as variables and the columns (time bins) as observations. The first two principal components (PCs) were kept, since 2 is the minimum number of components needed to embed nonlinear 1dimensional dynamics. Cells were sorted according to their loadings in PC1 and PC2, expecting that the relationship between these loadings would express the ordering in cell activation during the sequences.
The plane spanned by PC1 and PC2 was named the PC1–PC2 plane. In the PC1–PC2 plane, the loadings of each neuron (the components of the eigenvectors without being multiplied by the eigenvalues) defined a vector, for which we computed its angle \({\theta }_{i}={\rm{arctg}}\,\left(\frac{{l}_{{\rm{PC}}2}^{i}}{{l}_{{\rm{PC}}1}^{i}}\right)\in \left[{\rm{\pi }},{\rm{\pi }}\right)\), 1 ≤ i ≤ N, with respect to the axis of PC1, where \({l}_{{\rm{PC}}j}^{i}\) is the loading of cell i on PCj. Cells were sorted according to their angle θ in a descending manner.
Note that while we keep the first 2 principal components to sort the neurons, all principal components and the full matrices of calcium activity were used in the analyses (except for visualization purposes—for example, see ‘Manifold visualization for MEC sessions’). Finally, note that because in PCA a principal component is equivalent to −1 times the principal component, the sorting and an inversion of the sorting are equivalent. The sorting was chosen so that sequences would progress from the bottom to the top in the raster plot.
The PCA method was used throughout the paper for sorting the recorded cells unless otherwise stated.
Random sorting of cell identities
A random ordinal integer \(\in [1,N]\), where N is the total number of recorded cells with SNR > 4, was assigned to each neuron without repetition across cells. Neurons were sorted according to those assigned numbers (see example session in Extended Data Fig. 4d, top row).
Sorting of circularly shuffled data
A shuffled matrix of calcium activity was built by circularly shuffling the calcium activity of each cell separately. For each cell a random ordinal integer \(\in [1,T]\), where T is the total number of time bins (bin size = 129 ms), was chosen and the calcium activity was rigidly shifted by this integer using periodic boundary conditions. The assignment of random ordinal integers was made separately for each cell. The PCA method was then applied to the shuffled matrix of calcium activity (see example session in Extended Data Fig. 4d, second row).
Sorting of temporally shuffled data
Because circularly shuffling the data preserves the oscillations in the singlecell calcium activity, a second shuffling approach was considered (for singlecell data shuffling procedures see ‘Autocorrelations and spectral analysis of singlecell calcium activity’). A shuffled matrix of calcium activity was built by temporally shuffling the calcium activity of each cell separately. For each cell, each time bin of the calcium activity was assigned a random ordinal integer \(\in [1,T]\) without repetition across time bins, where T is the total number of time bins (bin size = 129 ms), and time bins were ordered according to their assigned number. The assignment of random ordinal integers was made separately for each cell, so that the obtained random orderings were not shared across cells. The PCA method was then applied to the shuffled matrix of calcium activity.
Sortings are preserved when different portions of data are used for obtaining the sortings
To determine whether using different portions of the session for sorting the neurons lead to different sortings, the PCA method was applied to: (i) all data within a session; (ii) the first half of the session; and (iii) the second half of the session. This procedure gave three sortings per session. Next, for each cell pair in a session the distance between the two cells in each of the three sortings was calculated. We illustrate this calculation with a toy example: if 5 neurons were recorded, and sorting (i) was: (1,4,5,2,3), the distance between cells 1 and 5 was 2, because those two cells were 2 positions apart in the sorting. The distance between cells 1 and 3 was 1 and not 4, however, because in the calculation of distances we took into account that the sorting mirrors the position of the cells in the ring, which has periodic boundary conditions.
We next calculated the correlation between the distances in: sorting (i) versus sorting (ii), sorting (i) versus sorting (iii) and sorting (ii) versus sorting (iii). If sortings obtained with different portions of data preserve the ordering of the neurons, we would expect high correlation values. We compared the obtained correlation values with the 95th percentile of a shuffled distribution obtained by assigning, to each cell, a random position in each of the sortings.

Sorting (i) versus sorting (ii): 15 of 15 oscillatory sessions (see ‘Oscillation score’) were above the cutoff of significance. Correlation values in experimental data ranged from 0.38 to 0.85. The 95th percentile of shuffled data ranged from 0.004 to 0.015 (n = 15 in both experimental and shuffled data).

Sorting (i) versus sorting (iii): 15 of 15 oscillatory sessions were above the cutoff of significance. Correlation values in experimental data ranged from 0.52 to 0.86. The 95th percentile of shuffled data ranged from 0.005 to 0.013 (n = 15 in both experimental and shuffled data).

Sorting (ii) versus sorting (iii): 15 of 15 oscillatory sessions were above the cutoff of significance. Correlation values in experimental data ranged from 0.17 to 0.53. The 95th percentile of shuffled data ranged from 0.005 to 0.013 (n = 15 in both experimental and shuffled data).
The high correlation values obtained provide support for what is illustrated in Extended Data Fig. 4e: using different portions of data for sorting the cells unveils the same dynamics.
Sorting methods based on nonlinear dimensionality reduction techniques
The PCA method for sorting cells relies on a twodimensional linear embedding. This linear embedding might not be optimal if the population vectors describe temporal trajectories that, despite being lowdimensional, lie on a curved surface. To take into account potential nonlinearities, four additional sorting methods were implemented, based on the following nonlinear dimensionality reduction techniques^{63}: tdistributed stochastic neighbour embedding (tSNE), LEM, Isomap and uniform manifold approximation and projection (UMAP)^{64} (see parameters below). First, to express in the sortings the ordering of the cells during the slow temporal progression of the sequences, the four methods used a resampled matrix of calcium activity as input. To compute this matrix, for each session, we downsampled each calcium activity by a factor 4 by calculating its mean in bins of 0.52 s. The calcium activity of all cells was then smoothed by convolving them with a gaussian kernel whose width was given by the oscillation bin size (see ‘Oscillation bin size’). After applying tSNE, LEM, Isomap or UMAP to the resampled matrix of calcium activity, we kept the first two dimensions obtained with each method, for the same reasons as presented for the PCA sorting method. To obtain the sorting, the following procedure was applied: We let Dim1 and Dim2 be the first two dimensions obtained with the chosen dimensionality reduction technique that we had applied to the resampled matrix. In analogy with the PCA method, the Dim1–Dim2 plane was spanned by Dim1 and Dim2 and for each cell the components on those dimensions defined a vector in this plane for which the angle \(\theta \in \left[{\rm{\pi }},{\rm{\pi }}\right)\) with respect to the axis of Dim1 was computed. Cells were then sorted according to their angles in a descending manner.
To apply tSNE to the population activity we used a perplexity value of 50. First, we applied PCA to the resampled matrix of calcium activity, and then we used the projection of the neural activity onto the first 50 principal components as input to tSNE. To apply LEM to the population activity, we used as hyperparameters k = 15 and σ = 2. Similarly, we used k = 15 for running isomap. Finally, we used n_neighbors=30, min_dist=0.3 and correlation as metric for running UMAP.
We used the MATLAB implementation of UMAP^{65} and the Matlab Toolbox for Dimensionality Reduction (https://lvdmaaten.github.io/drtoolbox/). Finally, when displaying the raster plots that resulted from the different sortings, the first cell (located at the bottom of the raster plot) was always the same. This was accomplished by circularly shifting the cells in the different sortings such that the initial cell in all sortings coincided with the initial cell of the sorting obtained with the PCA method.
Manifold visualization for MEC sessions
Sorting the cells and visualizing their combined neural activity through raster plots revealed the presence of oscillatory sequences of neural activity in the recorded data. To visualize the topology of the manifold underlying the oscillatory sequences of activity, both PCA and LEM were used.
PCA was applied to the matrix of calcium activity, which first had each row convolved with a gaussian kernel of width equal to four times the oscillation bin size (see ‘Oscillation bin size’). The manifold was visualized by plotting the neural activity projected onto the embedding defined by PC1 and PC2. In Fig. 2c (left) the neural activity of the entire session was projected onto the lowdimensional embedding. In Extended Data Fig. 4c, the neural activity corresponding to the concatenated epochs of uninterrupted oscillatory sequences was projected onto the embedding.
For the LEM approach, first PCA was applied to the matrix of calcium activity, which was previously resampled to bins of 0.52 s as in ‘Sorting methods based on nonlinear dimensionality reduction techniques’, and the first five principal components were kept. Next LEM was applied to the matrix composed of the 5 principal components, using as parameters k = 15 and σ = 2. We decided to keep 5 principal components prior to applying LEM to denoise the data, for which we leveraged the fact that sequences of activity constitute lowdimensional dynamics with intrinsic dimensionality equal to 1, and therefore truncating the data to the first 5 principal components should preserve the sequential activity. The manifold was visualized by plotting the neural activity projected onto the embedding defined by the first two LEM dimensions. In Fig. 2c (right) the neural activity of the entire session was projected onto the embedding.
Both approaches revealed a ringshaped manifold along which the population activity propagated repeatedly with periodic boundary conditions. One sequence was equivalent to one full turn of the population activity along the ringshaped manifold. Finally, we note that when using PCA for visualizing the manifold, in some sessions the ring was less evident (Extended Data Fig. 4c). This is because the population activity had more variations from sequence to sequence, which resulted on the rings that corresponded to each sequence not completely overlapping in the PC1 versus PC2 plane. While recovering rings with PCA is challenging due to PCA being a linear method, using a nonlinear method would have helped in visualizing the ring (as in Fig. 2c, right), but we decided not to do this for all quantifications because nonlinear methods require more fine tuning and are usually harder to interpret.
Phase of the oscillation
To track the progression of the population activity over time, we leveraged the low dimensionality of the ringshaped manifold and the circular nature of the population activity, and parametrized the population activity with a single timedependent parameter, which we called the phase of the oscillation. Hence, the phase of the oscillation varied as a function of time (bin size = 129 ms) and tracked the progression of the neural population activity during the oscillatory sequences. The neural activity was projected onto a twodimensional plane using PCA. The use of PCA avoided the selection of hyperparameters, which is required in all nonlinear dimensionality reduction techniques including LEM. Let \({\rm{PC}}{i}_{t}\left(t\right)\) be the projection of the neural population activity onto principal component i (PCi). The neural population activity at time point t projected onto the plane defined by PC1 and PC2 is then given by (\({\rm{PC}}{1}_{t}\left(t\right),{\rm{PC}}{2}_{t}\left(t\right)\)), which defines a vector in this plane. The phase of the oscillation is defined as the angle of this vector with respect to the PC1 axis and is given by
$$\varphi \left(t\right)={\rm{arctg}}\,\left(\frac{{\rm{PC}}{2}_{t}\left(t\right)}{{\rm{PC}}{1}_{t}\left(t\right)}\right).$$
(1)
During one sequence, the phase of the oscillation continuously traversed the range \([{\rm{\pi }},{\rm{\pi }})\) rad, which was consistent with the population activity propagating through the network and describing one turn along the ringshaped manifold. The repetitive and almost linear dependence between the phase of the oscillation and time illustrates how stereotyped the sequences were (Fig. 2d).
We note that the quantity \(\varphi \left(t\right)\) is always defined, regardless of whether the session is or is not classified as oscillatory. In the case of the oscillatory sessions, \(\varphi \left(t\right)\) tracks the progression of the oscillatory sequences.
Joint distribution of crosscorrelation time lag and angular distance in the PCA sorting
To further characterize the sequential activation in the MEC neural population and to introduce a score that would determine the extent to which a session exhibited oscillatory sequences (see ‘Oscillation score’), we determined the relationship between the time lags that maximized the crosscorrelation between the calcium activity of two cells (τ) and their angular distances in the PCA sorting (d). In the plane generated by PC1 and PC2, the loadings of each neuron defined a vector, for which we computed the angle \({\theta }_{i}={\rm{arctg}}\,\left(\frac{{l}_{{\rm{PC}}2}^{i}}{{l}_{{\rm{PC}}1}^{i}}\right)\in \left[{\rm{\pi }},{\rm{\pi }}\right)\), 1 ≤ i ≤ N, with respect to the axis of PC1, where \({l}_{{\rm{PC}}j}^{i}\) is the loading of cell i on PCj and N is the total number of recorded neurons (see ‘Correlation and PCA sorting methods’). The angular distance d between any two cells in the PCA sorting was calculated as the difference between their angles wrapped in the interval \(\left[{\rm{\pi }},{\rm{\pi }}\right)\) (see Extended Data Fig. 5b, left),
$${d}_{i,j}={(\theta }_{i}{\theta }_{j}),$$
(2)
where \(1\le i\le N,1\le j\le N\). The Matlab function angdiff was used for computing this distance. Note that the angular distance maps how far apart two cells are in the raster plot when cells are sorted according to the PCA method.
To estimate the joint distribution of crosscorrelation time lags and angular distances in the PCA sorting, the cross correlations between all pairs of cells were calculated using a maximum time lag of 248 s. For each cell pair the time lag at which the crosscorrelation peaked (τ) and the angular distance in the PCA sorting (d) were calculated. A discrete representation was used for these two variables: in all analyses, and unless stated otherwise, the range of possible τ values—that is, [−248,248] s—was discretized into 96 bins of size \(\varDelta \tau =\frac{496\,{\rm{s}}}{96} \sim 5\,{\rm{s}}\) and the range of possible d values—that is, [−π, π) rad—was discretized into 11 bins of size \(\varDelta d=\frac{2{\rm{\pi }}}{11} \sim 0.57\,{\rm{rad}}\). Using those bins, the joint distribution of τ and d was expressed as a twodimensional histogram that counted the number of cell pairs observed for every combination of τ bins and d bins, normalized by the total number of cell pairs.
An example of joint distribution of crosscorrelation time lags and angular distances in the PCA sorting is presented in Extended Data Fig. 5b, right, built on the example session shown in Fig. 2a. In sessions with clear periodic sequences, the time lag τ increased with the distance d. This dependence was observed a discrete number of times in each session, which indicated that cells were active periodically and at a fixed frequency or at an integer multiple of it (see Extended Data Fig. 5c, top for another example with a different time scale). In sessions without detectable periodic sequences such structure was not observed (Extended Data Fig. 5c, bottom).
Oscillation score
While striking oscillatory sequences were observed in multiple sessions and mice, the population activity exhibited considerable variability, ranging from nonpatterned activity to highly stereotypic and periodic sequences (Extended Data Fig. 5a). This variability prompted us to quantify, for each session, the extent to which the population activity was oscillatory, which we did by computing an oscillation score. For each session, we first calculated the phase of the oscillation \(\varphi \left(t\right)\) (bin size = 129 ms, equation (1)), which tracks the progression of the population activity in the presence of oscillatory sequences (see ‘Phase of the oscillation’ and Fig. 2d). Next the PSD of \({\rm{si}}{\rm{n}}\left(\varphi \left(t\right)\right)\) was calculated using Welch’s method with Hamming windows of 17.6 min (8,192 bins of 129 ms in each window) and 50% of overlap between consecutive windows (pwelch Matlab function, see ‘Autocorrelations and spectral analysis of singlecell calcium activity’). If the PSD peaked at 0 Hz and the PSD was strictly decreasing, the phase of the oscillation was not oscillatory and hence the population activity was not periodic in the analysed session. In this case the oscillation score was set to zero. Otherwise, prominent peaks in the PSD at a frequency larger than 0 Hz were identified. In order to disentangle largeamplitude peaks from small fluctuations in the PSD, a peak at frequency f_{max} was considered prominent and indicative of periodic activity if its amplitude was larger than (1) 9 times the mean of the tail of the PSD (that is, <PSD(f > f_{max})>, where <x> indicates the average over frequencies x) and (2) 9 times the minimum of the PSD between 0 Hz and f_{max} (that is, min(PSD(f < f_{max}))). If no peak in the PSD met these criteria the oscillation score was set to zero. Otherwise, the presence of a prominent peak in the PSD calculated on \({\rm{si}}{\rm{n}}\left(\varphi \left(t\right)\right)\) was considered indicative of periodic activity at the population level. Yet a crucial component for observing oscillatory sequences is that cells fire periodically and that the time lag that maximizes the cross correlations between the calcium activity of pairs of cells that are located at a fixed distance in the sequence comes in integer multiples of a minimum time lag, which ensures that cells oscillate at a fixed frequency and that the calcium activity of one cell is temporally shifted with respect to the other. To quantify the extent to which these features were present in the data, we computed the joint distribution of time lags and angular distance in the PCA sorting (τ was discretized into 240 bins and d was discretized into 11 bins, see ‘Joint distribution of crosscorrelation time lag and angular distance in the PCA sorting’). Next for each bin i of d, \(1\le i\le 11\), we calculated the PSD of the distribution of τ conditioned on the distance bin i (Welch’s methods, Hamming windows of 128 τ bins with 50% overlap between consecutive windows, pwelch Matlab function). The presence of a peak in this signal indicated that for bin i of d, the time lag that maximizes the cross correlations between cells was oscillatory (that is, it peaked at multiples of one specific time lag), as expected when cells are active periodically with an approximately fixed frequency and also with harmonics of the primary frequency (see example joint distribution in Extended Data Fig. 5b, right). The presence (or absence) of a peak that satisfied the condition of being larger than (1) 10 times the mean of the tail of the PSD (same definition as above), and (2) 4.5 times larger than the minimum between 0 Hz and the frequency at which the PSD peaked, was identified (same definition as above, the parameters are different from the ones used above because the signals are very different). The oscillation score was then calculated as the fraction of angular distance bins for which a peak was identified.
Based on the bimodal distribution of oscillation scores obtained in the calcium imaging data from MEC (Extended Data Fig. 5d), a session was considered to express oscillatory sequences if the oscillation score was ≥0.72. This cutoff (0.72) corresponded to the smallest oscillation score within the group with high scores (shown in green in Extended Data Fig. 5d). Note that because the distribution of oscillation scores was bimodal any other choice of threshold between 0.27 and 0.72 would have led to the same results. Using as cutoff 0.72 was also equivalent to asking that at least 8 out of the 11 distributions of τ conditioned on bin i of d, \(1\le i\le 11\), had a significant peak in their PSD, which accounted for the fact that for distances in the PCA sorting that are close to zero, cells exhibit instantaneous coactivity rather than coactivity shifted by some specific time lag, which makes the conditional probability not oscillatory. After applying the cutoff, 15 of 27 calcium imaging sessions in MEC in 5 mice were classified as oscillatory (Extended Data Fig. 5d, shown in green), and among those 15 sessions, 10 were recorded with synchronized behavioural tracking (see ‘Selfpaced running behaviour under sensoryminimized conditions’). The number of recorded cells in the calcium imaging oscillatory sessions ranged from 207 to 520. In the rest of the calcium imaging data, 0 of 25 PaS sessions in 4 mice were classified as oscillatory, and 0 of 19 VIS sessions in 3 mice were classified as oscillatory.
Oscillation bin size
The oscillatory sequences progressed at frequencies <0.1 Hz that varied from session to session. The oscillation bin size was a temporal bin size representative of the time scale of the oscillatory sequences in each session. It was used to quantify singlecell and neural population dynamics, for which describing the neural activity at the right time scale was fundamental (for example, see ‘Transition probabilities’). For each oscillatory session the period of the oscillatory sequences, denoted by P_{osc}, was calculated as the inverse of the frequency f_{max} at which the PSD of the signal \(\sin \left(\varphi \left(t\right)\right)\) peaked (see equation (1) and ‘Oscillation score’), that is, \({P}_{{\rm{osc}}}={f}_{\max }^{1}\). Note that this estimate of the period was reliable when during most of the session the network engaged in the oscillatory sequences, in which case the estimate was equivalent to the length of the session divided by the total number of sequences. However, it became less reliable the more interrupted the oscillatory sequences were.
The oscillation bin size \({T}_{{\rm{osc}}}\) was computed as the period of the oscillatory sequences divided by 10,
$${T}_{{\rm{osc}}}=\frac{{P}_{{\rm{osc}}}}{10}=\frac{1}{10\times {f}_{\max }}.$$
(3)
This choice of bin size was made so that each sequence would progress across ∼10 time points. Across 15 oscillatory sessions, the oscillation bin size ranged from 3 to 17 s (see Extended Data Fig. 9d).
In sessions without oscillatory sequences, there was not a welldefined peak in the PSD of \(\sin \left(\varphi \left(t\right)\right)\), and therefore the oscillation bin size was not possible or meaningful to calculate. Yet, to perform the quantifications of network dynamics at temporal scales similar to the ones investigated in oscillatory sessions, the mean oscillation bin size computed across all oscillatory sessions was used (mean oscillation bin size = 8.5 s).
Unless otherwise indicated, the utilized bin size was 129 ms.
Identification of individual sequences
The characterization of the oscillatory sequences required multiple analyses that relied on identifying individual sequences, for example to quantify the duration of the sequences and their variability. The procedure for identifying individual sequences was based on finding the time points at which each sequence began (visualized typically at the bottom of the raster plot) and ended (visualized typically at the top of the raster plot, see Extended Data Fig. 6a). Note that the beginning and the end of the sequence are arbitrary because of the periodic boundary conditions in the sequence progression, and therefore a different pair of phases that are 2π apart could have been used for defining the beginning and the end of the sequence.
One sequence was equivalent to one full turn of the population activity around the ringshaped manifold—that is, during one sequence the phase of the oscillation traversed 2π (see ‘Phase of the oscillation’). To calculate the phase of the oscillation and determine the time epochs during which it traversed 2π, we smoothed the calcium activity of all cells (bin size = 129 ms) using a gaussian kernel of width equal to the oscillation bin size. Next, the phase of the oscillation was calculated and discretized into 10 bins (that is, the range \([{\rm{\pi }},{\rm{\pi }})\) was discretized into 10 bins). Time points at which the phase of the oscillation belonged to a bin that was 3 or more bins away from the bin in the previous time point were considered as discontinuity points and were used to define the beginning and the end of putative sequences. Putative sequences were classified as sequences if the phase of the oscillation smoothly traversed the range \(\left[{\rm{\pi }},{\rm{\pi }}\right)\) rad in an ascending manner. To account for variability, decrements of up to 1 bin of the phase of the oscillation were allowed. This means that there could be fluctuations of up to 0.6 rad in the phase within one individual sequence, and still be considered a sequence. Points of sustained activity were disregarded. Segments of sequences in which the phase of the oscillation covered at least 5 bins (that is, 50% or more of the range \(\left[{\rm{\pi }},{\rm{\pi }}\right)\) rad) were also identified.
Sequence duration, sequence frequency and ISI
The duration of individual sequences was defined as the amount of time that it takes the phase of the oscillation to cover the range \(\left[{\rm{\pi }},{\rm{\pi }}\right)\) in a smooth and increasing manner, which is consistent with the population activity completing one full turn along the ringshaped manifold. To calculate the sequence duration, the time interval between the beginning and the end of the sequence was determined (see ‘Identification of individual sequences’).
To quantify the variability in sequence duration within and between sessions, two approaches were adopted. In approach 1 (Extended Data Fig. 6f left), the s.d. of sequence durations was computed for each oscillatory session. To estimate significance, in each of 500 iterations all sequences across 15 oscillatory sessions were pooled (421 sequences in total) and randomly assigned to each session while keeping the original number of sequences per session unchanged. For each iteration the s.d. of the sequence durations randomly assigned to each session was calculated. In approach 2 (Extended Data Fig. 6f, right), for each session i, 1 ≤ i ≤ 15, where 15 is the total number of oscillatory sessions, we considered all pairs of sequences within session i (within session group) or alternatively all pairs of sequences such that one sequence belongs to session i and the other sequence to session \(j,j\ne i\) (between session group). For each sequence pair in each group, the ratio between the shortest sequence duration and the longest sequence duration was calculated. The mean was computed over pairs of sequences in each group for each session separately. Notice that the larger this ratio the more similar the sequence durations are.
The sequence frequency was calculated as the total number of identified individual sequences in a session, divided by the total amount of time the network engaged in the oscillatory sequences during the session, which was computed as the length of the temporal window of concatenated sequences.
The ISI was defined as the length of the epoch from the termination of one sequence and the beginning of the next one. In other words, the ISI was calculated as the amount of time that elapsed between the time point at which the phase of the oscillation reached π (after completing one turn along the ringshaped manifold), and the time point at which it is equal to −π (prior to initiating the next turn along the ring).
Mean event rate during segments of the sequences
To determine how population activity varied during individual sequences (Extended Data Fig. 6c), the following approach was adopted. For each oscillatory session (see ‘Oscillation score’) all individual sequences were identified (see ‘Identification of individual sequences’). Each sequence was divided into ten segments of equal length. For each sequence segment, the mean event rate was calculated as the total number of calcium events across cells divided by sequence segment duration and number of cells. For each session the mean event rate per segment was calculated over sequences. Across sessions we found that the percentage rate change from the segment with the minimum event rate to the segment with the maximum rate was no more than 18% (Extended Data Fig. 6c).
Analysis of Neuropixels data
Neuropixels data was different from the calcium imaging data in that it consisted of spike times and not calcium traces. Despite this fundamental difference, for most of the analyses we applied the same methods to both datasets. When this was not possible (see below), we tried to minimize the differences between the two analyses pipelines.
Spike matrices
In order to create arrays that were similar to the matrices of calcium activity, for each recorded unit a spike train was built using a bin size of 120 ms (similar to the bin size used in calcium imaging data, 129 ms). Each time bin contained the number of spikes produced by the recorded unit in that bin. Spike matrices were built by stacking the spike trains of all recorded units (469 units in the example session presented in Fig. 2f, 410 units in the example session shown in Extended Data Fig. 4g).
Calcium traces are temporally correlated due to the slow dynamics of the calcium indicator. In addition, the observed periodic sequences unfolded over a time scale of minutes. To take these two factors into account, we smoothed the spike train of each recorded unit with a Gaussian kernel of width equal to 5 s.
Both the original spike matrix and the smoothed spike matrix were then binarized using, for each spike train, a threshold equal to the mean plus either 1 or 1.5 times the s.d. (1 for smoothed matrices; 1.5 for nonsmoothed matrices; as a reference, the threshold for binarization used in calcium data was the mean plus 1.5 times the s.d.; see ‘Binary deconvolved calcium activity and matrix of calcium activity’).
In the calcium imaging experiment, it took approximately 5 min to initiate the recording after the mouse was positioned on the wheel (mainly due to the time that was needed to find the imaging planes). In the Neuropixels data there was no such delay between positioning the mice on the wheel and starting the data acquisition. In order to make both datasets as comparable as possible, and in order to remove any effects due to arousal, the first 5 min of the Neuropixels sessions were discarded.
Autocorrelation and spectral analysis
The autocorrelations were calculated on the spike trains (without smoothing), and the PSD was calculated on the autocorrelations. Methods and parameters used for calculating the autocorrelation and PSDs were the same as in calcium imaging data (‘Autocorrelations and spectral analysis of singlecell calcium activity’).
Calculation of oscillation score
As in the calcium imaging data, in order to quantify the amount of oscillatory activity in the Neuropixels sessions, an oscillation score was computed. Because in the Neuropixels recordings (unlike in the calcium imaging data) there were some long periods of nonsequence activity between bouts of periodic sequences, possibly due to small differences in training protocol, we computed the oscillation score not on the full spike matrix but on the matrix of concatenated sequences (built by identifying all individual sequences in the smoothed spike matrix and concatenating them as described for the calcium imaging data in ‘Identification of individual sequences’ and ‘Sequence duration, sequence frequency and ISI’ above).
Sorting calculation and raster plot visualization
Neural population activity was visualized by means of raster plots, for which units were sorted using the PCA method (‘Correlation and PCA sorting methods’). The sorting was calculated on the smoothed spike matrix (Fig. 2f and Extended Data Fig. 4g, top), and the obtained sorting was applied also to the nonsmoothed spike matrices (Extended Data Fig. 4f,g, bottom).
While the sorting and visualization of neural population activity were performed as we did in calcium imaging data, there was one difference in how the two datasets were analysed. Because in the Neuropixels data the periodic sequences were more salient in some subsets of the sessions than others, for visualization purposes we calculated the sorting on a subset of the smoothed transition matrices. Those subsets are given by [1,200, 1,700] s for the example session of mouse no. 104368 (Fig. 2f) and [1,100, 1,400] s for the example session of mouse no. 102335 (Extended Data Fig. 4g). Note, however, that sequences were identified outside these session subsets too, indicating that the sorting unveils stereotyped sequences also outside the used subsets of data (see ‘Sortings are preserved when different portions of data are used for obtaining the sortings’).
Locking to the phase of the oscillation
To calculate the extent to which individual cells in the calcium imaging experiments were tuned to the oscillatory sequences, two quantities were used: the locking degree and the mutual information between the calcium event counts and the phase of the oscillation. For each oscillatory session, the phase of the oscillation \(\varphi \left(t\right)\) was computed (see equation (1)) and individual sequences were identified (see ‘Identification of individual sequences’). Next, the time points that corresponded to all individual sequences in one session were concatenated, which generated a new signal with the phase of the oscillation for all consecutive sequences, and a new matrix of calcium activity in which the network engaged in the oscillatory sequences uninterruptedly.
The locking degree was computed for each cell as the mean resultant vector length over the phases of the oscillatory sequences at which the calcium events occurred (bin size = 129 ms, function circ_r from the Circular Statistics Toolbox for Matlab^{66}). The locking degree has a lower bound of 0 and upper bound of 1. It is equal to 1 if all oscillation phases at which the calcium events occurred are the same (that is, perfect locking), and equal to zero if all phases at which the calcium events occurred are evenly distributed (total absence of locking). To estimate significance, for each cell a null distribution of locking degrees was built by temporally shuffling the calcium activity of that cell 1,000 times while the phase of the oscillation remained unchanged, and by computing, for each shuffle realization, the locking degree (shuffling was performed as in ‘Sorting of temporally shuffled data’). The 99th percentile of the estimated null distribution was used as a threshold for significance.
In order to assess the robustness of the locking degree, the obtained results were compared with a second measure based on information theory^{67}: the mutual information between the counts of calcium events (event counts) and the phase of the oscillation (bin size = 0.52 s). To estimate the reduction in uncertainty about the phase of the oscillation (P) given the event counts of the calcium activity (S), Shannon’s mutual information was computed as follows^{68}:
$${\rm{MI}}(S,P)=\sum _{p,s}{\rm{Prob}}(p,s){\log }_{2}\frac{{\rm{Prob}}(p,s)}{{\rm{Prob}}(p){\rm{Prob}}(s)},$$
where \({\rm{Prob}}\left(p,s\right)\) is the joint probability of observing a phase of the oscillation p and an event count s, \({\rm{Prob}}\left(s\right)\) is the marginal probability of event counts and \({\rm{Prob}}\left(p\right)\) is the marginal probability of the phase of the oscillation. All probability distributions were estimated from the data using discrete representations of the phase of the oscillation and the event counts. The event counts were partitioned into s_{max} + 1 bins to account for the absence of event counts as well as all possible event counts, where s_{max} is the maximum number of event counts per cell in a 0.52 s bin, and the phase of the oscillation was discretized into 10 bins of size \(\frac{2{\rm{\pi }}}{10}\).
The mutual information is a nonnegative quantity that is equal to zero only when the two variables are independent—that is, when the joint probability is equal to the product of the marginals \({\rm{Prob}}\left(p,s\right)={\rm{Prob}}\left(p\right){\rm{Prob}}\left(s\right)\). However, limited sampling can lead to an overestimation in the mutual information in the form of a bias^{69}. In order to correct for this bias, the calcium activity was temporally shuffled (as in ‘Sorting of temporally shuffled data’) and the mutual information between the event counts of the shuffled calcium activity and the phase of the oscillation, which remained unchanged, was calculated. This procedure, which destroyed the pairing between event counts and phase of the oscillation, was repeated 1,000 times and the average mutual information across the 1,000 iterations was computed and used as an estimation of the bias in the mutual information calculation. In the right panel of Fig. 3a, we report both the mutual information and the bias. In Extended Data Fig. 7a, the corrected mutual information was reported (MI_{c}), where the bias (⟨MI_{sh}⟩_{iterations}) was subtracted out from the Shannon’s mutual information (MI): MI_{c} = MI − ⟨MI_{sh}⟩_{iterations}.
Note that the locking degree and the mutual information between the event counts and the phase of the oscillation yielded consistent results (see Fig. 3a and Extended Data Fig. 7a).
Tuning of single cells to the phase of the oscillation
The selectivity of each cell to the phase of the oscillation in the calcium imaging data was visualized through tuning curves and quantified through their preferred phase. As in the analysis of ‘Locking to the phase of the oscillation’, the phase of the oscillation \(\varphi \left(t\right)\) was computed, individual sequences were identified, and the time points of the phase of the oscillation and the matrix of calcium activity that corresponded to all individual sequences in one session were concatenated.
Tuning curves
The range of phases \(\left[{\rm{\pi }},{\rm{\pi }}\right)\) rad was partitioned into 40 bins of size \(\frac{2{\rm{\pi }}}{40}\) rad. For each cell the tuning curve in the phase bin j, j = 0,…,39, was calculated as the total number of event counts that occurred at phases within the range \(\left[{\rm{\pi }}+j\frac{2{\rm{\pi }}}{40},{\rm{\pi }}+(\,j+1)\frac{2{\rm{\pi }}}{40}\right)\) divided by the total number of event counts during the concatenated oscillatory sequences.
Preferred phases
The preferred phase of each cell was calculated as the circular mean over the oscillation phases at which the calcium events occurred (function circ_mean from the Circular Statistics Toolbox for Matlab^{66}). In most of the analysis the preferred phase was calculated, for each cell, after concatenating all sequences. However, in a subset of analyses (see ‘Anatomical distribution of preferred phases’), the preferred phase was also calculated for individual sequences, as the circular mean over the oscillation phases at which the calcium events occurred in each sequence.
Unless otherwise stated, the preferred phase refers to the calculation performed on concatenated sequences (and not on individual sequences).
Distribution of preferred phases
To determine the extent to which the preferred phases across locked cells were uniformly distributed in one recorded session, the distribution of the cells’ preferred phases, that we shall denote Q, was estimated by discretizing the preferred phases into 10 bins of size \(\frac{2{\rm{\pi }}}{10}\) rad. The entropy of this distribution \({H}_{Q}={\sum }_{x=1}^{10}Q\left(x\right){\log }_{2}\left(Q\left(x\right)\right)\) was calculated and used to compute the entropy ratio H_{ratio} which quantifies how much Q departs from a flat distribution:
$${H}_{{\rm{ratio}}}=\frac{{H}_{Q}}{{H}_{{\rm{flat}}}}$$
(5)
where \({H}_{{\rm{flat}}}\) is the entropy of a flat distribution using 10 bins—that is, \({H}_{{\rm{flat}}}=3.32\) bits. The closer \({H}_{{\rm{ratio}}}\) is to 1 the flatter Q is, and therefore all preferred phases tend to be equally represented. The smaller \({H}_{{\rm{ratio}}}\) is, the more uneven Q is and some preferred phases tend to be more represented than others.
To estimate significance, for each session the procedure for calculating \({H}_{{\rm{ratio}}}\) was repeated for 1,000 iterations of a shuffling procedure where the preferred phase of the cells was calculated after the values of the phase of the oscillation were temporally shuffled. In Extended Data Fig. 7c, both panels, for each session the 1,000 shuffle realizations were averaged.
Participation index
The Participation Index (PI) quantifies the extent to which a cell’s calcium events were distributed across all sequences, or rather concentrated in a few sequences. For neurons that were active only in a few sequences the participation index was small (participation index ∼ 0), and for neurons that were reliably active during most of the sequences the participation index was high (participation index ∼ 1; Extended Data Fig. 7g shows three example neurons of the session in Fig. 2a).
The participation index was calculated for each cell separately as the fraction of sequences needed to account for 90% of the total number of calcium events. To compute the participation, individual sequences were identified (see ‘Identification of individual sequences’), and for each cell the number of calcium events per sequence was calculated and normalized by the total number of calcium events across all concatenated sequences, which yields the fraction of calcium events per sequence. This quantity was sorted in an ascending manner and its cumulative sum was calculated. The participation index is the minimum fraction of the total number of sequences for which the cumulative sum of the fraction of calcium events per sequence ≥0.9 (results remain unchanged when the cumulative sum is required to be ≥0.95).
Relationship between tuning to the phase of the oscillation and singlecell oscillatory frequency
To determine whether the frequency of oscillation of singlecell calcium activity was correlated with the extent to which the cell was locked and participated in the oscillatory sequences, for each cell the ratio between its oscillatory frequency (see ‘Autocorrelations and spectral analysis of singlecell calcium activity’) and the sequence frequency (see ‘Sequence duration, sequence frequency and ISI’) was calculated and denoted relative frequency. Next, for each session cells were divided into two groups: one group had cells with relative frequency ~1 (cells whose oscillatory frequencies were most similar to the sequence frequency), and the other group had cells with relative frequency ≠1 (cells whose oscillatory frequencies were most different from the sequence frequency). The size of each group was the same and was given by a percentage α of the total number of recorded cells in a session. For each group the locking degree (see ‘Locking to the phase of the oscillation’) and the participation index (see ‘Participation index’) were compared. For the quantification across all 15 oscillatory sessions, the mean locking degree and participation index were calculated for each group separately and for each session separately, and all 15 sessions were pooled. α varied from 5% to 50%.
Anatomical distribution of preferred phases
To determine whether the entorhinal oscillatory sequences resembled travelling waves, during which neural population activity moves progressively across anatomical space^{20,21,70,71,72,73,74}, we took three complimentary approaches.
Correlation between differences in preferred phase and anatomical distance
Preferred phases calculated using data from the entire session (after concatenating individual sequences)
For each of the 15 oscillatory sessions (across 5 mice) the Pearson correlation between the anatomical distance between cells in the FOV and the difference in their preferred phases (see ‘Tuning of single cells to the phase of the oscillation’) was calculated. In order not to count the same data twice, each correlation value was calculated using N × (N − 1)/2 samples (each sample was a cell pair), where N was the total number of cells recorded in the session. In the presence of travelling waves, a significant correlation between differences in preferred phase and anatomical distance between cells within the FOV is to be expected. To determine statistical significance the cells’ preferred phase were shuffled within the FOV 100 times, and for each shuffled realization the correlation values were calculated. Because we were interested in significant correlations, regardless of whether they were positive or negative, both in experimental and shuffled data we took the absolute value of the correlations. Next, the 95th percentile of the shuffled distribution (100 shuffled realizations per session) was used as cutoff for significance and compared with the correlation value in experimental data.
In order to rule out that the small correlation values observed in experimental data could be masking a dependency such that for larger distances the differences in preferred phase increased in absolute value, the same calculations were repeated but now taking the absolute value of the difference in preferred phase. Statistical significance was determined as in the previous paragraph.
Preferred phases calculated using data from individual sequences
Travelling waves could still be present if they move in different directions from sequence to sequence. To test for the presence of travelling waves without assuming similar wave directions across successive sequences, the quantification of correlation between the difference in preferred phase as a function of pairwise anatomical distance was repeated for each sequence separately. To calculate the preferred phase of each cell in each sequence (see ‘Tuning of single cells to the phase of the oscillation’), the mean phase at which the calcium events occurred in that individual sequence was computed. In each sequence, only cells that had at least 5 calcium events were included in the analysis. This analysis was performed separately on 421 sequences across 15 oscillatory sessions. Similarly to the analysis described above, when sequences were concatenated within a session, the calculations were repeated after taking the absolute value of differences in preferred phase.
Results are presented in Fig. 3f,g. In Fig. 3f, the correlation value was also nonsignificant when calculated using the absolute value of the differences in preferred phase (correlation = 0.0028, cutoff for significance of the correlation = 0.0146). In Fig. 3g, in the experimental data the absolute value of the correlations ranged from 6.4 × 10^{−6} to 0.147 (n = 421). Out of 421 sequences, 27 were classified as significant when compared to the 95th percentile of a shuffled distribution (cutoffs ranged from 0.007 to 0.237, n = 421). The fraction 27/421 was slightly above a chance level of 0.05 (0.05 × 421 = 21 sequences), yet for those 27 sequences the correlation values were very low, ranging from 0.008 to 0.137.
Calculation of local gradients of preferred phase
Previous studies have investigated the presence of travelling waves by computing local anatomical gradients of the phase of the oscillation, when the phase is calculated through the Hilbert transform applied to the activity of each electrode (for example, ref. ^{75}, Ecog data). In order to perform a similar analysis but applied to each sequence separately, two different approaches were taken.
Similarity of preferred phases in spatial bins of the FOV
First, the similarity in preferred phases of all cells within spatial bins of the FOV was used as a proxy for local gradients. The similarity in preferred phases was calculated as the mean vector length (MVL) of the distribution of preferred phases within each bin of the FOV. The analysis was performed for individual sequences separately.
For each of the 15 oscillatory sessions (over 5 mice), the FOV was divided into spatial bins of 100 μm x 100 μm (6 × 6 bins in 10 sessions, 10 × 10 bins in 5 sessions), or 200 μm x 200 μm (3 × 3 bins in 10 sessions, 5 × 5 bins in 5 sessions) (note that for 10 of the 15 oscillatory sessions the FOV was 600 μm x 600 μm, mice no. 60355, no. 60584, no. 60585; while for 5 of the 15 oscillatory sessions the FOV was 1,000 μm × 1,000 μm, mouse no. 59914; mouse no. 59911 did not show the oscillatory sequences). Next, the preferred phase of each cell per sequence was calculated (as we did in ‘Correlation between differences in preferred phase and anatomical distance’) and for each sequence and every spatial bin of the FOV the MVL was computed (only spatial bins with 10 or more cells were considered). If the MVL was 0, then all preferred phases in that bin were different and homogeneously distributed between −π and π, whereas if the MVL was 1 then all preferred phases were the same. In the presence of a travelling wave, each bin should have a high MVL value compared to chance levels. Statistical significance was determined by repeating the same MVL calculation after shuffling the cells’ preferred phases within the FOV 200 times, and using, for each spatial bin, a cutoff for significant of 95th percentile of the shuffled distribution. A nonsignificant fraction of spatial bins had a MVL value above the cutoff for significance.
Differences in preferred phase among pairs of cells in small neighbourhoods of the spatial domain
The analysis presented above is focused on the degree of similarity between preferred phases in spatial bins. In order to avoid small cell sample effects, and effects of adding a threshold number of cells for bins to be included when calculating similarity with the MVL measure above, we decided to also calculate the difference in preferred phases for all pairs of cells that were located within small neighbourhoods in the FOV, expecting that in the presence of travelling waves the differences in preferred phases of cell pairs within small neighbourhoods would be smaller than expected by chance. For each cell in the FOV, all other cells that were located within a circular neighbourhood of radius 50, 100 or 200 μm were identified and the differences in preferred phase between cell pairs within those areas were calculated. Next, for each sequence and each radius separately all phase differences were pooled, and the mean and the median of the obtained distributions were calculated. To determine significance, the preferred phases across all cells were shuffled 200 times and for each shuffled realization a distribution of differences in preferred phase was obtained and used to calculate the mean and median. Because in the presence of travelling waves smaller differences in preferred phases than in the shuffled data were expected, the mean and median calculated on experimental data were compared with the 5th percentile of the distribution of means and medians obtained from shuffled data. This comparison was performed for each sequence and each radius separately.
Centreofmass calculation of the population activity
To determine whether the population calcium activity was anatomically localized, as expected in the presence of travelling waves, we calculated its centre of mass (COM). First, all individual sequences were identified and the neural data was averaged in time bins of 5 s. We chose bins of 5 s because the sequences are very slow, however, results remain unchanged if bins of 1 s or 2 s are used instead. For each time point (bin size = 5 s) and for each sequence separately the COM of the population activity was calculated as:
$${\rm{COM}}=\frac{1}{M}\mathop{\sum }\limits_{i=1}^{N}{m}_{i}{{\bf{r}}}_{i},$$
where N is the total number of recorded cells in the session, r_{i} is the position of neuron i in the FOV, m_{i} is the total number of calcium events of neuron i within the 5 s time bin, and \(M={\sum }_{i=1}^{N}{m}_{i}\). The COM was visualized for one example sequence both in experimental data, and after randomly shuffling the position of the cells within the FOV (Extended Data Fig. 8d). To quantify the temporal trajectory of the COM across individual sequences, we calculated the cumulative distance travelled by the COM as the sum of the distances travelled by the COM between consecutive time points (bin size = 5 s). The cumulative distance travelled calculated on experimental data was compared with the 5th and 95th percentile of a distribution built by shuffling the positions of the cells in the FOV 500 times.
Procedure for merging steps
In order to average out the variability observed in single cells at the level of locking degree and participation index while preserving the temporal properties of the oscillatory sequences, an iterative process that defines new variables from combining the calcium activity of cells was implemented for each session separately (Extended Data Fig. 9a). This process is similar to a coarsegraining approach^{76}.
First, the N recorded cells in one session were sorted according to the PCA method. In the first iteration of the procedure, named merging step one, the calcium activity (see ‘Binary deconvolved calcium activity and matrix of calcium activity’) of pairs of cells that were positioned next to each other in the PCA sorting were added up (merging step 1 in Extended Data Fig. 9a). This resulted in \(\frac{N}{2}\) new variables, which in merging step 2 were grouped together in pairs of adjacent variables by adding up their activity, which yielded \(\frac{N}{4}\) new variables. Note that because in the PCA sorting cells whose activity is synchronous are positioned adjacent to each other, the new variables consist of groups of coactive cells.
In general, merging step j generates \(\frac{N}{{2}^{j}}\) variables by adding up the activity of pairs of \(\frac{N}{{2}^{j1}}\) variables from merging step j − 1, j > 1, with each new variable defined as:
$${\mathop{\sigma }\limits^{ \sim }}_{i}=\frac{{\sigma }_{2i1}+{\sigma }_{2i}}{2}\,\,\,\,\,i=1,\ldots ,\frac{N}{{2}^{j}}$$
where \({\widetilde{\sigma }}_{i}\) is the ith new variable that results from adding \({\sigma }_{2i1}\) and \({\sigma }_{2i}\), which were computed in the previous merging step, \(j1\). In merging step 1, \({\sigma }_{2i1}\) and \({\sigma }_{2i}\) are the calcium activity of cells in the position \(2i1\) and \(2i\), \(1\le i\le N\), in the sorting obtained with the PCA method.
This procedure was repeated 6 times until ~10 variables were obtained in each session (the exact number of variables depended on the number of recorded cells, N, in each session). If N was an odd number, the last cell in the sorting obtained with the PCA method was discarded and the procedure was applied to the first N − 1 cells in the sorting. In every merging step the participation index (see ‘Participation index’) of each new variable was calculated (see Extended Data Fig. 9b).
Division of cells into ensembles
After 5 merging steps (and for approximately 10 variables), the participation index reached a plateau (Extended Data Fig. 9b). This motivated the decision to split the recorded cells into 10 variables, which we later used to quantify the population dynamics (see ‘Analysis of population dynamics using ensembles of coactive cells’). From now on we will refer to those variables as ensembles, to highlight the fact that cells in each ensemble are coactive. The same number of ensembles was used in sessions that did not exhibit oscillatory sequences.
To distribute cells into 10 ensembles, cells were sorted according to the PCA method. If \(\frac{N}{10}\) is an integer, where N is the total number of cells in one session, then each ensemble contains \(\frac{N}{10}\) cells and the set of cells that belong to ensemble i, 1 ≤ i ≤ 10, is \(\left\{\left(i1\right)\times \frac{N}{10}+1,\left(i1\right)\times \frac{N}{10}+2,\ldots ,i\times \frac{N}{10}\right\}\). If \(\frac{N}{10}\) is not an integer then ensembles 1 to 9 contain \(\lfloor \frac{N}{10}\rfloor \) cells and ensemble 10 contains \(N9\times \lfloor \frac{N}{10}\rfloor \) cells, where \(\lfloor x\rfloor =max\{m\in {\mathbb{N}}/m\le x\}\) and \({\mathbb{N}}\) is the set of natural numbers. In this case the set of cells that belongs to each ensemble is:
$$\left\{\begin{array}{l}\left\{(i1)\times \lfloor \frac{N}{10}\rfloor +1,(i1)\times \lfloor \frac{N}{10}\rfloor +2,\ldots ,i\times \lfloor \frac{N}{10}\rfloor \right\},\,1\le {\rm{e}}{\rm{n}}{\rm{s}}{\rm{e}}{\rm{m}}{\rm{b}}{\rm{l}}{\rm{e}}\le 9\\ \left\{9\times \lfloor \frac{N}{10}\rfloor +1,9\times \lfloor \frac{N}{10}\rfloor +2,\ldots ,N\right\},\,{\rm{e}}{\rm{n}}{\rm{s}}{\rm{e}}{\rm{m}}{\rm{b}}{\rm{l}}{\rm{e}}=10\end{array}\right.$$
Note that each cell was assigned to only one ensemble.
After each cell was assigned to one of the ten ensembles, the activity of each ensemble as a function of time was calculated as the mean calcium activity across cells in that ensemble.
Finally, to calculate the oscillation frequency of ensemble activity, the PSD was calculated (Welch’s methods, 8.8 min Hamming window with 50% overlap between consecutive windows, pwelch Matlab function). The oscillation frequency was estimated as the frequency at which the PSD peaked. For each session, the oscillation frequency of the activity of the ensembles was compared to the sequence frequency, which was computed as the total number of sequences in the session divided by the amount of time the network engaged in the oscillatory sequences. The latter was calculated as the length of the temporal window of concatenated sequences (see ‘Identification of individual sequences’).
Analysis of population dynamics using ensembles of coactive cells
We adopted an ensemble approach to quantify the population dynamics (see ‘Procedure for merging steps’ and ‘Division of cells into ensembles’). With a total of 10 ensembles this approach averaged out the variability observed in singlecell locking degree and participation index while keeping the temporal progression of the oscillatory sequences (Extended Data Fig. 9f). In sessions with oscillatory sequences, all individual sequences were identified (see ‘Identification of individual sequences’) and the corresponding time bins were concatenated, which yielded a new matrix of calcium activity in which the oscillatory sequences were uninterrupted. Next, cells were divided into ensembles (see ‘Division of cells into ensembles’) and ensemble activity was downsampled using as bin size the oscillation bin size of the session (see ‘Oscillation bin size’). This procedure yielded a matrix, the ensemble matrix, with the activity of each ensemble corresponding to a single row (10 rows in total), and as many columns as time points sampled at the oscillation bin size. In nonoscillatory sessions, the full matrix of calcium activity was used and the temporal downsampling was conducted at the mean oscillation bin size computed across all 15 oscillatory sessions; that is, bin size = 8.5 s (see ‘Oscillation bin size’ for a description of the bin size used in nonoscillatory sessions). For both types of sessions (with and without oscillations), the activity of the 10 ensembles was described through a vector expressing, at each time point, the ensemble number with the highest activity at that time point (see Extended Data Fig. 9e,f). This vector was used to perform the following analyses: transition probabilities, probability of sequential activation of ensembles, and sequence score.
Transition probabilities
The transition probability from ensemble i to ensemble j was quantified as the number of times the transition \(i\to j\) was observed in the data of one session, normalized by the total number of transitions in one session. Transitions were identified from the vector that contained the ensemble number with maximum activity at each time point (transitions to the same ensemble between consecutive time points were disregarded). Transitions were allocated in a matrix of transition probabilities T of size 10 × 10, since 10 ensembles were used. In this matrix, the component \(\left(i,j\right)\) expressed the transition probability from ensemble i to ensemble j.
To establish statistical significance of the transition probabilities, the data was shuffled 500 times. In each shuffle realization, each row of the matrix of calcium activity (with concatenated sequences in the case of oscillatory sessions) was temporally shuffled (as in ‘Sorting of temporally shuffled data’), and the procedure for calculating the ensemble matrix and transition probabilities was applied to the shuffled data. For each transition, \(i\to j\) the 95th percentile of the shuffled distribution was used to define a cutoff.
Probability of sequential activation of ensembles
We calculated the probability of sequential ensemble activation according to the following procedure. From the vector expressing the ensemble number with the highest activity at each time point (sampled at the oscillation bin size), strictly increasing sequences of all possible lengths (from 2 to 10 ensembles) were identified. The number of ensembles in each sequence was the number of ensembles that were active in consecutive time points (epochs of sustained activity were disregarded). While the sequences had to be strictly increasing, they did not have to be continuous. Sequences could skip ensembles, in which case the maximum number of ensembles in one sequence was less than 10. The probability of the sequential activation of k ensembles, k = 2,…,10, was next estimated as the number of times a sequence of k ensembles was found, normalized by the total number of identified sequences. Note that all subsequences were also included in this estimation. For example, if the ensembles 1, 2 and 3 were active in consecutive time points, a sequence of three ensembles was identified, as well as three subsequences of two ensembles each: 1, 2, as well as 2, 3 and 1, 3.
In order to test for significance, the shuffled data from ‘Transition probabilities’ was used. The procedure to compute the probability of sequential activation of ensembles was applied to each of the 500 shuffle realizations performed per session. Shuffled data was compared with recorded data.
Sequence score
The sequence score measures how sequential the ensemble activity is. It is calculated from the probability of sequential activation of ensembles as the probability of observing sequences of three or more ensembles. The sequence score was calculated for each session of the dataset separately. To determine if the obtained scores were significant, for each session the 500 shuffle realizations used in ‘Probability of sequential activation of ensembles’ for assessing significance of the probability of sequential activation of ensembles were used to calculate the sequence score on shuffled data. Those values were used to build a shuffled distribution, and the 99^{th} percentile of this distribution was chosen as the threshold for significance.
Estimation of number of completed laps on the wheel, speed and acceleration
Features of the mouse’s behaviour were used to determine whether the MEC oscillatory sequences were modulated by running.
The wheel had a radius of 8.54 cm (see ‘Selfpaced running behaviour under sensoryminimized conditions’) and a perimeter of 53.66 cm. Therefore mice had to run for ∼53.7 cm to complete one lap on the wheel. For each session, we estimated the number of completed laps on the wheel from the position on the wheel recorded as a function of time. The number of completed laps during one sequence (see ‘Identification of individual sequences’) was calculated as the total distance run during the sequence divided by 53.7 cm.
The speed of the mouse was numerically calculated as the first derivative of the position on the wheel as a function of time (the sampling frequency of the position was 40 Hz for mice 60355 (MEC), 60353, 60354 and 60356 (PaS). The sampling frequency was 50 Hz for mice 60584 and 60585 (MEC), 60961, 92227 and 92229 (VIS). For mice 59911, 59914 (MEC) and 59912 (PaS), the wheel tracking was not synchronized to the ongoing image acquisition; see ‘Selfpaced running behaviour under sensoryminimized conditions’. The obtained speed signal from the former two groups of mice was interpolated so that the speed values matched the downsampled imaging time points (sampling frequency = 7.73 Hz), and smoothed using a square kernel of 2 s width. A threshold was applied such that all speed values that were smaller than 2 cm s^{−1} were set to zero and all speed values larger than 2 cm s^{−1} remained unchanged. We decided to threshold for immobility at a nonzero speed value (2 cm s^{−1}) in order to avoid classifying as running behaviour frames that only had minor movements of the wheel (‘twitches’), which were detected when mice slightly moved on the wheel but did not fully engage in locomotion. The threshold that we used is consistent with the one used in other studies, as in ref. ^{16}.
The speed signal obtained after applying the threshold was used to define immobility (running) bouts as the set of consecutive time points (bin size = 129 ms) for which the speed was equal to (larger than) zero (a similar approach was used in ref. ^{16}). We found that the median of velocities was 0 cm s^{−1} when all velocity values across the 10 MEC oscillatory sessions (over 3 mice) for which we had imaging data synchronized with behavioural data were pooled. This is because for some of the sessions the mice were immobile for most of the session.
When the threshold for immobility (2 cm s^{−1}, see above) was discarded (that is, set to 0 cm s^{−1}), the median was 1.3 cm s^{−1}—that is, still very low. In the absence of a threshold, our main result, which is that the oscillatory sequences traverse epochs of running and immobility, remained the same (median of probability of sequences during running = 0.85; median of probability of sequences during immobility = 0.65; two sample Wilcoxon signedrank test on the probability of sequences for running versus immobility, n = 10 oscillatory sessions over the 3 mice that had the tracking synchronized to imaging, P = 0.002, W = 55).
The acceleration was numerically calculated as the first derivative of the speed signal. Notice that in this case no interpolation was needed.
Because the available data did not have enough statistical power, it was not possible to compare the behaviour of the mice, for example in terms of its running speed and acceleration, between periods with and without ongoing oscillatory sequences.
Finally, mice that were imaged from the PaS or VIS performed the same minimalistic selfpaced running task as the mice that were imaged from the MEC recordings. The range of speed values in PaS or VIS mice across sessions = 0–58.6 cm s^{−1} (PaS) or 0–60.3 cm s^{−1} (VIS); median number of completed laps on rotating wheel in PaS or VIS mice across sessions = 145 (PaS) or 104 (VIS); maximum number of completed laps on rotating wheel in PaS or VIS mice across sessions = 502 (PaS) or 1,743 (VIS). These values are reported for MEC mice in the legend of Extended Data Fig. 2a.
Estimation of the probability of observing oscillatory sequences
To determine whether the MEC oscillatory sequences were observed during different behavioural states, the probability of observing the oscillatory sequences was calculated conditioned on whether the mouse was running or immobile. For each oscillatory session with behavioural tracking synchronized to the imaging data (10 sessions over 3 mice, see ‘Selfpaced running behaviour under sensoryminimized conditions’ and ‘Oscillation score’), all individual sequences were identified (see ‘Identification of individual sequences’). The subset of time bins that belonged to individual sequences were extracted and labelled as oscillation (bin size = 129 ms). The fraction of bins labelled as oscillation bins was 0.73 ± 0.07 (mean ± s.e.m., n = 10 sessions). Next, a second label was assigned to the time bins depending on whether they occurred during running or immobility bouts (bins labelled ‘running’ or ‘immobility’, respectively, see ‘Estimation of number of completed laps on the wheel, speed and acceleration’). The fraction of bins labelled as running = 0.43 ± 0.09, mean ± s.e.m., n = 10 sessions. After applying this procedure, each time bin had two labels, one indicating the running behaviour, and one indicating the presence (or absence) of oscillatory sequences. To estimate the probability of observing the oscillatory sequences conditioned on the mouse’s running behaviour, all bins labelled as running or immobility were identified and from each subset, the fraction of bins labelled as oscillation was calculated. These probabilities were computed for each session separately.
Sequences during immobility bouts of different lengths
The oscillatory sequences occurred both during running and immobility bouts. To quantify the extent to which individual sequences progressed during different lengths of immobility bouts, the following procedure was adopted. First, for each session, all immobility bouts were identified and assigned to bins of different lengths (see ‘Estimation of number of completed laps on the wheel, speed and acceleration’; length bins = 0–3 s, 3–5 s, 5–10 s, 10–15 s, 15–20 s, >25 s). Second, all individual sequences were identified (see ‘Identification of individual sequences’). Third, for each session and each length bin, the fraction of immobility bouts that were fully occupied by uninterrupted sequences was calculated. To estimate significance, for each session the time bins that belonged to all individual sequences were temporally shuffled. The third step of the procedure described above was performed for 500 shuffle iterations per session. In Fig. 4c, the recorded data has 10 data points per length bin, and the shuffled data has 5,000 data points per length bin, since 500 shuffled realizations per session were pooled.
Analysis of speed and sequence onset
To determine whether the onset of the MEC oscillatory sequences was modulated by the mouse’s running speed, changes in speed before and after sequence onset were investigated. For each session all individual sequences were identified (see ‘Identification of individual sequences’) and for each sequence the mean speed over windows of 10 s before and after sequence onset was calculated. Because no differences in the mean speed were observed before and after onset (Extended Data Fig. 2f left panel), we next determined whether changes in speed were correlated with the onset of sequence epochs, which were defined as epochs with uninterrupted sequences—that is, epochs with recurring sequences. The same analysis described above was repeated but only for the subset of sequences that were 10 s or more apart—that is, for sequences that belonged to different epochs.
The obtained results remained unchanged when the analysis was performed for 2 s windows before and after sequence onset.
We complemented this analysis by investigating whether new epochs of sequences were more likely to be initiated during running bouts. In each of the 10 oscillatory sessions we first identified all running and immobility bouts that were 20 s long, or longer. We then counted the number of times that a sequence onset occurred in each behavioural state. For this analysis we only considered sequences that were not preceded by other sequences (sequences that were 10 s apart or more). Results were upheld with running and immobility bouts of 40 s or longer, in which case sequence onset was 2.8 times more frequent during running.
Manifold visualization for example session in VIS and PaS
To visualize whether the topology of the manifold underlying the population activity in example sessions recorded in VIS and PaS was also a ring, PCA was used and a similar procedure to the one described in ‘Manifold visualization for MEC sessions’ was adopted.
For each example session, one corresponding to VIS and one corresponding to PaS (Fig. 5e,f), PCA was applied to the matrix of calcium activity, which first had each row convolved with a gaussian kernel of width equal to four times 8.5 s, which is the mean oscillation bin size computed across oscillatory sessions (see ‘Oscillation bin size’). Neural activity was projected onto the embedding generated by PC1 and PC2. Extended Data Fig. 11d,e shows the absence of a ringshaped manifold in VIS and PaS example sessions.
Coactivity and synchronization in PaS and VIS sessions
Sessions recorded in PaS and VIS did not exhibit oscillatory sequences. To further characterize their population activity, synchronization and neural coactivity were calculated.
Synchronization
Neural synchronization was calculated as the absolute value of the Pearson correlation between the calcium activity of pairs of cells (bin size = 129 ms). For each session, the Pearson correlation was calculated for all pairs of calcium activity (correlations with the same calcium activity were not considered) and used to build a distribution of synchronization values. In Extended Data Fig. 11j, these distributions were averaged across sessions for each brain area separately.
Coactivity
For each time bin in a session (bin size = 129 ms) the coactivity was calculated as the number of cells that had simultaneous calcium events divided by the total number of recorded cells in the session. This number represented the fraction of cells that was active in individual time bins. Using all time bins of the session, a distribution of coactivity values was calculated. In Extended Data Fig. 11k, the distributions were averaged across sessions for each brain area separately.
Model
To determine whether long sequences act as a template for the formation of given activity patterns in a neural population, we built a simple perceptron model in which 500 units were connected to an output unit (Extended Data Fig. 12a). There was a total of 500 weights in the network, one per input unit. The total simulation time was 120 s, with 3,588 simulation steps and a time step of 33.44 ms (original time step was 129 ms, to mimic the bin size used in calcium data, rescaled so that the length of one of the input sequences was 120 s, similar to the length of the sequences in Fig. 2b). The response of the output unit was given by R = WX, where W was the vector of weights, and X the matrix of input activity (each column is a time step, each row is the activity of one input unit). The weights were trained such that the output unit performed one of two target responses (see below). For each target, we trained the model using as input periodic sequences with 5 different lengths (one length per training), covering the range from very slow to very fast as compared to the characteristic time scale of the targets (100 s).
Inputs
The activity of input unit i was represented by a Gaussian: \({x}_{i}(t)={{\rm{e}}}^{\frac{{(t{\mu }_{i})}^{2}}{2{{\sigma }_{i}}^{2}}}\), \(1\le i\le 500,\,0\le t\le 240\,{\rm{s}}\), \({\sigma }_{i}=\sigma =7.6\,{\rm{s}}\), \(\forall i\). Across input units, the means of the Gaussians \({\mu }_{i}\) were temporally displaced such that, all together: (1) units fired in a sequence, and (2) the distance between the means of two consecutive cells in the sequence was the same for all pairs of consecutive cells.
This sequence was the slowest of the 5 sequence lengths we considered. Using this sequence as template, in order to build slower and periodic sequences we compressed the template and repeated it periodically by a factor of 2, 3, 4 and 8, to generate faster and periodic sequences of lengths 120, 60, 40 and 30 s respectively.
Targets
Two target responses were considered: ramp and Ornstein–Uhlenbeck process.
Ramp
The output neuron linearly increased its activity such that it was equal to 0 at time step = 0 (0 s), and to 1 at time step = 2,990 (100 s).
$${F}_{{\rm{R}}}(t)=\frac{t}{100}$$
Ornstein–Uhlenbeck process
Unlike the first target, which was deterministic, the second target was stochastic and generated by an Ornstein–Uhlenbeck process.
$$\frac{{\rm{d}}{F}_{{\rm{OU}}}}{{\rm{d}}t}=\frac{{\mu }_{{\rm{OU}}}{F}_{{\rm{OU}}}(t)}{\tau }+{\sigma }_{{\rm{OU}}}\xi (t)$$
where μ_{OU} = 1 denotes the longterm mean, ξ is a white noise of zero mean and variance σ_{OU} = 0.005, and τ = 25.6 s denotes the correlation time.
Training of weights
The weights between the inputs and the output unit were trained such that the output unit performed one of the two target responses explained above. At the end each of the 1,000 learning iterations, the weights were updated through the perceptron learning rule \(\varDelta {w}_{i}=\eta e{x}_{i}\), where x_{i} was the input from neuron i, \(1\le i\le 500\), and η = 1 was the learning rate. In each learning iteration, the error e was calculated as the sum over time steps t of the difference between the target response and the output response—that is, \(e=\sum _{t}T(t)WX(t),\) where T(t) is the target response (either the ramp or the Ornstein–Uhlenbeck process) at time point t, and X(t) is the vector of input activity at time point t. The mean total error plotted in Extended Data Fig. 12d was calculated as the mean error over the last 100 learning iterations.
Data analysis and statistical analysis
Data analyses were performed with customwritten scripts in Python and Matlab (R2021b). Results were expressed as the mean ± s.e.m. unless indicated otherwise. Statistical analysis was performed using MATLAB and P values are indicated in the figure legends and figures (NS: P > 0.05; *P < 0.05, **P < 0.01, ***P < 0.001). For data that displayed no Gaussian distribution and that was unpaired, the Wilcoxon ranksum test was used. For paired data or onesampled data, the Wilcoxon signedrank test was used. Twotailed tests were used unless otherwise indicated. Correlations were determined using Pearson or Spearman correlations. Friedman tests were used for analyses between groups. The Bonferroni correction was used when multiple comparisons were performed.
Power analysis was not used to determine sample sizes. The study did not involve any experimental subject groups; therefore, random allocation and experimenter blinding did not apply and were not performed.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Dr. Thomas Hughes is a UKbased scientist and science communicator who makes complex topics accessible to readers. His articles explore breakthroughs in various scientific disciplines, from space exploration to cuttingedge research.