Friday, January 21, 2022
HomeNatureToroidal topology of inhabitants exercise in grid cells

Toroidal topology of inhabitants exercise in grid cells


Rats

The information have been collected from three experimentally naive male Lengthy Evans rats (Rats Q, R and S, 300–500 g at time of implantation). The rats have been group-housed with three to eight of their male littermates earlier than surgical procedure and have been singly housed in massive Plexiglas cages (45 × 44 × 30 cm) thereafter. They have been saved on a 12-h mild–12-h darkish schedule, with strict management of humidity and temperature. All procedures have 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. Protocols have been authorized by the Norwegian Meals Security Authority (FOTS ID 18011 and 18013).

Electrode implantation and surgical procedure

The rats have been implanted with Neuropixels silicon probes25,26 concentrating on the MEC–parasubiculum (PaS) area. Two rats have been implanted bilaterally with prototype Neuropixels ‘section 3A’ single-shank probes and with one probe concentrating on MEC–PaS in every hemisphere; the third rat was implanted with a prototype Neuropixels 2.0 multi-shank probe within the left hemisphere. Probes have been inserted at an angle of 25° from posterior to anterior within the sagittal aircraft. Implantation coordinates have been AP 0.05–0.3 mm anterior to the sinus and 4.2–4.7 mm lateral to the midline. The probes have been inserted to a depth of 4,200–6,000 µm. The implant was secured with dental cement. The detailed implantation process has been described elsewhere6,26. After surgical procedure, the rats have been left to get well for roughly 3 h earlier than starting recording. Postoperative analgesia (meloxicam and buprenorphine) was administered in the course of the surgical restoration interval.

Recording procedures

The main points of the Neuropixels {hardware} system and the procedures for freely transferring recordings have been described beforehand. In short, electrophysiological alerts have been amplified with a achieve of 500 (for section 3A probes) or 80 (for two.0 probes), low-pass-filtered at 300 Hz (section 3A) or 0.5 Hz (2.0), high-pass-filtered at 10 kHz, after which digitized at 30 kHz (all steps carried out by the probe’s on-board circuitry). The digitized alerts have been multiplexed by an implant-mounted ‘headstage’ circuit board and have been transmitted alongside a light-weight 5-m tether cable, made utilizing both micro-coaxial (section 3A) or twisted pair (2.0) wiring.

Three-dimensional movement seize (OptiTrack Flex 13 cameras and Motive recording software program) was used to trace the rat’s head place and orientation, by attaching a set of 5 retroreflective markers to implant throughout recordings. The 3D marker positions have been projected onto the horizontal aircraft to yield the rat’s 2D place and head route. An Arduino microcontroller was used to generate digital pulses, which have been despatched to the Neuropixels acquisition system (by way of direct TTL enter) and the OptiTrack system (by way of infra-red LEDs), to allow exact temporal alignment of the recorded information streams.

Behavioural procedures

Information have been obtained from 4 recording classes carried out inside the first 72 h after restoration from surgical procedure. The recordings have been carried out whereas the rats engaged in three behavioural paradigms, every in a special enviornment inside the similar room. Plentiful distal visible and sonic cues have been obtainable to the rat. On every day of recording, the rat remained constantly related to the recording equipment throughout the assorted behavioural classes that have been carried out. Sometimes it was essential to take away twists that had gathered within the Neuropixels tether cable. In such instances, the continuing behavioural activity was paused whereas the experimenter gently turned the rat to take away the twists. Throughout pre-surgical coaching, the rats have been food-restricted, sustaining their weight at a minimal of 90% of their free-feeding physique weight. Meals was typically eliminated 12–18 h earlier than every coaching session. Meals restriction was not used on the time of recording.

Open-field foraging activity

Rats foraged for randomly scattered meals crumbs (corn puffs) in a sq. open-field (OF) enviornment of dimension 1.5 × 1.5 m, with black flooring and enclosed by partitions of peak 50 cm. A big white cue card was affixed to one of many enviornment partitions (peak similar because the wall; width 41 cm; horizontal placement on the center of the wall). On the time of the surgical procedure, every rat was extremely accustomed to the surroundings and activity (10–20 coaching classes lasting not less than 20 min every).

Wagon-wheel monitor foraging activity

The wagon-wheel (WW) monitor activity was designed to operate as a 1D model of the 2D OF foraging activity. The monitor’s geometry comprised an elevated round monitor with two perpendicular cross-linking arms spanning the circle’s diameter. The monitor was 10 cm large and was bounded on either side by a 1-cm-high lip. Every part of the monitor was fitted with a reward level, positioned midway between the 2 nearest junctions, within the centre of the part. Every reward level consisted of an elevated properly that may very well be remotely full of chocolate milk by way of connected tubing. To encourage foraging behaviour, a pseudorandom subset of the wells (between one and 4 of the eight wells) was stuffed at a given time, and the rat was allowed to discover the total maze freely and constantly. Wells have been refilled as crucial when the rat consumed rewards. Every rat was educated to excessive efficiency on the foraging activity earlier than the surgical procedure (accumulating not less than 30 rewards inside a 30-minute session). Coaching to this stage of efficiency took 5–10 half-hour classes.

Pure sleep

For sleep classes, the rat was positioned in a black acrylic ‘sleep field’ with a 40 × 40-cm sq. base and 80-cm-high partitions. The black coating of partitions was clear to infrared, which allowed the 3D movement seize to trace the rat by means of the partitions. The underside of the sleep field was lined with towels, and the rat had free entry to water. Throughout recording classes within the sleep field, the primary room mild was switched on and pink noise was performed by means of the pc audio system to attenuate disturbing background sounds. Sleep classes usually lasted 2–3 h, however have been aborted prematurely if the rat appeared extremely alert and unlikely to sleep.

Spike sorting and single-unit choice

Spike sorting was carried out with KiloSort 2.526. In short, the algorithm consists of three principal levels: (1) a raw-data alignment process that detects and corrects for shifts within the vertical place of the Neuropixels probe shank relative to the encompassing tissue; (2) an iterative template-matching process that makes use of low-rank, variable-amplitude waveform templates to extract and classify single-unit spikes; and (3) a curation process which detects acceptable template merging and splitting operations based mostly on spike practice auto- and cross-correlograms. Some customizations have been made to the usual KiloSort 2.5 methodology to enhance its efficiency on recordings from the MEC–PaS area, the place there’s a notably excessive spatiotemporal overlap of spike waveforms owing to the excessive density of cells. Subsequently, the utmost variety of spikes extracted per batch in step 1 above was elevated, as was the variety of template-matching iterations in step 2. To enhance the separation between cells with very similar-looking waveforms, the higher restrict on template similarity was raised from 0.9 to 0.975 in step 2 and to 1.0 on step 3, whereas supervising manually all merge and break up operations from step 3, utilizing a custom-made GUI working in MATLAB. The handbook supervision ensured that Kilosort 2.5 didn’t mechanically merge pairs of models with a dip within the cross-correlogram, which in our information was typically as a consequence of out-of-phase spatial tuning. The merge and break up operations have been repeated a number of occasions to make sure the perfect separation between single models.

Single models have been discarded if greater than 1% of their interspike interval distribution consisted of intervals lower than 2 ms. In additions, models have been excluded if they’d a imply spike charge of lower than 0.05 Hz or larger than 10 Hz (calculated throughout the total recording period).

Single-unit spike waveforms

Throughout spike sorting, Kilosort assigned every unit with a 2 ms spike waveform template on every recording channel. To calculate a consultant single waveform for every unit, the peak-to-peak amplitude of the template was calculated on each channel, and the templates from the three highest-amplitude channels have been averaged to generate the consultant spike waveform. To calculate spike width, a unit’s consultant waveform was finely interpolated (from 61 to 1,000 factors) utilizing a cubic spline operate. Spike width was outlined because the time distinction between the waveform’s adverse peak (to which the waveform was aligned by Kilosort), and the next constructive peak.

Spatial place and route tuning

Throughout awake foraging classes within the OF enviornment or wagon-wheel monitor, solely time epochs wherein the rat was transferring at a velocity above 2.5 cm s−1 have been used for spatial or toroidal analyses. To generate 2D charge maps for the OF enviornment, place estimates have been binned right into a sq. grid of three × 3-cm bins. The spike charge in every place bin was calculated because the variety of spikes recorded within the bin, divided by the point the rat spent within the bin. To interpolate the values of unvisited bins, two auxiliary matrices have been used, M1 and M2, setting visited bins equal to the worth of the unique charge map in M1 and to 1 in M2, and setting unvisited bins to zero in each. One iteration of the image-processing ‘closing’ operation was then carried out (binary dilation adopted by erosion, filling out a subset of the non-visited bins) on M2, utilizing a disk-shaped structuring component, first padding the matrix border by one bin. Each matrices have been then spatially smoothed with a Gaussian kernel of smoothing width 2.75 bins. Lastly, the speed map was obtained by dividing M1 by M2. Fee-map spatial autocorrelograms and grid scores have been calculated as described beforehand28. The selectivity of every cell’s place tuning was quantified by computing its spatial info content material42, measured in bits per spike (see ‘Info content material’).

Head-direction tuning curves have been calculated by binning the head-direction estimates into 6° bins. The spike charge in every angular bin was calculated because the variety of spikes recorded within the bin divided by the point that the rat spent within the bin. The resultant tuning curve was smoothed with a Gaussian kernel with σ = 2 bins, with the ends of the tuning curve wrapped collectively. The selectivity of head-direction tuning was quantified utilizing the imply vector size (MVL) of the tuning curve. This was calculated in line with:

$${rm{M}}{rm{V}}{rm{L}}=,frac{|{sum }_{j=1}^{M}{{bf{f}}}_{j}exp (i{{boldsymbol{alpha }}}_{j})|}{{sum }_{j=1}^{M}{{bf{f}}}_{j}},$$

the place vector f represents the tuning curve values (firing charges), vector α represents the corresponding angles, M is the variety of tuning curve values, and |∙| represents absolutely the worth of the enclosed time period.

Grid module classification

A novel methodology was carried out to detect populations of cells equivalent to grid modules by discovering clusters of cells that expressed related spatially periodic exercise within the open discipline (Prolonged Information Fig. 2). Opposite to earlier clustering-based strategies for grid modules3, this method makes no assumptions in regards to the particular geometry of the grid sample, thus making it much less prone to the detrimental results of geometric distortions corresponding to ellipticity3,30.

For every MEC–PaS cell in a given recording, a coarse-resolution charge map of the OF session was constructed, utilizing a grid of 10 × 10-cm bins, with no smoothing throughout bins. The 2D autocorrelogram of this charge map was calculated, and the central peak was eliminated by excluding all bins situated lower than 30 cm from the autocorrelogram centre. Bins situated greater than 100 cm from the autocorrelogram centre have been additionally excluded. The autocorrelograms for all cells have been subsequently transformed into column vectors, z-standardized, then concatenated to type a matrix with spatial bins as rows and cells as columns. The nonlinear dimensionality discount algorithm UMAP43,44 was then utilized to this matrix, yielding a two-dimensional level cloud wherein every information level represented the autocorrelogram of 1 cell (Prolonged Information Fig. 2a–d; UMAP hyperparameters: ‘metric’=‘manhattan’, ‘n_neighbors’=5, ‘min_dist’=0.05, ‘init’=‘spectral’). Within the resultant 2D level cloud, cells with small absolute variations between their autocorrelogram values have been situated close to to 1 one other. The purpose cloud was partitioned into clusters utilizing the DBSCAN clustering algorithm (MATLAB operate ‘dbscan’, minimal 30 factors per cluster, eta = 0.6–1.0). In each recording, the biggest cluster was primarily composed of cells that both lacked sturdy spatial selectivity or have been spatially selective however with out clear periodicity. All remaining clusters contained cells with excessive grid scores, and with related grid spacing and orientation (Prolonged Information Fig. 2a–d); cluster membership was subsequently used as the idea for grid module classification. In a single recording (rat ‘R’ day 1), two clusters have been recognized that had related common grid spacing and orientation (labelled as ‘R1a’ and ‘R1b’ in Prolonged Information Fig. 2a–d), suggesting that they represented the identical grid module. R1b appeared to comprise cells with increased variability within the within-field firing charges of the spatial charge maps, accompanied by extra irregularities within the autocorrelograms. These two clusters have been merged collectively in subsequent evaluation (wherein the resultant cluster is named ‘R1’).

A subset of the cells that have been assigned to grid module clusters by the above process have been tuned to each location and head route (conjunctive grid × route cells). These cells, which have been outlined as having a head-direction tuning curve with imply vector size above 0.3, have been discarded from additional evaluation.

Classification of sleep states

SWS and REM durations have been recognized on the idea of a mix of behavioural and neural exercise, following beforehand described approaches6,45,46. First, sleep durations have been outlined as durations of sustained immobility (longer than 120 s with a locomotion velocity of lower than 1 cm s−1 and head angular velocity of lower than 6° s−1). Qualifying durations have been then subclassified into SWS and REM on the idea of the amplitude of delta- and theta-band rhythmic exercise within the recorded MEC–PaS cells. Spike occasions for every cell have been binned at a decision of 10 ms and the resultant spike counts have been binarized, such that ‘0’ indicated the absence of spikes and ‘1’ indicated a number of spikes. The binarized spike counts have been then summed throughout all cells (Prolonged Information Fig. 9A). The rhythmicity of this aggregated firing charge with respect to delta (1–4 Hz) and theta (5–10 Hz) frequency bands was quantified by making use of a zero-phase, fourth-order Butterworth band-pass filter, then calculating the amplitude from absolutely the worth of the Hilbert remodel of the filtered sign, which was smoothed utilizing a Gaussian kernel with σ = 5 s after which standardized (‘z-scored’). The ratio of the amplitudes of theta and delta exercise was therefore calculated (theta/delta ratio, ‘TDR’). Intervals wherein TDR remained above 5.0 for not less than 20 s have been categorized as REM; durations wherein TDR remained beneath 2.0 for not less than 20 s have been categorized as SWS (Prolonged Information Fig. 9B).

Spectral evaluation was carried out on 10-ms-binned multi-unit exercise utilizing the multi-tapered Fourier remodel, carried out by the Chronux toolbox (http://chronux.org/, operate ‘mtspectrumsegc’). Non-overlapping 5-second home windows have been used, with a frequency bandwidth of 0.5 Hz and the utmost variety of tapers.

Visualization of toroidal manifold

For every module of grid cells, spike occasions of co-recorded cells within the OF have been binned for every cell at a decision of 10 ms, and the binned spike counts have been convolved with a Gaussian filter with σ = 50 ms. Time bins wherein the rat’s velocity was beneath 2.5 cm s−1 have been then discarded. To account for variability of common firing charges throughout cells, the smoothed firing charge of every cell was z-scored. For computational causes, the time bins have been downsampled, taking each twenty fifth time bin (equating to 250-ms intervals between chosen samples). Collectively, the downsampled firing charges of the total inhabitants of cells fashioned a matrix with time bins in rows and cells in columns. PCA was utilized to this matrix (treating time bins as observations and cells as variables), and the primary six principal parts have been retained (Prolonged Information Figs. 3Aa–c, 4a–d). UMAP43,47 was then run on these six principal parts (with time bins as observations and principal parts as variables). The hyperparameters for UMAP have been: ‘n_dims’=3, ‘metric’=‘cosine’, ‘n_neighbours’=5000, ‘min_dist’=0.8 and ‘init’=‘spectral’.

For visualizing the toroidal manifold throughout WW, smoothed firing charges have been first calculated by the identical process described above for OF. Subsequently, to permit comparability of the toroidal manifold between OF and WW, the identical PCA and UMAP transformations calculated for the OF information have been re-applied to the WW information, by supplying the fitted OF UMAP transformation because the argument ‘template_file’ to the ‘run_umap’ operate within the MATLAB implementation47.

Preprocessing of inhabitants exercise

Every topological evaluation was based mostly on the exercise of a single module of grid cells, throughout a single experimental situation in a single recording session. Topological evaluation of multi-module and conjunctive grid × route cell exercise was not thought-about as we count on such information to exhibit higher-dimensional topological construction requiring the next variety of cells27. The experimental circumstances have been: open-field foraging (OF), wagon-wheel monitor foraging (WW), slow-wave sleep (SWS), and fast eye-movement sleep (REM). Sleep epochs of the identical sort have been collected from throughout the recording and concatenated for evaluation functions. Equally, in a single case (rat ‘S’), two WW activity classes have been concatenated to extend the pattern dimension.

In complete there have been 27 mixtures of module (Q1, Q2, R1, R2, R3, S1) and experimental situation (OF day 1, OF day 2, WW, REM, SWS).

Preprocessing of spike trains started by computing delta features centred on the spike occasions (valued 1 at time of firing; 0 in any other case), and convolving these temporally with a Gaussian kernel with σ = 50 ms (OF, WW and REM) or 25 ms (SWS). Samples of the smoothed firing charges of all cells (‘inhabitants exercise vectors’) have been then computed at 50-ms intervals. The awake states have been additional refined by excluding vectors which originated from time durations when the rat’s velocity was beneath 2.5 cm s−1.

Computing the persistent cohomology of a degree cloud is computationally costly and could also be delicate to outliers (for instance, spurious factors breaking the topology of nearly all of factors within the level cloud). Because of this, it is not uncommon to preprocess the info by downsampling and dimension-reducing the purpose cloud. The identical preprocessing process was used for all datasets within the current research.

First, the info factors have been downsampled by retaining the 15,000 most energetic inhabitants exercise vectors (as measured by the imply inhabitants firing charge). Throughout SWS, this choice criterion had the consequence of mechanically discarding inhabitants exercise vectors throughout down-states, when neural exercise is near-silent. As noise is inherently extra prevalent and cosine distances much less dependable in high-dimensional areas (“the curse of dimensionality”)48, dimensionality-reduction and a normalization of distances have been subsequently carried out. The diminished level cloud was z-scored and projected to its six first principal parts, thus lowering noise whereas retaining a lot of the variance (see Prolonged Information Fig. 4a). This was supported by the shortage of grid construction and the clear drop in defined deviance after six parts (see Prolonged Information Fig. 4b, c). The defined deviance was computed by becoming a GLM mannequin to every element individually, utilizing the spatial coordinates as covariate, suggesting that the upper parts are much less spatially modulated and presumably higher described by different (unknown) covariates. In step with this, the toroidal construction was most clearly detected within the barcodes when evaluating the ratio of the lifetimes of the 2 most persistent H1 bars versus the third longest-lived H1 bar for the barcodes obtained when utilizing completely different numbers of parts within the evaluation (see Prolonged Information Fig. 4d). These analyses each indicated that dimensionality discount was required to firmly exhibit the toroidal topology within the grid cells. The empirical findings are supported theoretically; see ‘Theoretical rationalization of the six-dimensionality proposed by PCA’ in Supplementary Strategies.

To additional simplify the low-dimensional level cloud, a special downsampling approach was launched, based mostly on a point-cloud density technique motivated by a topological denoising approach launched beforehand49 and a fuzzy topological illustration utilized in UMAP43,50. Elements of the open-source implementation of the latter have been copied on this computation. This method consisted of assigning, for every level, a neighbourhood energy to its okay nearest neighbours, and subsequently sampling factors that symbolize essentially the most tight-knit neighbourhoods of the purpose cloud in an iterative method. First, we outlined ({m}_{i,{i}_{j}}^{{prime} }=exp left(-frac{{d}_{i,{i}_{j}}}{{sigma }_{i}}proper),) the place ({d}_{i,{i}_{j}}) is the cosine distance between level xi and its j-th nearest neighbour and σi is chosen to make (mathop{sum }limits_{j=1}^{okay}{m}_{i,{i}_{j}}^{{prime} }={log }_{2}okay,)utilizing okay = 1,500. The neighbourhood energy was then obtained by symmetrizing: ({m}_{i,{i}_{j}}={m}_{i,{i}_{j}}^{{prime} }+,{m}_{j,{j}_{i}}^{{prime} }-,{m}_{i,{i}_{j}}^{{prime} }cdot {m}_{j,{j}_{i}}^{{prime} }). Lastly, the purpose cloud was diminished to 1,200 factors by iteratively drawing the i-th level as: (mathop{max }limits_{{x}_{i}}sum _{jin tilde{I}}{m}_{j,{j}_{i}},) the place Ĩ denotes the indices of the factors not already sampled. In different phrases, for every iteration, the sampled level is the one with the strongest common membership of the neighbourhoods of the remaining factors.

To compute the persistent cohomology of the downsampled level cloud, the neighbourhood strengths have been first computed for the diminished level cloud (utilizing okay = 800) and its adverse logarithm was taken, acquiring a distance matrix. This matrix was then given as enter to the Ripser implementation51,52 of persistent cohomology, returning a barcode. Briefly, the barcode gave an estimate of the topology of the fuzzy topological illustration of the six principal parts of the grid-cell inhabitants exercise. Thus, in essence, step one of UMAP was utilized earlier than describing the ensuing illustration with persistent cohomology, as a substitute of utilizing it to challenge every level of the purpose cloud to a illustration of user-specified dimensionality for visualization (Prolonged Information Fig. 3Ad, e). This provides a extra direct and steady quantification of the worldwide information construction, with out having to decide on an initialization53 or optimize a lower-dimensional illustration.

Persistent cohomology

Persistent cohomology, a instrument in topological information evaluation, was used to characterize the manifold assumed to underlie the info. This has clear ties with persistent homology and the primary consequence (the barcode) is similar, thus the 2 phrases are sometimes used interchangeably. Persistent cohomology was chosen as a result of the computation is (to our information) sooner and is required to acquire cocycle representatives, that are essential to carry out decoding (see ‘Cohomological decoding’). Persistent (co-)homology has beforehand been profitable in analysing neural information, describing the ring topology of head route cell exercise22,23,24, the spherical illustration of inhabitants exercise in major visible cortex54, and the exercise of place cells55,56,57,58.

The final define of the algorithm is as follows. Every level within the cloud is changed by a ball of infinitesimal radius, and the balls are regularly expanded in unison. Taking the union of balls at a given radius leads to an area with holes of various dimensions. The vary of radii for which every gap is detected is tracked; that is known as the ‘lifetime’ of the outlet and is represented by the size of a bar. The totality of bars is known as the barcode.

The software program package deal Ripser51,52 was used for all computations of persistent cohomology. Ripser computes the persistent cohomology of ‘Vietoris-Rips complexes’ (which approximate the union of balls for various radii), constructed based mostly on the enter distance matrix and a selection of coefficients (in our case, 47-coefficients), and outputs the barcode and cocycle representatives for all bars. The prime 47 was chosen as homology and cohomology coincide on this case and as it’s unlikely that this divides the torsion subgroup of the homology of the area. Torsion could point out, for instance, orientability of a manifold and in selecting 47 as our prime, we disregard all however 47-torsion. Testing with different primes (for instance, 43) gave related outcomes (information not proven) and the Betti numbers stayed the identical no matter selection of prime.

To confirm that the lifetimes of outstanding bars within the barcodes have been past probability, shuffled distributions have been generated for the persistence lifetimes in every dimension. In every shuffling, the spike practice of every cell was shifted independently in time by rolling the firing charge arrays a random size between 0 and the size of the session. The identical preprocessing and persistence evaluation have been then carried out on the shifted spike trains as for the unshuffled information. This was carried out 1,000 occasions, and every time a barcode was obtained. The barcodes have been concatenated for all shuffles and the utmost lifetime was discovered for every dimension. This lifetime served as a significance criterion for the bar lifetimes. It’s famous, nevertheless, that this can be a heuristic and that statistics of barcodes are nonetheless not properly established.

Cohomological decoding

As there are different areas with related barcodes as for a torus, the outcomes recognized by the barcode have been additional investigated, utilizing the ‘cohomological decoding’ process launched beforehand59 to calculate a toroidal parametrization of the purpose clouds of inhabitants exercise. This assigns to every level corresponding positions on every of the 2 round options recognized by the 1D bars with the longest lifetime, leading to coordinates that additional characterize the underlying form of the info.

Cohomological decoding is motivated by the statement that the 1D cohomology (with integer coefficients) of a topological area X is equal to the set of homotopy-equivalent lessons of steady maps from X to the circle (S1)60; that’s:

$${H}^{1}(X;{mathbb{Z}})cong [X,{S}^{1}].$$

This subsequently signifies that for every 1D bar current at a given radius, there exists a corresponding steady map from the Vietoris-Rips complicated of that radius to the circle. Thus, we could first use persistent cohomology to detect which components symbolize significant (long-lived) options of the info and select a radius for which these options exist. Because the vertices of the Vietoris-Rips complicated are factors within the level cloud, the round values of the corresponding maps on the vertices describe round coordinates of the info.

Within the current case, persistent cohomology was first utilized to the grid-cell inhabitants exercise and X was recognized because the Vietoris-Rips complicated for which the 2 longest-lived one-dimensional bars within the barcode (representing every of the 2 circles of the torus) existed. To outline the specified toroidal coordinates on a website that was as massive as potential, we selected the complicated given on the scale of the delivery plus 0.99 occasions the lifetime of the second longest-lived one-dimensional bar within the barcode22,59,61. Subsequent, the cocycle representatives (given by the persistent cohomology implementation of Ripser51,52) of every of the chosen 1D bars outlined 47-values for every of the sides within the complicated. These edge values have been then lifted to integer coefficients and subsequently smoothed by minimizing the sum over all edges (utilizing the scipy implementation ‘lsmr’). The values on the vertices (factors) of every edge adopted from the sting values and gave the round parametrizations of the purpose cloud. The product of the 2 parametrizations thus offered a mapping from the neural exercise to the two-dimensional torus—that’s, giving a toroidal coordinatization (decoding) of the info.

As persistent cohomology was computed for a diminished dataset of 1,200 factors and subsequently round parametrizations have been obtained just for this level cloud, every parametrization was interpolated to the inhabitants exercise from the remainder of the session(s). First, the 1,200 toroidal coordinates have been weighted by the normalized (‘z-scored’) firing charges of the cells at these time factors, acquiring a distribution of the coordinates for every grid cell. The decoded toroidal coordinates have been then computed by discovering the mass centre of the summed distributions, weighted by the inhabitants exercise vector to be decoded. These exercise vectors have been calculated by first making use of a Gaussian smoothing kernel of 15-ms normal deviation to delta features centred on spike occasions, sampling at 10-ms intervals after which z-scoring the exercise of every cell independently. Time intervals that contained no spikes from any cell have been subsequently excluded. When decoding was used to evaluate or examine the tuning properties of single cells (for instance, comparability of toroidal versus spatial description), the coordinates have been computed utilizing the weighted sum of the distributions of the opposite cells; that’s, the contribution of the cell to be assessed or in contrast was eliminated. When evaluating preservation of toroidal tuning throughout two classes, coordinates have been interpolated both utilizing the toroidal parametrization in every session independently (‘Separate’) or utilizing the identical toroidal parametrization in each classes (‘Frequent’).

Toroidal charge map visualization

For visualization, toroidal firing charge maps have been calculated in the identical means because the bodily area covariate (see ‘Spatial place and route tuning’), first binning the toroidal floor right into a sq. grid of seven.2° × 7.2° bins and computing the common spike charge in every place bin. Nonetheless, for toroidal maps, it was crucial to handle the 60° angle between the toroidal axes earlier than smoothing. After binning the toroidal coordinates, the speed map was ‘straightened’ by shifting the bins alongside the x axis (‘horizontally’) the size of (y mod 2)/2 bins, the place y is the vertical enumeration of the given bin. Copies of the speed map have been then tiled in a three-by-three sq. (much like Prolonged Information Fig. 5d), earlier than making use of the closing and smoothing operations as for the spatial firing charge map. The one toroidal charge map was lastly recovered by reducing out the centre tile, rotating it 90° and defining 15° shear angles alongside each the x and the y axis to right for the 60° offset between them.

Comparability of spatial periodicity

Variations in grid periodicity between OF and WW environments have been quantified for a given cell by evaluating the grid scores within the two behavioural circumstances. Two various strategies have been used to generate the spatial autocorrelograms for this comparability: (1) evaluating the autocorrelograms for OF and WW instantly; and (2) evaluating autocorrelograms for OF and WW after first equalizing the spatial protection between the 2 circumstances.

For methodology (1), charge maps have been calculated as specified within the above part ‘Spatial place and route tuning’, utilizing the identical grid of three × 3-cm bins for each environments. This set of bins spanned the whole lot of the OF enviornment and coated many of the WW monitor aside from some small areas on the outer extrema, which have been discarded for the aim of this evaluation. For every of the 2 charge maps, the autocorrelogram was computed and the grid rating was calculated.

Technique (2) was much like methodology (1), besides that the cell’s OF charge map was transformed right into a ‘masked OF’ charge map, by eradicating all bins that have been unvisited by the rat within the WW session. This successfully equalized the place protection between the 2 circumstances, and thus allowed for a extra legitimate comparability.

Toroidal versus spatial description

The explanatory significance of the toroidal description was evaluated by evaluating statistical measures of how properly the toroidal coordinates defined neural exercise on the torus and in bodily area. For a good comparability, it was necessary to keep away from overfitting, which could happen if a toroidal parametrization of a degree cloud is used to explain that very same set of information factors. Two precautions have been taken to keep away from such overfitting: first, the info have been decoded utilizing the toroidal parametrization from a special situation (an OF session for a WW recording and a WW session for an OF recording), and second, the cell for which the statistical measurement was made was omitted from the decoding.

The comparability of toroidal and environmental representations additionally accounted for monitoring error within the bodily place estimate, which primarily resulted from the roughly 4 cm vertical offset of the monitoring gadget above the rat’s head. This causes a discrepancy when the angle α between the animal’s zenith and the axis of gravitation is completely different from 0°, measured as 4 tan(α) cm. The imply discrepancy within the recorded place information was measured to 1.5 cm. To account for this error of the place estimate, proportional Gaussian noise was added to the toroidal coordinates, utilizing an ordinary deviation of 1.5 cm/Ω, the place Ω denotes the grid spacing of the actual grid-cell module, estimated from the imply interval of the fitted cosine waves of the toroidal coordinates within the open discipline (see ‘Toroidal alignment’).

Info content material

The knowledge content material (I) was calculated as beforehand described42, to quantify and examine the quantity of knowledge carried by single-cell exercise in regards to the location on the torus and bodily area per spike. Each covariates have been binned in a M = 15 × 15 grid of sq. bins. For every bin  j, the common firing charge fj (given in spikes per second), and the occupancy ratio, pj, have been computed. The knowledge content material for every grid cell was then given as:

$$I=frac{1}{bar{{bf{f}}}}mathop{sum }limits_{j=1}^{M}{{bf{f}}}_{j}{log }_{2}frac{{{bf{f}}}_{j}}{bar{{bf{f}}}}{{bf{p}}}_{j},$$

the place (bar{{bf{f}}}) is the imply firing charge of the cell throughout your complete session.

Be aware that though the speed maps for bodily area have a number of firing fields, whereas the toroidal charge maps have single firing fields, we count on the spatial info to be comparable, because the measure primarily will depend on the ratio of bins with excessive firing exercise. This quantity ought to be comparable because the firing discipline dimension (in bins) will likely be inversely associated to the variety of fields within the charge map, assuming that the discretization of the map captures the related firing charge variations. For instance, given an identical binning of area, a bigger OF surroundings will embrace extra fields, however the variety of bins per discipline will lower correspondingly. The binning used ought to be adequate to resolve the smallest fields, as the identical discretization was utilized in classifying the grid cells within the recorded inhabitants.

Deviance defined

Deviance defined was computed to measure how properly a Poisson GLM mannequin fitted to the spike rely was at representing the info, utilizing both the toroidal coordinates or the tracked place as regressors. An identical set-up was used to that of a earlier research62, with a smoothness prior for the GLM to keep away from overfitting.

Each the toroidal and spatial coordinates have been binned right into a 15 × 15 grid of bins, and GLM design matrices have been constructed with entries Xi(t) = 1 if the covariate at time t fell within the i-th bin and Xi(t) = 0 in any other case.

The Poisson likelihood of recording okay spikes in time bin t is:

$$P(okay|mu (t),beta )=exp (-mu (t)){frac{mu (t)}{okay!}}^{okay},$$

the place (mu (t)=exp (sum _{i}{beta }_{i}{X}_{i}(t))) is the anticipated firing charge in time bin t. The parameters β of the Poisson GLM have been optimized for every covariate by minimizing the price operate:

$$L(beta |mu (t),gamma ,okay)=-sum _{t},{rm{l}}{rm{n}}(P(okay(t)|mu (t),beta ))+frac{1}{2}gamma sum _{i,jin N}{({beta }_{i}-{beta }_{j})}^{2},$$

the place N is the set of neighbour pairs. The primary time period is the adverse log-likelihood of the spike rely within the given time bin, whereas the second time period places a penalty on massive variations in neighbouring parameters, imposing smoothness within the covariate response of the expected spike rely.

The parameters, β, have been initialized to zero after which modified to reduce the loss operate by first working two iterations of gradient descent, earlier than optimizing utilizing the ‘l-bfgs-b’-algorithm (as carried out within the ‘scipy.optimize’-module) with ‘gtol’=1e-5 because the cut-off threshold, and eventually working two extra iterations of gradient descent. A 3-fold cross validation process was used, repeatedly becoming the mannequin to two-thirds of the info and testing on the held-out final third.

The smoothness hyperparameter γ was optimized a priori on every grid-cell module based mostly on the summed probability, testing γ (1, √10, 10, √1,000), and located to be both 1 or √10 in all instances.

Equally, after becoming a null mannequin (utilizing solely the intercept time period) and the saturated mannequin (completely becoming every spike rely), the deviance defined may very well be computed as:

$$1-frac{{{rm{ll}}}_{{rm{s}}}-{{rm{ll}}}_{{rm{p}}}}{{{rm{ll}}}_{{rm{s}}}-{{rm{ll}}}_{0}},$$

the place llp, ll0 and lls denote the cross-validated log probability of the fitted mannequin, the null mannequin and the saturated mannequin, respectively. This gives a normalized comparability describing the distinction between the fitted mannequin and the idealized mannequin.

Toroidal alignment

To deduce a geometrical interpretation of the tori, as characterised by way of the cohomological decoding, and examine the toroidal parametrizations throughout modules and circumstances, two cosine waves of the shape cos(ωt + okay) have been fitted to the OF mappings of the decoded round coordinates (Prolonged Information Fig. 5a), the place t is the centre 1002-bins of a 540° × 540°-valued 1502-bin grid rotated θ levels. The parameters (ω, okay, θ) have been optimized by minimizing the sq. distinction between the cosine waves and the cosine of the imply of the round coordinates in 1002 bins of the bodily surroundings (smoothed utilizing a Gaussian kernel with 1-bin normal deviation). Estimates have been first obtained by discovering the minimal when testing all mixtures within the following intervals, every discretized in 10 steps: ω [1,6], ϕ [0, 360) and θ [0,180). The parameters of the cosine waves were further optimized using the ‘slsqp’minimization algorithm (as implemented in the ‘scipy.optimize’-module using default hyperparameters). The period of each cosine wave was computed as 1.5 m/ω, giving a spatial scale estimate of the grid-cell modules.

As circular coordinates have arbitrary origin and orientation (that is, clockwise or counterclockwise evolution) we needed to realign the directions of the circular coordinates to compare these across modules and sessions (see Extended Data Fig. 4b). The clockwise orientation of each circular coordinate was first determined by noting whether (ωt + k) or 360° − (ωt + k) best fit the spatial mapping of the circular means of the toroidal coordinates, and subsequently reoriented to obtain the same orientation for both coordinates. The coordinate for which cos(θ) was largest (intuitively, the ‘x axis’) was then defined as the first coordinate (denoted ϕ1, with parameters (ω1, k1, θ1)) and the other as the second coordinate (ϕ2). Although (ϕ1, ϕ2) fully describe the toroidal location, the hexagonal torus allows for three axes, and the two axes obtained are thus oriented at either 60° or 120° relative to each other (see Extended Data Fig. 5b). The difference in directions was given by θ1θ2 and if this difference was greater than 90°, ϕ2 was replaced with ϕ2 + 60° ϕ1. Finally, the origin of the coordinates was aligned to a fixed reference, by subtracting the mean angular difference between the decoded coordinates and the corresponding coordinates obtained when using the toroidal parametrization of the reference OF session.

For visualization (Extended Data Fig. 5), it was furthermore necessary, in some cases, to rotate both vectors of the rhombi 30 degrees depending on whether one of the axes was directed outside of the box.

Preservation of toroidal tuning

Centre-to-centre distance and Pearson correlation were computed between toroidal tuning maps of different sessions to measure the degree of preservation between the toroidal descriptions.

First, the preferred toroidal firing location for each cell was computed as the centre of mass of the toroidal firing distribution:

$${T}_{{rm{c}}}=arctan ,2left(frac{{sum }_{i},sin ,{theta }_{i}cdot {{bf{y}}}_{i}}{{sum }_{i}{{bf{y}}}_{i}},frac{{sum }_{i},cos ,{theta }_{i}cdot {{bf{y}}}_{i}}{{sum }_{i}{{bf{y}}}_{i}},right),$$

where yi denotes the mean spike count of the given cell in the i-th bin whose binned toroidal coordinates are given by θi. The distance between mass centres found in two sessions (‘S1’ and “S2”) was then defined as:

$$d=Vert arctan ,2{(sin ({T}_{{rm{c}}}^{{S}_{2}}-{T}_{{rm{c}}}^{{S}_{1}}),cos ({T}_{{rm{c}}}^{{S}_{2}}-{T}_{{rm{c}}}^{{S}_{1}}))}_{2}Vert $$

where || ||2 refers to the L2-norm.

Pearson correlation between two tuning maps was computed by flattening the smoothed 2D rate maps to 1D arrays and calculating the correlation coefficient, r, using the ‘pearsonr’-function given in the ‘scipy.stats’-library.

To determine how much the preservation of the toroidal representations across two sessions (measured with Pearson correlation and peak distance) differed from a random distribution, the indices of the cells in one of the sessions were randomly re-ordered before computing correlation and distance for the pair of conditions. This process was repeated 1,000 times, and the P value was calculated from the rank of the original r value or distance with respect to the shuffled distribution.

Classification of grid cells

Temporal autocorrelograms were computed, for each cell, by calculating a histogram of the temporal lags between every spike and all surrounding spikes within a 200 ms window, using 1 ms bins. The histogram was then divided by the value of the zero-lag bin, which was subsequently set to zero. The autocorrelogram was smoothed using a gaussian kernel with smoothing window 4 ms. Considering the autocorrelograms of all modules during OF foraging (day 2 for R1–3) as a point cloud, the cosine distances between all points were calculated, and hence each point’s 80 nearest neighbours were found. This defined a graph in which each point described a vertex and the neighbour pairs gave rise to edges. A density estimate was then calculated as the exponential of the negative distances summed over each neighbour for each point. The graph and the density estimate were given as the input to the Gudhi implementation63 of ToMATo64. ToMATo uses a hill-climbing procedure to find modes of the density function and uses persistence to determine stable clusters. In the present case, the algorithm finds three long-lived clusters.

Minimum number of cells for torus detection

To address the question of how many cells are minimally needed to expect to see toroidal structure, random samples of n = 10, 20, …, 140 cells were taken from R2 (n = 149 cells) during OF foraging, and the same topological analysis was repeated as for the whole population. The cells were resampled 1,000 times for each number of cells in the subsample. To determine whether toroidal structure was detected, a heuristic was introduced based on the circular parameterization given by the two most persistent 1D bars in the barcode mapped onto physical space. An estimate of the resulting planar representation of the torus was obtained by fitting planar cosine waves to each mapping (see ‘Toroidal alignment’). For the analysis to be determined ‘successful’ in detecting toroidal structure, we required: (i) the mean value of the least-squares fitting (across bins of the mapping) to be less than 0.25; (ii) the angle of the rhombus to be close to 60° (between 50° and 70°); and (iii) the side lengths to be within 25% of each other.

Toroidal peak detection

The number of peaks per toroidal rate map was detected to assert the number of grid cells whose toroidal rate map portrayed single fields. First, 1,000 points were sampled from the toroidal distribution given by the mean activity of each cell in 150 × 150 bins of the stacked toroidal surface (that is, as described in ‘Toroidal rate map visualization’, each 50 × 50-binned toroidal rate map is first ‘straightened’ and subsequently stacked in 3 × 3 to address the toroidal boundaries) and then spatially smoothed using a Gaussian kernel with smoothing widths 0, 1, 2, …, 10 bins with mode set to ‘constant’ in the ‘scipy.gaussian_filter’ function. Next, the points were clustered by computing a density estimate, using the Euclidean distance, and defining neighbours as points closer than 5 bins. Cluster labels were iteratively assigned to each point and all its neighbours in a downhill manner, instantiating a new cluster identity if the point was not already labelled. Finally, the centroids for each cluster were computed and counted as a peak depending on whether its position fell within the centre 50 × 50 bins of the stacked rate maps.

Simulated CAN models

To confirm the expected outcomes of topological analyses of grid cell CAN models, grid cells were simulated using two different, noiseless CAN models (Extended Data Fig. 7).

First, a 56 × 44 grid cell network was simulated based on the CAN model proposed previously9, but using solely lateral inhibition (for details see ref. 11) in the connectivity matrix, W. The animal movement was given as the first 1,000 s of the recorded trajectory of rat ‘R’ during OF session, originally sampled at 10 ms, and interpolated to 2-ms time steps. The speed, v(t), and head direction θ(t) of the animal was calculated as the (unsmoothed) displacement in position for every time step. The activity, s, was updated as:

$${{bf{s}}}_{i+1}={{bf{s}}}_{i}+frac{1}{tau }(-{{bf{s}}}_{i}+{(I+{{bf{s}}}_{i}cdot W+alpha v(t)cos (theta (t)-mathop{theta }limits^{ sim }))}_{+}),$$

where (…)+ is the Heaviside function and (mathop{theta }limits^{ sim }) is the population vector of preferred head directions. The following parameters were used: I = 1, α = 0.15, l = 2, W0 = −0.01, R = 20 and τ = 10, and let the activity pattern stabilize by first initializing to random and performing 2,000 updates, disregarding animal movement. For computational reasons, the activity was set to 0 if si < 0.0001. The simulation was subsequently downsampled keeping only every 5th time frame.

Next, a 20 × 20 grid-cell network was simulated, for a synthetically generated OF trajectory (‘random walk’), based on the twisted torus model formulated in a previous study10. The parameter values and the code for computing both the grid cell network (choosing a single grid scale by defining the parameter ‘grid_gain’ = 0.04) and the random navigation (using 5,000 time steps) were given by the implementation by Santos Pata65.

Idealized torus models

To compare the results of both the original and simulated grid cell networks with point clouds where the topology is known, a priori, to be toroidal, points were sampled from a square and a hexagonal torus. First, a 50 × 50 (angle) mesh grid (θ1, θ2) was created in the square [0,2π)×[0,2π) and slight Gaussian noise (ϵ = 0.1N(0,1)) was added to each angle. The square torus was then constructed via the 4D Clifford torus parametrization: (cos(θ1), sin(θ1), cos(θ2), sin(θ2)). The hexagonal torus was constructed using the 6D embedding: (cos(θ1), sin(θ1), cos(a1θ1 + θ2), sin(a1θ1 + θ2), cos(a2θ1 + θ2), sin(a2θ1 + θ2)), where a1=1/√3 and a2 = −1/√3.

Histology and recording locations

Rats were given an overdose of sodium pentobarbital and were perfused intracardially with saline followed by 4% formaldehyde. The extracted brains were stored in formaldehyde and a cryostat was used to cut 30-µm sagittal sections, which were then Nissl-stained with cresyl violet. The probe shank traces were identified in photomicrographs, and a map of the probe shank was aligned to the histology by using two reference points that had known locations in both reference frames: (1) the tip of the probe shank; and (2) the intersection of the shank with the brain surface. In all cases, the shank traces were near-parallel to the cutting plane, therefore it was deemed sufficient to perform a flat 2D alignment in a single section where most of the shank trace was visible. The aligned shank map was then used to calculate the anatomical locations of individual electrodes (Extended Data Fig. 1).

Data analysis and statistics

Data analyses were performed with custom-written scripts in Python and MATLAB. Open-source Python packages used were: umap (version 0.3.10), ripser (0.4.1), numba (0.48.0), scipy (1.4.1), numpy (1.18.1), scikit-learn (0.22.1), matplotlib (3.1.3), h5py (2.10.0) and gudhi (3.4.1.post1). Samples included all available cells that matched the classification criteria for the relevant cell type. 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. All statistical tests were one-sided.

The most intensive computations were performed on resources provided by the NTNU IDUN/EPIC computing cluster66.

Additional discussion

The demonstration that populations of grid cells operate on a toroidal manifold, which is preserved across environments and behavioural states, confirms a central prediction of CAN models. The present observations provide the first—to our knowledge—population-level visualization of a two-dimensional CAN manifold, though there is accumulating evidence for one-dimensional CANs in a number of neural systems. The most powerful support for the latter has been obtained in fruit flies, in which CAN-like dynamics can be visualized in a ring of serially connected orientation-tuned cells of the central complex67,68,69. In mammals, analysis of data from dozens of simultaneously recorded head direction cells has shown that population activity in these cells faithfully traverses a conceptual ring22,23,24, in accordance with ring-attractor models17,18,19. Dynamics along low-dimensional manifolds with line, ring, or sheet topologies is also thought to underlie a wide range of other mammalian brain functions that operate on continuous scales, spanning from visual orientation tuning14 to neural operations underlying place-cell formation70,71,72, as well as motor control73, decision making and action selection74,75,76, and certain forms of memory39,77,78,79,80. The present analyses provide a visualization of 2D CAN dynamics in pure grid cells within a module and, together with the previous work, point to a widespread implementation of CAN dynamics in the brain. The existence of CAN structure to constrain activity to low-dimensional manifolds does not preclude additional mechanisms for pattern formation, however. Grid cell patterns may emerge also by feedforward mechanisms12,38,81,82,83,84,85,86. Such mechanisms may operate in parallel with recurrent networks87 and may even be the primary mechanism for grid-like firing at early stages of development, before the full maturation of recurrent connectivity11,88,89,90.

Reporting summary

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

RELATED ARTICLES

LEAVE A REPLY

Please enter your comment!
Please enter your name here

Most Popular

Recent Comments