Subjects

All procedures were approved by the Animal Care Committee of the Montreal Neurological Institute at McGill University in accordance with Canadian Council on Animal Care guidelines (protocol MNI-7839). The subjects (n = 43 for electrophysiology and n = 6 for anatomy) were adult (>8 weeks old) male mice bred by crossing wild-type females on C57BL/6J background (Jackson Laboratories, 000664) with either homozygous male VGAT-IRES-Cre mice (Jackson Laboratories, 028862; n = 38), PV-IRES-Cre mice (Jackson Laboratories, 017320; n = 2) or SST-IRES-Cre mice (Jackson Laboratories, 013044; n = 2). One additional mouse implanted with a Neuropixel probe (Fig. 1a–c) was a cross-bred C57BL/6J and FVB (Jackson Laboratory, 001800). Mice were kept in standard conditions (room temperature and 50% humidity) on a 12-h light/12-h dark cycle and were housed in group cages (2–5 mice per cage) before electrode implantation surgery and individually afterwards.

Electrode implantation

Mice were implanted under isofluorane anesthesia. Silicon probes were mounted on in-house built movable microdrives and implanted through a small craniotomy. Probes were implanted either vertically above left PoSub (from Bregma: AP, −4.24 mm; ML, 2.05 mm; DV, −1.00 mm)60 or at a 26° angle pointing away from the midline into left pRSC (from Bregma: anterior-posterior (AP), −4.24 mm; medial-lateral (ML), 1.70 mm; dorsal-ventral (DV), −1.00 mm). A mesh wire cap was then built around the implanted microdrive and was reinforced with UV-cured adhesive. Mice were allowed to recover for at least 1 d before electrophysiological recordings.

Each mouse was implanted with one probe. The probes were either a Neuropixel 1.0 probe (384 active sites arranged in a dense checkerboard layout, that is two columns, 20 μm between each row), a single shank with 64 recording sites (H5; Cambridge NeuroTech) or four shanks with 8 recording sites each (Buzsaki32; NeuroNexus). In all experiments, both ground and reference wires were soldered to a single 100 μm silver wire, which was then implanted above the cerebellum.

Behavioral procedures

Before the implant surgery, mice were habituated over several days to forage for small pieces of Honey Cheerios cereal (Nestle) in the open field. For most recordings, the recording chamber consisted of a metal frame (90 × 90 × 180 cm) supporting a plastic platform with removable walls (width, 80 cm; height, 50 cm) that could be arranged into either a square or triangular open field. Recordings of two mice implanted in ADN were conducted in a circular open field (diameter, 85 cm; height, 50 cm). The recording protocol consisted of a sleep session in the home cage, followed by open-field exploration in a square arena and another sleep session. A subset of animals then explored a triangular arena. A white rectangular cue card on one of the walls served as a salient landmark.

Recording procedures

During the recording sessions with 32- or 64-channel silicon probes, neurophysiological signals were acquired continuously at 20 kHz on a 256-channel RHD USB interface board (Intan Technologies) and captured with Intan RHX software (Intan Technologies). For the Neuropixel recording, neurophysiological signals were acquired continuously at 20 kHz on a control board (Native Instruments) and captured with SpikeGLX 3.0 (https://billkarsh.github.io/SpikeGLX/). The wide-band signal was downsampled to 1.25 kHz and used as a local field potential(LFP) signal. The recording cables were tethered to a motorized electrical rotary joint (AERJ; Doric Lenses). Ahead of the main recording session, the microdrive was lowered over several hours in small (35–70 μm) increments until the whole shank was positioned in PoSub or ADN. A short open-field session was then recorded to map the HD receptive fields of all neurons. For PoSub, the recording depth was adjusted so that sharply tuned HD cells (a hallmark of PoSub) were present along the whole length of the shank. Data collection did not commence until at least 2 h after the last depth adjustment. For animals implanted with single-shank linear probes, only one session per mouse was included in each analysis to prevent double-counting of cells. For animals implanted with four-shank Buzsaki32 probes, multiple sessions per mouse (obtained on separate days) were included in the analysis, ensuring that the probe was moved by at least 70 μm between the recording sessions.

Animal position and orientation were tracked in 3D using seven infrared cameras (Flex 13; Optitrack) placed above the enclosure and coupled to the Motive 2.0 motion capture system (Optitrack). Seven small tracking markers were attached to the headcap. Video recording was captured by an overhead camera (Flex 13; Optitrack) placed close to the rotary joint. Animal position and head orientation were sampled at 100 Hz and were synchronized with the electrophysiological recording via voltage pulses registered by the RHD USB interface board (Intan Technologies).

Optogenetic experiments

VGAT-IRES-Cre mice were injected with an adeno-associated virus vector (AAV 2/9 CAG-Flex-ArchT-EGFP or CAG-Flex-EGFP, titer: 4–5 × 1012 genomic copies (GC) per ml; Neurophotonics) into TRN (from Bregma: AP, −0.70 mm; ML, 1.25 cm; DV, −3.25 cm) under isofluorane anesthesia either unilaterally for ADN recordings or bilaterally for PoSub recordings. Injections (300–400 nl per injection site) were done with a microinjector (Harvard Apparatus) through a small craniotomy, at the speed of 100 nl s−1. The needle was left in place for 2–5 min after injection.

At least 3 weeks after the injection surgery, optic fiber implants (Doric Lenses; MFC_200/240-0.22_25mm_SM3) were implanted unilaterally (left hemisphere, ADN recordings) or bilaterally (PoSub recordings) above ADN at a 20° angle from the sagittal plane (from Bregma: AP, −0.82 mm; ML, 1.00 cm; DV, −2.25 cm). Mice were then implanted with a microdrive-mounted Buzsaki32 probe above left ADN (Bregma: AP, −0.82 mm; ML, 0.85 cm; DV, −2.00 cm) or either Buzsaki32 probe (Neuronexus) or H5 probe (Cambridge Neurotech) above left PoSub, as described above.

Laser light was delivered from a 520 nm fiber-coupled laser diode module (Doric Lenses) controlled with a laser diode module driver (Doric Lenses). Light power output at the tip of the fiber implant was measured before each implantation and an output curve was calculated individually for each implant. Light output was then set to 14–16 mW before each recording session. For this subset of mice, the second sleep session was followed by a second exploration session in the open field, during which a laser stimulation protocol was delivered via patch cords attached to the optic fiber implants. After 5 min of exploration, light pulses (1 s) were delivered at 0.2 Hz in groups of 60 (5 min total), each followed by 5 min of no stimulation. Four such epochs were delivered in total, resulting in 240 light pulses over a 45-min recording session.

Cue rotation experiment

A subset of animals implanted into PoSub with linear probes underwent a cue rotation experiment in a separate recording session. To this end, the frame of the recording chamber was fitted with black plastic insets that covered the floor (90 × 90 cm) and walls to obstruct any visual cues. Each of the four walls had an identical panel made of two light-emitting diode (LED) strips (yellow V-shape or blue X-shape) in the center. A small (diameter: 30 cm) elevated circular platform was placed in the center of the arena. Before each recording session, two adjacent LED panels were chosen as distal visual cues. The LED light intensity was titrated so that the panels were visible in the dark but provided minimal illumination of the surrounding area. The on/off cycle of the LED panels was controlled with an Arduino microcontroller. To reduce the specificity of olfactory and auditory cues, the whole recording chamber was thoroughly cleaned with antibacterial wipes before the experiment and white noise was emitted from speakers placed underneath the chamber.

The recording protocol consisted of a 1-h sleep session in the home cage followed by a 75-min cue rotation session. At the start of the cue rotation session, a single LED cue was switched on, and the mouse was placed directly on the circular platform and was left undisturbed for the duration of the session. After 15–20 min of exploration with a stable cue, the cue rotation protocol was initiated. The protocol consisted of the illuminated cue switching back and forth between two adjacent walls every 200 s for a total of 16 rotations. Additionally, to habituate the mouse to the cue disappearing from its field of view, the cue was switched off for 0.1 s every 20 s. To encourage the mouse to explore the platform for the whole duration of the session, the experiment was conducted in the middle of the dark phase of the light cycle and the recording chamber was sprayed with a new odor (scented air freshener) right before the cue rotation session.

Tissue processing and probe recovery

Following the termination of the experiments, animals were deeply anesthetized and perfused transcardially first with 0.9% saline solution followed by 4% paraformaldehyde solution. The microdrive was then retracted to remove the probe from the brain. Brains were sectioned with a freezing microtome coronally in 40 μm slices. Sections were washed, counterstained with DAPI and Green Neurotrace and mounted on glass slides with ProlongGold fluorescence antifade medium. Sections containing probe tracts were additionally incubated with a Cy3 anti-mouse secondary antibody (1:200 dilution; Cedarlane, 715-165-150) to help visualize the electrode tract. Widefield fluorescence microscope (Leica) was used to obtain images of sections and verify the tracks of silicon probe shanks, optic fiber position and virus expression.

Spike sorting and unit classification

Spike sorting was performed semi-automatically using Kilosort 2.0 (ref. 61) followed by manual curation of the waveform clusters using the software Klusters62 or Phy63. At this stage, any cluster without a clear waveform and clear dip in the spike train auto-correlogram at the 0–1 ms time bin was classified as noise and cluster pairs with similar waveforms and a clear dip in their spike train cross-correlograms at the 0–1 ms time bin were merged.

For PoSub recordings, viable units were first identified as units that (1) had an average firing rate of at least 0.5 Hz during open-field exploration and (2) had a waveform with negative deflection (criterion aiming to exclude spikes from fibers of passage). Next, putative excitatory cells and putative FS interneurons were classified on the basis of their mean firing rate during open-field exploration and the through-to-peak duration of their average waveforms (Extended Data Fig. 2a). Putative FS interneurons were defined as cells with short trough-to-peak duration (<0.4 ms) and high mean firing rates (>10 Hz). Conversely, cells with long trough-to-peak (>0.4 ms) and low mean firing rates (<10 Hz) were classified as putative excitatory cells.

HD tuning curves and tuning metrics

The animal’s HD was calculated as the horizontal orientation (yaw) of a polygon constructed in Optitrack tracking software by connecting the three-dimensional coordinates of all tracking markers on the animal’s headcap. The yaw of the polygon was measured in the global coordinates (that is, the axes of the environment, not the axes of the polygon), and these were constant across the whole study. HD tuning curves were then computed as the ratio between histograms of spike count and total time spent in each direction in bins of 1° and smoothed with a Gaussian kernel of 3° s.d. Tuning curves were computed from epochs when the animal’s speed exceeded 2 cm s−1 for all analyses except cue rotation and optogenetic experiments, where epoch duration (200 and 240 s, respectively) was too short to allow for further refinement.

To prevent any bias in HD cell population due to assumptions about the unipolar shape of HD cell tuning curves, we chose to define HD cells based on HD information contained in the tuning curves64, calculated for n angular bins as:

$$I=\sum _\lambda (_)_\left(\frac_)}\right)p(_)$$

where λ(Θi) is the firing rate of the cell for the ith angular bin, λ is the average firing rate of the neuron during exploration and p(Θi) is the occupancy (that is, normalized time spent) in direction Θi. This information rate (measured in bit s−1) was normalized by the cell’s average firing rate to provide an information content (measured in bit spike-1) independent of firing rate.

For each cell, we obtained the control tuning curve using a time-reversed HD angle—a method that preserves the dynamics of both the spike train and the HD angle but decouples the two from each other. We computed the time-reversed HD angle by reversing the order of HD angle values with respect to their timestamps within a particular epoch. We then classified HD cells as those with HD information scores higher than the 99th percentile of the null distribution (>0.2 bits per spike, 85% of putative excitatory cells). We did not apply any HD information threshold to PoSub-FS cells.

Cross-validated HD tuning curve auto-correlograms and cross-correlograms

We obtained HD tuning curve auto-correlogram and cross-correlogram by computing Pearson’s correlation coefficients between the reference tuning curve vector and the second tuning curve vector (from either the same or another cell), which was circularly shifted by 0–359 bins. To minimize the effect of non-HD factors on tuning curves computed from the same epoch, we used a cross-validation procedure whereby the two tuning curves were computed from separate halves of the epoch.

When computing HD tuning curves during the exploration of the triangular open field, we noticed that sometimes the cells’ receptive fields were rotated with respect to the prior square open-field exploration. To correct for this, we first calculated for each HD cell the degree of tuning curve rotation between the two environments via cross-correlation. We then used the average rotation per recording session to circularly shift all triangular open-field tuning curves in the opposite direction by the equivalent amount. This allowed us to compute the true tuning curve correlation between the two environments.

Detection of monosynaptic connections

Spike train cross-correlograms of ±50 ms binned in 0.2-ms windows were convolved with a Gaussian kernel of 4 ms s.d, resulting in a predictor of the baseline rate. At each time bin, the 99th percentile of the cumulative Poisson distribution (at the predicted rate) was used at the statistical threshold for significant detection of outliers from baseline. A putative connection was considered significant when at least two consecutive bins in the cross-correlograms passed the statistical test.

Analysis of HD cell realignment after cue rotation

To estimate the degree of realignment of the HD system following cue rotation, HD cell spike times were binned into population vectors (50 ms window, smoothed in 100 ms s.d. Gaussian windows). Based on cells’ tuning curves from the baseline period of exploration with a stable cue, the population vectors were converted into a Bayesian probabilistic map under the hypothesis that HD cells fire as a Poisson process. The instantaneous internal HDs were taken as the maxima of these probabilistic maps. These estimates faithfully tracked the head orientation of the mouse during the period preceding cue rotation. The degree of realignment of the HD system was calculated as the decoder error—angular difference between the real HD and the decoded HD at each time bin. Because not all cue rotations resulted in HD realignment, we excluded all cue rotation epochs that (1) resulted in less than 45° of mean decoder error in the following 200 s epoch and (2) occurred when the animal was stationary in the preceding epoch (average velocity <2 cm s−1).

To estimate the point at which the HD system remaps, we fitted a sigmoid curve to the decoder error values following each cue rotation using the sigm_fit function (https://www.mathworks.com/matlabcentral/fileexchange/42641-sigm_fit). We defined the beginning and end of realignment epochs as the timestamps corresponding to the values of 0.01 and 0.99 of the normalized sigmoids. We then, for each cell, calculated an HD tuning curve for the remainder of each cue epoch (from the end of realignment to the next cue rotation) and computed the cross-correlation (see above) between HD tuning curves from consecutive epochs. For each cell, the degree of realignment was defined as the tuning curve offset that results in the highest correlation coefficient. The difference in realignment between FS and HD cells was defined as, for each FS cell, the angular difference between its degree of realignment and the average realignment of HD cells in the same epoch.

Classification of sleep states

Sleep scoring was performed using the automated SleepScoreMaster algorithm and TheStateEditor software65,66 (Buzsaki Laboratory, https://github.com/buzsakilab). The wide-band signal was downsampled to 1.25 kHz and used as the LFP signal. Electromyograph (EMG) was computed from correlated high-frequency noise across several channels of the linear probe. Recording sessions with less than 100 s of REM sleep detected (n = 6) were excluded from the analyses involving REM sleep because of low number of samples.

Pairwise spike–rate coupling

Quantification of pairwise spike–rate coupling between cells was quantified using a GLM according to the method described in ref. 49. Spike trains were binned in 100 ms bins and smoothed in 100 ms s.d. Gaussian windows. The population firing rate was calculated by aggregating all spike times from all recorded units in a given recording and processing them in the same manner as single spike trains. Both binned trains were then restricted to either WAKE or REM epoch (see above).

The GLM was fitted using the MATLAB glmfit function. The binned spike train of cell A was modeled as a Poisson process, as a function of both the binned spike train of cell B and the binned population firing rate, using a log link function. The model produced a coefficient of coupling between the spike trains of cells A and B (‘β’), as well as a coefficient for the coupling of cell A to the population firing rate (‘βPOP’). The procedure was repeated by offsetting the spike train of cell A by ±10 s in 100 ms intervals, to yield the equivalent of the spike train cross-correlogram that discounts the coupling of cell A to the local population rate. Because this procedure, unlike Pearson’s correlation, is not symmetric, it was repeated by swapping cell A and cell B and averaging the coupling coefficient values at equivalent offset intervals. A cell pair was removed from the analysis in rare cases when the glmfit function identified the model as ill-fitted or reached the iteration limit.

Visualization and analysis of the ring manifold

For the visualization of the HD manifold during WAKE, HD cell spike times from the whole epoch were binned into population vectors in 200 ms bins, smoothed in 400 ms s.d. Gaussian windows and a square root of each value was taken. Then, nonlinear dimensionality reduction was performed using the Isometric Feature Mapping (Isomap) algorithm39 implemented in the MATLAB Toolbox for Dimensionality Reduction (version 0.8.1b, https://lvdmaaten.github.io/drtoolbox/). The parameters were set to 12 nearest neighbors and three dimensions—the latter to inspect if there is no higher dimensional manifold in the data. Shuffled Isomap embeddings were computed by shifting each cell’s binned spike train in time by a random number of bins.

Internal HD at each time bin was then calculated as a four-quadrant arctangent of the first two Isomap dimensions (range: −180° to 180°). Notably, the ‘virtual HD’ generated this way has arbitrary directionality (clockwise/counterclockwise) and an arbitrary point of origin. To align it, Isomap directionality was established by computing the Isomap error as the difference between real HD values and (1) virtual HD values and (2) additive inverse of virtual HD values, and selecting the directionality for which the distribution of angular differences between real and virtual HD had smaller circular variance. Internal HD tuning curves were then computed as the ratio between histograms of spike count and total time spent in each virtual HD bin (1°) and smoothed with a Gaussian kernel of 3° s.d. Real HD tuning curves were computed by downsampling the real HD into 200 ms bins and applying the same procedure as above. To correct for the arbitrary point of origin, HD tuning curve cross-correlations (see above) were computed between the real and virtual HD tuning curve of each HD cell, and the mean offset of maximum correlation was then used to circularly shift all tuning curve vectors by the equivalent number of bins. Although this procedure (as well as Isomap mapping in general) was dependent on the real HD tuning of HD cells, it was independent of FS cell tuning. Control HD tuning curves were computed by time-reversing the Isomap angle (see above).

For the comparison between Isomap HD tuning curves during WAKE and REM, population vectors across WAKE and REM were computed in the same manner as above. WAKE population vectors were then randomly downsampled to equal in number to the REM population vectors. Isomap algorithm was then run on both WAKE and REM population vectors together. Virtual HD was computed in the same way as above. HD tuning curves were computed in bins of 6° and smoothed with a Gaussian kernel of 6° s.d for both real and internal HD. Larger bin size was chosen due to sometimes uneven occupancy of virtual HD during REM.

Fourier analysis of HD tuning curves

To decompose tuning curves into Fourier series, we projected the tuning curves onto a basis of sine and cosine functions whose frequencies were the harmonic of the unit circle, that is, from the fundamental frequency (period of 360°) to the highest possible frequency (2°, the inverse of the Nyquist frequency as tuning curves were computed in 1° bins). The power, or Fourier coefficient, at a particular frequency was defined as the root mean square of the projection values onto the sine and cosine basis at that frequency. Similarly, the phase was defined as the arctangent of the projections. The validity of the projection was verified by checking that the sum of squared Fourier coefficients is equal to the variance of the tuning curves (Parseval’s identity), which was indeed the case (Extended Data Fig. 3e). Because higher Fourier components likely represent noise fluctuations in the tuning curves, we focused our analysis on the relative power of the first ten Fourier components, normalizing their individual power values to the sum of their power.

Kullback–Leibler (KL) divergence was used as a measure to assess the similarity between the individual Fourier spectra and the population means. While KL divergence is regularly used to compare probability distributions, we deemed it appropriate to apply it to normalized Fourier spectra as they were mathematically indistinguishable from probability distributions. We thus computed the KL divergence between the spectrum of an individual neuron σi (k) (with \(k[1.10]\) the angular frequencies) and the average Fourier spectrum of a population Si (k) as follows:

$$_}(_||S)=\sum __(k)\log \left(\frac_(k)}\right)$$

While KL divergence is not symmetrical, that is, KL divergence (DKL) (σi||S) does not equal DKL (S||σi), we always applied it in the same direction, that is, DKL (σi||S).

To further validate our results, we repeated the analysis after applying a higher velocity threshold (10 cm s−1) as well as after the exclusion of HD cells with multiple receptive fields or noisy unit clusters. To detect multipeaked cells, multipeak score was computed for each cell, defined as the mean firing rate of a normalized tuning curve at the angular values outside of its primary peak (±width of the curve from the angle of its maximum firing rate). To identify noisy unit clusters, for each cell, a spike train auto-correlogram was computed at ±100 ms (1 ms bins) and the cluster contamination score was defined as the ratio between the value at 0.1 ms and the highest value. Noisy clusters were defined as those with cluster contamination scores above 0.2 (11% of all cells in the dataset).

Isomap analysis of HD tuning curve auto-correlograms

Cross-validated auto-correlograms of each cell’s HD tuning curve were computed as described above. Then, nonlinear dimensionality reduction was performed using the Isomap algorithm39. The parameters were set to 12 nearest neighbors and two dimensions. When mapping the first three Fourier components onto the resulting embedding, we normalized their power to the total sum of their powers.

Anatomical tract tracing

VGAT-IRES-Cre mice were injected with an adeno-associated virus vector (AAV 2/9 CAG-Flex-EGFP, titer: 4 to 5 × 1012 GC per ml; Neurophotonics; 500 μl per injection site) bilaterally into the TRN as described above. Four weeks after injections, animals were perfused transcradially with 4% PFA in phosphate-buffered saline and their brains were then cut coronally in 40 μm sections with a freezing microtome. The sections were counterstained with blue NeuroTrace (Thermo Fisher Scientific; 1:200 dilution) and mounted on slides with coronal z-stacks of the sections containing the rostral thalamus that were taken with a Leica SP-8 confocal microscope at ×10 magnification, using the same settings for all sections. GFP signal was acquired using the 473 nm excitation laser line. Z-projections of each stack were then obtained using ImageJ (version 1.52). Quantification of anterograde tracing was done in ImageJ. The images were converted to grayscale, and rectangular regions of interest (ROI) were defined within each thalamic nucleus. Average pixel intensity per ROI was then calculated using the ‘Measure’ function.

Gain modulation analysis: experiment

Epochs of optogenetic stimulation (light ON) consisted of time periods when the laser was switched on (240 pulses of 1-s duration, 240 s per session). Control epochs (light OFF) were defined as time periods in between light pulses (240 periods of 4-s duration, 960 s per session). Light ON and light OFF tuning curves were computed from these periods. For analysis of tuning curve width, HD tuning curves were then computed in bins of 1° and smoothed with a Gaussian kernel of 3° s.d.

For analysis of additive and multiplicative gain, to preserve the independence of individual angular bins, HD tuning curves were instead computed in bins of 6° with no Gaussian smoothing applied. Light ON (that is, high gain) and light OFF (that is, baseline) tuning curves for each cell were then normalized by dividing them by the maximum value of the baseline tuning curve. To calculate the additive and multiplicative factors for each cell, a linear fit between high gain and baseline tuning curve vectors was then obtained using the MATLAB polyfit function. The slope of the resulting linear fit and its Y intercept were then taken as multiplicative and additive factors, respectively.

Because estimation of gain factors is not possible for cells that are not gain-modulated and/or their baseline and high gain tuning curves are not significantly correlated, we excluded cells in the ArchT group that did not show statistically significant modulation by light and did not show a statistically significant correlation between baseline and high gain tuning curves. Significant modulation by light was computed by, for each cell, calculating the average firing rate during each baseline and each high gain epoch, and subsequently comparing the two sets of values with a Wilcoxon signed-rank test with a significance threshold of 0.05. Significance of the tuning curve correlation was determined by the Pearson correlation coefficient (significance threshold of 0.05). This led to the exclusion of 1 of 47 PoSub-FS cells and no cells from other populations.

Gain modulation analysis: simulations

The relationship between tuning curve correlations and additive/multiplicative factors was first explored in simple simulations involving monotonically increasing vectors (correlation coefficient (r) between vectors = 1). Noise was added separately to each vector by offsetting each point by a random value drawn from a Gaussian distribution. We varied the s.d. of the noise distribution to alter the value of the Pearson correlation coefficient between the two vectors. We then performed a linear regression between the two vectors using the MATLAB polyfit function and computed the additive and multiplicative factors as outlined above. The procedure was performed 1,000 times for each value of s.d. of the noise distribution. This enabled us to quantify how the apparent gain factor estimations change as a function of vector correlation.

On the basis of above simple simulations, it became apparent that our methods might have underestimated the amount of multiplicative gain in PoSub-FS cells. Consequently, we next sought to determine whether considerably low tuning curve correlations in the PoSub-FS cells preclude the detection of multiplicative gain via linear regression. To that end, we applied a method analogous to the one described above to real PoSub-FS cell tuning curves. We used the baseline tuning curves from the ArchT group and computed simulated high-gain tuning curves by applying Gaussian noise independently to each tuning curve until we obtained a similar distribution of tuning curve correlations to the one observed experimentally. We subsequently applied either multiplicative gain (multiplying each point on the tuning curve by 1.2) or additive gain (by adding 20% of the maximum tuning curve value to each point on the tuning curve) before applying Gaussian noise. These simulations allowed us to establish that multiplicative gain is indeed detectable in tuning curves with correlation values similar to those observed experimentally in PoSub-FS cells and that our experimental results in PoSub-FS cells are best explained by significant additive gain.

Model fitting and theory

First, we derived network parameters based on the statistics of observed PoSub-FS cell tuning curves. Specifically, we characterized the connectivity from ADN- and PoSub-HD cells to PoSub-FS cells using mean and s.d. of connection weights. To differentiate the contributions of ADN and PoSub inputs, we optimized these parameters. The goal was to closely replicate the firing rate statistics and Fourier tuning profiles of the PoSub-FS cells. The optimization process is governed by the following three constraints, formulated as a loss function minimized through a gradient descent procedure: (1) matching the firing rates between simulated and recorded PoSub-FS cells, (2) equating the variance of the firing rates and (3) ensuring similarity in tuning curve shapes by aligning the variances of the actual Fourier coefficients.

Next, we theoretically derived the asymptotic behavior of the proportion of simulated cells with a certain Fourier power. We showed that this proportion becomes independent of the weight mean and s.d. in the limit of low or high variance relative to the mean.

Finally, we aimed to determine the actual distribution of output neuron tuning curves. We show that tuning curves of randomly connected linear neurons, with normally distributed weights W, are constrained to lie in the subspace spanned by the singular components of their input tuning curves. Furthermore, they are distributed in this subspace according to a multivariate normal distribution with diagonal covariance. As a result, Fourier components of an output tuning curve are independent of each other, as observed in the population of PoSub-FS cells.

Supplementary Methods provide further details of the model fitting procedure, the theoretical distribution of FS tuning curves and a generalization to other synaptic weight distributions.

Data analysis and statistics

Analyses presented in Fig. 2 and Extended Data Fig. 5 were conducted using software custom-written in Python 3.9.7 with the following libraries: Scipy (version 1.7.3), Matplotlib (3.5.1), Numpy (1.20.3), Uncertainties (3.1.17) and Hdf5storage (0.1.18). All other analyses were conducted using MATLAB R2020b (Mathworks) with TStoolbox 2.0 and Circular Statistics Toolbox version 1.21 (ref. 67). Statistical comparisons were performed with nonparametric tests (Mann–Whitney U test, Wilcoxon signed-rank test) or analysis of variance (ANOVA) with multiple comparisons, where applicable. For ANOVA, data distribution was assumed to be normal but this was not formally tested. All statistical tests were two-tailed. No statistical methods were used to predetermine sample sizes, but our sample sizes are similar to those reported in previous publications34,48,49. Data collection and analysis were not performed blind to the conditions of the experiments. Data collection was not randomized as individual experiments were performed in sequence.

Reporting summary

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