Irregular distribution of grid cell firing fields in rats exploring a 3D volumetric space

We investigated how entorhinal grid cells encode volumetric space. On a horizontal surface, grid cells usually produce multiple, spatially focal, approximately circular firing fields that are evenly sized and spaced to form a regular, close-packed, hexagonal array. This spatial regularity has been suggested to underlie navigational computations. In three dimensions, theoretically the equivalent firing pattern would be a regular, hexagonal close packing of evenly sized spherical fields. In the present study, we report that, in rats foraging in a cubic lattice, grid cells maintained normal temporal firing characteristics and produced spatially stable firing fields. However, although most grid fields were ellipsoid, they were sparser, larger, more variably sized and irregularly arranged, even when only fields abutting the lower surface (equivalent to the floor) were considered. Thus, grid self-organization is shaped by the environment’s structure and/or movement affordances, and grids may not need to be regular to support spatial computations.


Fig. S1. Arena and lattice firing rate maps for all analysed grid cells.
Colorbars can be seen in Fig. 1. Two grid cells are shown per row. For each grid cell: Left) XY projected firing rate maps for the 1st and 2nd arena sessions respectively (top row) and autocorrelations of these (bottom row). Firing rate maps are scaled between 0Hz and the maximum firing rate value is given beneath the map. Autocorrelations are scaled between -0.2 and 1. Middle left) the volumetric firing rate map in oblique view, voxels with a firing rate less than 10% of the maximum are transparent. Lines show the outermost edges of the lattice maze. Text above the map gives the rat number, date, tetrode and cluster of the cell. Text below the map gives the maximum firing rate displayed. Middle right) projected planar firing rate maps shown in their planes of projection. These are scaled between 0Hz and the maximum value in all 3 maps. Lines show the outermost edges of the lattice maze. Right) autocorrelations of the projected firing rate maps shown in their planes of projection. These are scaled between -0.2 and 1. Lines show the outermost edges of the lattice maze. It is unknown why cell 922 21022018 t11 c2 did not fire in the lattice maze but was recorded in both arena sessions.  Fig. 1F. Histological Nissl-stained sections. For each rat, the sagittal brain section best demonstrating the position of the electrodes is given, text gives the animal identification and estimated layer of the mEC where recordings were made. Close-ups of the mEC are provided for every rat except rat 'ONE' where the whole brain section is given.  Figure 2D. a, Grid fields were distributed throughout the lattice. Proportion of fields in every lattice layer for all grid cells (n=47, black lines) in the X, Y and Z axes (top, middle and bottom). See Methods: Grid field distribution in the lattice. b, Grid fields were more variable in the lattice. Filled markers represent grid cells, open circles denote mean, error bars denote SEM. The consistency (coefficient of variation) of field volume within grid cells (n=47) and in shuffled data (n=470, see Methods: Grid field size). Lower values indicate that a cell's firing fields have consistent volumes. In the arena grid field volume was more consistent than chance (F(1,328) = 4.9, p = .0280, η 2 = 0.0146) as expected in grid cells. But this was not the case in the lattice, although this test is associated with a marginal probability (F(1,383) = 3.1, p = .0815, η 2 = 0.0079, one-way ANOVA).

Fig. S4
. Schematics and verification of the field radius and grid spacing analyses a, Two-dimensional schematic showing the field size (radius) estimation procedure. The autocorrelation (top left) was thresholded at 0.25 and field size was estimated as the radius of the remaining central peak. b, This analysis was validated on simulated FCC field arrangements with varying field sizes (n=100 arrangements with field sigmas between 1 and 10 bins and grid spacing in mm equal to 120σ; see Methods: Simulated field arrangements). Note that if fields overlap this analysis will tend to overestimate field size. LLS: linear least squares line fit. Text gives the result of a Pearson linear correlation and the slope of the LLS. c, Two-dimensional schematic showing the field spacing (scale or wavelength) estimation procedure. The median correlation found at different distances from the central peak were mapped and grid spacing was estimated as the location of the first peak after the central one. d, This analysis was validated on HCP, FCC and columnar (COL) field arrangements (n=100 each, random fields excluded as they do not have a reference spacing to validate against). Note that the size of our mazes limited field spacing to a maximum of approximately 120cm. LLS: linear least squares line fit. Text gives the result of a Pearson linear correlation and the slope of the LLS. e, The field spacing estimation results from real grid cells (mean and SEM across cells). For each cell the spacing was estimated as in c, shown here is the median correlation at every radius normalized by the first peak (the estimated spacing). A clear peak can be seen in both the arena (n=47, top) and lattice (n=27, bottom). Cells for which no spacing could be estimated are not shown.

Fig. S5
. Schematic demonstrating the planar symmetry analysis. a, Example simulated field arrangements in oblique and side views. b, Taking the autocorrelation of the FCC arrangement shown in a, it is clear there are four planes that can be cut through the centre (top row) which slice through fields arranged in a planar hexagonal pattern (bottom row). c, We extract all possible planes and calculate the hexagonal grid score for each one. Mapping these according to their pitch/azimuth yields different patterns for the different field arrangements which can be used to differentiate them. After extracting all possible planes through a cell's spatial autocorrelation we calculated the hexagonal and square grid score for each one. Mapping these according to their pitch/azimuth yields different patterns for the different field arrangements which can be used to differentiate them. In each example the grid score maps are shown as cylindrical projections and have already been corrected so that the 'best plane' is horizontal (placing the highest hexagonal grid scores at the top and bottom of the maps). Left column demonstrates the process for the χFCC score: first, the three values at 72° from the best plane and separated by 120° between them with the maximum total sum were found. In an FCC arrangement, low square grid scores are found at 50° to the best plane at the same azimuthal angles and high square grid scores can be found offset 60° from them in azimuth. Thus, χFCC was calculated as the difference between these points. Middle column demonstrates the process for the χHCP score: again, the three values at 72° from the best plane and separated by 120° between them with the maximum total sum were found. In an HCP arrangement, high square grid scores are found at 50° to the best plane and offset 60° from them in azimuth, while low square grid scores are found at 72° to the best plane and at the same azimuth angles. Thus, χHCP was calculated as the difference between these points. Right column demonstrates the process for the χCOL score: in a columnar arrangement there is a large field of pitch angles that intersect the columns and exhibit a high hexagonal grid score, however, there are virtually no high grid scores perpendicular to this plane. Thus, χCOL was calculated as the difference between the median grid score found at the angles 0 -60° & 120 -180° and the median score found at the angles 60 to 120°. These angles were chosen as they divide the spherical domain into two regions with equal surface area (the poles and an equatorial 'belt').

Fig. S7. Explanation and validation of the planar symmetry analysis.
For panels b & d: open circles denote mean, error bars denote SEM. a, Autocorrelations were randomly rotated before analysis by the planar symmetry method to avoid assumptions of the structure's orientation. This schematic shows the correct identification of the 'best plane' or the plane with the highest grid score and successful reorientation of the autocorrelation based on this. b, Filled markers represent simulated cells. Left) for every simulated field arrangement (n=100 simulated cells per arrangement) we calculated the three symmetry scores (χFCC, χHCP and χCOL), in each case the correct score is highest for its corresponding field arrangement while the random ('RND') and shuffled field ('Field') arrangements are low in all three scores. c, If simulated cells are categorized based on these scores (category = configuration with max score) they are grouped with 100% accuracy, confirming that the planar symmetry analysis can correctly identify field arrangements. d, Filled markers represent grid cells. Left) the three symmetry scores found for the real grid cells (n=46), all scores are much lower than expected if the arrangements were present and there was no significant difference between them (F(2,138) = 2.9, p = .0584, η 2 = 0.040, one-way ANOVA). Right) if grid cells are categorized as in c an equivalent proportion of cells fall into each group meaning that no one field arrangement dominates the data. e, For each grid cell (n=47 & 46) we calculated the maximum possible planar grid score (the maximum value found in the hexagonal grid score maps in Fig. S6; blue markers and distributions) and compared these to shuffled arrangements of firing fields (n=2 per grid cell; Methods: Grid field shuffle; grey markers and distributions). As expected, arena grid scores were higher than chance (t(138) = 2.7, p = .0072, Cohen's d = 0.47). But lattice scores were not different from chance (t(135) = 0.9, p = .3668, Cohen's d = 0.16) suggesting that grid cells did not exhibit hexagonal arrangements of fields on any arbitrary plane through the lattice.

Fig. S8. Planar symmetry analysis yields similar results using a field shuffle a,
We repeated the planar symmetry analysis after correcting grid cell autocorrelations for anisotropy (Methods: Autocorrelation correction for anisotropy). An example of a grid cell ratemap is shown (left), the autocorrelation of this ratemap (middle) and the same autocorrelation after correcting for field anisotropy (right). Note that the corrected autocorrelation is compressed along the vertical axis so that the central peak is more spherical than in the uncorrected version. b, The configuration results when using these autocorrelations did not differ from uncorrected ones. n=46 cells. c, In addition to the FCC, HCP, columnar and random field arrangements described elsewhere we also used a grid field shuffle which randomly shuffles field positions while maintaining their local firing activity as best as possible. Owing to the computational demands of this analysis we performed the shuffle only twice per grid cell. An example of a grid cell ratemap is shown on the left and two versions of the same ratemap with shuffled firing fields on the right. d, The ratemaps in this shuffled dataset exhibit some higher FCC and HCP scores (n=92 shuffles, red polygon) than the real grid cells (n=46, black markers). However, on average the scores are all centred on zero, like the real grid cells (values shown in Fig. S7B, 'Field' group).

Fig. S9. Schematic demonstrating the field orientation analysis.
An example firing rate map (left) is thresholded at 20% of the peak firing rate and regions which passed our criteria were considered place fields (2nd plot shows an example field). We can visualize these regions (3rd plot) and extract features such as their principal axis (red line). To visualize the orientation of multiple fields we project these axes onto a unit sphere and generate a spherical Von-Mises kernel smoothed density map (right plot), where hot colours denote that many fields 'pointed' their principal axis in this orientation.

Fig. S10. Very little change in speed cell activity and speed coding by grid cells in the lattice.
For panels b & c: open circles denote mean, error bars denote SEM. a, The mean and SEM proportion of time animals spent moving at various running speeds, averaged across sessions. b, Speed scores of speed modulated non-grid cells (n=222) dropped in the lattice (F(2,663) = 12.6, p = 4.13 × 10 -6 , η 2 = 0.0367, one-way ANOVA), however, the slope and y-intercept of speedfiring rate relationships were unaffected (p = .8047 & p = .1418, respective one-way ANOVAs). Right plot shows the mean and SEM running speed vs firing rate relationship for all speed cells. c, Speed scores, slope and y-intercepts of speed-firing rate relationships were also unaffected in grid cells (n=47, 46 & 47; p = .0921, p = .0895 & p = .6482, respective one-way ANOVAs). Right plot shows the mean and SEM running speed vs firing rate relationship for all grid cells. Averages and errors computed across sessions. a, Left) the mean PSD in each maze, averaged across sessions. Red lines denote the frequency band associated with theta rhythm. Middle) maximum power associated with the theta band in each session, these did not differ (F(2,123) = 0.1, p = .9180, η 2 = 0.0014, one-way ANOVA). Right) frequency associated with max theta power in each maze session, these did not differ (F(2,123) = 2.5, p = .0828, η 2 = 0.0397). b, Mean and SEM theta power (amplitude of the Hilbert transform) observed at various running speeds, averaged across sessions. Below this are boxplots where markers represent sessions. These show the Pearson's correlation between power and speed which did not differ between mazes (left; F(2,123) = 1.2, p = .2957, η 2 = 0.0196, one-way ANOVA) and the y-intercept of the power and speed relationship which also did not differ (right; F(2,123) = 2.9, p = .0602, η 2 = 0.0447, one-way ANOVA). c, Same as b but for theta frequency. The correlation between theta frequency and speed decreased in the second arena session with respect to the lattice (F(2,123) = 5.6, p = .0046, η 2 = 0.0839, one-way ANOVA) while the y-intercept did not differ (right; F(2,123) = 2.9, p = .8367, η 2 = 0.0029, one-way ANOVA). For each maze the left plot shows the 400ms spike autocorrelogram (1ms bins) and the right plot shows the spike-theta phase histogram. Red lines trace the amplitude of a scaled theta wave at every phase. Theta modulation and phase preference were unaffected in the lattice. b, Scatter plots showing the frequency of intrinsic theta (extracted from the spike autocorrelograms) vs extrinsic theta (extracted from the power spectral density estimate of the LFP). Markers represent grid cells, numbers represent counts that are above or below the diagonal. Markers above the dashed line represent cells firing at a rate faster than global theta, an indicator of phase precession. Cells were not modulated at a higher frequency than extrinsic theta in the arena (n=47, ߯ ଶ (1) = 0.02, p = .8840) or lattice (n=45, ߯ ଶ (1) = 2.69, p = .1011, Chi-square tests of expected proportions). c, Histograms showing the preferred theta phase (value of circular mean) of all grid cells. Red lines trace the amplitude of a scaled theta wave at every phase, blue vertical lines show the circular mean of the preferred phases. Cells showed an equal preference for one theta phase in both mazes (F(1,92) = 0.1, p = .7067, η 2 = 0.0016, one-way ANOVA comparing spike phase Rayleigh vector lengths) and these preferred phases did not differ between the mazes (k = 423, p > .99, two-sample Kuiper test). Shown below each plot is the firing probability of all grid cells relative to theta phase. Each row represents a cell and rows are sorted from top to bottom by the circular mean of the cell's phase angles from early to late. d, Filled markers represent cells, open circles denote mean, error bars denote SEM. General firing properties such as the theta index extracted from each grid cell's (n=47, 46 & 47) spike autocorrelation (top; F(2,135) = 0.9, p = .3905, η 2 = 0.0138, one-way ANOVA), burst index (middle; F(2,137) = 0.1, p = .9025, η 2 = 0.0015, one-way ANOVA) and the full width at half maximum of the interspike interval (fwhm of the ISI: bottom; F(2,130) = 0.1, p = .8616, η 2 = 0.0023, one-way ANOVA) were unchanged in the lattice maze. -Session and cell totals for each rat. Histological verification of the mEC layers can be seen in Fig. S2. *Note that two rats did not yield any grid cells but contributed to the behaviour and other electrophysiological data. a Rats 873 and 895 exhibited horizontal planar grid firing patterns.