Abstract
As animals navigate on a twodimensional surface, neurons in the medial entorhinal cortex (MEC) known as grid cells are activated when the animal passes through multiple locations (firing fields) arranged in a hexagonal lattice that tiles the locomotion surface^{1}. However, although our world is threedimensional, it is unclear how the MEC represents 3D space^{2}. Here we recorded from MEC cells in freely flying bats and identified several classes of spatial neurons, including 3D border cells, 3D headdirection cells, and neurons with multiple 3D firing fields. Many of these multifield neurons were 3D grid cells, whose neighbouring fields were separated by a characteristic distance—forming a local order—but lacked any global lattice arrangement of the fields. Thus, whereas 2D grid cells form a global lattice—characterized by both local and global order—3D grid cells exhibited only local order, creating a locally ordered metric for space. We modelled grid cells as emerging from pairwise interactions between fields, which yielded a hexagonal lattice in 2D and local order in 3D, thereby describing both 2D and 3D grid cells using one unifying model. Together, these data and model illuminate the fundamental differences and similarities between neural codes for 3D and 2D space in the mammalian brain.
Access options
Subscribe to Journal
Get full journal access for 1 year
$199.00
only $3.90 per issue
All prices are NET prices.
VAT will be added later in the checkout.
Tax calculation will be finalised during checkout.
Rent or Buy article
Get time limited or full article access on ReadCube.
from$8.99
All prices are NET prices.
Data availability
The data generated and analysed in the current study are available from the corresponding author on reasonable request. Source data are provided with this paper.
Code availability
The code generated for the current study is available from the corresponding author on reasonable request.
References
 1.
Hafting, T., Fyhn, M., Molden, S., Moser, M.B. & Moser, E. I. Microstructure of a spatial map in the entorhinal cortex. Nature 436, 801–806 (2005).
 2.
Finkelstein, A., Las, L. & Ulanovsky, N. 3D maps and compasses in the brain. Annu. Rev. Neurosci. 39, 171–196 (2016).
 3.
Krupic, J., Burgess, N. & O’Keefe, J. Neural representations of location composed of spatially periodic bands. Science 337, 853–857 (2012).
 4.
Stensola, T., Stensola, H., Moser, M.B. & Moser, E. I. Shearinginduced asymmetry in entorhinal grid cells. Nature 518, 207–212 (2015).
 5.
Hayman, R., Verriotis, M. A., Jovalekic, A., Fenton, A. A. & Jeffery, K. J. Anisotropic encoding of threedimensional space by place cells and grid cells. Nat. Neurosci. 14, 1182–1188 (2011).
 6.
Hayman, R. M., Casali, G., Wilson, J. J. & Jeffery, K. J. Grid cells on steeply sloping terrain: evidence for planar rather than volumetric encoding. Front. Psychol. 6, 925 (2015).
 7.
Casali, G., Bush, D. & Jeffery, K. Altered neural odometry in the vertical dimension. Proc. Natl Acad. Sci. USA 116, 4631–4636 (2019).
 8.
Yartsev, M. M., Witter, M. P. & Ulanovsky, N. Grid cells without theta oscillations in the entorhinal cortex of bats. Nature 479, 103–107 (2011).
 9.
Taube, J. S., Muller, R. U. & Ranck, J. B. Jr Headdirection cells recorded from the postsubiculum in freely moving rats. I. Description and quantitative analysis. J. Neurosci. 10, 420–435 (1990).
 10.
Sargolini, F. et al. Conjunctive representation of position, direction, and velocity in entorhinal cortex. Science 312, 758–762 (2006).
 11.
Finkelstein, A. et al. Threedimensional headdirection coding in the bat brain. Nature 517, 159–164 (2015).
 12.
Solstad, T., Boccara, C. N., Kropff, E., Moser, M.B. & Moser, E. I. Representation of geometric borders in the entorhinal cortex. Science 322, 1865–1868 (2008).
 13.
Savelli, F., Yoganarasimha, D. & Knierim, J. J. Influence of boundary removal on the spatial representations of the medial entorhinal cortex. Hippocampus 18, 1270–1282 (2008).
 14.
Yartsev, M. M. & Ulanovsky, N. Representation of threedimensional space in the hippocampus of flying bats. Science 340, 367–372 (2013).
 15.
Hales, T. C. A proof of the Kepler conjecture. Ann. Math. 162, 1065–1185 (2005).
 16.
Stella, F. & Treves, A. The selforganization of grid cells in 3D. eLife 4, e05913 (2015).
 17.
Mathis, A., Stemmler, M. B. & Herz, A. V. M. Probable nature of higherdimensional symmetries underlying mammalian gridcell activity patterns. eLife 4, e05979 (2015).
 18.
Horiuchi, T. K. & Moss, C. F. Grid cells in 3D: reconciling data and models. Hippocampus 25, 1489–1500 (2015).
 19.
Boccara, C. N. et al. Grid cells in pre and parasubiculum. Nat. Neurosci. 13, 987–994 (2010).
 20.
Krupic, J., Bauza, M., Burton, S. & O’Keefe, J. Local transformations of the hippocampal cognitive map. Science 359, 1143–1146 (2018).
 21.
Boccara, C. N., Nardin, M., Stella, F., O’Neill, J. & Csicsvari, J. The entorhinal cognitive map is attracted to goals. Science 363, 1443–1447 (2019).
 22.
SanguinettiScheck, J. I. & Brecht, M. Home, head direction stability, and grid cell distortion. J. Neurophysiol. 123, 1392–1406 (2020).
 23.
Krupic, J., Bauza, M., Burton, S., Lever, C. & O’Keefe, J. How environment geometry affects grid cell symmetry and what we can learn from it. Phil. Trans. R. Soc. Lond. B 369, 20130188 (2013).
 24.
Stensola, H. et al. The entorhinal grid map is discretized. Nature 492, 72–78 (2012).
 25.
Kropff, E. & Treves, A. The emergence of grid cells: intelligent design or just adaptation? Hippocampus 18, 1256–1269 (2008).
 26.
Burak, Y. & Fiete, I. R. Accurate path integration in continuous attractor network models of grid cells. PLoS Comput. Biol. 5, e1000291 (2009).
 27.
McNaughton, B. L., Battaglia, F. P., Jensen, O., Moser, E. I. & Moser, M.B. Path integration and the neural basis of the ‘cognitive map’. Nat. Rev. Neurosci. 7, 663–678 (2006).
 28.
Fiete, I. R., Burak, Y. & Brookings, T. What grid cells convey about rat location. J. Neurosci. 28, 6858–6871 (2008).
 29.
Mathis, A., Herz, A. V. M. & Stemmler, M. B. Resolution of nested neuronal representations can be exponential in the number of neurons. Phys. Rev. Lett. 109, 018103 (2012).
 30.
Stemmler, M., Mathis, A. & Herz, A. V. M. Connecting multiple spatial scales to decode the population activity of grid cells. Sci. Adv. 1, e1500816 (2015).
 31.
Sarel, A., Finkelstein, A., Las, L. & Ulanovsky, N. Vectorial representation of spatial goals in the hippocampus of bats. Science 355, 176–180 (2017).
 32.
Omer, D. B., Maimon, S. R., Las, L. & Ulanovsky, N. Social placecells in the bat hippocampus. Science 359, 218–224 (2018).
 33.
Ulanovsky, N. & Moss, C. F. Hippocampal cellular and network activity in freely moving echolocating bats. Nat. Neurosci. 10, 224–233 (2007).
 34.
Yovel, Y., Falk, B., Moss, C. F. & Ulanovsky, N. Optimal localization by pointing off axis. Science 327, 701–704 (2010).
 35.
Skaggs, W. E., McNaughton, B. L., Wilson, M. A. & Barnes, C. A. Theta phase precession in hippocampal neuronal populations and the compression of temporal sequences. Hippocampus 6, 149–172 (1996).
 36.
Henriksen, E. J. et al. Spatial representation along the proximodistal axis of CA1. Neuron 68, 127–137 (2010).
 37.
Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO—the Open Visualization Tool. Model. Simul. Mater. Sci. Eng. 18, 015012 (2010).
 38.
Larsen, P. M., Schmidt, S. & Schiøtz, J. Robust structural identification via polyhedral template matching. Model. Simul. Mater. Sci. Eng. 24, 055007 (2016).
 39.
Brandon, M. P. et al. Reduction of theta rhythm dissociates grid cell spatial periodicity from directional tuning. Science 332, 595–599 (2011).
 40.
Koenig, J., Linder, A. N., Leutgeb, J. K. & Leutgeb, S. The spatial periodicity of grid cells is not sustained during reduced theta oscillations. Science 332, 592–595 (2011).
 41.
Hansen, J.P. & Verlet, L. Phase transitions of the LennardJones system. Phys. Rev. 184, 151–161 (1969).
 42.
Langston, R. F. et al. Development of the spatial representation system in the rat. Science 328, 1576–1580 (2010).
 43.
Wills, T. J., Cacucci, F., Burgess, N. & O’Keefe, J. Development of the hippocampal cognitive map in preweanling rats. Science 328, 1573–1576 (2010).
 44.
Eliav, T. et al. Nonoscillatory phase coding and synchronization in the bat hippocampal formation. Cell 175, 1119–1130 (2018).
 45.
Derdikman, D. et al. Fragmentation of grid cell maps in a multicompartment environment. Nat. Neurosci. 12, 1325–1332 (2009).
 46.
Torquato, S. & Stillinger, F. H. Jammed hardparticle packings: from Kepler to Bernal and beyond. Rev. Mod. Phys. 82, 2633 (2010).
 47.
D’Albis, T. & Kempter, R. A singlecell spiking model for the origin of gridcell patterns. PLoS Comput. Biol. 13, e1005782 (2017).
 48.
MonsalveMercado, M. M. & Leibold, C. Hippocampal spiketiming correlations lead to hexagonal grid fields. Phys. Rev. Lett. 119, 038101 (2017).
 49.
Weber, S. N. & Sprekeler, H. Learning place cells, grid cells and invariances with excitatory and inhibitory plasticity. eLife 7, e34560 (2018).
 50.
Yoon, K. et al. Specific evidence of lowdimensional continuous attractor dynamics in grid cells. Nat. Neurosci. 16, 1077–1084 (2013).
 51.
Guanella, A., Kiper, D. & Verschure, P. A model of grid cells based on a twisted torus topology. Int. J. Neural Syst. 17, 231–240 (2007).
 52.
Fuhs, M. C. & Touretzky, D. S. A spin glass model of path integration in rat medial entorhinal cortex. J. Neurosci. 26, 4266–4276 (2006).
 53.
Klukas, M., Lewis, M. & Fiete, I. Efficient and flexible representation of higherdimensional cognitive variables with grid cells. PLoS Comput. Biol. 16, e1007796 (2020).
 54.
Burak, Y. & Fiete, I. Do we understand the emergent dynamics of grid cell activity? J. Neurosci. 26, 9352–9354 (2006).
 55.
Rowland, D. C., Roudi, Y., Moser, M.B. & Moser, E. I. Ten years of grid cells. Annu. Rev. Neurosci. 39, 19–40 (2016).
Acknowledgements
We thank A. Treves for extensive discussions, A. V. M. Herz and E. Katzav for helpful suggestions; F. Stella, E. D. Karpas, T. Tamir, A. Rubin, A. Finkelstein, S. R. Maimon, S. Ray, D. Omer, S. Palgi and A. Sarel for comments on the manuscript; J.M. Fellous for contributing to the initial stages of this project; S. Kaufman, O. Gobi and I. Shulman for bat training; A. Tuval for veterinary support; C. Ra’anan and R. Eilam for histology; M. P. Witter for advice on histological delineation of MEC borders; B. Pasmantirer and G. Ankaoua for mechanical designs; and G. Brodsky and N. David for graphics. We thank the Methods in Computational Neuroscience course at MBL, Woods Hole, where the modelling work was initiated. N.U. is the incumbent of the Barbara and Morris Levinson Professorial Chair in Brain Research. Y.B. is the incumbent of the William N. Skirball Chair in Neurophysics. This study was supported by research grants to N.U. from the European Research Council (ERCCoG – NATURAL_BAT_NAV), Israel Science Foundation (ISF 1319/13), and Minerva Foundation, and by the André Deloro Prize for Scientific Research and the Kimmel Award for Innovative Investigation to N.U. Y.B. acknowledges support from the Israel Science Foundation (ISF 1745/18), the European Research Council (ERCSyG – KiloNeurons), and the Gatsby Charitable Foundation. H.S. acknowledges support from the Gatsby Charitable Foundation.
Author information
Affiliations
Contributions
G.G, L.L. and N.U. conceived and designed the experiments. G.G. conducted the experiments, with contributions from L.L. G.G. analysed the experimental data. L.L. and N.U. guided the data analysis. J.A., Y.B. and H.S. conducted the theoretical modelling, and J.A. analysed the model results. G.G. and N.U. wrote the first draft of the manuscript, with major input from L.L. All authors participated in writing and editing of the manuscript. N.U. supervised the project.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Peer review information Nature thanks Mark Brandon, Torkel Hafting and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Extended data figures and tables
Extended Data Fig. 1 Expected properties of 3D grid cells, and coverage of 3D space by the bats’ behaviour.
a, Five basic properties of 2D grid cells. b, Expected properties of 3D grids under five hypothetical possibilities, ranging from highly ordered to random arrangement of fields. c, d, Bat behaviour covers 3D space. c, Examples of bat trajectories in the flight room. Rows: one example session from each of the four bats included in this study. Columns show different viewing perspectives: 3D view; top view of the room (xy); and side views of the room (xz and yz). Note that for bat 4 only, part of the room (grey area) was blocked off by a seethrough nylon mesh (see Methods). d, Histograms of distance flown per session, plotted separately for each of the four bats; we included here only sessions in which we recorded at least one of the 125 wellisolated neurons that passed the inclusion criteria for analysis.
Extended Data Fig. 2 A series of histological sections from one bat, showing two tetrode tracks in MEC layer 5.
Two tetrode tracks that enter from postrhinal cortex (POR) and progress along layer 5 of MEC in bat 2. Four sagittal Nisslstained sections are arranged from medial (top) to lateral (bottom). Left, wide view (scale bar, 1,000 μm); right, zoomin onto the region of the tetrode tracks (scale bar, 500 μm). The postrhinal (dorsal) border of MEC is marked by a black line. The combination of the angled tetrode penetration (see Methods) and the sagittal slicing of the brain resulted in tetrode tracks that are recognizable as a small hole (sometimes surrounded by glial scar) that proceeds over many consecutive sections. The distance of the hole (track) from the postrhinal border of MEC is indicated on the right. Coloured circles mark the tracks of the two tetrodes (purple, tetrode 1 (TT1); green, tetrode 4 (TT4)). The lesion at the end of the track of tetrode 1 (purple) is visible in section 31c (section numbers are indicated). The lesion at the end of the track of tetrode 4 (green) is visible in section 29b.
Extended Data Fig. 3 Firing fields, and additional examples of cells.
a, Peak detection: firingrate map of an example cell, with overlaid black dots marking the positions of the detected field peaks (see Methods). Note that only the peaks of the visible fields are displayed here (not plotting the peaks for fields ‘buried’ deep inside the 3D volume of the map). b, Threshold for peak detection shuffle (one of three criteria for a field to be valid; these three criteria narrowed the n = 125 cells to n = 78 cells): an identified peak was included in the analysis if the firing rate in its peak voxel was higher than the firing rate in the same voxel in 75% of the spike shuffles (see Methods). Shown are the shuffle percentiles for all fields (both above and below this threshold); the 75% threshold was set at a natural kink in the distribution. c, Number of fields per neuron: a histogram for all the neurons with valid fields (n = 78 cells), showing the number of fields per cell. A threshold of 10 fields was set for a neuron to be considered as ‘multifield’ (n = 66 cells, black bars; grey bars show cells below the threshold); see Methods for the rationale of setting this 10field threshold. Example cells with varying numbers of fields are shown in d. d, Examples of cells with varying numbers of fields, including <10 fields (two left cells, not considered as multifield cells), and both grid and nongrid multifield cells with ≥10 fields. Plotted as in Fig. 2a, b. Top: 3D firingrate maps. Bottom: box plots showing the median field sizes in the x, y and z directions. Horizontal line, median; box limits, 25th to 75th percentiles; whiskers,10th to 90th percentiles. Shown are the results of Wilcoxon ranksum tests comparing field dimensions in x vs. y, y vs. z and x vs. z, Bonferronicorrected for three pairwise comparisons; n.s., nonsignificant; the number of fields for each cell (that is, the n for each test) is written above each firingrate map. e, Two example neurons with significant field elongations; plotted as in d (overall, 9 out of 66 neurons showed significant field elongations). Note that these two neurons differed in their elongation direction. Wilcoxon ranksum test comparing field dimensions in x vs. y, y vs. z and x vs. z, Bonferronicorrected for three pairwise comparisons: left cell: P_{yz} = 0.03, P_{xz} = 0.003; right cell: P_{xy} = 0.006 (Bonferronicorrected for three comparisons); the number of fields for each cell (that is, the n for each test) is written above each firingrate map. In the 9 cells with significant field elongation, fields did not resemble columns and elongation was weak and not systematic, that is, the elongation direction was heterogeneous and differed across these 9 elongated neurons.
Extended Data Fig. 4 Nonstereotypy of flights, and diversity of passes through fields: part I.
a, Passes per field: histograms showing the number of flight passes creating each firing field (passes per field and passes with spikes per field serve as two of three criteria for a detected field to be considered valid). The four histograms correspond to all passes within the fields (top) and all passes with spikes within the field (bottom), computed for two different values of the radius from the field’s peak: 30 cm (left) and 50 cm (right). b, Flight nonstereotypy: a histogram showing the percentage of flights of each cell that passed through the same sequence of fields, within two different values of the distance (radius) from the field’s peak: 30 cm (left) and 50 cm (right). Flights that pass through the same sequence of fields correspond to a stereotypic trajectory of the bat: for example, if the same flight sequence repeats four times, it would appear here in the fourth bin, that is, with a value of 4 for ‘No. of flights with same sequence’. The results here show that most flights are in fact unique and do not repeat more than once for the same cell (see large bin with value 1 for ‘No. of flights with same sequence’), which suggests that the bats’ flights were nonstereotypic. c, The xy projections of the trajectories that passed through 20 example fields. Shown are trajectories within a 0.5m radius around the field’s peak (this is the typical radius of grid fields in our data). We extracted the trajectories for all the 1,113 fields of the 66 multifield neurons, as follows: for each field we took a 3D sphere with 0.5m radius around the field’s peak, then extracted the 3D trajectories through this field, and then projected these 3D trajectories onto xy (2D projections). The Rayleigh vector length (RV) of the xyprojected trajectories was then computed for each field. Here we plotted one example field from approximately every fifth percentile of the RV distribution—a total of 20 example fields, ordered from low to high RV. The RV value and percentile (in parentheses) are indicated: top left example is from the first percentile of the RV = lowest RV (most isotropic trajectories); bottom right example is from the 99th percentile of the RV = highest RV (most unidirectional trajectories). Scale bar, 0.5 m. d, The RV distribution across all the fields (n = 1,113), calculated on the trajectories forming each field (see examples in c). Median RV of this distribution = 0.25: this is a relatively low value, indicating rather isotropic flights through most fields.
Extended Data Fig. 5 Nonstereotypy of flights, and diversity of passes through fields: part II.
a, Two examples of individual cells (columns), depicting the firing rates in each of the fields based on the trajectories coming from specific takeoff balls (top matrices) or specific landing balls (bottom matrices); colourcoded from 0 (blue) to maximum firing rate (red; value indicated). In each matrix, the variety of firing rates in each row of the matrix (corresponding to each field) shows that the firing of each field was created by trajectories originating from a variety of takeoff balls and ending on a variety of landing balls. b, Percentage of fields per cell that are created from diverse trajectories. Box plots show the percentages of fields per cell whose behavioural trajectories involved at least 50% of the takeoff balls (left) or landing balls (right). This analysis was done for the three bats for which we used landing balls. Horizontal line, median; box limits, 25th to 75th percentiles; whiskers,10th to 90th percentiles. The high percentage shown here means that across cells, the large majority of fields per neuron involved most of the takeoff balls and landing balls—that is, the bats’ flights were diverse and nonstereotypic. c, Histograms showing the number of fields that were formed by trajectories that originated from various numbers of different takeoff balls (top) and ending at various numbers of different landing balls (bottom). This analysis was done for the three bats for which we used landingballs. d–f, As in a–c, but here instead of takeoff and landing balls we examined the firingrate at different directions of the trajectories passing through the field. This analysis was done on all the four bats. d, Examples of two individual cells. e, Percentage of fields in which the trajectories with spikes that occurred close to the fieldpeak (within 50 cm of the field peak) came from ≥50% of the direction bins. f, Distribution of the percentage of direction bins for trajectories with spikes that passed through each field. For example, 60% on the xaxis means that the spikes of that field came from 60% of all the possible direction bins; that is, the flight trajectories with spikes that passed through that field spanned 60% of all the possible direction bins (60% of 360°; we used here 36° bins)—a high diversity of directions.
Extended Data Fig. 6 Searching for global and local order.
a, Two lattices were proposed as theoretically probable structures of 3D grids^{16,17}: these lattices exhibit a hexagonal arrangement of spheres, and were proved to be the optimal arrangement for densest sphere packing in infinite 3D space^{15}. These lattice arrangements—face centred cubic (FCC) and hexagonal close pack (HCP)—are widely found in chemistry. Both of these lattices are comprised of layers of spheres arranged in a hexagonal pattern (denoted here as layers A, B and C); the two lattices differ in the phase of these layers. Each layer is placed in one of two optional sets of holes created by the lower layer: the HCP lattice has an A–B–A pattern of layers (left); the FCC lattice has an A–B–C pattern of layers (right). b, Illustration of the procedure of fitting data to hexagonal lattices (see Methods). Cartoon in 2D depicts the fitting that we used in 3D to fit the recorded 3D field arrangements to perfect FCC and HCP lattices. To find the bestfitting lattice for each of our recorded cells, we compared each cell with 3 million variations of FCC and HCP perfect lattices; this was done separately for FCC and HCP lattices (1.5 million variations for FCC and 1.5 million for HCP). In each iteration, we found for each of the recorded fields (purple spheres) the closest field in the perfect lattice (the simulated perfect lattice fields are depicted by dotted blue circles; red lines depict the distances between a recorded field and its nearest simulated lattice field). We then computed the squared distance from each recorded field to its closest simulated synthetic field; we summed these squared distances and computed the root, to yield the root mean square (RMS) error of each cell for each of the 1.5 million FCCs or 1.5 million HCPs. The cell’s lowest RMS error across all the 1.5 million FCC lattices was termed the ‘FCC fit error’, and was then compared to the FCC fit errors that were analogously computed for all the shuffles (that is, for each shuffle we found the minimal FCC fit error over the 1.5 million synthetic FCC’s). The same was done for HCP. c, Distribution of local angles. There was no anglebias of local 3D angles, suggesting no local orientational (angular) order in the field arrangements. Top, local angles were computed between triplets of nearestneighbouring fields; to avoid edge effects, we computed the angles only for fields that were within the convex hull of field arrangement. Grey, all recorded multifield cells; black line, shuffled data (using spike shuffling). Middle and bottom, distributions of local angles plotted separately for 3D grid cells (middle: red) and nongrid multifield cells (bottom: blue), together with shuffle data (black lines). The distributions for recorded cells did not differ significantly from the spikeshuffled data (black; indicated are the P values of Kolmogorov–Smirnov tests between data and shuffles). d, e, The calculation of the CV of distances to nearest neighbouring fields. d, Cartoon illustrating the three distances between each field and its three nearest neighbouring fields. These distances were then pooled across all fields (while avoiding doublecounts), and the CV over these distances was calculated to quantify how fixed is the distance scale of the neuron (see Methods). e, Rationale for why the CV was computed over three nearest neighbours, and not 1, 2 or >3 nearest neighbours. Right, a perfect 2D lattice arrangement will yield low CV values regardless of how many nearest neighbours are considered: three (bottom) or two (top). Left, a periodic 1D structure (‘fields along a string’)—which is not a 2D lattice—may still yield a low CV when considering only two neighbours (top left), but a high CV when considering three or more nearest neighbours (bottom left). Thus, three is the minimal number of nearest neighbours that can differentiate these possibilities, and unequivocally characterize the 3D arrangement of fields. We did not use more than three neighbours, because using a large number of neighbours may accentuate possible edge effects. f, Comparison of the CV measure and gridness score for cells in 2D. Left, scatter plot of the CV of nearestneighbour distances versus the classical 2D gridness score, for medial entorhinal cortex neurons recorded in rats running in 2D in the Moser laboratory (courtesy of M.B. Moser and E. I. Moser, from https://www.ntnu.edu/kavli/research/gridcelldata; n = 259 cells; we excluded from analysis putative interneurons (mean firing rate >5 Hz), nonspatial neurons (cells with spatial information <0.2 bits per spike) and border cells (manually excluded)). The Pearson correlation of this scatter and its P value are given in the topright corner. Right, example firingrate maps of two cells. Top, a spatial nongrid cell with low gridness score and high CV; bottom, a grid cell with high gridness score and low CV. The locations of these two cells on the scatter plot are indicated. Here the gridness score was computed as described in the Methods, and the CV was computed using the distances of the two nearest neighbours of each field, rather than three nearest neighbours as done for 3D cells, because in 2D many cells did not have enough fields. g, h, The lack of global order could not be explained by proximity to room borders. g, Local angles throughout the room: shown is a scatter plot of the local angles between triplets of neighbouring fields as a function of their distance from the nearest wall. Points were slightly jittered in the xaxis, for display purposes only. Pearson correlation and significance are shown. The nonsignificant correlation suggests that the walls did not affect the 3D arrangement of fields. h, CV of nearestneighbour distances throughout the room: shown is a histogram of the differences between CVs of fields at the edges of the room and the CVs of fields in the centre of the room (CV_{edges} − CV_{centre}). The border between ‘centre’ and ‘edges’ of the room was defined as the median distance of all the fields from their closest wall. Grey bars, histogram for all multifield cells (n = 66 cells); black line, histogram for shuffles done per neuron, using spike shuffling. P value of Kolmogorov–Smirnov test of data versus shuffles is indicated. The nonsignificant difference suggests that the walls did not affect the 3D arrangement of fields.
Extended Data Fig. 7 The firing of 3D grid cells and nongrid multifield cells cannot be explained via a global lattice that undergoes distortions (scaling, shearing, or barrelling).
a–c, Mean CVs of distances between neighbouring fields, calculated after undistorting the field positions of multifield cells (n = 66); shown separately for each individual bat (columns 1–4), and for all the cells from all bats pooled together (column 5). We did three types of undistortion: unscaling, unshearing, and unbarrelling (see Methods). a, Unscaling: field positions were scaled by a factor α along the specified axis. α < 1 corresponds to shrinking towards the middle of the room along the specified axis, α > 1 means expansion, α = 1 corresponds to no change. The three rows correspond to unscaling along three pairs of axes: xy, yz, and xz. Each heatmap describes the mean CV for all possible combinations of unscaling of two axes: α_{1} describes the first axis (in xy (first row) it corresponds to x), and α_{2} describes the second axis (in xy it corresponds to y). CV is colourcoded from minimal value (blue) to maximal value (red), and these values are written above each matrix. White squares indicate no scaling (α_{1} = α_{2} = 1, that is, original data). b, Unshearing: the fields of each cell underwent shearing/unshearing as described^{4} (see Methods), where shearing was done by two walls simultaneously (yielding 15 possible combinations of shearing, all shown here (rows)); we computed the two possible shearing directions (positive and negative γfactor, which can also be thought of as shearing and unshearing). The writing on the left indicates the two walls that create the shearing and the axes in which the shearing is occurring (for example, xz shear x means that the xz wall is creating a shearing along the xaxis). Factor γ_{1} corresponds to the first line (for example, xz shear x in the top row), and γ_{2} corresponds to the second line (for example, yz shear y in the top row). Heatmap shows mean values of CV, colourcoded from minimum (blue) to maximum (red; values indicated above the panel). White squares indicate no shearing (γ_{1} = γ_{2} = 0, that is, original data). c, Unbarrelling: a barrelling transformation (also known as fisheye lens distortion) and unbarrelling transformation (also known as pincushion distortion) was performed on the field arrangement of each cell. Plots show mean CV values as a function of the barrel and unbarrel factor β, for the different 2D planes along which the barrelling/unbarrelling was done (xy, yz, xz (rows)). All the transformations in a–c (scaling, shearing, barrelling) showed that distorting or undistorting away from the original field arrangement increases the CV—which argues against the notion that the original field arrangements constitute a global lattice that underwent distortions.
Extended Data Fig. 8 Both 3D grid cells and nongrid multifield cells exhibit an anatomical gradient along the dorsoventral axis of MEC, as seen in rodent 2D grid cells.
a, Examples of recording sites and firingrate maps of cells along the interfield distance gradient, here using the spikeshuffling method as in Fig. 3b. Central scatterplot: both 3D grid cells (red) and nongrid multifield cells (blue) exhibited an increase in interfield distance as a function of the recording depth along the dorsoventral axis of MEC—akin to what is known from rodent 2D grid cells^{1}. 3D grid cells were identified using spike shuffling; data points were pooled here across tetrodes and bats, and include only cells with highquality histology that allowed precise assessment of recording depth within MEC (n = 8 grid cells, n = 28 nongrid multifield cells). Side panels: two histology examples, and their associated 3D firingrate maps. Left, 3D grid cell recorded in layer 3; right, nongrid multifield cell recorded in layer 5. b, The increase in interfield distance along the dorsoventral axis of MEC is also evident in individual tetrode tracks (TT). Shown are examples of three tetrode tracks, which proceeded mostly along a single cortical layer (indicated). Pearson correlation coefficients and their P values are indicated. c, Dorsoventral gradient, plotted as in Fig. 3b, but here classifying grid cells based on the fieldshuffling method (rather than the spikeshuffling method as in Fig. 3b). We identified 20 neurons as 3D grid cells using field shuffling; included in these scatterplots were only cells with highquality histology that allowed precise assessment of recording depth within MEC (n = 11 grid cells, n = 25 nongrid multifield cells). Left, gradient of interfield distances between neighbouring fields as a function of recording depth in MEC (distance from the postrhinal border of MEC); right, gradient of number of fields per cell as a function of recording depth. Note that the gradient is exhibited by both 3D grid cells (red) and nongrid multifield cells (blue), suggesting they are not two separate populations of neurons, but a continuum. Shown are Pearson correlation coefficients and their P values.
Extended Data Fig. 9 Border cells and firing characteristics of the various cell classes.
a–c, Border cells. a, Definition of border cells in 3D arenas. We considered only neurons for which a 2D firing field was detected in each of the cell’s 2D planes (see Methods). 2D firingrate maps were plotted in the three principal 2D views of the room geometry: the top view xy, and the two side views xz and yz (green, orange and blue, respectively). We defined an elongation index for each field in each 2D plane as the field size along the more elongated axis/field size along the less elongated axis; the elongation index was computed separately for all three 2D planes (xy, xz, yz). A border field was defined as a field that was significantly elongated in exactly two out of the three views (for determining the significance we used the elongation index for each field and compared it to the elongation indices of spike shuffles; 99% criterion). Significant 1D border cells (left) were defined as ‘sausageshaped’ neurons in which: (i) two of the 2D planes (2D views) exhibited firing fields with an elongation index > 2, which was also significantly higher than shuffled data (>99% of the shuffles); and (ii) the third plane (2D view) exhibited a firing field that was nonelongated, and was smaller than in the two elongated planes (left, blue square), thus defining a 1D ‘sausage’. Significant 2D planar border cells (right) were defined as ‘pancakeshaped’ neurons in which: (i) two of the 2D planes (2D views) exhibited firing fields with elongation index >2, which was also significantly higher than shuffled data (>99% of the shuffles); and (ii) the third plane (2D view) exhibited a firing field that was nonelongated, and was larger than in the two elongated planes (right, green square), thus defining a 2D ‘pancake’. b, c, A full depiction of the example border cells plotted in Fig. 4b, here shown also with trajectories and spikes, as well as an additional example. b, Example of a border cell that exhibited firing along a 1D border of space—in this case, along the wall edge. Left, a top view of the room (xy); middle, right, side views (xz and yz). Top, flight trajectory (grey) with spikes superimposed (red dots); bottom, firing rate maps, colourcoded from zero (blue) to the maximal firing rate across the three views of each cell (red). White line, edges of the detected firing field (enlarged by one bin for display purposes); indicated are the peak firingrate and elongation index (EI; see Methods). Grey rectangle, area of the room that was made inaccessible for one of the four bats (see Methods). c, Examples of border cells that exhibited firing along a 2D planar border of space—in this case, the floor. Plotted as in b. Top cell also shown in Fig. 4b. d, e, Characteristics of the various cell classes. d, Firing characteristics for the different functional cell classes in MEC during 3D flight (mean ± s.d. across cells; the column of headdirection cells is based on significant azimuthal headdirection cells). Peak firing rate (FR): the highest firing rate in the 3D firing rate map. Mean firing rate: the mean of the 3D firing rate map (second row) and the ratio of total spikes in flight/total time in flight (third row); these values were very similar because the 3D behaviour was quite uniform. Fourth row: Theta index^{8,19,39,40} (TI) computed as in ref. ^{19}, where a criterion of TI > 5 was used as threshold for defining thetamodulated cells. All the cells in our study, including all the grid cells, exhibited TI < 5, indicating a lack of theta rhythmicity in the MEC of flying bats, consistent with previous findings on bat MEC cells in 2D^{8,44}. e, Examples of spiketrain autocorrelations (top) and their power spectra (bottom) for the various cell classes. The theta index is indicated above the power spectra. Note that for the four grid cells, the two leftmost examples have theta index below the population mean (see d), whereas the two rightmost examples have theta index above the population mean; none of these four examples exhibited clear oscillations.
Extended Data Fig. 10 Distribution of bestfitting temperatures in the Lennard–Jones pairwise interactions model, and comparison to alternative pairwiseinteractions models.
a–c, Simulations of highly structured cells with low CV (CV < 0.25, n = 29 neurons). a, Distribution of the bestfitting temperatures (T) in the Lennard–Jones (LJ) pairwise interactions model (red, 3D grid cells identified by spike shuffling; blue, nongrid multifield cells). The bestfitting temperature for the model (log_{10}T = 0.5; black dashed line) was computed as the average logarithm of the temperatures across all these cells; this temperature was consistent with the individual temperatures for almost all of the cells. The small variation in the fitted temperatures allowed us to use a single temperature for all the cells (log_{10}T = 0.5, dashed line). b, c, Comparison of multiple models of shortrange repulsion and longrange attraction (Lennard–Jones, modified Lennard–Jones, and differenceofGaussians (DOG) potentials). b, Quantifying the match between model and data. Three rightmost bars: distances between the RDF of the data and the RDFs of the pairwise interactions models: Lennard–Jones, modified Lennard–Jones, and DOG (log_{10}T = 0.5 for all models). Left bar, distance between the RDF of the data and the RDF of the random Poisson model (in which the fields were randomly distributed in the 3D region visited by the bat). N = 1,000 simulated RDFs for each model; shown mean ± s.e.m. The distance between the RDFs of the model simulations and the RDF of the data was substantially smaller for the three pairwise interactions models than for the random Poisson model (ttest: P < 10^{−300} for all three comparisons). c, Bestfitting temperatures (in log scale, mean ± s.d.) for neurons with CV < 0.25, plotted for the three pairwise interactions models. All these potentials yielded similar bestfitting temperatures: the Lennard–Jones and modified Lennard–Jones potentials yielded very similar bestfitting temperatures, and the DOG potential yielded a slightly lower temperature. The vertical dashed line denotes the temperature that was used for the bulk of the simulations in this study with the Lennard–Jones potential (log_{10}T = 0.5). We note that for the modified Lennard–Jones and DOG potentials, the match between the RDF of the data and the RDF of the model was also very good for log_{10}T = 0.5 (see b, where for all three pairwise interactions models we used log_{10}T = 0.5). d, e, Simulations for highCV cells (CV > 0.25, neurons not exhibiting an ordered arrangement of fields; n = 37 neurons). d, Distribution of the bestfitting temperatures in the Lennard–Jones pairwise interactions model (red, 3D grid cells; blue, nongrid multifield cells). The average logarithm of the bestfitting temperatures across these neurons was ~ log_{10}T = 3.5 (black dashed line). e, Quantifying the match between model and data; ttest, P values indicated. Three rightmost bars: distances between the RDF of the data and the RDFs of the three pairwise interactions models (temperatures used here: Lennard–Jones, log_{10}T = 3.5; modified Lennard–Jones, log_{10}T = 1.5; DOG, log_{10}T = 1; these were the bestfitting temperatures for each of the models; note that these temperatures correspond to a much higher noise level than in the pairwise interactions models in a–c, where we considered neurons with CV < 0.25 and fitted a temperature of log_{10}T = 0.5). Left bar, distance between the RDF of the data and the RDF of the random Poisson model. N = 1,000 simulated RDFs for each model; shown mean ± s.e.m.
Extended Data Fig. 11 Modelling: index definitions, 2D grids, and effective density.
a, Definition of the regularity index (RI), computed as the relative height of the maximum versus the minimum in the RDF (excluding its first peak); that is, the RI reflects the modulation depth of the second (or third) peak of the RDF (see Methods). b, Definition of effective density ρ_{eff}, which is defined as r_{0}/mean interfield distance (see Methods). Note that when the fields are at their equilibrium position, at distances equalling the minimal energy of the Lennard–Jones potential, the effective density is ρ_{eff} = 0.9; that is, ρ_{eff} = 0.9 is the ‘natural’ effective density. c, Two examples of 2D simulated cells, showing field positions (bottom) and autocorrelograms (top, zoomed in on the centre). Numbers above autocorrelograms indicate the classical 2D gridness score (g). The left example is also shown in Fig. 5d. d, Distribution of gridness scores in all the pairwiseinteraction 2D simulations (room size: 2 × 2 m). Many of the 2D simulations resulted in a high gridness score, indicating a highly ordered hexagonal organization. e, Plot of the classical 2D gridness score^{8,10,19} versus 2D room size, for our 2D model simulations (simulations in a room of size a × a × 10cm height; the gridness score was computed as described above (see Methods) and is shown as mean ± s.e.m.; minimum 500 simulations per room size a (the number of simulations varied because they were chosen post hoc to match the nominal ρ_{eff})). This plot is shown for fixed effective density ρ_{eff} = 0.9, which is the ‘natural’ value for ρ_{eff}, where the local distances equal the minimal energy distance r_{min} of the Lennard–Jones potential (see b; note that ρ_{eff} is defined as: ρ_{eff} = r_{0}/d_{NN}, where d_{NN} is the mean nearest neighbour distance (see Methods); because r_{min} = r_{0} × 1.12 in the Lennard–Jones function, this yields: ρ_{eff} = (r_{min}/1.12)/d_{NN} ≈ 0.9 × r_{min}/d_{NN}. That is, when d_{NN} = r_{min}, this yields ρ_{eff} = 0.9). In other words, ρ_{eff} = 0.9 is the effective density that yields interfield distances that match the minimum of the Lennard–Jones potential. We note that this effective density is also well within the empirically measured range of effective densities for our experimental data (see f). f, Distribution of effective density of fields ρ_{eff} in our experimental data, computed using the bestfitted scale parameter r_{0} for each neuron (shown here are the n = 29 multifield cells with CV < 0.25). g, Regularity index as a function of room height. Simulations were done in 10cm increments of room height, going from effectively 2D arena (0.1 m ceiling) to a full 3D room (2.5 m ceiling), as in Fig. 5e. Each curve was computed for a fixed effective density ρ_{eff} (see b); different curve shading corresponds to different ρ_{eff} values, from light to dark grey (ρ_{eff} = 0.75 to ρ_{eff} = 1). Thick curve (shown in Fig. 5e): ρ_{eff} = 0.9, a value within the range seen in our data (see f), and which is the ‘natural’ effective density corresponding to the equilibrium positions of the fields (see b). Left, computed for parameters corresponding to the 3D data, using r_{0} = 1 m and a 5 × 5 m floor (with varying ceiling height). Right: computed for parameters corresponding to 2D data, using r_{0} = 0.5 m and a 2 × 2 m floor (with varying ceiling height). Dashed line marks the value of r_{0} used. Note the drop in regularity values when going from 2D to 3D; at this transition point the room height is slightly smaller than the distance parameter r_{0} of the Lennard–Jones potential (and thus the room height allowed a second layer of fields).
Extended Data Fig. 12 3D grid cells are better explained by the pairwise interactions model than by jittered FCC or jittered HCP models; and pairwise interactions lead to a local but not global order at low temperatures of the model.
a, RDF plots for three models (top row, green: Lennard–Jones pairwise interactions model; middle and bottom rows, dark and light purple: jittered FCC and jittered HCP models, where the RDFs were computed for jittered FCC and jittered HCP models with the same nearestneighbour mean distance (d_{NN}) and CV of distances (CV_{NN}) as used in the pairwise interactions model). Each model was simulated using three temperatures (T) of the model (left, zero temperature T = 0, a nonbiological scenario with no noise; middle, low temperature (log_{10}T = −1) with little noise; right, the higher temperature that fitted well our experimental 3D data (log_{10}T = 0.5)). For the temperature fitted to the data (log_{10}T = 0.5), note the comparison of the RDF for the pairwise interactions model (top right, green) to the RDF for a random Poisson model (grey): the peak of the RDF is much sharper around r = 1 in the pairwise interactions model (green), indicating preserved local distances between nearby fields. Note that for very low temperatures with little or no jitter (left, middle), the jittered FCC and jittered HCP models showed very repetitive and distinct peaks in the RDF as a function of distance r, indicating a nearperfect lattice—as expected from FCC and HCP with no jitter (T = 0) or very small jitter (low temperature log_{10}T = −1). By contrast, for the pairwise interactions model, at low temperatures (both at log_{10}T = −1 and T = 0) the RDF was nonzero at all distances (see insets), which is not consistent with a perfect global lattice. At T = 0, the RDF seems more consistent with a distorted FCC lattice, but with preserved local distances. In other words, even at very low temperatures, the pairwise interactions model exhibited local distances without a global lattice. b, Distributions of local angles for the different models, following the same organization as a. Even at T = 0, the pairwise interactions model yielded field arrangements with a continuous angle distribution, with all possible angles (except very small angles that cannot occur for nearestneighbour triangles)—in contrast to the jittered FCC and jittered HCP arrangements where angles concentrated solely in discrete values, which indicates again that even at T = 0, the pairwise interactions model does not exhibit a perfect global lattice. Here again, at the temperature fitted to the 3D data (log_{10}T = 0.5), the pairwise interactions model yielded a more ordered arrangement than the random Poisson (top right, green curve for Lennard–Jones: see the mild peak at an angle of 60°, as compared to the grey curve for random Poisson: no peak). This means that the pairwise interactions model yielded stronger local orientational order than randomly arranged fields. c, d, Energy convergence of the Lennard–Jones pairwise interactions model at the temperature that best fitted the 3D data (log_{10}T = 0.5, orange) and at T = 0 (light blue). c, Energy per pair of fields versus the step number of the simulation. d, Left, energy change versus the step number of the simulation, plotted in linear–log scale for both temperatures. Right, energy change of the simulations for T = 0, plotted in log–log scale. Shading in c, d depicts s.e.m. Note that c, d show that for the temperature which was used to fit our model to the data (log_{10}T = 0.5, orange), the system converged to an asymptote very rapidly, long before we finished running the simulations—that is, the system rapidly reached equilibrium; therefore, the simulations reported in this study correspond to a system in equilibrium.
Extended Data Fig. 13 Comparing the pairwise interactions model for 3D versus 2D, for six combinations of model parameters T and r_{0} (varying the number of fields).
a–f, Regularity phase planes for 2D and 3D (see Extended Data Fig. 11a for illustration of the regularity index)—plotted as a function of room size, a, and number of fields, N. Room dimensions: 2D, a × a × 10cm height; 3D, roomsize a × a × a/2. The three rows (a, b and c, d and e, f) correspond to a changing temperature parameter T (c, d show the temperature that was best fitted to the experimental data: log_{10}T = 0.5). The two columns (a, c, e and b, d, f) correspond to a changing distance parameter r_{0} of the Lennard–Jones potential (values of T and r_{0} are indicated). For each panel, we compared simulations in 2D (left) and in 3D (right). The top graphs in each panel show the regularity index (grey scale) for various combinations of a and N (we conducted simulations for N ≥ 10, corresponding to multifield cells, which have ≥10 fields). Note that all these simulations exhibited much higher regularity of field arrangements in 2D than in 3D. Red hatched area denotes a nonplausible region, where the experimentally observed field diameters (~20 cm for published 2D rat grid cells^{24} and ~1 m in our 3D bat grid cells) do not allow fitting N such fields into a room of size a (see Methods). The bottom graphs show the RDF as a function of r (normalized distance) for two specific values of a, N (marked by small squares in the top graphs). Left, 2D simulation for a = 2 m, N = 30 fields (the approximate arena size and number of fields for a typical rat 2D grid cell from dorsal MEC^{24}); right, 3D simulation for a = 5 m, N = 20 fields (corresponding to our typical 3D data from bats). RDFs plotted as in Fig. 5d. Note that for all six parameter combinations (T, r_{0}; a–f), the 2D simulations showed a high regularity index (owing to a prominent second peak emerging in the RDF in 2D) for a large portion of the phase plane, and in particular in the a, N combination corresponding to a typical 2D grid cell from rat dorsalMEC^{24} (black pixel in the left panels)—whereas all 3D simulations showed a low regularity index for a, N corresponding to our data (white pixel in the right panels). A similar difference between 2D and 3D was also found when plotting phase planes of room size a versus the linear density of fields, ρ_{0} (Extended Data Fig. 14).
Extended Data Fig. 14 Comparing the pairwise interactions model for 3D versus 2D, for six combinations of model parameters T and r_{0} (varying the linear field density).
a–f, Similar phase plane plots as in Extended Data Fig. 13, comparing 2D and 3D, where now the regularity index was plotted as a function of room size a and linear field density ρ_{0} (with ρ_{0} defined as the square root of the field density in 2D and cubic root of the field density in 3D). We used r_{0} = 1 m in the right column (b, d, f) and r_{0} = 0.5 m in the left column (a, c, e), for three temperatures T (rows; c, d show the temperature that was best fitted to the experimental data: log_{10}T = 0.5). Similarly to the plots of the regularity index as a function of room size a and particle number N (Extended Data Fig. 13), these plots show the following: (i) arrangements in 2D are more ordered than in 3D (for matched density); (ii) decreasing the temperature increases the amount of order (compare the three rows); and (iii) increasing the density increases the amount of order (see monotonic dependence of the regularity index on ρ_{0} within each panel).
Supplementary information
Rights and permissions
About this article
Cite this article
Ginosar, G., Aljadeff, J., Burak, Y. et al. Locally ordered representation of 3D space in the entorhinal cortex. Nature 596, 404–409 (2021). https://doi.org/10.1038/s4158602103783x
Received:
Accepted:
Published:
Issue Date:
Further reading

Neuroscientists go wireless
Nature Methods (2021)

Disordered grids in the third dimension
Nature Neuroscience (2021)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.