All experiments had been carried out in accordance with the Norwegian Animal Welfare Act and the European Conference for the Safety of Vertebrate Animals used for Experimental and Different Scientific Functions, Allow numbers 18011 and 29893.
Topics
Male C57/Bl6 mice had been housed in social teams of two–6 people per cage (calcium imaging experiments) or individually (electrophysiology experiments, after implantation). The mice had entry to nesting materials and a planar working wheel and had been stored on a 12 h gentle/12 h darkness schedule in a temperature and humiditycontrolled vivarium. Meals and water had been supplied advert libitum. Twophoton calcium imaging knowledge had been collected from a cohort of 12 mice (5 implanted in MEC, 4 in PaS, and 3 in VIS). Electrophysiological knowledge from the MEC had been collected from 2 mice.
Surgical procedures
For all surgical procedures, anaesthesia was induced by putting the themes in a plexiglass chamber crammed with isoflurane vapour (5% isoflurane in medical air, circulation of 1 l min^{−1}). Surgical procedure was carried out on a heated surgical procedure desk (38 °C). Air circulation was stored at 1 l min^{−1} with 1–3% isoflurane as decided from physiological monitoring of respiration and heartbeat. The mice had been allowed to get better from surgical procedure in a heated chamber (33 °C) till they regained full mobility and application. Postoperative analgesia was given within the type of subcutaneous injections of Metacam (5 mg kg^{−1}) 24 and 48 h after the primary Metacam injection so long as was deemed needed. Moreover, the mice got 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 primary 36 h after the primary Temgesic injection.
Surgical procedures for calcium imaging
Surgical procedures had been carried out in line with a twostep protocol. Through the first process, new child pups or grownup mice had been injected in MEC or PaS, or grownup mice had been injected in VIS with a virus carrying a assemble 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, College of Pennsylvania, USA) was diluted 1:1 in sterile DPBS (1× Dulbecco’s Phosphate Buffered Saline, Gibco, ThermoFisher). Through the second process, two weeks later, a microprism was implanted to realize optical entry to contaminated neurons situated in MEC and PaS, or a glass window was inserted to acquire comparable entry in VIS.
Virus injection and microprism implantation in MEC and PaS
Within the first surgical process, new child pups obtained injections of AAV1SynGCaMP6m in the future after delivery^{51}. An analgesic was supplied instantly earlier than the surgical procedure (Rymadil, Pfizer, 5 mg kg^{−1}). Preheated ultrasound gel (39 °C, Aquasonic 100, Parker) was generously utilized on the pup’s head with a purpose to create a big medium for the transmission of ultrasound waves. Actualtime ultrasound imaging (Vevo 1100 System, Fujifilm Visualsonics) allowed for focused supply of the viral combination to particular areas of the mind. Throughout ultrasound imaging, the pup was immobilized by means of a custommade mouth adapter. The ultrasound probe (MS550S) was lowered to be in shut contact with the gel and thus the pup’s head to permit visualization of the focused constructions. The probe was stored in place for the entire period of the process through the VEVO injection mount (VEVO Imaging Station. Imaging in BMode, frequency: 40 MHz; energy: 100%; acquire: 29 dB; dynamic vary: 60 dB). Goal areas had been recognized by structural landmarks: the MEC or PaS had been recognized within the antero–posterior and medio–lateral axis by the looks of the aqueduct of Sylvius and the lateral sinus. The goal space for injection was similar to a coronal part at ∼−4.7 mm from bregma within the grownup mouse. The answer containing the virus (250 ± 50 nl per injection) was injected within the goal areas through beveled glass micropipettes (Origio, {custom} made; outer tip opening: 200 μm; interior tip opening: 50 μm) utilizing a pressurepulse system (Visualsonics, 5 pulses, 50 nl per pulse). The pipette tip was pushed by means of the mind with none incision on the pores and skin, or a craniotomy, and, to cut back the period of the process, retracted instantly after depositing the virus within the goal space. The anatomical specificity of the an infection was verified by imaging serial sections of the contaminated hemispheres after experiment completion (see ‘Histology of calcium imaging mice and reconstruction of fieldofview location’).
Two weeks after the viral injection, we carried out a second process, through which a microprism was implanted within the left hemisphere to realize optical entry to the superficial layers of MEC and PaS^{52}. The implanted microprism was a rightangle prism with 2 mm aspect size and reflective enhanced aluminium coating on the hypotenuse (Tower Optical). The prism was glued to a 4mmdiameter (CS4R, thickness no. 1) spherical coverslip with UVcurable adhesive (Norland). On the day of surgical procedure, mice had been anaesthetized with isoflurane (IsoFlo, Zoetis, 5% isoflurane vapourised in medical air delivered at 0.8–1 l min^{−1}) after which two analgesics had been supplied by means of 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 native analgesic was utilized beneath the pores and skin overlaying the cranium (Marcain, Aspen, 1–3 mg kg^{−1}). Their scalp was eliminated with surgical scissors and the floor of the bone was dried earlier than being generously coated with optibond (Kerr). To extend the thickness and stability of the cranium and total preparation, a skinny layer of dental cement (Charisma, Kulzer) was utilized on the uncovered cranium, besides within the location above the implant, the place a 4mmwide round craniotomy was made. The craniotomy was positioned over the dorsal floor of the cortex and cerebellum, with the centre positioned ∼ 4 mm lateral from the centre of the medial sinus, and above the transverse sinus simply above the MEC and PaS. After the dura was eliminated above the cerebellum, the decrease fringe of the prism was slowly pushed within the empty area between the forebrain and the cerebellum, simply posterior to the transverse sinus. The sides of the coverslip had been secured to the encircling cranium with UVcurable dental cement (Venus Diamond Move, Kulzer). A customdesigned metal headbar was connected to the dorsal floor of the cranium, centred upon and positioned parallel to the highest face of the microprism. All uncovered areas of the cranium, together with the headbar, had been lastly coated with dental cement (Paladur, Kulzer) and made opaque by including carbon powder (Sigma Aldrich) till the dental cement powder turned darkish gray.
Virus injection and glass window implantation in VIS
In a special cohort of mice than these used for MEC/PaS imaging, we induced the expression of GCaMP6m in neurons of the grownup VIS for subsequent imaging. We focused the injection of the identical AAV1SynGCaMP6m viral answer used within the creating MEC and PaS to the first visible cortex. On the day of surgical procedure, 3 to 5monthold mice had been anaesthetized with isoflurane (IsoFlo, Zoetis, 5 % isoflurane vapourized in medical air delivered at 0.8–1 l min^{−1}) after which two analgesics had been supplied by means of 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 native anaesthetic was utilized beneath the pores and skin overlaying the cranium (Marcain, Aspen, 1–3 mg kg^{−1}). The virus was injected at three areas in VIS, all of which had been inside the following anatomical ranges in the suitable hemisphere: 2.3–2.5 mm lateral from the midline, 0.9–1.3 mm anterior from lambda^{53}. At every injection web site, 50 nl of the virus was injected 0.5 mm beneath the dura and the pipette was left in place for 3–4 min to allow the virus to diffuse. The pipette was then delivered to 0.3 mm beneath the dura and one other 50 nl was injected. The pipette was then left in place for five–10 min earlier than retracting it utterly. The pace of the injections was 5 nl s^{−1}.
Two weeks after the viral injection, a surgical procedure to chronically implant a glass window over VIS was carried out. The mice had been dealt with as beforehand described for the prism surgical procedure in MEC/PaS, together with anaesthesia, supply of analgesics, and scalp removing. Optibond was utilized to the uncovered cranium besides within the location of the craniotomy. A 4mmwide craniotomy was made, centred on the virus injection coordinates, and a 4mm glass window was positioned beneath the cranium edges of the craniotomy. The glass was barely bigger than the craniotomy, so after it was manoeuvred in place, the upward strain exerted by the mind secured it in place in opposition to the cranium, thereby minimizing the presence of empty gaps that may favour tissue and bone regrowth. The sides of the window had been secured with UVcurable dental cement and superglue earlier than the positioning of the headbar as described for the MEC–PaS implantation. All uncovered areas of the cranium, together with the headbar, had been lastly coated with dental cement (Paladur, Kulzer) that was made opaque by including carbon powder (Sigma Aldrich) till the dental cement powder turned darkish gray.
Neuropixels probe implants
Two grownup mice (4 to five months previous) had been implanted with fourshank Neuropixels 2.0 silicon probes^{54} focusing on the superficial layers of MEC within the left hemisphere. Previous to the surgical procedure, the mice got basic analgesics (Metacam, Boehringer Ingelheim, 5 mg kg^{−1} and Temgesic, Indivior, 0.05–0.1 mg kg^{−1}) subcutaneously and one native anaesthetic was utilized beneath the pores and skin overlaying the cranium (Marcain, Aspen, 1–3 mg kg^{−1}). After incision, a gap was drilled over the cerebellum for an anchor screw linked to a floor wire. Craniotomies had been then drilled. Probes focusing on the MEC had been lowered from the floor to depths between 2.5 mm and a couple of.7 mm relative to the dura mater. They had been implanted with essentially the most medial shank positioned on the mind floor 3.2 mm lateral to the midline and 0.4 mm anterior to the transverse sinus edge. The 4 shanks had been oriented with the electrode websites on the posterior aspect. In one of many two mice (no. 104638), the probe was first rotated 7° within the horizontal airplane (angle on the subject of the coronal airplane), with essentially the most lateral shank in essentially the most posterior place such that the shanks had been parallel to the transverse sinus. The 4 shanks had been then lowered vertically from this place.
The Neuropixels probe of the second mouse (no. 102335) was not rotated within the horizontal airplane—that’s, all shanks had the identical anterior–posterior coordinates. The electrode shanks of this mouse had been lowered from the floor with a 2° angle relative to the coronal airplane, such that the shank ideas had been essentially the most posterior. The shanks remained inside the similar sagittal airplane as they had been lowered. This second mouse was additionally implanted with a probe focusing on the CA1 area in the suitable hemisphere, 1.225–1.975 mm relative to the midline, at a depth of three mm relative to dura mater, with all shanks 2.1 mm posterior to bregma. The hippocampal knowledge weren’t used within the current research. The probes had been secured to the cranium utilizing an adhesive (OptiBond, Kerr), UVcurable dental cement (Venus Diamond Move, Kulzer), and dental cement (Meliodent, Kulzer). A headbar was connected as described above for the calcium imaging research.
Selfpaced working behaviour below sensoryminimized circumstances
Coaching of mice started 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) had been headrestrained by a headbar with their limbs resting on a freely rotating styrofoam wheel with a steel shaft fastened by means of the centre. The radius of the wheel was ∼85 mm and the width 70 mm. Low friction ball bearings (HK 0608, Kulelager) had been affixed to the ends of the steel shaft and held in place on the optical desk utilizing a {custom} mount. This association allowed the mice to selfregulate their motion. The place of the mouse on the rotating wheel was measured utilizing a rotary encoder (E6B2CWZ3E, YUMO) connected to its centre axis. Step values of the encoder (4,096 per full revolution, ∼130 µm decision) had been digitized by a microcontroller (Teensy 3.5, PJRC) and recorded utilizing {custom} Python scripts at 40–50 Hz. Wheel monitoring was triggered at the beginning of imaging and synchronized to the continued picture acquisition by means of a digital enter 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 exact synchronization was not out there to us and these knowledge had been therefore not used for comparability of motion and imaging knowledge. A Tslot picture interrupter (EESX672, Omron) served as a lap (full revolution) counter. Design and code of the wheel are publicly out there below https://github.com/kavlintnu/wheel_tracker.
The opposite Neuropixels probeimplanted mouse (no. 102335) was headrestrained by a headbar whereas resting on a round disc coated with rubber spray. The radius of this wheel was ∼85 mm. The mouse was allowed selfpaced motion on the wheel. Threedimensional movement seize (OptiTrack Flex 6 cameras and Motive recording software program) was used to trace the rotation of the wheel by monitoring retroreflective markers positioned on the wheel edge. Digital pulses had been generated utilizing an Arduino microcontroller which had been used to align the Neuropixels acquisition system and the OptiTrack system through direct TTL enter and infrared LEDs.
In all mice, the selfpaced job was carried out below circumstances of minimal sensory stimulation, in darkness, and with no rewards to sign elapsed time or distance run^{16,17}. Previous to the imaging classes, the calcium imaging mice had been accustomed to the setup by means of day by day exposures over the course of between 5 and 15 classes, one session per day. Neuropixelsimplanted mice had been habituated to the setup by step by step rising the time spent on the wheel over 4 days. In every session, after the mice had been positioned on the wheel, they had been gently headrestrained and free to run or relaxation^{55,56} for 30, 45 or 60 min.
Recording classes of Neuropixelsimplanted mice additionally consisted of trials the place the mice had been freely foraging in a 80 cm × 80 cm open area area for 30 min. These open area trials preceded the selfpaced wheel trials and weren’t used within the current research.
Twophoton imaging in headfixed mice
A custombuilt 2photon benchtop microscope (Femtonics, Hungary) was used for 2photon imaging of the goal areas (that’s, 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 because the excitation supply. Common laser energy on the pattern (after the target) was 50–120 mW. Emitted GCaMP6m fluorescence was routed to a GaAsP detector by means of a 600 nm dichroic beamsplitter plate and 490–550 nm bandpass filter. Gentle was transmitted by means of a 16×/0.8 NA waterimmersion goal (MRP07220, Nikon) fastidiously lowered in shut contact to the coverslip glued to the microprism (for MEC–PaS imaging) or above the coverslip involved with the mind floor (for VIS imaging). For the microprismimplanted mice, the target lens was aligned to the ventro–lateral nook of the prism, to constantly establish the place of MEC and PaS throughout mice. Ultrasound gel (Aquasonic 100, Parker) or water was used to fill the hole between the target lens and the glass coverslips. The software program MESc (v 3.3 and three.5, Femtonics, Hungary) was used for microscope management and knowledge acquisition. Imaging time collection of both ∼30 min or ∼60 min had been acquired at 512 × 512 pixels (sampling frequency: 30.95 Hz, body period: ∼32 ms; pixel measurement: both 1.78 × 1.78 µm^{2} or 1.18 × 1.18 µm^{2}). Time collection acquisition was initiated arbitrarily after the mouse was headrestrained on the setup.
Neuropixels recordings in headfixed mice
Indicators had been recorded utilizing a Neuropixels acquisition system as described beforehand^{25,57}. Briefly, the electrophysiological sign was amplified with a acquire of 80, lowpassfiltered at 0.5 Hz, highpassfiltered at 10 kHz, and digitized at 30 kHz on the probe circuit board. The digitized sign was then multiplexed by the ‘headstage’ circuit board and transmitted alongside a 5 m tether cable utilizing twisted pair wiring to a Neuropixels PXIe acquisition module. The info was visualized and recorded utilizing SpikeGLX model 20201103 software program (https://billkarsh.github.io/SpikeGLX).
Histology
Histology of calcium imaging mice and reconstruction of fieldofview location
On the final day of imaging, after the imaging session, the mice had been anaesthetized with isoflurane (IsoFlo, Zoetis) after which obtained an overdose of sodium pentobarbital earlier than transcardial perfusion with freshly ready PFA (4% in PBS). After perfusion, the mind was extracted from the cranium and stored in 4% PFA in a single day for postfixation. The PFA was then exchanged with 30% sucrose to cryoprotect the tissue.
To confirm the anatomical location of the imaged FOVs within the microprismimplanted mice, we used small, custommade pins, derived from a skinny piano wire coated with an answer of 1,1′dioctadecyl3,3,3′,3′tetramethylindocarbocyanine perchlorate (DiI; DiIC18(3)) (ThermoFischer), to mark the placement of the imaged tissue in relation to the prism footprint. A DiIcoated pin was inserted into the mind tissue on the location left empty by the prism footprint, and particularly focused to the ventro–lateral nook of the footprint (see ‘Surgical procedures’). The pin was left in place to favour switch of DiI from the steel pin to the mind tissue, and to depart a fluorescent mark on the placement of the imaged FOV. After 30 to 60 s, the pin was eliminated and the mind was sliced on a cryostat in 30–50 µm thick sagittal sections. All slices had been collected sequentially in a 24well plate crammed with PBS, earlier than being mounted of their applicable anatomical order on a glass slide in custommade mounting medium. For confocal imaging, a Zeiss LSM 880 microscope (Carl Zeiss) was used to scan by means of the entire collection of slices and find the place of the DiI fluorescent mark. Pictures had been then acquired utilizing an EC PlanNeofluar 20×/0.8 NA air immersion, 40×/1.3 oil immersion, or 63×/1.4 oil immersion goal (Zeiss, laser energy: 2–15%; optical slice: 1.28–1.35 ethereal items, step measurement: 2 µm). Earlier than acquisition, acquire and digital offset had been established to optimize the dynamic vary of acquisition to the dynamic vary of the GCaMP6m and DiI alerts. Settings had been stored fixed throughout acquisition throughout brains. Based mostly on the placement of the purple fluorescent mark, we may infer the place, on the medio–lateral and dorso–ventral extent of the mind, the ventro–lateral nook of the microprism (and therefore the 2photon FOV aligned to it) was situated.
We used the Paxinos mouse mind atlas^{53} to provide 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 inside the anatomical boundaries of both the MEC and adjoining PaS, and allowed for a standardized comparability throughout mice. For every imaged mouse, we mapped the dorso–ventral and medio–lateral location of the DiI mark on the refence flat map (Prolonged Information Fig. 1c). Mice had been assigned to ‘MEC imaging’ or ‘PaS imaging’ teams relying on the placement of the FOV: a mouse could be additional analysed as being a part of the MEC imaging group if greater than 50% of the world of the FOV occupied by GCaMP6mexpressing cells may very well be situated within the MEC.
To confirm the anatomical location of the FOVs in VIS within the glass window implanted mice, we sliced the mind till we reached the anatomical coordinates at which the virus was infused (see ‘Surgical procedures’). Coronally reduce slices of fifty µm thickness had been collected sequentially in a 24well plate, and instantly mounted of their applicable anatomical order on a glass slide in custommade mounting medium. For confocal imaging, a Zeiss LSM 880 microscope (Carl Zeiss) was used in line with the identical specification as described above for MEC/PaS.
Histology and reconstruction of Neuropixels probe placement
After the tip of experiments, the mice had been anaesthetized and obtained an overdose of isoflurane (IsoFlo, Zoetis) earlier than transcardial perfusion with saline adopted by 4% formaldehyde. The mind was both extracted after perfusion or stored in a single day in 4% formaldehyde for postfixation earlier than extraction. The brains had been then saved in 4% formaldehyde. Frozen 30 µm thick sagittal sections had been reduce on a cryostat, mounted on glass, and stained with Cresyl violet (Nissl). To estimate the shank areas, we used an Axio Scan.Z1 (Carl Zeiss) slide scanner microscope for brightfield detection at 20x magnification. We used Paxinos mouse mind atlas^{53} and the Allen Mouse Mind Frequent Coordinate Framework^{58} model 3 by means of the siibraexplorer (Forschungszentrum Juelich, https://atlases.ebrains.eu/viewer/) to estimate anatomical location of recording websites. A map of the probe shank was aligned to the histology assuming that the chopping airplane was nearparallel to the sagittal airplane. When doable, the anatomical areas had been calculated utilizing the tip of the probe shanks and the intersection of the shank with the mind floor as reference frames. When this was not doable, the profile of a closeby mind area (for instance, the hippocampus) was used to estimate the MEC implant web site. We noticed thetarhythmicity of neural exercise on all recorded shanks, as anticipated for recording areas within the MEC.
Evaluation of imaging time collection
Imaging time collection knowledge had been analysed utilizing the Suite2p^{59} Python library (https://github.com/MouseLand/suite2p). We used its builtin routines for movement correction, area of pursuits (ROI) extraction, neuropil sign estimation, and spike deconvolution. Nonrigid movement correction was chosen to align every body iteratively to a template. High quality was assessed by visible inspection of the corrected stacks and builtin movement correction metrics. The Suite2p GUI was used to manually subselect putative neurons based mostly on anatomical and sign traits and to discard apparent artefacts that collected in the course of the evaluation—for instance, ROIs with footprints spanning giant areas of the FOV, ROIs that didn’t have clearly delineated circumferences within the generated most depth projection, or ROIs that had been extracted robotically however confirmed no seen calcium transients.
Uncooked fluorescence calcium traces of every ROI had been neuropilcorrected to create a fluorescence calcium sign F_{corr} by subtracting 0.7 occasions the neuropil sign from the uncooked fluorescence traces. We used the Suite2p builtin model of nonnegative deconvolution^{60} with tau = 1 s to deconvolve F_{corr}, yielding the idea for the binarized sequences that we consult with because the calcium exercise (see ‘Binary deconvolved calcium exercise and matrix of calcium exercise’). As a result of absence of floor fact knowledge for our mixture of indicator, area, and imaging circumstances, we used a decay tau that was on the decrease finish of biologically believable values (tau = 1 s), which allowed even quick and low amplitude spiking responses to be picked up by the evaluation and due to this fact didn’t bias our evaluation in direction of largeamplitude calcium transients (presumed bursting responses). To estimate the signaltonoise ratio (SNR) of every cell individually, we additional thresholded the calcium exercise (with out binarization) at 1 s.d. over the imply, yielding filtered calcium exercise, and labeled the remaining exercise as noise. We moreover ensured that noise was temporally properly segregated from filtered calcium exercise by requiring knowledge factors labeled as noise to be separated by no less than one second earlier than and ten seconds after filtered calcium exercise. The SNR of the cell was then estimated because the ratio of the imply amplitude of F_{corr} throughout episodes of filtered calcium exercise over the s.d. of F_{corr} throughout episodes of noise. If no knowledge factors remained after the filtering of calcium exercise, the cell was assigned a SNR of zero.
Binary deconvolved calcium exercise and matrix of calcium exercise
As a way to denoise the recorded fluorescence calcium alerts and have good temporal decision, all analyses within the research had been carried out utilizing the deconvolved calcium exercise of the recorded cells. For every cell whose SNR was bigger than 4, the deconvolved calcium exercise (see ‘Evaluation of imaging time collection’) was downsampled by an element of 4 by calculating the imply over time home windows of ∼129 ms (unique sampling frequency = 30.95 Hz, sampling frequency used within the analyses = 7.73 Hz). As a result of the ultraslow oscillations and periodic sequences unfolded on the time scales of seconds to minutes, this downsampling step gave temporal decision for all quantifications whereas permitting us to work with smaller arrays (ultraslow oscillations and the oscillatory sequences had been additionally detectable when utilizing the unique sampling frequency), which in a few of the analyses decreased the computing time. Subsequent, the downsampled deconvolved calcium exercise was averaged over time and its s.d. was calculated. A threshold equal to this common plus 1.5 occasions the s.d. was used to transform the deconvolved calcium exercise right into a binary deconvolved calcium exercise, such that every one values above the edge had been set to 1 (calcium occasions), and all values beneath or equal to that threshold had been set to 0. Except acknowledged in any other case, for all analyses all through the research we used the deconvolved and binary calcium exercise, to which for simplicity we consult with as ‘deconvolved calcium exercise’ or just ‘calcium exercise’. The calcium exercise of all cells in a session with SNR > 4 was stacked to assemble a binary matrix of calcium exercise which had as many rows as neurons, and as many columns as time bins sampled at 7.73 Hz. The inhabitants vectors are the columns of the matrix of calcium exercise.
Word that the recorded calcium alerts probably mirror a mix of teams of single spikes and higherfrequency bursts, though it was not doable to differentiate between the 2 varieties of firing. The sensitivity of the calcium indicator was probably not excessive sufficient to detect subthreshold potentials.
Spike Sorting and singleunit choice
Spike sorting of Neuropixels knowledge was carried out utilizing a model of KiloSort 2.5 (ref. ^{54}) with some customizations to enhance efficiency on recordings from the MEC area as described beforehand^{25}. All trials in a session had been clustered collectively. Single items had been discarded from evaluation based mostly on a < 20% estimated contamination fee with spikes from different neurons. These items had been robotically labelled by the KiloSort 2.5 algorithm as ‘good’ items. Within the instance session from mouse no. 104638 solely good items had been thoughtabout. Within the instance session of mouse no. 102335, as a result of the variety of good items was decrease (<250), we additionally used multiunit exercise (MUA).
Autocorrelations and spectral evaluation of singlecell calcium exercise
To find out if the calcium exercise of single cells shows ultraslow oscillations, for every neuron the PSD was calculated on the autocorrelation of its calcium exercise. The PSD was computed utilizing Welch’s technique (pwelch, builtin Matlab perform), with Hamming home windows of 17.6 min (8,192 bins of 129 ms in every window) and 50% of overlap between consecutive home windows. Word that when calculating the PSD a big window was wanted to establish oscillation frequencies ≪0.1 Hz.
To visualise whether or not particular oscillatory patterns at fastened frequencies had been current within the neural inhabitants, all autocorrelations from one session had been sorted and stacked right into a matrix, the place rows are cells and columns are time lags. The sorting of autocorrelations was carried out in line with the utmost energy of every PSD in a descending method. The frequency at which the PSD peaked was used as an estimate of the oscillatory frequency of the cell’s calcium exercise.
As a way to decide significance for the height of the PSD, we thoughtabout two excessive and reverse shuffling procedures: On the one hand, provided that circularly shuffling the information preserves all inter calcium occasions (Prolonged Information Fig. 3c,d), taking this strategy would protect the form and the place of the height within the PSD calculated on experimental knowledge. Then again, destroying the inter calcium occasion intervals by assigning a random place to every calcium occasion within the time collection would result in a flat PSD (Prolonged Information Fig. 3c,d). Within the latter strategy, all cells could be labeled as oscillatory. To bridge these two approaches we developed a brand new shuffling process. For every cell we divided its calcium exercise vector into n epochs of size W, with (n=T/(W,bullet ,{rm{SF}})), the place T is the overall variety of time bins sampled at a frequency SF = 7.73 Hz (that’s, bin measurement = 129 ms). We subsequent shuffled these epochs (and preserved the ordering of the time bins inside every epoch). This technique preserved the inter calcium occasion interval, however on the similar time disrupted the periodicity. Within the restrict the place W = 129 ms, this technique coincides with shuffling all calcium occasions with out preserving the inter calcium occasion intervals; within the restrict the place W = T/SF, this technique is equal to circularly shuffling the information. For every of the 200 shuffled realizations we calculated the PSD and the fraction of cells for which the height of the PSD in experimental knowledge was above the ninety fifth percentile of a shuffled distribution constructed with the values of the PSDs calculated on shuffled knowledge (and on the frequency at which the PSD computed on experimental knowledge peaked). Right here we current the outcomes for five completely 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 beneath the standard period of the sequences (W < 50 s), the nice majority of cells are labeled as having a peak within the PSD. As anticipated, when W is just like the period of the sequences (W ≥ 50 s), the fraction of oscillatory cells rapidly drops. This fraction is now not considerably above an opportunity degree of 5%.
This strategy was used for figuring out the fraction of oscillatory cells each in calcium imaging and in Neuropixels knowledge. In the primary textual content we current the outcomes similar to W = 20 s.
Lastly, we be aware that there was some variability within the frequency at which the PSD peaked throughout cells inside a session. For instance, within the instance session proven in Fig. 1b–d and Fig. 2a, some singlecell PSDs peaked at a frequency of 0.0066 Hz, whereas others did so at a frequency of 0.0075 Hz. Nonetheless, in lots of instances the PSDs had been broad sufficient to exhibit excessive energy in neighbouring frequencies too, offering assist to the frequencies being quite clustered amongst a subset of values, with some slight variability round these values. When all cells had been analysed (n = 6,231 cells pooled throughout 15 oscillatory classes, 5 mice), in roughly half of the MEC knowledge the oscillatory frequency on the singlecell degree was similar to the frequency on the inhabitants degree (Prolonged Information Fig. 7e). This discovering factors to a small variability within the frequency of singlecell exercise in MEC, as anticipated within the presence of recurring sequences.
Correlation and PCA sorting strategies
To find out whether or not neural inhabitants exercise displays temporal construction we visualized the inhabitants exercise by the use of raster plots through which we sorted all cells in line with completely different strategies.
Correlation technique
This technique kinds cells such that these which might be close by within the sorting are extra synchronized than these which might be additional away. First, every calcium exercise was downsampled by an element 4 by calculating the imply over counts of calcium occasions in bins of 0.52 s. The obtained calcium exercise was then smoothed by convolving it with a gaussian kernel of width equal to 4 occasions the oscillation bin measurement, a bin measurement that was consultant of the temporal scale of the inhabitants dynamics (see ‘Oscillation bin measurement’). The cross correlations between all pairs of cells had been calculated utilizing time bins as knowledge factors, and a most time lag of 10 time factors, equal to ∼ 5 s. This small time lag allowed us to establish close to instantaneous correlation whereas retaining details about the temporal order of exercise between cell pairs. The utmost worth of the crosscorrelation between cell i and cell j was saved within the entry (i,j) of the correlation matrix C, which was a sq. matrix of N rows and N columns, the place N was the overall variety of recorded neurons within the session with SNR > 4. If the crosscorrelation peaked at a destructive time lag the worth within the entry (i,j) was multiplied by −1. The entry with the very best crosscorrelation worth was recognized and its row, denoted by i_{max,} was used because the ‘seed’ cell for the sorting process and chosen to be the primary cell within the sorting. Cells had been then sorted in line with the values within the entries ({(i}_{max },j)), (j=mathrm{1,2},ldots ,N), j ≠ i_{max}, that’s, their correlations with the seed cell, in a descending method.
PCA technique
Computing correlations from the calcium exercise or the calcium alerts may be noisy attributable to wonderful tuning of hyperparameters (for instance, the dimensions of the kernel used to easy the calcium exercise of all cells). To keep away from this, we leveraged the truth that the periodic sequences of neural exercise represent lowdimensional dynamics with intrinsic dimensionality equal to 1, and sorted the cells based mostly on an unsupervised dimensionality discount^{61} strategy (an analogous strategy was utilized in ref. ^{62}). For every recording session, PCA was utilized to the matrix of calcium exercise (bin measurement = 129 ms; utilizing Matlab’s builtin pca perform), together with all epochs of motion and immobility and utilizing the rows (neurons) as variables and the columns (time bins) as observations. The primary two principal parts (PCs) had been stored, since 2 is the minimal variety of parts wanted to embed nonlinear 1dimensional dynamics. Cells had been sorted in line with their loadings in PC1 and PC2, anticipating that the connection between these loadings would specific the ordering in cell activation in the course of the sequences.
The airplane spanned by PC1 and PC2 was named the PC1–PC2 airplane. Within the PC1–PC2 airplane, the loadings of every neuron (the parts of the eigenvectors with out being multiplied by the eigenvalues) outlined a vector, for which we computed its angle ({theta }_{i}={rm{arctg}},left(frac{{l}_{{rm{PC}}2}^{i}}{{l}_{{rm{PC}}1}^{i}}proper)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]), the place N is the overall variety of recorded cells with SNR > 4, was assigned to every neuron with out repetition throughout cells. Neurons had been sorted in line with these assigned numbers (see instance session in Prolonged Information Fig. 4d, prime row).
Sorting of circularly shuffled knowledge
A shuffled matrix of calcium exercise was constructed by circularly shuffling the calcium exercise of every cell individually. For every cell a random ordinal integer (in [1,T]), the place T is the overall variety of time bins (bin measurement = 129 ms), was chosen and the calcium exercise was rigidly shifted by this integer utilizing periodic boundary circumstances. The task of random ordinal integers was made individually for every cell. The PCA technique was then utilized to the shuffled matrix of calcium exercise (see instance session in Prolonged Information Fig. 4d, second row).
Sorting of temporally shuffled knowledge
As a result of circularly shuffling the information preserves the oscillations within the singlecell calcium exercise, a second shuffling strategy was thoughtabout (for singlecell knowledge shuffling procedures see ‘Autocorrelations and spectral evaluation of singlecell calcium exercise’). A shuffled matrix of calcium exercise was constructed by temporally shuffling the calcium exercise of every cell individually. For every cell, every time bin of the calcium exercise was assigned a random ordinal integer (in [1,T]) with out repetition throughout time bins, the place T is the overall variety of time bins (bin measurement = 129 ms), and time bins had been ordered in line with their assigned quantity. The task of random ordinal integers was made individually for every cell, in order that the obtained random orderings weren’t shared throughout cells. The PCA technique was then utilized to the shuffled matrix of calcium exercise.
Sortings are preserved when completely different parts of information are used for acquiring the sortings
To find out whether or not utilizing completely different parts of the session for sorting the neurons result in completely different sortings, the PCA technique was utilized to: (i) all knowledge inside a session; (ii) the primary half of the session; and (iii) the second half of the session. This process gave three sortings per session. Subsequent, for every cell pair in a session the space between the 2 cells in every of the three sortings was calculated. We illustrate this calculation with a toy instance: if 5 neurons had been recorded, and sorting (i) was: (1,4,5,2,3), the space between cells 1 and 5 was 2, as a result of these two cells had been 2 positions aside within the sorting. The gap between cells 1 and three was 1 and never 4, nonetheless, as a result of within the calculation of distances we took under consideration that the sorting mirrors the place of the cells within the ring, which has periodic boundary circumstances.
We subsequent 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 completely different parts of information protect the ordering of the neurons, we might anticipate excessive correlation values. We in contrast the obtained correlation values with the ninety fifth percentile of a shuffled distribution obtained by assigning, to every cell, a random place in every of the sortings.

Sorting (i) versus sorting (ii): 15 of 15 oscillatory classes (see ‘Oscillation rating’) had been above the cutoff of significance. Correlation values in experimental knowledge ranged from 0.38 to 0.85. The ninety fifth percentile of shuffled knowledge ranged from 0.004 to 0.015 (n = 15 in each experimental and shuffled knowledge).

Sorting (i) versus sorting (iii): 15 of 15 oscillatory classes had been above the cutoff of significance. Correlation values in experimental knowledge ranged from 0.52 to 0.86. The ninety fifth percentile of shuffled knowledge ranged from 0.005 to 0.013 (n = 15 in each experimental and shuffled knowledge).

Sorting (ii) versus sorting (iii): 15 of 15 oscillatory classes had been above the cutoff of significance. Correlation values in experimental knowledge ranged from 0.17 to 0.53. The ninety fifth percentile of shuffled knowledge ranged from 0.005 to 0.013 (n = 15 in each experimental and shuffled knowledge).
The excessive correlation values obtained present assist for what’s illustrated in Prolonged Information Fig. 4e: utilizing completely different parts of information for sorting the cells unveils the identical dynamics.
Sorting strategies based mostly on nonlinear dimensionality discount methods
The PCA technique for sorting cells depends on a twodimensional linear embedding. This linear embedding won’t be optimum if the inhabitants vectors describe temporal trajectories that, regardless of being lowdimensional, lie on a curved floor. To consider potential nonlinearities, 4 further sorting strategies had been carried out, based mostly on the next nonlinear dimensionality discount methods^{63}: tdistributed stochastic neighbour embedding (tSNE), LEM, Isomap and uniform manifold approximation and projection (UMAP)^{64} (see parameters beneath). First, to specific within the sortings the ordering of the cells in the course of the gradual temporal development of the sequences, the 4 strategies used a resampled matrix of calcium exercise as enter. To compute this matrix, for every session, we downsampled every calcium exercise by an element 4 by calculating its imply in bins of 0.52 s. The calcium exercise of all cells was then smoothed by convolving them with a gaussian kernel whose width was given by the oscillation bin measurement (see ‘Oscillation bin measurement’). After making use of tSNE, LEM, Isomap or UMAP to the resampled matrix of calcium exercise, we stored the primary two dimensions obtained with every technique, for a similar causes as introduced for the PCA sorting technique. To acquire the sorting, the next process was utilized: We let Dim1 and Dim2 be the primary two dimensions obtained with the chosen dimensionality discount approach that we had utilized to the resampled matrix. In analogy with the PCA technique, the Dim1–Dim2 airplane was spanned by Dim1 and Dim2 and for every cell the parts on these dimensions outlined a vector on this airplane 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(tright)) 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(tright),{rm{PC}}{2}_{t}left(tright))), 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(tright)={rm{arctg}},left(frac{{rm{PC}}{2}_{t}left(tright)}{{rm{PC}}{1}_{t}left(tright)}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(tright)) is always defined, regardless of whether the session is or is not classified as oscillatory. In the case of the oscillatory sessions, (varphi left(tright)) 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 (1le ile N,1le jle 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 measurement (varDelta tau =frac{496,{rm{s}}}{96} sim 5,{rm{s}}) and the vary of doable d values—that’s, [−π, π) 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(tright)) (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(tright)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(tright)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, (1le ile 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, (1le ile 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(tright)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}{10times {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(tright)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,jne 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 instance session of mouse no. 104368 (Fig. 2f) and [1,100, 1,400] s for the instance session of mouse no. 102335 (Prolonged Information Fig. 4g). Word, nonetheless, that sequences had been recognized exterior these session subsets too, indicating that the sorting unveils stereotyped sequences additionally exterior the used subsets of information (see ‘Sortings are preserved when completely different parts of information are used for acquiring the sortings’).
Locking to the section of the oscillation
To calculate the extent to which particular person cells within the calcium imaging experiments had been tuned to the oscillatory sequences, two portions had been used: the locking diploma and the mutual data between the calcium occasion counts and the section of the oscillation. For every oscillatory session, the section of the oscillation (varphi left(tright)) was computed (see equation (1)) and particular person sequences had been recognized (see ‘Identification of particular person sequences’). Subsequent, the time factors that corresponded to all particular person sequences in a single session had been concatenated, which generated a brand new sign with the section of the oscillation for all consecutive sequences, and a brand new matrix of calcium exercise through which the community engaged within the oscillatory sequences uninterruptedly.
The locking diploma was computed for every cell because the imply resultant vector size over the phases of the oscillatory sequences at which the calcium occasions occurred (bin measurement = 129 ms, perform circ_r from the Round Statistics Toolbox for Matlab^{66}). The locking diploma has a decrease sure of 0 and higher sure of 1. It is the same as 1 if all oscillation phases at which the calcium occasions occurred are the identical (that’s, excellent locking), and equal to zero if all phases at which the calcium occasions occurred are evenly distributed (complete absence of locking). To estimate significance, for every cell a null distribution of locking levels was constructed by temporally shuffling the calcium exercise of that cell 1,000 occasions whereas the section of the oscillation remained unchanged, and by computing, for every shuffle realization, the locking diploma (shuffling was carried out as in ‘Sorting of temporally shuffled knowledge’). The 99th percentile of the estimated null distribution was used as a threshold for significance.
As a way to assess the robustness of the locking diploma, the obtained outcomes had been in contrast with a second measure based mostly on data principle^{67}: the mutual data between the counts of calcium occasions (occasion counts) and the section of the oscillation (bin measurement = 0.52 s). To estimate the discount in uncertainty concerning the section of the oscillation (P) given the occasion counts of the calcium exercise (S), Shannon’s mutual data 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)},$$
the place ({rm{Prob}}left(p,sright)) is the joint chance of observing a section of the oscillation p and an occasion depend s, ({rm{Prob}}left(sright)) is the marginal chance of occasion counts and ({rm{Prob}}left(pright)) is the marginal chance of the section of the oscillation. All chance distributions had been estimated from the information utilizing discrete representations of the section of the oscillation and the occasion counts. The occasion counts had been partitioned into s_{max} + 1 bins to account for the absence of occasion counts in addition to all doable occasion counts, the place s_{max} is the utmost variety of occasion counts per cell in a 0.52 s bin, and the section of the oscillation was discretized into 10 bins of measurement (frac{2{rm{pi }}}{10}).
The mutual data is a nonnegative amount that is the same as zero solely when the 2 variables are impartial—that’s, when the joint chance is the same as the product of the marginals ({rm{Prob}}left(p,sright)={rm{Prob}}left(pright){rm{Prob}}left(sright)). Nonetheless, restricted sampling can result in an overestimation within the mutual data within the type of a bias^{69}. As a way to appropriate for this bias, the calcium exercise was temporally shuffled (as in ‘Sorting of temporally shuffled knowledge’) and the mutual data between the occasion counts of the shuffled calcium exercise and the section of the oscillation, which remained unchanged, was calculated. This process, which destroyed the pairing between occasion counts and section of the oscillation, was repeated 1,000 occasions and the common mutual data throughout the 1,000 iterations was computed and used as an estimation of the bias within the mutual data calculation. In the suitable panel of Fig. 3a, we report each the mutual data and the bias. In Prolonged Information Fig. 7a, the corrected mutual data was reported (MI_{c}), the place the bias (⟨MI_{sh}⟩_{iterations}) was subtracted out from the Shannon’s mutual data (MI): MI_{c} = MI − ⟨MI_{sh}⟩_{iterations}.
Word that the locking diploma and the mutual data between the occasion counts and the section of the oscillation yielded constant outcomes (see Fig. 3a and Prolonged Information Fig. 7a).
Tuning of single cells to the section of the oscillation
The selectivity of every cell to the section of the oscillation within the calcium imaging knowledge was visualized by means of tuning curves and quantified by means of their most wellliked section. As within the evaluation of ‘Locking to the section of the oscillation’, the section of the oscillation (varphi left(tright)) was computed, particular person sequences had been recognized, and the time factors of the section of the oscillation and the matrix of calcium exercise that corresponded to all particular person sequences in a single session had been concatenated.
Tuning curves
The vary 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 }}+jfrac{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}Qleft(xright){log }_{2}left(Qleft(xright)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), (1le ile 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(i1right)times frac{N}{10}+1,left(i1right)times frac{N}{10}+2,ldots ,itimes 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 (N9times lfloor frac{N}{10}rfloor ) cells, where (lfloor xrfloor =max{min {mathbb{N}}/mle 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 ,itimes lfloor frac{N}{10}rfloor right},,1le {rm{e}}{rm{n}}{rm{s}}{rm{e}}{rm{m}}{rm{b}}{rm{l}}{rm{e}}le 9 left{9times lfloor frac{N}{10}rfloor +1,9times lfloor frac{N}{10}rfloor +2,ldots ,Nright},,{rm{e}}{rm{n}}{rm{s}}{rm{e}}{rm{m}}{rm{b}}{rm{l}}{rm{e}}=10end{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 (ito 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,jright)) 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, (ito 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}}}), (1le ile 500,,0le tle 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, (1le ile 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.