Animals

All animal procedures were performed in accordance with national and international guidelines and were approved by the local health authority (Das Landesamt für Natur, Umwelt und Verbraucherschutz). For this study, 10–25-week-old Vgat-ires-cre knock-in mice (The Jackson Laboratory) and C57BL/6 male mice were used, except for studies using MoSeq, which involved female mice. Mice were housed under standard conditions (air temperature 20–24 °C, relative humidity 45–65%) in the animal facility and kept on a 12 h light–dark cycle. Before all experiments, mice were handled by the experimenter and habituated to the experimental enclosure for 3–5 days63. This habituation procedure is important for minimizing the potential influence of unfamiliar experimental procedures or enclosure’s novelty on innate behaviors, for example, to ensure that animals consume food pellets in the experimental enclosure. Before the experiments with the optogenetic manipulations contingent on feeding, food was taken out from home cages for about 1 h; mice received water ad libitum.

Viral injections

Viral injections in the LH, LPO and mPFC were performed according to previously described protocols10,64. Mice were treated with buprenorphine (0.1 mg kg−1), anesthetized with isoflurane and placed in the stereotactic apparatus (David Kopf Instruments). A small hole was drilled in the skull with a dental drill for each virus injection site according to the stereotactic coordinates. A sterile glass pipette made using a micropipette puller (Sutter Instruments) was mounted on a syringe (Hamilton CS-Chromatographie Service) to infuse viruses at a rate of 100 nl min−1; injection volume and speed were controlled with a micro-pump (Harvard Apparatus, Hugo Sachs Elektronik). After the injection, the injection pipette remained in the injection area for about 10 min and then was slowly lifted before the incision was sutured. Optogenetic constructs were purchased from the University of North Carolina (UNC) Vector Core or provided by K. Deisseroth. For manipulation of LH and LPO Vgat cells, Vgat-cre mice were injected bilaterally into the LPO (anterior-posterior (AP) 0 mm, mediolateral (ML) ± 1 mm, dorsal-ventral (DV) −5 mm and −5.25 mm) with 0.3 μl per injection site of AAV8-Ef1a-DIO-ChRmine-mScarlet (provided by K. Deisseroth, titer 5 × 1012 vg ml−1) or 0.3 μl per injection site of AAV8-Ef1a-DIO-mScarlet (provided by K. Deisseroth, titer 5 × 1012 vg ml−1). In the LH (AP −1.7 mm, ML ± 1 mm, DV −5 mm and −5.25 mm), 0.3 μl per injection site of AAV2-Ef1a-DIO-ChETA-eYFP (UNC Vector Core, titer 3.5 × 1012 vg ml−1) or 0.3 μl per injection site of AAV2-EF1a-DIO-eYFP (UNC Vector Core, titer 4.5 × 1012 vg ml−1) were injected bilaterally. For manipulations of the mPFC–LH or mPFC–LPO projections, C57BL/6 mice were injected bilaterally in the mPFC (AP 1.7 mm, ML ± 0.3 mm, DV −2.4 mm and −2.8 mm) with 0.2 μl per injection site of eNPAC2.0 (provided by K. Deisseroth, titer 1.84 × 1013 vg ml−1) or AAV5-hSyn-eYFP (UNC Vector Core, titer 3.3 × 1012 vg ml−1).

Implantation of optic fibers and electrodes

Optic fiber implants were manufactured from 100-μm diameter multimode optic fiber (numerical aperture 0.22) and zirconia ferrules (Thorlabs). For optogenetic manipulations of signaling between the LH and LPO, mice were implanted with optic fibers in the LH (AP −1.7 mm, ML 1 mm at a 21.8° angle, ML −1 mm, DV −4.7 mm) and in the LPO (AP 0 mm, ML −1 mm at a 21.8° angle, ML 1 mm, DV −4.7 mm). For the optogenetic manipulations of the mPFC–LH projections, mice were implanted bilaterally with optic fibers in the LH (AP −1.7 mm, ML 1 mm at a 21.8° angle, ML −1 mm, DV −4.7 mm). For the optogenetic manipulations of the mPFC–LPO projections, mice were implanted bilaterally with optic fibers in the LPO (AP 0 mm, ML 1 mm at a 21.8° angle, ML −1 mm, DV −4.7 mm). For the extracellular neuronal and LFP recordings, silicon probes (B32, NeuroNexus Technologies) were mounted on custom-made microdrives and implanted as described previously10,64,65. For the mPFC and LH simultaneous recordings, the following implantation coordinates were used: mPFC (AP 1.7 mm, ML 0.2 mm, medial shank, DV −2.4 mm) and LH (AP −1.58 mm, ML 0.8 mm, medial shank, DV −4.9 mm). For the LH and LPO recordings, the following implantation coordinates were used: LH (AP −1.58 mm, ML 0.8 mm, medial shank, DV −5 mm) and LPO (AP 0 mm, ML 0.5 mm, medial shank, DV −5 mm) combined with optic fibers implanted at a 21.8° angle in the LH (AP −1.7 mm, ML 1 mm, DV −4.7 mm) and the LPO (AP 0 mm, ML 1 mm, DV −4.7 mm). For the VTA recordings, the following implantation coordinates were used: VTA (AP −3.1 mm, rostral shank, ML 0.4 mm, DV −4.2 mm). For simultaneous LFP recordings from the LH, LPO and VTA, a custom stationary probe (four recording sites × eight shanks, NeuroNexus Technologies; Fig. 7a) was implanted along the line defined by the following coordinates: first shank (AP 0 mm, ML 1 mm, DV 5.4 mm) and last shank (AP −3.8 mm, ML 0.5 mm, DV 4.8 mm).

Data acquisition

The recording setup was a custom-made enclosure10 (length/width/height 50 × 30 × 20 cm) with two interconnected compartments (25 × 30 × 20 cm each). Water presented in a water cup, food provided in a food cup and a new object from Lego or similar toy sets were placed in three corners of the enclosure. Mice were freely moving in the enclosure during the recordings. During the recordings, silicon probes were connected to a preamplifier (NeuraLynx) to eliminate cable movement artifacts. Signals were differentially amplified and band-pass-filtered (1 Hz–8 kHz) and acquired continuously at 32 kHz (Digital Lynx, NeuraLynx). Synchronization with the acquisition of electrophysiological data recording of the animals’ behavior was performed from different angles by four cameras at 25 Hz (Motif, Loopbio). A light-emitting diode was attached to the headset to track the animal’s position at 25 Hz using a top-mounted camera. For pose estimation using DeepLabCut, the behavior of pairs of mice in the enclosure was recorded at 15 Hz. For behavioral motion segmentation (MoSeq21), female mice were recorded for 20 min while they were freely exploring an arena (length/width/height 45 × 25 × 40 cm) with a female conspecific and either a new object (piece of Lego) or (high-fat) food behind a mesh (length/width/height 8 × 8 × 7 cm) on the left and right sides, respectively. Behavior in these experiments was captured at 30 Hz with a depth camera (Kinect for Windows v.2, Microsoft) positioned 65 cm above the floor of the arena.

Optogenetic stimulation

All mice used in the behavioral assays were allowed to recover after the fiber implantation for at least 1 week. Mice were randomly assigned to the experimental conditions. For optogenetic manipulation, 473 nm and 589 nm diode-pumped solid-state lasers (Laserglow Technologies) were used. For the stimulation of projections of LPO cells expressing ChRmine, a light delivery from a 589-nm laser was controlled by a shutter (Doric Lenses). Stimulation protocols were implemented using a stimulus generator (Multi Channel Systems). One side of the patch cord was connected to the implanted optical fiber with a zirconia sleeve (components from Thorlabs) and the other side was connected to the laser with an FC/PC adapter. The optogenetic experiments were performed in the test enclosure described above, once for each type of the stimulation and 14 ± 7 times for different types of optogenetic experiments lasting approximately 20 min each. During the electrophysiological recordings (38 ± 5 sessions per mouse), time stamps of laser pulses were acquired synchronously with neuronal signals and video frames. The behavior of mice was recorded from different angles by four cameras at 25 Hz (Motif, Loopbio).

Optogenetic manipulations of the LH–LPO circuit

For a closed-loop optogenetic manipulation of the LH–LPO circuit in mice expressing ChRmine-mScarlet in LPO Vgat cells and ChETA-eYFP in LH Vgat cells, beta out-of-phase, nonrhythmic or beta in-phase stimulation was applied unilaterally (Extended Data Fig. 6a–e). A separate control group of mice expressing mScarlet in LPO Vgat cells and eYFP in LH Vgat cells, without optogenetic actuators, also received closed-loop stimulation with the beta out-of-phase protocol. Beta out-of-phase stimulation consisted of 5-ms 589-nm light pulses in the LH and 5-ms 473-nm light pulses in the LPO at 20 Hz with a 25-ms offset between brain regions. During nonrhythmic stimulation, the amount of light irradiation was matched to a 10% duty cycle of the beta out-of-phase and in-phase protocols: 20 pulses, 5 ms each, were randomly assigned times outside the beta band (mean interpulse interval 5 ms) during 200-ms epochs of each 1-s window. These 200-ms epochs of the 589-nm stimulation in the LH and 473 nm in the LPO did not overlap. The light power output was 1–4 mW during the light pulses measured at the tip of each of the two patch cords using an optical power meter (Thorlabs). Beta out-of-phase, nonrhythmic or beta in-phase stimulation was started when an animal spontaneously initiated F, S or E, in separate experiments, and lasted for 5 or 10 s for each manipulated behavioral episode during four corresponding 5-min parts of a 20-min session (Extended Data Fig. 6i). The stimulation was repeated each time when an animal engaged in the investigated behavior. The time elapsed from stimulation onset to the end of the behavioral episode was defined as the latency to behavioral transition.

Noncontingent on the animals’ behavior, either beta out-of-phase or unidirectional LPO–LH stimulation (Extended Data Fig. 6d) was applied for 20 min in repeated blocks of 10-s stimulation alternating with 20-s breaks.

Optogenetic manipulations of mPFC–LH and mPFC–LPO projections

In mice expressing eNPAC2.0-eYFP or eYFP in the mPFC, stimulation was performed in the LH or in the LPO. Each experimental session lasted for 20 min. Pulses of 473-nm light for 5 ms at 20 Hz, a light power output of 10–15 mW from the tip of the patch cord or 589-nm light for 10 s continuously and a light power output of 20 mW from the tip of the patch cord were applied in the optogenetic excitation or inhibition experiments, respectively. The closed-loop optogenetic stimulation was triggered by the spontaneous initiation of behaviors and lasted for 10 s for each manipulated behavioral episode. For the social behavior tests, stimulation at theta frequency (9 Hz: 11-ms light-on, 100-ms light-off phases for 10 s in each stimulation episode) and nonrhythmic stimulation with light intensity matched to the beta frequency protocol (as described for the LH–LPO circuit stimulation) were performed.

Brain dissection and imaging

After completion of the experiments, mice were deeply anesthetized and electrolytic lesions at selected recording sites were performed to visualize the locations of the recording electrodes. Mice were perfused with 4% paraformaldehyde in PBS. Brains were fixed overnight in paraformaldehyde, placed for cryoprotection in 30% sucrose at 4 °C for 24 h and then coronally or sagittally sectioned into 40-μm slices on a cryostat (CM1900, Leica Biosytems). To visualize the electrolytic lesions and silicon probe tracks, brain sections were imaged using a widefield Axio Imager M2 microscope (ZEISS). To visualize the projections and control viral expression, sections were imaged using a confocal microscope (Leica SP8, Leica Biosytems).

Behavioral analysisBehavioral scoring

Ethograms were obtained using a frame-by-frame scoring of behaviors using Adobe Premiere Pro (v.2020) (Adobe) in multiangle synchronized video recordings66. Frames when a resident mouse (implanted with electrodes or optic fibers) was consuming food pellets were scored as feeding. Social contact was defined as sniffing or following an intruder mouse. During the stimulation of the LH–LPO circuit, the latter behavior evolved into a prolonged chasing, defined as uninterrupted pursuing of an intruder for longer than 2 s. New object exploration was defined as sniffing, gnawing, touching or climbing a new object.

DeepLabCut

Markerless pose estimation was performed with the DeepLabCut toolbox (v.2.2.0.2)22. First, k-means clustering and manual selection were performed to select frames from each video across behaviors. Six key points (snout, left ear, right ear, left side (middle-left part of body), right side (middle-right part of body) and tail base) of each animal were localized on each frame. In total, 610 labeled frames were selected across eight video recordings and used to train a multiscale deep learning model DLCR-Net_ms5. Randomly assigned 95% of the data were used for training and the rest for testing. The network was trained for 120,000 iterations until cross-entropy loss plateaued. The estimated coordinates of key points in each frame were used to define behaviors67: animal in the food zone (the distance between the food zone center and the snout or an ear was less than the radius of food zone, 3.6 cm); feeding (snout and both ears in the food zone for at least 1.3 s); mouse in the water zone (the distance between the water zone center and the snout or an ear was less than the radius of the water zone, 3 cm); drinking (snout in the water zone for at least 1.5 s); social contact (snout or an ear inside or on the edge of the polygon area defined by the six key points of an intruder mouse); new object exploration (the distance between the new object zone and the snout or an ear was less than 1.3 cm); rearing (the snout was at least 2.2 cm over the enclosure wall or the snout was at least 0.7 cm over the middle separator wall and the length of a vector between the snout and the tail base was less than 10.3 cm); and immobility (the coordinates of each of the six key points changed less than 0.5 cm s−1). Frames with pose patterns not meeting any of these criteria were classified as behaviorally undefined.

Behavioral motion segmentation

Using custom Python scripts (adapted from MoSeq (v.1)21 by R. Ung from the G. Stuber’s laboratory), depth images and frame time stamps were converted into a binary format for further analysis. Region-of-interest polygons delimiting the boundaries of the arena and two-dimensional images to inspect the behavioral syllables after analyses were acquired simultaneously.

MoSeq was performed in a Debian GNU/Linux 8 virtual environment running on a Linux (Ubuntu 16.94.3 LTS) compute cluster. Behavior was classified using MoSeq v.1 (ref. 21). Briefly, depth mouse images were cropped along the arena boundaries, extracted from the arena background, parallax-corrected and orientated along the spine axis. Time series data were subjected to wavelet transformation and dimensionally compressed using principal component analysis. To classify the behavioral syllables, an autoregressive hidden Markov model was applied to the first ten principal components. A template-matching procedure ensured that only repeated principal component trajectories (that is, meaningful ones) were selected as syllables. One of the model parameters, the self-transition bias kappa, was set to match the median syllable duration with the median approximate change point of each dataset identified using a filtered derivative algorithm (κ = 5). To qualitatively verify behavioral syllables, we manually assessed the visualizations of each syllable using the two-dimensional recordings.

Analyses of syllable usage were performed in Python v.2.7. A narrow zone of 10 × 5 cm before the mesh with food, conspecific or new object was defined as a contact zone. The 4-s periods just before entering the zone (with the center of the head), with a minimum zone visit duration of 333 ms and a minimum zone visit interval of 2 s, was divided into eight 0.5-s bins (hence, for example, ‘−2 s’ corresponds to the period from 2 s to 1.5 s before social contact). Frames corresponding to the previous zone visit were excluded from the transition periods. For the calculation of syllable usage during contact, all frames inside the contact zone were included. The probability of syllable usage during randomly selected 2-s epochs and repeated 1,000 times, equal or greater than the usage during the transition epochs, was computed and normalized across bins.

Electrophysiological data analysisSpike sorting and unit characterization

Electrophysiological signals were preprocessed using NDManager (http://neurosuite.sourceforge.net/)68 and analyzed using custom-written MATLAB v.2014b algorithms (MathWorks) as described previously64. Action potentials (spikes) were detected in a high-pass filtered signal and spike waveforms were represented by the first three principal components and by the amplitudes of the action potentials. Spike sorting was performed automatically69 (https://github.com/klusta-team/klustakwik) followed by manual cluster adjustment based on auto-correlations and cross-correlations of spikes trains, the Mahalanobis distance between pairs of clusters and the visual comparison of waveform profiles across channels68 (Extended Data Fig. 1b). Isolation distance69 was computed for the sorted units: LH = 62.5 ± 0.9, n = 2,417 cells; LPO = 63.8 ± 2.0, n = 415 cells; mPFC = 66.9 ± 1.1, n = 2,374 cells; VTA = 82.0 ± 3.2, n = 308 putative dopamine cells.

For individual behaviors, we computed the firing rate of cells. A surrogate distribution of 1,000 firing rate values was obtained for each cell by 2–4-min offsets of the behavior time stamps. For each behavior, the match score was calculated as the percentile of the firing rate during the behavior in the surrogate distribution. Multimodal LH cells were defined as units with a firing rate preference in the upper quartile for each of the three behaviors.

LFP analysis

The LFP was obtained by downsampling the wide band signal to 1,250 Hz using NDManager68. High-resolution time frequency analysis was performed using a continuous Morlet wavelet transform. The multitaper method (NW = 3, window length of 1,024) was used to compute power spectral density and coherence according to the ethogram times. Beta oscillations were detected in the 15–30 Hz band-pass-filtered, rectified and smoothed signal. Events with amplitudes exceeding 2 s.d. above the noise mean for at least 80 ms were detected. The beginning and the end of the oscillatory epochs were designated at times when the amplitude fell below 1 s.d.

Discharge phases

Spikes fired during the detected oscillation episodes were assigned beta oscillation phases, computed using the Hilbert transform of the 15–30 Hz filtered signal. Histograms of spike counts in 20 phase bins were convolved with a Gaussian kernel (size = 0.65 s.d.) and normalized by the total number of spikes in the histogram65. This approach was also used to compute the discharge phases during the gamma oscillations (30–60 Hz, minimum duration of 25 ms (ref. 10).

To examine the timing of neuronal discharge in the LH and LPO during beta out-of-phase stimulation, we used a linear approximation of the 20-Hz sinewave as a reference for the spike phase assignment. The times of blue light pulses stimulating projections of LH cells in the LPO and of red light pulses stimulating projections of LPO cells in the LH defined the period of the stimulation rhythm for the assignment of phases to the spikes of LH and LPO cells, respectively. As we optogenetically stimulated the inhibitory inputs from the LH to the LPO and from the LPO to the LH, we evaluated the proportion of units inhibited by the optogenetic stimulation. For this purpose, we summed the normalized binned firing probability within the first 7 ms after pulse termination and the normalized binned firing probability within the following 7 ms. We calculated the ratio of these sums and detected any outlier units defined as more than three scaled median absolute deviations away from the median. Units falling below the 30th percentile in the ratio distribution were defined as inhibited units, the population firing probability of which was summarized in stimulation phase histograms. The onset of the first stimulation pulse was assigned as phase π and the onset of every second pulse was assigned as phase −π. Every mid-interpulse interval was assigned as phase 0 radian. Other phases were linearly interpolated at 20 kHz (the sampling rate of spike trains). Each spike was assigned a corresponding beta phase in the stimulation cycle. The obtained phases were offset by π for the LH spike strains (that is, stimulation of the LH cell projections at 0° and 360°) and by 3 π for the LPO spike trains (that is, stimulation of the LPO cell projections at 180° and 520°) according to the out-of-phase timing of the pulses in this stimulation protocol. The spike phase distribution of each unit was binned into 20 bins per beta cycle.

Machine learning modelingPhase signatures

Firing probability versus beta oscillation phase histograms were computed for individual LH and VTA cells using the spikes fired during 2 s before transition to F, S and E separately for each behavior for control or behavior epochs of the same duration. Control epochs excluded transitions to the aforementioned behaviors. In a separate analysis, transition and control epochs were additionally selected for the same behavioral state, locomotion and posture change.

Cells with histograms containing at least 20 (168 ± 15) spikes for the LH cells and at least 10 (29 ± 3) spikes for the presumed dopamine cells (VTA cells with a spike width greater than 0.3 ms (ref. 29) and firing rate lower than 10 Hz (ref. 47) were used for the subsequent decoding of behavioral transitions. For each of the 20 phase bins, the population distribution of firing probabilities was estimated; cells with the firing probability in the upper quartile of the distribution, that is, ‘highly active cells’ at a given phase, were selected. A phase signature was defined through an asymmetry of individual behavior match scores’ distribution in a population of highly active cells in each phase bin as:

where \(c=|\|\) and D is a set of match scores m of N cells.

Phase signatures were computed based on match scores for F, S and E during transitions to these three behaviors resulting in three behavior-specific phase signatures for each of the three types of transitions or, if specified, for a combination of different types of transitions. To account for the variability of the phase signatures in each phase bin, the distribution of match scores in the set of highly active cells at a given phase was bootstrapped with replacement 1,000 times to derive the datasets for the modeling70. To generate the control sets, the time stamps of the transitions were randomly offset, excluding overlaps with the transition epochs from native ethograms. The first 1,000 offset trials with the number of spikes that were sufficient for the estimation of firing probability and the phase histograms were selected.

Support vector classifications

SVM models were implemented using the Python package Scikit-learn (v.1.2.2)71. Phase signatures from either individual or multiple phase bins were the inputs to one SVM. The classes and input datasets for the SVMs are described below (see also the design description of the SVM model in Supplementary Information).

Models 1 and 2 in Fig. 2a–c aimed to classify transition (2-s epochs) versus control epochs (2 s, excluding transitions to the three behaviors) within individual phase bins (eight bins) in the peak neighborhood (peak ± 72°). To do so, an SVM was trained and tested with tenfold cross-validation. Model 1 was trained and tested on phase signatures related to individual behaviors (F, S and E), while model 2 used the phase signatures of all three behaviors in one SVM.

Model 3 in Fig. 2d assessed the phase signatures (related to F, S and E) in transition versus control epochs across phases near the peaks of beta oscillations. To generate the phase-shuffled datasets, the phase of each spike was jittered by a random offset from a uniform distribution. Then all spikes of all cells were additionally offset by the same random phase between 2.5 radian (the width of the peak neighborhood) and π using different random offsets for control and transition epochs. Separate SVMs were computed on phase signatures (related to F, S and E) from the original and phase-shuffled datasets to classify transition versus control epochs. Training was done on the data from the phase bin with the highest (in the entire oscillation cycle) difference of phase signature amplitude between transition and control. Training was performed this way in the original and in the phase-shuffled data. Subsequently, the SVM trained on the original dataset was tested on individual phase bins in the peak neighborhood in the original dataset, excluding the bin used for training. Testing in that bin was performed using a separate SVM, trained on the phase bin with the second highest (in the entire oscillation cycle) amplitude phase signature. Testing of the SVM trained on phase-shuffled data was performed on the phase-shuffled data from the bin, in which the original dataset could be decoded with the highest accuracy. The resulting decoding accuracy in the phase-shuffled data was close to the chance level, which was typical for phase-shuffled data also in other phase bins.

Model 4 in the statistical information in Supplementary Information for Fig. 3c and Extended Data Fig. 5b classified three upcoming behaviors (F, S and E) using in one SVM phase signatures from all eight phase bins in the peak neighborhood during transitions. The SVMs were trained and tested on the phase signatures of individual behaviors (F, S and E) with tenfold cross-validation.

Model 5 in the statistical information in Supplementary Information for Fig. 3c and Extended Data Fig. 5b was similar to model 4 except that, instead of upcoming behaviors, it classified three current behaviors (F, S and E) using the phase signatures during random 2-s epochs of behaviors.

Models 6 and 7 in Fig. 3d were similar to models 1 and 2, respectively, but they classified three upcoming behaviors (F, S and E) using the phase signatures in individual phase bins in the peak neighborhood during transitions.

Models 8 and 9 in Fig. 3e were similar to models 1 and 2, respectively, but they classified three current behaviors (F, S and E) using the phase signatures in individual phase bins in the peak neighborhood. These SVMs were trained and tested on phase signatures during random 2-s epochs of behaviors (F, S and E).

Model 10 in Fig. 7i classified transition (2-s epochs) versus control epochs (2 s, excluding transitions to the above behaviors) using in one SVM the phase signatures from all phase bins in the entire cycle to account for the phase offset of dopamine cell discharge in relation to the LH (LH in Fig. 4a; VTAdopamine in Extended Data Fig. 10d).

Models 11 and 12 in Extended Data Fig. 4d,e were similar to model 3 except that they classified pooled transitions (2-s epochs preceding any of the three behaviors, that is, F, S and E) versus control epochs (2 s, excluding transitions to the three behaviors), either using the data from all mice pooled (model 11) or separately from individual mice (model 12).

In the LH recordings, as described above, two-class and three-class (one-versus-rest multiclass classification72) models were computed using a nonlinear radial basis or linear kernel (depending on the dimensionality of the feature space). Linear SVMs were used to classify the population activity of putative dopamine cells. To minimize overfitting, training and testing were done on different data subsets. Except for models 3, 11 and 12, which were designed to be trained and tested on different phase bins, a stratified tenfold cross-validation procedure was used: each training set was randomly divided into ten subsamples with the same proportion of samples from each class as in the complete set. One subsample was then retained for testing the model, while the other nine subsamples were used for training, with this procedure repeated using all ten subsamples so that each subsample was used only once to evaluate the performance of the model. Decoding accuracies were initially averaged within tenfold cross-validation trials and the resulting accuracies were averaged across 1,000 repeated cross-validations.

The significance of classifications was assessed using permutation tests. We randomly permuted the labels and then used the same decoding approach as for decoding the original labels, except for using the stratified tenfold cross-validation once (instead of 1,000 times) for each model. The permutation of labels was repeated 1,000 times to assess the chance performance of a classifier computed as the average of accuracies across permuted sets. The performance of a classifier was considered significant when it fell in the 5% upper tail of its chance performance distribution.

Statistical analyses

Statistical analyses were performed using MATLAB v.2014b (MathWorks), Python v.3 (https://www.python.org/) or Prism 9 (GraphPad Software). The level of significance and the number of neurons and mice are indicated in the figure legends. A likelihood ratio test73 was used to compare bivariate circular distributions (see the statistical information related to Fig. 5d,e). All statistical tests were two-tailed unless indicated otherwise; permutation and randomization tests were right-tailed. Two-group comparisons were performed using a t-test, Mann–Whitney U-test or Wilcoxon matched-pairs test depending on the normality of a distribution. Multiple group comparisons were performed using an ANOVA or multiple two-group tests with α correction, adjusting for multiple comparisons. The Grubbs’ test was used to exclude outlier points from behavioral datasets. A median absolute deviation outlier test was used to exclude outlier points from the analysis of optogenetic entrainment. No further data points or animals were excluded. Sample size was determined according to the accepted practice for the applied assays. No statistical methods were used to predetermine sample sizes; sample sizes are similar to those reported in previous publications38,49,54,66. Data analysis was performed blindly using automatic selection of data from a database. The full description of the statistical analyses corresponding to each dataset is provided in the statistical information in Supplementary Information. Unless specified otherwise, descriptive statistics are reported as the mean ± s.e.m.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.