Quantitative modeling of the emergence of macroscopic grid-like representations

Abstract

When subjects navigate through spatial environments, grid cells exhibit firing fields that are arranged in a triangular grid pattern. Direct recordings of grid cells from the human brain are rare. Hence, functional magnetic resonance imaging (fMRI) studies proposed an indirect measure of entorhinal grid-cell activity, quantified as hexadirectional modulation of fMRI activity as a function of the subject's movement direction. However, it remains unclear how the activity of a population of grid cells may exhibit hexadirectional modulation. Here, we use numerical simulations and analytical calculations to suggest that this hexadirectional modulation is best explained by head-direction tuning aligned to the grid axes, whereas it is not clearly supported by a bias of grid cells toward a particular phase offset. Firing-rate adaptation can result in hexadirectional modulation, but the available cellular data is insufficient to clearly support or refute this option. The magnitude of hexadirectional modulation furthermore depends considerably on the subject's navigation pattern, indicating that future fMRI studies could be designed to test which hypothesis most likely accounts for the fMRI measure of grid cells. Our findings also underline the importance of quantifying the properties of human grid cells to further elucidate how hexadirectional modulations of fMRI activity may emerge.

Keywords: entorhinal cortex; functional magnetic resonance imaging; grid cells; human; navigation; neuroscience; rodents.

Figure 1.. Qualitative hypotheses on the emergence of macroscopic grid-like representations in the human entorhinal cortex.

(A) Grid-cell properties comprise grid orientation, grid spacing, and grid phase offset. (B) Conjunctive grid by head-direction cell hypothesis. Macroscopic grid-like representations (right) may emerge from the firing of conjunctive grid by head-direction cells (left and middle) that exhibit increased firing when the subject moves aligned as compared to misaligned with the grid axes (right) (Doeller et al., 2010). (C) Repetition suppression hypothesis. In a grid cell population with similar grid orientations and grid spacings but with distributed phases (left; colored circles represent firing fields of different grid cells), aligned movements (horizontal pink arrow) lead to more frequent activation (shorter distance between firing fields) of a smaller number of different grid cells, whereas misaligned movements (diagonal gray arrow) lead to less frequent activation (larger distance between firing fields) of a higher number of different grid cells (Doeller et al., 2010). Thus, aligned movements may lead to more pronounced repetition suppression as compared to misaligned movements (middle), resulting in a hexadirectional modulation of population spiking activity and thus in the emergence of grid-like representations (right). (D) Structure-function mapping hypothesis. Because anatomically adjacent grid cells exhibit similar grid phase offsets (in addition to similar grid orientations and grid spacings) (Gu et al., 2018), recordings from a limited number of grid cells with a non-random distribution of phase offsets may lead to macroscopic grid-like representations. The left panel shows the grid phase offset of four different grid cells, whose anatomical locations are illustrated in the middle panel. Depending on the subject’s starting location relative to the phase offset of the grid fields, movements aligned or misaligned with the grid axes lead to higher sum grid cell activity as compared to misaligned or aligned movements (right panel). Furthermore, the orientation of hexadirectional modulation may shift when recording from neighboring voxels in anatomical space due to a shift in the clustered phase offsets. D, dorsal; L, lateral; M, medial; V, ventral. Figure 1 has been adapted from Figure 3 from Kunz et al., 2019.

Figure 1—figure supplement 1.. Effect of tortuosity on the hexasymmetry for the repetition suppression hypothesis with random walks (Δt=0.01s s, T=9000s, τr=3s, wr=1).

(A) The hexasymmetry ${\displaystyle |{\tilde{A}}_6|}$ (blue) and the scaled path hexasymmetry ${\displaystyle |{\tilde{T}}_6|\cdot {\tilde{A}}_0}$ (orange) as a function of the movement tortuosity ${\displaystyle {\sigma }_{\theta }}$ of a random walk. The shaded areas represent the standard error when averaging over 100 realizations of the trajectory, each with an initial direction sampled from a random uniform distribution on the interval ${\displaystyle [0,2\pi )}$. (B) The mean difference between ${\displaystyle |{\tilde{A}}_6|}$ and ${\displaystyle |{\tilde{T}}_6|\cdot {\tilde{A}}_0}$ as a function of the movement tortuosity ${\displaystyle {\sigma }_{\theta }}$. Note that ${\displaystyle |{\tilde{A}}_6|\ge |{\tilde{T}}_6|\cdot {\tilde{A}}_0}$ for ${\displaystyle {\sigma }_{\theta }\le 0.4}$. For ${\displaystyle {\sigma }_{\theta }>0.4}$ the two curves begin to overlap: repetition suppression ceases to have a significant effect and the hexasymmetry ${\displaystyle |{\tilde{A}}_6|}$ is primarily dictated by the path hexasymmetry. The dip in the mean difference near ${\displaystyle {\sigma }_{\theta }\sim 0.1}$ is thought to be due to numerical noise. (C) The mean firing rate ${\displaystyle \widetilde{A_0}}$ does not depend on movement tortuosity ${\displaystyle {\sigma }_{\theta }}$, indicating that any effect of repetition suppression on the population firing rate is small. (D) Path hexasymmetry ${\displaystyle |{\tilde{T}}_6|}$ as a function of the movement tortuosity ${\displaystyle {\sigma }_{\theta }}$. The solid line depicts the mean path hexasymmetry averaged over 100 realizations of a trajectory, while the dashed line plots the analytical result in Equation 65. (E) Leftmost panel: Five examples (colored line segments) of a piecewise linear (‘p-l’) trajectory; each linear path segment has a length of 300 cm. Three rightmost panels: Five example random-walk (‘rand’) trajectories (colored curves) for three different values of the movement tortuosity ${\displaystyle {\sigma }_{\theta }}$ with total simulation time ${\displaystyle T=60}$ s and a total path length of 600 cm. For illustration purposes, the initial directions of example random-walk trajectories were chosen such that they are regularly distributed on the interval ${\displaystyle [0,2\pi )}$. When viewed on a scale within the range of ±600 cm, increasing the tortuosity from 0.1 to 0.5 results in increasingly more curved trajectories. The bright dots in all panels of (E) show the grid fields with grid spacing 30 cm, which is small when viewed on this scale. In the simulations for Figures 2—5, the random-walk trajectories have total lengths of 90,000 cm or longer and use a movement tortuosity of ${\displaystyle {\sigma }_{\theta }=0.5{\text{rad/s}}^{1/2}}$.

Figure 1—figure supplement 2.. Effect of size and shape of finite environments on hexasymmetry.

(A) Examples of bounded random-walk trajectories (${\displaystyle \mathrm{\Delta }t=0.01}$ s, ${\displaystyle T=9000}$ s, ${\displaystyle v=10}$ cm/s, ${\displaystyle {\sigma }_{\theta }=0.5}$ rad/s^1/2; see Section ‘Environments with boundaries’ for details on the generation of bounded trajectories) within square boundaries (bottom) or circular boundaries (top) of varying sizes increasing from left to right (numbers at top labeled by ‘size’ indicate ‘half of side length’ for square boundaries or ‘radius’ for circular boundaries). Trajectories are overlaid on the firing field of an example grid cell. The lengths of the depicted trajectories are modified for illustration purposes, and do not reflect the full extent of trajectories used to calculate hexasymmetries. (B) Hexasymmetry for the three hypotheses for different sizes of the boundaries; blue, circular boundary; orange, square boundary. For all three hypotheses, the ‘ideal’ parameter sets were used; for values of the parameters, see Table 1. For the conjunctive grid by head-direction hypothesis and the repetition suppression hypothesis, the obtained values of the hexasymmetry have a weak (if any) dependence on boundary shape and size (apart from fluctuations due to noise in different realizations), and the obtained values are similar in magnitude to those obtained in infinite environments: Figure 2F for ‘conjunctive’, Figure 3E for ‘repetition suppression’. In the case of the structure-function mapping hypothesis, the strong clustering of grid phase offsets combined with the tendency of the simulated agent to move along the walls of the boundary results in a dependence of the hexasymmetry on the environment size. This effect is largest for smaller environment sizes. As the size of the environment is increased, the agent more uniformly samples the space, and the hexasymmetry approaches the magnitude of that obtained in infinite environments (Figure 4C). (C) Path hexasymmetry for different sizes of the boundaries. In (B) and (C), lines represent the mean and shaded areas represent the standard error as obtained from 50 trajectories.

Figure 1—figure supplement 3.. Effect of rotation of finite environments on hexasymmetry.

(A) Examples of bounded random-walk trajectories (${\displaystyle \mathrm{\Delta }t=0.01}$ s, ${\displaystyle T=9000}$ s, ${\displaystyle v=10}$ cm/s, ${\displaystyle {\sigma }_{\theta }=0.5}$ rad/s^1/2; see Section ‘Environments with boundaries’ for details on the generation of bounded trajectories) within circular (top) or square (bottom) boundaries for different rotation angles (numbers at top) relative to the grid orientation. Trajectories are overlaid on the firing field of an example grid cell. The lengths of the depicted trajectories are modified for illustration purposes, and do not reflect the full extent of trajectories used to calculate hexasymmetries. (B) Hexasymmetry for the three hypotheses for different rotation angles of the trajectory; blue: circular boundary; orange: square boundary. For the conjunctive grid by head-direction hypothesis and the repetition suppression hypothesis, the obtained values of the hexasymmetry have a weak (if any) dependence on boundary shape and orientation (apart from fluctuations due to noise in different realizations), and the obtained values are similar in magnitude to those obtained in infinite environments: Figure 2F for ‘conjunctive’, Figure 3E for ‘repetition suppression’. In the case of the structure-function mapping hypothesis, the strong clustering of grid phase offsets combined with the tendency of the simulated agent to move along the walls of the boundary results in a periodic dependence on the boundary orientation in square environments when the sides of the square align with the strongly clustered grid firing fields. For all three hypotheses, the ‘ideal’ parameter sets were used; for values of the parameters, see Table 1. (C) Path hexasymmetry does not depend on the rotation angle of the trajectory. In (B) and (C), lines represent the mean and shaded areas represent the standard error as obtained from 50 random-walk trajectories with the same parameters as in Table 1.

Figure 1—figure supplement 4.. Hexasymmetry of random-walk trajectories.

The (horizontal) orange dashed line shows the offset (or maximum) of the path hexasymmetry, ${\displaystyle {\tilde{T}}_0\equiv 1}$. The (diagonal) red dashed line shows the average path hexasymmetry ${\displaystyle |{\tilde{T}}_6|=1/\sqrt{M}}$ associated with randomly sampling a movement direction at each time step from a uniform distribution on the interval ${\displaystyle [0,2\pi )}$. Solid colored curves show an upper bound for the mean path hexasymmetry of a random walk as a function of the number of time steps (square root of Equation 57) for five different movement tortuosities ${\displaystyle {\sigma }_{\theta }}$ and a simulation time step size ${\displaystyle \mathrm{\Delta }t=10}$ ms. The corresponding five colored dashed lines show an approximation (${\displaystyle |{\tilde{T}}_6|\propto 1/\sqrt{M}}$, Equation 65) to the solid curves; the approximation is excellent if the number ${\displaystyle M}$ of time steps is large enough. Colored dots show the respective mean path hexasymmetries obtained from numerical simulations (Equation 1). The black arrow shows the multiplicatory shift in the number of time steps that is necessary to obtain the same hexasymmetry for a trajectory with ${\displaystyle {\sigma }_{\theta }=0.5{\text{rad/s}}^{1/2}}$ as for a trajectory with ${\displaystyle {\sigma }_{\theta }=1{\text{rad/s}}^{1/2}}$, as derived in Equation 70.

Figure 2.. Conjunctive grid by head-direction (HD) cell hypothesis.

(A) Left: The preferred head direction (red arrow) of a single conjunctive cell is aligned to one of the grid axes (Doeller et al., 2010). Two factors add noise to this relation: the HD tuning has a certain width (${\displaystyle 1/\sqrt{{\kappa }_c}}$) and the alignment of grid orientation to HD tuning angle is jittered (${\displaystyle {\sigma }_c}$). Middle: Distribution of the tuning width ${\displaystyle 1/\sqrt{{\kappa }_c}}$ around all possible grid axes for two example values. Right: Convolution (gray) between the distributions of the jitter ${\displaystyle {\sigma }_c=5^{\circ }}$ and tuning width ${\displaystyle 1/\sqrt{{\kappa }_c}\approx {8.1}^{\circ }}$ (${\displaystyle {\kappa }_c=50{\text{rad}}^{-2}}$) to obtain the effective HD tuning of a single conjunctive cell around one grid axis. (B–D) Simulation of the conjunctive hypothesis using ‘ideal’ parameters of ${\displaystyle {\kappa }_c=50{\text{rad}}^{-2}}$ and ${\displaystyle {\sigma }_c=0^{\circ }}$. The scale bars (red) represent a distance of 120 cm. (B) Left: Illustration of a star-like walk (path segments are cut for illustration purposes), overlaid onto the firing pattern of a single grid cell. Right: Population firing rate as a function of the subject’s movement direction (which is identical with heading direction in our simulations) for star-like runs with mean firing rate ${\displaystyle {\tilde{A}}_0=1279.7}$ spk/s (for 1024 cells) and path hexasymmetry ${\displaystyle |{\tilde{T}}_6|<{10}^{-10}}$; see Methods for definitions of ${\displaystyle {\tilde{A}}_0}$ and ${\displaystyle {\tilde{T}}_6}$. (C) Left: Illustration of a piecewise linear walk (cut for illustration purposes), overlaid onto the firing pattern of a single grid cell. Right: Population firing rate as a function of movement direction for piecewise linear walks with ${\displaystyle {\tilde{A}}_0=1279.3}$ spk/s and ${\displaystyle |{\tilde{T}}_6|<{10}^{-10}}$. (D) Left: Illustration of a random walk (cut for illustration purposes), overlaid onto the firing pattern of a single grid cell. Right: Population firing rate as a function of movement direction for random walks with ${\displaystyle {\tilde{A}}_0=1281.0}$ spk/s and ${\displaystyle |{\tilde{T}}_6|=6.7\cdot {10}^{-3}}$. (E) Hexasymmetry (color coded) as a function of HD tuning width and alignment jitter for star-like walk trajectories at an offset of ${\displaystyle (0,0)}$. Higher hexasymmetry values are achieved for stronger HD tuning and tighter alignment of the preferred head directions to the grid axes. The red symbols correspond to the three parameter combinations used in subplots (B–D) and (F–H) for further illustration. Large tuning widths (${\displaystyle {\kappa }_c\to 0}$) correspond to cosine tuning for which the hexasymmetry approaches 0. (F) Population firing rate as a function of movement direction for a random-walk trajectory with jitter ${\displaystyle {\sigma }_c=0}$ and concentration parameter ${\displaystyle {\kappa }_c=50{\text{rad}}^{-2}}$ (tuning width ${\displaystyle \approx {8.1}^{\circ }}$). (G) Population firing rate as a function of movement direction for a random walk with jitter ${\displaystyle {\sigma }_c={1.5}^{\circ }}$ and concentration parameter ${\displaystyle {\kappa }_c=10{\text{rad}}^{-2}}$ (tuning width ${\displaystyle \approx {18.1}^{\circ }}$). (H) Population firing rate as a function of movement direction for a random walk with jitter ${\displaystyle {\sigma }_c=3^{\circ }}$ and concentration parameter ${\displaystyle {\kappa }_c=4{\text{rad}}^{-2}}$ (tuning width ${\displaystyle \approx {28.6}^{\circ }}$). The hexasymmetry for the case of ${\displaystyle {\mathrm{p}}_{\mathrm{c}}=33}$% is stated in brackets. All other simulations presented in this figure use ${\displaystyle {\mathrm{p}}_{\mathrm{c}}=100}$% conjunctive (${\displaystyle N=1024}$) cells, which is higher than in empirical studies (Sargolini et al., 2006; Boccara et al., 2010). In subplots (B–D) and (F–H), the black solid lines and magenta dashed lines represent the results from the numerical simulations of Equation 8 and the analytical derivation in Equation 40, respectively. The radial subplots in (F–H), right, illustrate the effective HD tuning width around a single grid axis analogous to A, right. H, neural hexasymmetry; spk/s, spikes per second.

Figure 2—figure supplement 1.. Comparison between the ‘ideal’ tuning width and ‘realistic’ tuning widths found in existing literature for the conjunctive grid by head-direction (HD) cell hypothesis.

(A) Illustration of tuning curves for different tuning widths parameterized using the large-${\displaystyle {\kappa }_c}$ approximation. (B) Comparison between tuning widths using the large-${\displaystyle {\kappa }_c}$ approximation (magenta line) and the tuning widths from Sargolini et al., 2006, which use the mean vector approach (numerical: black line, analytical: blue line). The star corresponds to an angular standard deviation of ${\displaystyle {8.1}^{\circ }}{\displaystyle ({\kappa }_c=50)}$, which is the tuning width we use for the ‘ideal’ case in the text. The plus denotes an angular standard deviation of ${\displaystyle {29}^{\circ }({\kappa }_c=4)}$, which is the tuning width we use for the ‘realistic’ case in the text. The triangle marks an angular standard deviation of ${\displaystyle {55}^{\circ }}{\displaystyle ({\kappa }_c=1.7)}$, which is in the range of tuning widths reported by Sargolini et al., 2006 (the large-${\displaystyle {\kappa }_c}$ approximation does not hold in this case, and ${\displaystyle {55}^{\circ }}$ angular standard deviation would erroneously correspond to ${\displaystyle {\kappa }_c=1.1}$). ‘Mean vector’ in the legend indicates the angular standard deviation calculated using the vector summation of the mean firing rate over angular bins of ${\displaystyle 6^{\circ }}$. ‘Bessel functions’ denote the angular standard deviation calculated by using a ratio of Bessel functions instead of the resultant vector. See section ‘Tuning widths and the large-${\displaystyle \kappa }$ approximation’ for details on the calculation of the angular standard deviation. (C) Median hexasymmetry (color coded) as a function of HD tuning width (using the large-${\displaystyle {\kappa }_c}$ approximation) and alignment jitter for star-like walk trajectories at an offset of ${\displaystyle (0,0)}$. Higher hexasymmetry values are achieved for stronger HD tuning and tighter alignment of the preferred head directions to the grid axes. The star corresponds to a tuning width of ${\displaystyle {8.1}^{\circ }({\kappa }_c=50)}$. The plus represents a tuning width of ${\displaystyle {29}^{\circ }({\kappa }_c=4)}$. The triangle marks a Sargolini-like tuning width of ${\displaystyle {44}^{\circ }({\kappa }_c=1.7)}$, which corresponds to the triangle in panel B. To transform the tuning width from ${\displaystyle {55}^{\circ }}$ when using the mean vector method to a tuning width of ${\displaystyle {44}^{\circ }}$ when using the large-${\displaystyle {\kappa }_c}$ approximation, the value of ${\displaystyle {\kappa }_c}$ for a tuning width of ${\displaystyle {55}^{\circ }}$ was found numerically using the ‘mean vector’ curve in panel B. This value of ${\displaystyle {\kappa }_c}$ was then used in the expression ${\displaystyle 1/\sqrt{{\kappa }_c}}$ to calculate the tuning width in the large-${\displaystyle {\kappa }_c}$ approximation. The units of ${\displaystyle 1/\sqrt{{\kappa }_c}}$ have been transformed from radians to degrees by multiplication with a factor ${\displaystyle 180/\pi }$. The hexasymmetry was taken as a median over four trials.

Figure 2—figure supplement 2.. Dependence of the hexasymmetry on the proportion of conjunctive grid by head-direction cells.

The hexasymmetry exhibits a linear dependence on the proportion of conjunctive cells above a noise floor. A black dashed line proportional to ${\displaystyle c{p}_c}$ is plotted for comparison, where ${\displaystyle c=17}$ is an offset parameter chosen to illustrate the linear dependence of the hexasymmetry on the proportion of conjunctive cells for ${\displaystyle {\mathrm{p}}_{\mathrm{c}}>0.1}$. In random walks, the noise floor is the contribution of the path hexasymmetry to the neural hexasymmetry. For star-like walks and piecewise linear walks, the noise floor is due to sampling a finite number of grid phase offsets from a uniform distribution on the unit rhombus. The ‘realistic’ parameters were used for each trajectory (${\displaystyle {\mathrm{p}}_{\mathrm{c}}=0.33}$, ${\displaystyle {\kappa }_c=4{\text{rad}}^{-2}}$, ${\displaystyle {\sigma }_c=3^{\circ }}$). The shaded area represents the standard error of the hexasymmetry.

Figure 3.. Repetition suppression hypothesis.

(A) Left: Tuning of an example grid cell with aligned (pink arrow) and misaligned (gray arrow) movement directions. Right: Examples of firing-rate adaptation (dashed lines, adaptation weight ${\displaystyle w_r=1}$ and time constant ${\displaystyle {\tau }_r=3}$ s) for an aligned run (pink) and a misaligned run (gray). More firing-rate adaptation (i.e. stronger repetition suppression) occurs along the aligned run compared to the misaligned run (attenuation of peaks: 24% and 18% respectively). For both runs, firing rates are reduced compared to the case without adaptation (solid lines). (B) Simulations of hexasymmetry as a function of the weight (${\displaystyle w_r}$) and time constant (${\displaystyle {\tau }_r}$) of adaptation for star-like runs. The red dot marks the optimal parameters (${\displaystyle w_r=1}$ and ${\displaystyle {\tau }_r=3}$ s) used in panels A and C–E. (C) Population firing rate as a function of the subject’s movement direction for a star-like run at an offset of ${\displaystyle (0,0)}$. The solid line represents firing rates for single runs where adaptation does not carry over when sampling different movement directions, i.e., the ‘teleportation’ between path segments resets the repetition suppression effects, and movement directions are sampled consecutively from ${\displaystyle 0^{\circ }}$ to ${\displaystyle {359}^{\circ }}$ in steps of ${\displaystyle 1^{\circ }}$ (mean firing rate ${\displaystyle {\tilde{A}}_0=866.4}$ spk/s, path hexasymmetry ${\displaystyle |{\tilde{T}}_6|<{10}^{-10}}$). The dashed line represents firing rates for single runs with adaptation carry-over and randomly sampled movement directions without replacement (${\displaystyle {\tilde{A}}_0=839.8}$ spk/s, path hexasymmetry ${\displaystyle {\displaystyle |{\tilde{T}}_6|<{10}^{-10}}}$, the corresponding hexasymmetry is shown in brackets). (D) Population firing rate as a function of movement direction for a piecewise linear walk (${\displaystyle {\tilde{A}}_0=839.0}$ spk/s, ${\displaystyle |{\tilde{T}}_6|<{10}^{-10}}$). (E) Population firing rate as a function of movement direction for a random walk (${\displaystyle {\tilde{A}}_0=839.7}$ spk/s, ${\displaystyle |{\tilde{T}}_6|=6.7\cdot {10}^{-3}}$). spk/s, spikes per second. For all repetition suppression simulations, the grid phase offsets of 1024 grid cells were sampled randomly from a uniform distribution across the unit rhombus, and the hexasymmetries were averaged over 20 realizations. The scale bars (red) in panels C–E represent a distance of 120 cm.

Figure 4.. Structure-function mapping hypothesis.

(A–C) Left: Short example trajectories (white) overlaid onto the firing-rate pattern of an example grid cell (colored). Shown trajectories are for illustration purposes only, and do not reflect the full length of the simulation. The scale bars (red) represent a distance of 120 cm. Right: Population firing rates of 1024 grid cells as a function of the movement direction. Black and light gray lines represent the results from the numerical simulations of Equation 8 and the analytical derivation in Equation 32, respectively. Grid phase offsets cells are strongly clustered at ${\displaystyle (0,0)}$ with ${\displaystyle {\kappa }_s=10}$ (left insets). (A) Left: Subsegment of a star-like walk. Right: Population firing rate as a function of the movement direction for star-like runs originating at phase offset ${\displaystyle (0,0)}$ with mean firing rate ${\displaystyle {\tilde{A}}_0=1362.4}$ spk/s (spikes/second) and path hexasymmetry ${\displaystyle |{\tilde{T}}_6|<{10}^{-10}}$. (B) Left: Subsegment of a piecewise linear walk. Right: Population firing rate as a function of the movement direction for a piecewise linear trajectory. ${\displaystyle {\tilde{A}}_0=1244.3}$ spk/s, ${\displaystyle |{\tilde{T}}_6|<{10}^{-10}}$. (C) Left: Subsegment of a random walk. Right: Population firing rate as a function of movement direction for a random walk. ${\displaystyle {\tilde{A}}_0=1305.4}$ spk/s, ${\displaystyle |{\tilde{T}}_6|=6.7\cdot {10}^{-3}}$. (D) Left: Hexasymmetry as a function of the subject’s starting location (relative to grid phase offset). Right: Movement direction associated with the highest sum grid-cell activity, i.e., the phase of peaks in (A, right), has a bimodal distribution (0 or 30°). (E) Same as (D) but for piecewise linear walks. (F) Same as (D) but for random walks. (G) Example of a two-dimensional slice of a three-dimensional (3D) random-field simulation with a spatial resolution of 15 μm. The simulated volume of 3 × 3 × 3 mm³ represents the approximate spatial extent of a voxel in functional magnetic resonance imaging (fMRI) experiments (for details, see Methods). (H) The pairwise phase distance in the rhombus is shown as a function of the pairwise anatomical distance for all pairs of simulated cells from (G). Since no spatial correlation structure is induced, the pairwise phase distance between grid cells remains constant when varying the pairwise anatomical distance between them. (I) The resulting phase clustering for 200³ simulated grid cells in a 3 × 3 × 3 mm³ voxel. Brighter colors indicate a higher prevalence of a particular grid phase offset. The distribution of grid phases appears to be homogeneous with a clustering concentration parameter of ${\displaystyle {\kappa }_s=0.001}$. (J–M) Same as (G–I) but for a convolution of the 3D random field with the 3D Gaussian correlation kernel of width 30 μm shown in (K). Grid cells located next to each other in anatomical space (${\displaystyle \lesssim 30\text{µm}}$) exhibit similar grid phase offsets (L). No clear clustering is visible, and the fitted clustering concentration parameter was ${\displaystyle {\kappa }_s=0.014}$ (M). (N–Q) Same as (J–M) but for the grid-like correlation kernel with projections shown in (O), which adds some longer-range spatial autocorrelation to the grid phase offsets of different grid cells. The red bars in (O) measures a distance of 300 μm, which corresponds to the separation between two peaks of the correlation kernel lying on the same plane. The pairwise phase distance (P) exhibits a dip around 300 μm. Note in (Q) that the prevalence of particular grid phase offsets is more biased than in (M), with a fitted clustering concentration parameter of ${\displaystyle {\kappa }_s=0.052}$. (R) Dependence of the clustering concentration parameter ${\displaystyle {\kappa }_s}$ on the number ${\displaystyle N}$ of grid cells in a voxel. A random distribution of grid cells in anatomical space was obtained by subsampling from the 200³ grid cells simulated in (G–Q) over 300 realizations. We found ${\displaystyle {\kappa }_s\approx 0.05}$ for 10³ grid cells in a voxel, and that ${\displaystyle {\kappa }_s}$ decreases monotonically as ${\displaystyle N}$ is increased. Convolution of grid phase offsets with a correlation kernel (as in J–Q) leads to saturation of ${\displaystyle {\kappa }_s}$ for large ${\displaystyle N}$. Note that the range of values of ${\displaystyle {\kappa }_s}$ here is three orders of magnitude smaller than the strong clustering considered in (A–C). The dashed-dotted line depicts the line ${\displaystyle 1/\sqrt{N}}$ for comparison. The vertical dotted line at ${\displaystyle N={20}^3}$ corresponds to the empirically estimated count of grid cells in a 3 × 3 × 3 mm³ fMRI voxel. The secondary lower horizontal axis shows the average distance between regularly distributed grid cells in anatomical space. (S) The distribution of the clustering concentration parameter ${\displaystyle {\kappa }_s}$ when using either no kernel, a Gaussian kernel, or a grid-like kernel for 20³ grid cells subsampled from the 200³ grid cells simulated in (G–Q) over 300 realizations. The grid-like kernel results in a larger maximum value of the clustering concentration parameter, and a large number of realizations results in relatively low clustering. ***, p<0.001. (T) The dependence of the hexasymmetry ${\displaystyle |{\tilde{A}}_6|}$ on the number of angles sampled when unwrapping the star for the piecewise linear walk, averaged over 20 trajectories for each data point. Each additional angle sampled adds 300 cm to the total length of the path. A line proportional to ${\displaystyle c/\sqrt{M}}$ is plotted for comparison, where ${\displaystyle c}$ is an offset parameter (${\displaystyle c=4000}$ spk/s) chosen such that the slope of ${\displaystyle |{\tilde{A}}_6|}$ can be compared to the slope of ${\displaystyle 1/\sqrt{M}}$, which is proportional to the path hexasymmetry ${\displaystyle |{\tilde{T}}_6|}$ (Equation 67). The close fit between the solid and dashed lines indicates that the neural hexasymmetry ${\displaystyle |{\tilde{A}}_6|}$ is highly correlated with the path hexasymmetry ${\displaystyle |{\tilde{T}}_6|}$. The vertical dotted line at 360 sampled angles corresponds to the number of angles sampled in the star-like run and the piecewise linear run for all main figures. The value of the clustering concentration parameter used was ${\displaystyle {\kappa }_s=10}$.

Figure 4—figure supplement 1.. Comparison of hexasymmetry resulting from the standard grid cell model (‘standard’) and the structure-function mapping hypothesis (‘Gu’, identical to Figure 5).

Both hypotheses are implemented for three different types of trajectories: star-like walks (‘star’), piecewise linear walks (‘p-l’), and random walks with small step size (‘rand’). For each setting, we show the average neural hexasymmetry ${\displaystyle |{\tilde{A}}_6|}$ for 1024 cells (dark blue bars) and the average contribution of the trajectory ${\displaystyle |{\tilde{T}}_6|\cdot {\tilde{A}}_0}$ (light blue bars) where ${\displaystyle {\tilde{A}}_0}$ is the (variable) average population activity and ${\displaystyle |{\tilde{T}}_6|}$ is the hexasymmetry of the trajectory; see Methods for definitions of symbols. The violin plots depict the distributions for firing-rate hexasymmetry (red) and path hexasymmetry (orange). For the ‘standard’ grid cell model, we implemented the structure-function mapping hypothesis, but each realization of the grid phase offsets was sampled from a two-dimensional von Mises distribution with a clustering concentration parameter ${\displaystyle {\kappa }_s=0}$ (equivalent to a uniform distribution on the unit rhombus). In this scenario, random fluctuations give rise to a sample mean clustering concentration parameter of ${\displaystyle {\hat{\kappa }}_s=0.04}$. The sample mean clustering concentration parameter in this case depends on the number of grid cells in a voxel (Figure 4R). Distributing the grid phases regularly throughout the unit rhombus results in ${\displaystyle {\hat{\kappa }}_s=0}$ and reduces the mean neural hexasymmetry (star: ${\displaystyle 5.6\cdot {10}^{-11}}$ spk/s; p-l: ${\displaystyle 4.8\cdot {10}^{-11}}$ spk/s). Using a constant function (such that each grid cell has the same firing rate for any position of the agent in space) in place of the grid-like firing fields (Equation 2) while sampling the grid phase offsets from a uniform distribution on the unit rhombus similarly reduces the mean neural hexasymmetry (star: ${\displaystyle 3.6\cdot {10}^{-11}}$ spk/s; p-l: ${\displaystyle 3.4\cdot {10}^{-11}}$ spk/s). Either regularly distributing the grid phase offsets or using a constant grid function results in a hexasymmetry that is similar in magnitude to the path hexasymmetry. Each hypothesis condition was simulated for 300 realizations. For ‘Gu’ we used ${\displaystyle {\kappa }_s=0.1}$; ***, ${\displaystyle \mathrm{p}<0.001}$; n.s., not significant.

Figure 5.. Comparison of hexasymmetry resulting from the three hypotheses.

Each of the three hypotheses is implemented for three different types of navigation trajectories: star-like walks (‘star’), piecewise linear walks (‘p-l’), and random walks with small step size (‘rand’). For each setting, we show the average neural hexasymmetry ${\displaystyle |{\tilde{A}}_6|}$ for 1024 cells (dark blue bars) and the average contribution of the trajectory ${\displaystyle |{\tilde{T}}_6|\cdot {\tilde{A}}_0}$ (light blue bars) where ${\displaystyle {\tilde{A}}_0}$ is the (variable) average population activity and ${\displaystyle |{\tilde{T}}_6|}$ is the hexasymmetry of the trajectory; see Methods for definitions of symbols. The violin plots depict the distributions for firing-rate hexasymmetry (red) and path hexasymmetry (orange). For each hypothesis, we calculate the hexasymmetry for ‘ideal’ parameters (conjunctive: ${\displaystyle {\mathrm{p}}_{\mathrm{c}}=1}$, ${\displaystyle {\kappa }_c=50{\text{rad}}^{-2}}$, tuning width ${\displaystyle 1/\sqrt{{\kappa }_c}={8.1}^{\circ }}$, ${\displaystyle {\sigma }_c=0}$; repetition suppression: ${\displaystyle {\tau }_r=3}$ s, ${\displaystyle w_r=1}$; clustering: ${\displaystyle {\kappa }_s=10}$) as well as more realistic parameters (conjunctive ‘Doeller’: ${\displaystyle {\mathrm{p}}_{\mathrm{c}}=0.33}$, ${\displaystyle {\kappa }_c=4{\text{rad}}^{-2}}$, tuning width ${\displaystyle 1/\sqrt{{\kappa }_c}={29}^{\circ }}$, ${\displaystyle {\sigma }_c=3^{\circ }}$ motivated by Doeller et al., 2010; Boccara et al., 2010; Sargolini et al., 2006; repetition suppression ‘weaker’: ${\displaystyle {\tau }_r=1.5}$ s, ${\displaystyle w_r=0.5}$; clustering ‘Gu’: ${\displaystyle {\kappa }_s=0.1}$ motivated by Gu et al., 2018). For clustering, we also simulated grid phase offsets randomly sampled from a uniform distribution with ${\displaystyle {\kappa }_s=0}$ and found that neural hexasymmetries were reduced (compared to ‘Gu’) by ${\displaystyle \sim 30\mathrm{\%}}$ (Figure 4—figure supplement 1). The parameters for the random-walk scenario are ${\displaystyle T=9000}$ s and ${\displaystyle \mathrm{\Delta }t=0.01}$ s; see Table 1 for a summary of parameters. Each hypothesis condition was simulated for 300 realizations. ***, p<0.001; n.s., not significant; (***), a seemingly significant result (p<0.001) that we think is spurious (see Results section). For pairwise comparisons of the hexasymmetry values from different trajectory types for each set of parameters, see Figure 5—figure supplement 1. For a quantification of path hexasymmetry as a function of the percentage of subsampled path segments, see Figure 5—figure supplement 2. For a comparison between our method and previously used methods for evaluating hexasymmetry, see Figure 5—figure supplement 4.

Figure 5—figure supplement 1.. Pairwise comparisons of the hexasymmetry values from different trajectory types for each set of parameters.

The compared trajectory types are star-like walks (‘star’), piecewise-linear walks (‘p-l’), and random walks (‘rand’). For each hypothesis, we calculate the hexasymmetry for ideal parameters (conjunctive: ${\displaystyle {\mathrm{p}}_{\mathrm{c}}=1}$, ${\displaystyle {\kappa }_c=50{\text{rad}}^{-2}}$; ${\displaystyle {\sigma }_c=0}$, repetition suppression: ${\displaystyle {\tau }_r=3}$ s, ${\displaystyle w_r=1}$; clustering: ${\displaystyle {\kappa }_s=10}$) as well as more realistic parameters (conjunctive: ${\displaystyle {\mathrm{p}}_{\mathrm{c}}=0.33}$, ${\displaystyle {\kappa }_c=4{\text{rad}}^{-2}}$; ${\displaystyle {\sigma }_c=3^{\circ }}$ repetition suppression: ${\displaystyle {\tau }_r=1.5}$ s, ${\displaystyle w_r=0.5}$; clustering: ${\displaystyle {\kappa }_s=0.1}$). In the case of the repetition suppression hypothesis, the solid blue bars show the star-like walk with a carry-over of the repetition suppression mechanism when teleporting between different path segments, while the transparent bars with the dotted blue borders represent the star-like walk with no carry-over of the repetition suppression effect across different path segments. For the star-like walk, the starting phase of the star is sampled from a uniform distribution across the unit rhombus between realizations, and remains constant within each realization of the star-like walk trajectory. The direction of movement for both the star-like walk and the piecewise linear walk is sampled randomly without replacement from the integer angles ${\displaystyle \{0,1,2,...,359{\}}^{\circ }}$. The parameters for the random-walk scenario are ${\displaystyle T=9000}$ s and ${\displaystyle \mathrm{\Delta }t=0.01}$ s. Each hypothesis condition was simulated for 300 realizations. Mann-Whitney U test; ***, p<0.001; **, p<0.01; n.s., not significant. Note that the scales of the horizontal axes are different across subpanels.

Figure 5—figure supplement 2.. Hexasymmetry resulting from piecewise linear walks (‘p-l’) and random walks (‘rand’) as a function of the percentage of subsampled path segments for the structure-function mapping hypothesis.

Values of ‘ideal’ parameters (left) and more realistic parameter values (‘Gu’; right) are identical to the ones used in Figure 5. The thinner gradient bars show the percentage of subsampled path segments required to produce the corresponding scaled path hexasymmetry ${\displaystyle |{\tilde{T}}_6|\cdot {\tilde{A}}_0}$. The thicker dark blue and light blue bars represent the hexasymmetry ${\displaystyle |{\tilde{A}}_6|}$ and the path hexasymmetry multiplied by the mean firing rate ${\displaystyle |{\tilde{T}}_6|\cdot {\tilde{A}}_0}$, respectively. The percentages of path segments were subsampled from 1% to 100% in steps of 1%. In the case of the random-walk trajectories, e.g., subsampling 100% of the path segments yields ${\displaystyle |{\tilde{T}}_6|\cdot {\tilde{A}}_0=|{\tilde{A}}_6|=7.92}$ spk/s. For each percentage value, the path hexasymmetry was averaged over 100 different realizations of the piecewise linear or random-walk trajectory.

Figure 5—figure supplement 3.. Comparison of the hexasymmetry resulting from the three hypotheses when using trajectories from rats and humans.

For each hypothesis condition, we used 97 rat trajectories (Sargolini et al., 2006) and 181 human trajectories across two datasets (Kunz et al., 2015, Kunz et al., 2021) to estimate neural hexasymmetry and path hexasymmetry. Example trajectories are overlaid on the firing field of a grid cell (bottom panels). The human trajectories were rescaled such that the subject traversed three to four grid firing fields (Hafting et al., 2005) when crossing the environment. For each dataset, the corresponding neural hexasymmetry and path hexasymmetry were evaluated. In the top panel we show for each setting the average neural hexasymmetry ${\displaystyle |{\tilde{A}}_6|}$ for 1024 cells (dark blue bars) and the average contribution of the trajectory ${\displaystyle |{\tilde{T}}_6|\cdot {\tilde{A}}_0}$ (light blue bars) where ${\displaystyle {\tilde{A}}_0}$ is the (variable) average population activity; see Methods for definitions of symbols. The violin plots depict the distributions for firing-rate hexasymmetry (red) and path hexasymmetry (orange). While the human trajectories do not result in significant neural hexasymmetry under the repetition suppression hypothesis in this data, we note that there is a variety of different virtual-reality navigation tasks in humans with varying degrees of straightness. For example, in Horner et al., 2016, the navigation paths very clearly resemble piecewise-linear walks. In such studies, repetition suppression effects may play a role in the emergence of grid-like representations. For each hypothesis, we calculated the hexasymmetry for ‘ideal’ parameters (conjunctive:${\displaystyle {\mathrm{p}}_{\mathrm{c}}=1}$, ${\displaystyle {\kappa }_c=50{\text{rad}}^{-2}}$; ${\displaystyle {\sigma }_c=0}$; repetition suppression: ${\displaystyle {\tau }_r=3}$ s, ${\displaystyle w_r=1}$; clustering: ${\displaystyle {\kappa }_s=10}$). ***, p<0.001; n.s., not significant.

Figure 5—figure supplement 4.. Comparison of hexasymmetry resulting from the three hypotheses when using different methods to calculate the hexasymmetry.

Analogous to Figure 5, each of the three hypotheses (‘conjunctive’, ‘repetition suppression’, ‘structure-function mapping’) is implemented for three different types of trajectories: star-like walks (‘star’), piecewise linear walks (‘p-l’), and random walks with small step size (‘rand’). For each trajectory type within each hypothesis condition, we show four violin plots with the color corresponding to hexasymmetry calculated using different methods: a generalized linear model (‘GLM’, Doeller et al., 2010), a generalized linear model with binning (‘GLM (binning)’, Kunz et al., 2015), our new method (‘New method’), and a circular-linear correlation (‘CL-corr’, Maidenbaum et al., 2018). The mean path hexasymmetry is indicated by a black square. The black star marks the mean hexasymmetry of a surrogate distribution obtained by circularly shifting the firing rates. In all cases, ‘ideal’ parameters were used (conjunctive: ${\displaystyle {\mathrm{p}}_{\mathrm{c}}=1}$, ${\displaystyle {\kappa }_c=50{\text{rad}}^{-2}}$; ${\displaystyle {\sigma }_c=0}$; repetition suppression: ${\displaystyle {\tau }_r=3}$ s, ${\displaystyle w_r=1}$; clustering: ${\displaystyle {\kappa }_s=10}$). In order to visualize the GLM methods, which can take on negative values, we use the ‘symlog’ scale provided by pyplot, with a linear threshold of ±10 spk/s. Each hypothesis condition for each method was simulated for 300 realizations. Mann-Whitney U test; ***, p<0.001; **, p<0.01; n.s., not significant; for further details, see Methods (Implementation of previously used metrics). For the conjunctive grid by head-direction cell hypothesis, all of the methods find significant neural hexasymmetry for all three types of trajectories, although the circular-linear correlation finds a statistically less significant effect (${\displaystyle \mathrm{p}<0.01}$ for circular-linear correlation, all ${\displaystyle \mathrm{p}<0.001}$ for other methods). In the case of the repetition suppression hypothesis, all methods find significant neural hexasymmetry for star-like walks and piecewise linear random walks (${\displaystyle \mathrm{p}<0.01}$ for circular-linear correlation, all ${\displaystyle \mathrm{p}<0.001}$ for the other methods), but find no significant neural hexasymmetry for random walks. With respect to the structure-function mapping hypothesis, all of the methods find significant neural hexasymmetry for star-like walks (${\displaystyle \mathrm{p}<0.01}$ for circular-linear correlation, all ${\displaystyle \mathrm{p}<0.001}$ for the other methods). For piecewise linear walks and random walks, the new methods detect a spurious significant neural hexasymmetry (labeled ‘(***)’). In section ‘Structure-function mapping hypothesis’, we explain that this spurious effect is due to an inhomogenous sampling of space. All other methods find no significant hexasymmetry under the structure-function mapping hypothesis for piecewise linear walks and random walks. On a side note, when there is significant hexasymmetry, the GLM (binning) method used in Kunz et al., 2015, shows an attenuation in hexasymmetry of roughly 60% due to the averaging effect of binning firing rates into aligned and misaligned movement directions, but this reduction does not change the statistical outcomes. Overall, these observations align with the results obtained through our new method, as depicted in Figure 5.

Author response image 1.. Variability of peak firing rates of grid cells when using a random walk trajectory in a finite environment.

(A) Left: The distribution of the peak firing rates for the two conditions (repetition-suppression hypothesis: orange; without repetition-suppression: blue) for 1024 grid cells across 300 simulations and a simulated duration of 9000 seconds, which is kept constant for all panels of (A). Right: The distribution for the coefficient of variation (CV) of the peak firing rates across all cells. The dashed lines represent distribution means. (B) The dependence of the mean CV on the simulated duration for 1024 grid cells across 300 simulations. As the simulated duration increases, the difference in CV between the two cases decreases. The shaded area represents the standard error. (C) Four examples of the rate map for different values of the simulated duration. In the rightmost panel, the trajectory is not shown so that the rate map can be seen clearly. We employed the 'ideal’ parameter set for repetition suppression: τr = 3 s, wr = 1; in all simulations, cells did not show any head direction-tuning, and there was no clustering of grid-phase offsets. The size of the environment was limited to a square with side lengths of 100 cm. Mann-Whitney U test; ***, P < 0.001; **, P < 0.01; *, P < 0.05; n.s., not significant.