Quantitative modeling of the emergence of macroscopic gridlike 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 gridcell 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 headdirection tuning aligned to the grid axes, whereas it is not clearly supported by a bias of grid cells toward a particular phase offset. Firingrate 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.
Editor's evaluation
This computational work represents a valuable and long overdue assessment of the potential mechanisms associating patterns of activity of entorhinal grid cells, recorded mostly in rodents, with the population property of hexasymmetry detected in noninvasive human studies. The methodic comparison of alternative hypotheses is compelling, and the conclusions are important for the future design of experiments assessing the neural correlates of human navigation across physical, virtual, or conceptual spaces.
https://doi.org/10.7554/eLife.85742.sa0Introduction
The neural basis of spatial navigation comprises multiple specialized cell types such as place cells (O’Keefe and Dostrovsky, 1971), headdirection (HD) cells (Taube et al., 1990), and grid cells (Hafting et al., 2005), whose activity profiles result from intricate mechanisms of microcircuits in the medial temporal lobes (Tukker et al., 2022). Grid cells are neurons that activate whenever an animal or human traverses the vertices of a triangular grid tiling the entire environment into equilateral triangles (Hafting et al., 2005; Jacobs et al., 2013). Grid cells may allow the navigating organism to perform vector computations and may thus constitute an essential neural substrate for different types of spatial navigation including path integration, though their exact functional role still remains unclear (Stemmler et al., 2015; Bush et al., 2015; Moser et al., 2017; Stangl et al., 2018; Gil et al., 2018; Banino et al., 2018; Bierbrauer et al., 2020; Ginosar et al., 2023).
In rodents, grid cells can be recorded using electrodes inserted into the medial entorhinal cortex (EC). In humans, measuring grid cells using invasive methods is only rarely possible, e.g., by recording singleneuron activity in epilepsy patients who are neurosurgically implanted with intracranial depth electrodes (Jacobs et al., 2013; Nadasdy et al., 2017). Hence, to enable the detection of grid cells in healthy humans, a functional magnetic resonance imaging (fMRI) method has been developed that tests for a hexadirectional modulation of the bloodoxygenleveldependent signal as a function of the subject’s movement direction through a virtual environment (Doeller et al., 2010). We here refer to this phenomenon of a hexadirectional modulation of the fMRI signal as ‘macroscopic gridlike representations’, which has been replicated repeatedly in recent years (e.g. Kunz et al., 2015; Bellmund et al., 2016; Horner et al., 2016; Constantinescu et al., 2016; Bierbrauer et al., 2020). The mechanisms underlying the emergence of such macroscopic gridlike representations remain still unclear, however.
To provide possible explanations for the emergence of macroscopic gridlike representations, previous studies presented several qualitatively different hypotheses on how the activity of single grid cells translates into a macroscopically visible hexadirectional fMRI signal (Doeller et al., 2010; Kunz et al., 2019). Three main hypotheses have been developed: (i) the ‘conjunctive grid by HD cell hypothesis’; (ii) the ‘repetition suppression hypothesis’; and (iii) the ‘structurefunction mapping hypothesis’ (Figure 1).
The conjunctive grid by HD cell hypothesis rests on two key findings: first, the existence of conjunctive grid by HD cells, which were found in the deeper layers of the entorhinal cortex and in pre and parasubiculum (Sargolini et al., 2006); second, the observation that the directional tuning of grid cells within the entorhinal cortex is aligned with the grid axes (Doeller et al., 2010), though further studies are needed to corroborate this observation. Assuming that the directional tuning width of these conjunctive grid by HD cells is not too broad, movements aligned with the grid axes (as compared to misaligned movements) result in increased spiking activity of the conjunctive grid by HD cell population. Given some correlation between population spiking activity and the fMRI signal, this systematic difference in the firing of conjunctive grid by HD cells when moving aligned versus misaligned with the grid axes may thus cause a macroscopically visible fMRI signal with hexadirectional modulation (Figure 1B). We note that the conjunctive grid by HD cell hypothesis does not necessarily depend on an alignment between the preferred head directions with the grid axes.
The repetition suppression hypothesis (Figure 1C) is based on the assumption that the phenomenon of repetition suppression—i.e., neural activity being reduced for repeated stimuli (GrillSpector et al., 2006)—also applies to grid cells (Doeller et al., 2010; Killian et al., 2012). Critical to this hypothesis is that relatively fewer different grid cells are activated more often during movements aligned with the grid axes, and relatively more different grid cells are activated less often during misaligned movements. Due to this systematic difference in how many grid cells are activated how often, a higher degree of repetition suppression at the level of spiking activity or the fMRI signal (i.e. fMRI adaptation) during aligned movements as compared with misaligned movements can emerge, again resulting in a hexadirectional modulation of fMRI activity as a function of the subject’s movement direction through the spatial environment.
Regarding the structurefunction mapping hypothesis (Figure 1D), studies in rodents have demonstrated that the firing fields of anatomically adjacent grid cells do not only have similar spacing and orientation (Stensola et al., 2012), but also a similar grid phase offset to a reference location (Heys et al., 2014; Gu et al., 2018). Studies in rodents observed that grid cells also cluster anatomically (Obenhaus et al., 2022; Naumann et al., 2018), which may lead to stronger gridlike representations in some voxels (with an increased percentage of grid cells) than in others (with a decreased percentage of grid cells). Therefore, recordings from a small area of the EC (e.g. a sufficiently small voxel of an fMRI scan) may sample grid cells with similar firing fields, which basically behave similarly. It has been suggested that such a grid cell population might show higher average firing rates during aligned movements (because more firing fields are traversed) versus misaligned movements, again resulting in macroscopically visible gridlike representations (Kunz et al., 2019).
In this study, we aimed at quantitatively evaluating the three hypotheses on the emergence of macroscopic gridlike representations using a modeling approach. Our results show that all three hypotheses can result in macroscopic gridlike representations under ideal and specific conditions, but that the strength of the hexadirectional modulation varies by orders of magnitude. Key findings are also that the subjects’ type of navigation paths through the spatial environments and the exact biological characteristics of grid cells determine to what extent a given hypothesis can explain a hexadirectional population signal in the EC. In this way, our results help understand how grid cells may have a specific correlate in fMRI, make predictions on how future fMRI studies could establish evidence in favor of or against one of the three hypotheses, and suggest that the biological properties of grid cells in humans should be investigated in greater detail in order to support or weaken the plausibility of either of the three hypotheses.
Results
In what follows, we first describe the navigation strategies that we use in our model, then define a new measure to quantify neural hexasymmetry and path hexasymmetry, and finally calculate neural hexasymmetries for the three hypotheses, which we model for idealized as well as more realistic choices of parameters.
Navigation strategies
To evaluate the different hypotheses on the emergence of gridlike representations, we considered three different types of navigation trajectories: ‘starlike walks’, ‘piecewise linear walks’, and ‘random walks’. We opted for this approach in order to test whether a subject’s navigation pattern—which in itself comes with a certain degree of hexasymmetry (‘path hexasymmetry’)—influences the emergence of hexadirectional sum signals of neuronal activity.
During each path segment of starlike walks, the simulated agent started from the same (x/y)coordinate and navigated along 1 of 360 predefined allocentric navigation directions (0° to 359° in steps of 1°; Figure 2B, left). This ensured that the navigation trajectory itself exhibited a hexasymmetry that was essentially zero. Each path had a length of 300 cm, which was 10 times the grid scale (see Table 1 for a summary of all model parameters). After each path segment was traversed, the agent was ‘teleported’ back to the initial (x/y)coordinate and completed the next path segment. For realworld experiments, this type of navigation including teleportation is unusual, but it can be implemented in virtualreality experiments (Vass et al., 2016; Deuker et al., 2016).
During piecewise linear walks, the subject also completed 360 path segments of 300 cm length along the same 360 predefined allocentric navigation directions, as in the starlike walks. In this case, however, the path segments were ‘unwrapped’ such that the starting location of a path segment was identical with the end location of the preceding path segment (Figure 2C, left). The sequence of allocentric navigation directions was randomly chosen. As for starlike walks, piecewise linear walks do not exhibit hexasymmetry a priori.
For random walks, we modeled navigation trajectories following Kropff and Treves, 2008; Si et al., 2012; D’Albis and Kempter, 2017 (for details, see Methods around Equation 1), which allowed us to vary the tortuosity of the paths. For a certain value of the tortuosity parameter (${\sigma}_{\theta}=0.5$ rad/s^{1/2}) and a time step $\mathrm{\Delta}t=0.01$ s, this led to navigation paths that we considered similar to those seen in rodent studies (Figure 2D, left; Figure 1—figure supplement 1E). Tortuosity describes how convoluted a navigation path is (the higher the tortuosity, the more convoluted the path). A value of 0.5 rad/s^{1/2} means that in 1 s of time the standard deviation of movement direction change is about 0.5 rad = 29° (see Figure 1—figure supplement 1E, and Equation 42 for $m\mathrm{\Delta}t=1$ s).
Apart from random walks in basically infinite environments, we also simulated random walks in finite enclosures (circles and squares) with different sizes and orientations; we found that these restrictions had a negligible effect on path hexasymmetry (Figure 1—figure supplement 2C and Figure 1—figure supplement 3C; for a definition of path hexasymmetry, see the paragraph after the next one). Because the allocentric navigation directions are not predefined for random walks, they exhibit varying degrees of path hexasymmetry. The longer the simulated random walks, the smaller the path hexasymmetry (Figure 1—figure supplement 4). We simulated random walks with a total length of typically 900 m ($M=9\cdot {10}^{5}$ steps).
As we describe below, the emergence of gridlike representations based on the conjunctive hypothesis is robust against the specific type of navigation strategy, whereas the other two hypotheses are sensitive to particular navigation strategies. Future studies on hexadirectional signals should thus consider the kind of navigation paths subjects will use during a given task.
Quantifying neural hexasymmetry generated by the three hypotheses
To test how the activity of grid cells could give rise to hexasymmetry of a macroscopic signal, we used a firingrate model of grid cell activity (Equation 2). Furthermore, we developed a new measure $H$ to quantify neural hexasymmetry (see Methods, Equation 12), which is the magnitude of the hexadirectional modulation of the summed activity of many cells. The hexasymmetry $H$ has a value $H=0$ if there is no hexadirectional modulation, e.g., if the population firing rate does not depend on the movement direction. Conversely, if the population firing rate has a hexadirectional sinusoidal modulation, half of its amplitude (in units of the firing rate) equals the value of the hexasymmetry $H$. Using the same approach, we quantified the hexasymmetry of a trajectory and called this path hexasymmetry (see Methods after Equation 14).
Conjunctive grid by HD cell hypothesis
The conjunctive grid by HD cell hypothesis (Doeller et al., 2010) suggests that hexadirectional activity in the EC emerges due to grid cells whose firing rate is additionally modulated by head direction, whereby the preferred head direction is closely aligned with one of the grid axes (Figure 2A). By modulating the activity of individual grid cells with an HD tuning term aligned with the grid axes (Methods, Equation 3), our simulations indeed showed that these properties resulted in a clear hexadirectional modulation of sum gridcell activity (Figure 2B–D). When considering different types of navigation trajectories, we found that they led to similar distributions of sum gridcell activity as a function of movement direction and, accordingly, to similar hexasymmetry values (Figure 2B–D). In all three cases, the directions of maximal activity were aligned with the grid axes.
These results were obtained using ideal values for the preciseness of the HD tuning (i.e. large concentration parameter of HD tuning, $\kappa}_{c$) and for the alignment of the preferred head directions to the grid axes (i.e. small alignment jitter of the HD tuning to the grid axes, $\sigma}_{c$; Figure 2A). We were thus curious how the hexasymmetry changed when using a wide range of parameter values that would also include biologically plausible values. We varied $\kappa}_{c$ between values corresponding to narrow tuning widths ($\kappa}_{c}=50\text{}{\text{rad}}^{2$, which corresponds to an angular variability of $1/\sqrt{{\kappa}_{c}}\approx {8}^{\circ}$) and wide tuning widths ($\kappa}_{c}=4\text{}{\text{rad}}^{2$, i.e. an angular variability of approximately $29}^{\circ$), and we varied $\sigma}_{c$ between values of no jitter (${\sigma}_{c}=0$) and significant jitter ($\sigma}_{c}={3}^{\circ$); see also Table 1 for a summary of parameters. We found that a combination of narrow HD tuning widths and no jitter resulted in the largest hexasymmetry $H$ (Figure 2E–H), while wider tuning widths with nonzero jitter resulted in smaller values for the hexasymmetry (for larger tuning widths and a comparison to experimental values, see Figure 2—figure supplement 1).
Repetition suppression hypothesis
Next, we performed simulations to understand whether the repetition suppression hypothesis (Doeller et al., 2010) results in a hexadirectional modulation of population gridcell activity. This hypothesis proposes that gridcell activity is subject to firingrate adaptation and thus leads to reduced gridcell activity when moving along the grid axes as compared to when moving along other directions than the grid axes (Figure 3A). This difference is due to the fact that when the subject moves along the grid axes the grid fields of fewer grid cells are traversed relatively more often (strong repetition suppression), whereas when moving not along the grid axes the grid fields of more grid cells are traversed relatively less often (weak repetition suppression) (Doeller et al., 2010).
In our model, the repetition suppression hypothesis depends on two adaptation parameters: the adaptation time constant $\tau}_{r$ and the adaptation weight $w}_{r$ (Equation 5). We explored a large range of adaptation time constants and found that the time constant that leads to the largest hexasymmetry is roughly the subject’s speed $v$ divided by the grid scale $s$ (Figure 3B). We constrained values of the adaptation weight to the full range of reasonable values ($0<{w}_{r}\le 1$) and found that the larger the value of the adaptation weight the larger is the hexasymmetry (Figure 3B). When examining how the different types of navigation trajectories affected hexadirectional modulations based on repetition suppression, we found that starlike and piecewise linear walks resulted in clear and significant hexasymmetry values (Figure 3C and D), which is driven by the long linear segments in these trajectory types. In contrast, random walks did not result in a significant hexadirectional modulation of sum gridcell activity because the tortuosity of the random walk that we typically used (${\sigma}_{\theta}=0.5$ rad/s^{1/2}) is too large (i.e. trajectories turn too much between two adjacent firing fields of a grid cell) to be able to exploit the movementdirection dependence of repetition suppression (Figure 3E). Examining in more detail which tortuosity values would still lead to some hexadirectional modulation due to repetition suppression, we found that ${\sigma}_{\theta}\lesssim 0.25$ rad/s^{1/2} was the upper bound. For smaller values of the tortuosity parameter trajectories are straight enough to allow for a hexadirectional modulation of sum gridcell activity (Figure 1—figure supplement 1).
A notable difference of the repetition suppression hypothesis compared to the conjunctive grid by HD cell hypothesis is that the apparent preferred grid orientation (i.e. the movement directions resulting in the highest sum gridcell activity) is shifted by $30}^{\circ$ and is thus exactly misaligned with the grid axes of the individual grid cells (Figure 3C and D). This is due to the fact that the adaptation mechanism suppresses gridcell activity more strongly when moving aligned with a grid axis as compared to when moving misaligned with a grid axis (Figure 3A).
Structurefunction mapping hypothesis
We next investigated the structurefunction mapping hypothesis, according to which a hexadirectional modulation of EC activity emerges in situations when a population of grid cells is recorded whose grid phase offsets are biased toward a particular offset (Kunz et al., 2019). In the ideal case, all grid phase offsets are identical, and thus all grid cells behave like a single grid cell. This hypothesis is called ‘structurefunction mapping hypothesis’ because of a direct mapping between the anatomical locations of the grid cells in the EC and their functional firing fields in space.
We found indeed that highly clustered grid phase offsets (‘ideal’ concentration parameter for clustering ${\kappa}_{s}=10$) resulted in significant hexadirectional modulations of sum gridcell activity when the subject performed starlike walks starting at a phase offset of $(0,0)$, i.e., the center of the cluster of firing fields of grid cells (Figure 4A). Interestingly, the hexasymmetry values during starlike walks were strongly dependent on the subject’s starting location relative to the locations of the grid fields as only particular starting phases (within the unit rhombus of grid phase offsets) such as $(0,0)$ or $(0.3,0.3)$ led to clear hexasymmetry whereas others, e.g., $(0.6,0)$, did not (Figure 4D, left). Additionally, the ‘apparent preferred grid orientation’ (i.e. the movement directions associated with the highest sum gridcell activity) was shifted by $30}^{\circ$ for certain offsets in the unit rhombus illustrating the subject’s starting locations relative to the firingfield locations (Figure 4D, right). It was furthermore notable that the summed gridcell activity as a function of movement direction exhibited relatively sharp peaks at multiples of $60}^{\circ$ with additional small peaks in between (Figure 4A, right). This pattern is clearly distinct from the more sinusoidal modulation of sum gridcell activity resulting from the conjunctive grid by HD cell hypothesis (Figure 2) and the more fullwave rectified sinusoidal modulation resulting from the repetition suppression hypothesis (Figure 3).
When examining the structurefunction mapping hypothesis for piecewise linear walks and random walks, hexasymmetry appeared to be considerably lower as compared to simulations with starlike walks (Figure 4B and C). For piecewise linear walks (but not for random walks), population gridcell activity as a function of movement direction again exhibited sharp peaks at multiples of $60}^{\circ$, in both positive and negative directions (Figure 4B). For both piecewise linear walks and random walks, the hexasymmetry values did not show a systematic dependency on the starting location of the subject’s navigation trajectory relative to the locations of the grid fields (Figure 4E and F, left), which is due to the fact that these navigation trajectories randomize the starting locations of all path segments. Accordingly, the apparent preferred grid orientations varied randomly as a function of the starting location of the very first path segment and did not exhibit systematic shifts (Figure 4E and F, right), which is in contrast to the clear shifts in the apparent preferred grid orientations for starlike walks (Figure 4D).
In the abovementioned simulations for the structurefunction mapping hypothesis, we chose high values for the clustering of the grid phase offsets from all grid cells. Specifically, we set the clustering concentration parameter to ${\kappa}_{s}=10$, due to which the centers of the firing fields of different grid cells were very close to each other (insets in Figure 4A–C, right). However, previous empirical studies (Gu et al., 2018; Heys et al., 2014) suggest that the clustering is much weaker. We thus performed numerical simulations to obtain a more realistic estimation of the clustering concentration parameter, $\kappa}_{s$. We observed that our numerically determined clustering concentration parameter that matched the empirical clustering was in fact closer to randomly distributed grid phase offsets (${\kappa}_{s}<0.1$; Figure 4G–S) than to strong, ideal clustering $({\kappa}_{s}=10)$. Specifically, we implemented a ‘simple’ clustering using a Gaussian kernel that led to a clustering concentration parameter ${\kappa}_{s}\approx 0.014$ (Figure 4M). The clustering concentration parameter remained low when we added an empirically guided ‘higherorder’ spatial clustering of grid phase offsets by adding some longerrange spatial autocorrelations to the gridphase offsets (Figure 4N–Q; Gu et al., 2018), which led to a clustering concentration parameter ${\kappa}_{s}\approx 0.052$ (Figure 4Q). Accordingly, we found that these more ‘realistic’ clustering concentration parameters resulted in substantially reduced hexasymmetries for the structurefunction mapping hypothesis (described in the next section and in Figure 5). Furthermore, the neuronal hexaysymmetries were much lower than the ones for the ‘realistic’ conjunctive grid by HD cell hypothesis and for the ‘realistic’ repetition suppression hypothesis.
Our abovementioned results for the structurefunction mapping hypothesis indicate that the navigation pattern has a strong influence on the neuronal hexasymmetry $H$ (Figure 4A–C). We thus wondered how much the hexasymmetries of the associated navigation paths contributed to the neuronal hexasymmetries. We first considered random walks and found that the path hexasymmetry was proportional to $1/\sqrt{M}$ for a large number $M$ of steps (Figure 1—figure supplement 4). For the random walk in Figure 4C with tortuosity ${\sigma}_{\theta}=0.5$ rad/s^{1/2}, we derived from Figure 1—figure supplement 4 that the expected path hexasymmetry for 9000 s simulation time ($M=0.9\cdot {10}^{6}$ steps for $\mathrm{\Delta}t=0.01$ s) is about 0.007, which results in a contribution to the neural hexasymmetry of $H\approx 7$ spk/s for a population rate of about 1024 spk/s. This number is close to the obtained neural hexasymmetry ($H=8.0$ spk/s) in Figure 4C. We thus believe that for random walks the hexasymmetry $H$ obtained in the structurefunction mapping case is mainly determined by the path hexasymmetry.
We then turned to piecewise linear walks and examined neural hexasymmetries $H$ across a broad range of total trajectory distances. Surprisingly, we observed that $H$ also decreased as $1/\sqrt{M}$ (Figure 4T) even though the path hexasymmetry of piecewise linear walks is basically zero by construction. This indicates that the apparent neural hexasymmetry $H$ of sum gridcell activity for piecewise linear walks (e.g. in Figure 4B) is driven by random subsamples of all path segments—specifically those path segments crossing through grid fields. These subsamples of path segments necessarily exhibit higher path hexasymmetries than the full set of path segments that has zero path hexasymmetry. We thus conclude that, for piecewise linear walks, hexasymmetry values were driven by a subsampling of the movement directions due to the sparsity of the gridfield locations.
Taken together, the structurefunction mapping hypothesis with very strong (‘ideal’) clustering produces neural hexasymmetry values that are larger than expected from path hexasymmetries only with respect to starlike walks, including a strong dependence on the subject’s starting location (Figure 4D, left). This range of hexasymmetry values is comparable to those of the repetition suppression hypothesis with ‘ideal’ parameters (Figure 3), but values are at least an order of magnitude smaller than in the ‘ideal’ conjunctive grid by HD cell case (Figure 2). When using realistic clustering concentration parameters, neural hexasymmetry values are further reduced as compared to simulations with ideal clustering (Figure 5).
Overall evaluation of the three hypotheses
To provide a systematic evaluation of the three hypotheses, we computed 300 realizations of each hypothesis (using the simulated activity of 1024 cells), separately for each type of navigation and for both ideal and more realistic parameter settings (see Table 1 for a summary of parameter values and section ‘Parameter estimation’ for their justification). This resulted in 18 different hypothesis conditions (Figure 5; see also Figure 5—figure supplement 1 for pairwise comparisons). For each hypothesis condition, we assessed its statistical significance by performing nonparametric MannWhitney U tests between the neural hexasymmetries ($H:={\stackrel{~}{A}}_{6}$; see also Equation 12) and the product ${\stackrel{~}{T}}_{6}\cdot {\stackrel{~}{A}}_{0}$ with the multipliers path hexasymmetry ${\stackrel{~}{T}}_{6}$ and average population activity $\stackrel{~}{A}}_{0$.
For the conjunctive hypothesis, we found that all three types of navigation led to significant neural hexasymmetries. This was the case for both ideal parameters (MannWhitney U tests, all U=0, all p<0.001) and for more realistic parameters (MannWhitney U tests, all U=0, all p<0.001). For starlike walks and piecewise linear walks, the path hexasymmetries multiplied with the average population activities were near zero because we designed these navigation paths to equally cover all different movement directions. For random walks, the path hexasymmetries multiplied with the average population activities were at values of about 7 spk/s and showed larger variability because the directions of the random walks were not predefined and could thus vary with regard to their hexasymmetries. If the conjunctive hypothesis was true, fMRI studies should thus see a hexadirectional modulation of entorhinal fMRI activity for all three types of the subjects’ navigation paths.
For the repetition suppression hypothesis, we found that starlike walks and piecewise linear walks resulted in significant neural hexasymmetries for both ideal and weaker parameters (MannWhitney U tests, all U=0, all p<0.001), whereas random walks did not (MannWhitney U tests, both U>1753, both p>0.405). For random walks, the path hexasymmetries multiplied with the average population activities were lower as compared to the other two hypotheses because the repetition suppression necessarily leads to lower average population activities. If the repetition suppression hypothesis was true, fMRI studies should thus observe significant neural hexasymmetries only for starlike walks and piecewise linear walks with long enough segments, whereas random walks (with large enough tortuosities, Figure 1—figure supplement 1) should not lead to significant neural hexasymmetries.
Regarding the structurefunction mapping hypothesis, the statistical tests showed that most of the hypothesis conditions resulted in significant neural hexasymmetries as compared to the path hexasymmetries multiplied with the average population activities. This was the case for both ideal parameters (MannWhitney U tests, all U<736, all p<0.001) and more realistic parameters (MannWhitney U tests for starlike walks and piecewise linear walks, all U=0, all p<0.001). Only the hypothesis condition with random walks and more realistic parameters exhibited no significant neural hexasymmetries (MannWhitney U test, U=1786, p=0.472). However, these results regarding the structurefunction mapping hypothesis should be treated with great caution. First, the neural hexasymmetries for starlike walks heavily depend on the starting location of the subject relative to the grid fields, and different starting locations lead to different apparent grid orientations (Figure 4D). Second, the significant results for the navigation conditions with piecewise linear walks and random walks actually result from an inhomogeneous sampling of movement directions through the grid fields and therefore do not reflect true neural hexasymmetries (Figure 4T; see also Figure 5—figure supplement 2 for path hexasymmetry as a function of the number of subsampled path segments). Only the structurefunction mapping hypothesis is susceptible to this effect because the grid fields are clustered at similar spatial locations (whereas the grid fields are homogeneously distributed in the case of the conjunctive hypothesis and the repetition suppression hypothesis). In simulations with infinitely long paths, the neural hexasymmetries (for the navigation types of piecewise linear walks and random walks) would not be significantly higher than the path hexasymmetries multiplied with the average population activities. In empirical studies, this effect can be detected by correlating the subjectspecific path distances with the subjectspecific neural hexasymmetries: if there is a generally negative relationship, this will hint at the fact that the neural hexasymmetries are basically due to relevant path hexasymmetries of path segments crossing the grid fields. In essence, we therefore believe that the structurefunction mapping hypothesis leads to true neural hexasymmetry only in the case of starlike walks.
Influence of other factors
Our simulations above were performed in an infinite spatial environment, which is different from empirical studies in which subjects navigate finite environments. We were thus curious whether the size and shape of finite environments could affect the strength of hexadirectional modulations of population gridcell activity.
Our simulations showed that for both circular and square environments, hexasymmetry strengths did not considerably depend on the size of the environment when the subject performed random walks (Figure 1—figure supplement 2). Similarly, rotating the navigation trajectories relative to the grid patterns did not affect the hexasymmetry strengths (Figure 1—figure supplement 3). These results thus suggest that experiments in animals and humans can use various types and sizes of the environments to investigate hexadirectional modulations of sum gridcell activity.
Discussion
We performed numerical simulations and analytical estimations to examine how the activity of grid cells could potentially lead to a neural population signal in the EC. Such a neural population signal has been observed in multiple fMRI studies including Doeller et al., 2010; Kunz et al., 2015; Constantinescu et al., 2016; Horner et al., 2016; Bellmund et al., 2016; Stangl et al., 2018; Nau et al., 2018; Julian et al., 2018; Bierbrauer et al., 2020; Julian and Doeller, 2021; Bongioanni et al., 2021; Moon et al., 2022, and in multiple iEEG/MEG studies including Maidenbaum et al., 2018; Chen et al., 2018; Staudigl et al., 2018; Chen et al., 2021; Wang and Wang, 2021; Convertino et al., 2023. The recorded population signals showed a hexadirectional modulation of the entorhinal fMRI/iEEG signal as a function of the subject’s movement direction through its spatial environment. We note though that some studies did not find evidence for neural hexasymmetry. For example, a surface EEG study with participants ‘navigating’ through an abstract vowel space did not observe hexasymmetry in the EEG signal as a function of the participants’ movement direction through vowel space (Kaya et al., 2020). Another fMRI study did not find evidence for gridlike representations in the ventromedial prefrontal cortex while participants performed valuebased decisionmaking (Lee et al., 2021). This raises the question whether the detection of macroscopic gridlike representations is limited to some recording techniques (e.g. fMRI and iEEG but not surface EEG) and to what extent they are present in different tasks.
Potential mechanisms underlying hexadirectional population signals in the EC
We examined three hypotheses that have been previously suggested as potential mechanisms underlying the emergence of the hexadirectional population signal in the EC (Doeller et al., 2010; Kunz et al., 2019) and found that all three hypotheses can—in principle and in ideal situations—lead to a hexadirectional modulation of EC population activity. However, none of the three hypotheses described here may be true and another mechanism may explain macroscopic gridlike representations. This includes the possibility that neural hexasymmetry is completely unrelated to gridcell activity, previously summarized as the ‘independence hypothesis’ (Kunz et al., 2019). For example, a population of HD cells whose preferred head directions occur at offsets of 60° from each other could result in neural hexasymmetry in the absence of grid cells. The conjunctive grid by HD cell hypothesis thus also works without grid cells, which may explain why gridlike representations have been observed (using fMRI) in regions outside the EC (Doeller et al., 2010; Constantinescu et al., 2016), where rodent studies have not yet identified grid cells except for some evidence of grid cells in the primary somatosensory cortex of foraging rats (Long and Zhang, 2021). In that case, however, another mechanism would be needed that could explain why the preferred head directions of different HD cells occur at multiples of 60°. Attractornetwork structures may be involved in such a mechanism, but this remains speculative at the current stage.
Navigation paths have a major influence on the hexadirectional population signal
A major observation of this study is that the way how subjects navigate through the environment has a major influence on whether a hexadirectional population signal can be observed. We distinguished three major types of navigation: navigation with randomwalk trajectories in which straight paths are quite short, which resembles the navigation pattern in rodents; navigation with piecewise linear trajectories in which the subject navigates along straight paths combined with random sharp turns between the straight segments; and navigation with starlike trajectories in which the subject starts each path from a fixed center location and navigates along a straight path with a predetermined allocentric direction and distance. Critically, we found that the conjunctive hypothesis leads to a hexadirectional population signal irrespective of the specific type of navigation; that the repetition suppression hypothesis leads to hexadirectional population signals only in the case of starlike trajectories and piecewise linear trajectories (but not for random trajectories); and that for the structurefunction mapping hypothesis ‘true’ hexadirectional population signals can only be observed for starlike trajectories.
The finding that the type of navigation paths influences whether a hypothesis can lead to hexadirectional population signals of the EC is informative to future fMRI/iEEG studies, which could empirically evaluate which of the three hypotheses is most likely to be true: By asking or requiring the subjects to navigate in different ways through the task environments (in a starlike fashion, in a piecewise linear fashion, and in a random fashion), these future fMRI/iEEG studies could demonstrate whether hexadirectional population signals are present in all three navigation conditions (speaking in favor of the conjunctive hypothesis); whether they are present only during starlike or piecewise linear trajectories (in favor of the repetition suppression hypothesis); or whether they are mainly visible for starlike trajectories and exhibit a systematic decrease with increasing total trajectory length for piecewise linear and random walks (in this case speaking in favor of the structurefunction mapping hypothesis). We note that in our simulations, each linear path segment of a piecewise linear walk has a length that is 10 times larger than the grid scale. While humans performing navigation tasks exhibit straighter trajectories than those of rats (Doeller et al., 2010; Kunz et al., 2015; Horner et al., 2016), they are still more similar in scale to random walks than piecewise linear walks. Thus, when using empirical navigation paths from previous studies, neural hexasymmetries were significant only for the conjunctive hypothesis (Figure 5—figure supplement 3). In contrast to this major effect of navigation type on the presence of hexadirectional signals, the size and shape of the environment did not influence the strength of the hexadirectional signals in a relevant manner (Figure 1—figure supplement 2 and Figure 1—figure supplement 3).
The structurefunction mapping hypothesis predicts that in fMRI studies involving starlike walks, there may be up to 30° shifts in the orientation of hexadirectional modulation between neighboring voxels (Figure 4D). This is in contrast with the other two hypotheses, where the preferred grid orientations are similar between voxels, and may differ only slightly when recording distinct grid modules with different grid orientations (Stensola et al., 2012).
A note on our choice of the values of model parameters
Another insight of this study is that the exact biological properties of grid cells play a major role regarding the question whether hexadirectional population signals can be observed. For example, the conjunctive hypothesis cannot lead to hexadirectional population signals if the tuning width of the conjunctive grid by HD cells is too broad (Figure 2E). Here, the percentage of strongly tuned conjunctive cells also plays a relevant role. Empirical studies in rodents found a wide range of tuning widths among grid cells ranging from broad to narrow (Doeller et al., 2010; Sargolini et al., 2006). The percentage of conjunctive cells in the EC with a sufficiently narrow tuning may thus be low. Such distributions (with a proportionally small amount of narrowly tuned conjunctive cells) lead to low values in the absolute hexasymmetry. The neural hexasymmetry in this case would be driven by the subset of cells with sufficiently narrow tuning widths. If this causes the neural hexasymmetry to drop below noise levels, the statistical evaluation of this hypothesis would change.
We also found that hexadirectional population signals emerge only if the preferred directions of the conjunctive cells are aligned precisely enough with the grid axes of the grid cells (Figure 2E). Furthermore, no hexadirectional population activity will emerge if the preferred head directions of conjunctive grid by HD cells exhibit other types of rotational symmetry such as 10fold rotational symmetry reported in recordings from rats (Keinath, 2016). Additionally, while we assumed that all conjunctive grid cells maintain the same preferred head direction between different firing fields, conjunctive grid cells have also been shown to exhibit different preferred head directions in different firing fields (Gerlei et al., 2020). This could lead to hexadirectional modulation if the different preferred head directions are offset by 60° from each other, but will not give rise to hexadirectional modulation if the preferred head directions are randomly distributed. To the best of our knowledge, the distribution of preferred head directions was not quantified by Gerlei et al., 2020, thus this remains an open question. Together, it currently seems possible that grid cells and conjunctive grid by HD cells meet the necessary biological properties for the conjunctive hypothesis (Doeller et al., 2010), but further studies on the characteristics of conjunctive grid by HD cells (also in humans) will be useful because they may find tuning widths which are too wide or tuning curves which are not aligned with the grid axes.
The exact biological properties of grid cells also play a major role for the structurefunction mapping hypothesis, which heavily relies on the property of neighboring grid cells to share a similar grid phase offset (i.e. a high spatial autocorrelation of grid phases) and whether there might be longerrange spatial autocorrelations between the grid phases (Heys et al., 2014; Gu et al., 2018). With the evidence at hand, it seems unlikely to us that the grid phases are clustered strongly enough to facilitate a hexadirectional population signal. Specifically, we found that the clustering concentration parameters of grid phase offsets gleaned from existing empirical studies (Gu et al., 2018) produced the smallest neural hexasymmetry when compared to the other conditions we considered, while the neural hexasymmetries obtained from simulations randomly sampling grid phase offsets from a uniform distribution on the unit rhombus (that one might refer to as the ‘standard grid cell model’) were only slightly smaller.
Regarding realistic parameters for the repetition suppression hypothesis, we are currently not aware of detailed measurements of the relevant gridcell properties (i.e. the adaptation time constant $\tau}_{r$ and the adaptation weight $w}_{r$) so that it remains unclear to us to what extent the repetition suppression hypothesis is biologically plausible. Future studies may thus quantify the relevant properties of (human) grid cells in greater depth to help clarify which hypothesis regarding the emergence of hexadirectional population signals may most likely be true.
Influence of grid scale on hexasymmetry
In this study, we assumed that all grid cells have the same grid spacing of 30 cm. For rats, this is at the lower end of grid spacing values, corresponding to the dorsal region of the medial EC (Stensola et al., 2012). How would the three hypotheses perform for larger values of grid spacing? The conjunctive hypothesis would remain unchanged as it does not depend on grid spacing. In contrast, the repetition suppression hypothesis depends on the grid scale. Here, the strongest hexasymmetry is achieved when the adaptation time constant is similar to the average time that it takes the subject to move between neighboring grid fields (Figure 3B). Altering the grid scale changes this traversal time between grid fields and thus hexasymmetry. In Figure 3, we chose the ‘ideal’ adaptation time constant to maximize hexasymmetry. A change in grid scale would thus reduce hexasymmetry for this particular value of the adaptation time constant. As there are different gridcell modules in the EC that exhibit different grid spacings (Stensola et al., 2012), repetition suppression effects may vary across modules (if the adaptation time constant is not tuned in accordance with the grid scale). Regarding the structurefunction mapping hypothesis and starlike walks, different grid spacings can have a major effect on neural hexasymmetry if the path segments are short compared to the grid scale, similar to effects of the starting location (Figure 4D); for longer path segments, altering the grid scale should have a smaller effect on neural hexasymmetry. For the structurefunction mapping hypothesis, random walks and piecewiselinear walks do not lead to real neural hexasymmetries, and different grid spacings are thus irrelevant. Overall, different grid spacings are thus most important in the context of the repetition suppression hypothesis.
Another factor that might complicate the picture is an interaction of different grid scales: rodent grid cells have been found to be organized into four to five discrete modules with different grid scales (Stensola et al., 2012) and this may also be true in humans. Together with the anatomical dorsoventral length of the human EC (≈2 cm, Behuet et al., 2021), this means that a module might cover roughly 4–5 mm. A typical voxel of 3 × 3 × 3 mm^{3} would thus either contain one or two grid modules. The case of a voxel containing a single module was already discussed. If a voxel contained two modules, there would be two subpopulations of grid cells with different grid spacings: For the structurefunction mapping hypothesis, significant neural hexasymmetry can be achieved only for starlike trajectories. In this scenario, the hexasymmetry depends on the starting position of the trajectories, i.e., the center position of the star (Figure 4D). Two subpopulations of grid cells from different modules essentially act like two individual grid cells in the context of the structurefunction mapping hypothesis. Thus, they can either add their hexasymmetry values (e.g. when the star center is in the middle of a grid field for both subpopulations) or suppress each other (when the star center is in a regime of 0° grid phase for one subpopulation and in a regime of ±30° grid phase for the other). For the repetition suppression hypothesis, two populations of grid cells with the same adaptation time constant but different grid spacing would result in two different values of hexasymmetry, depending on the interaction between gridfield traversal time and adaptation time constant. The joint hexasymmetry value would then depend on the two subpopulation sizes. We note that the adaptation time constant might be tuned to match the grid spacing. This would result in high hexasymmetry values for both subpopulations, and thus also for the full population.
Sources of noise
The neural hexasymmetry as defined by Equation 12 also contains contributions from the path hexasymmetry of the underlying trajectory. Starlike walks and piecewise linear walks have a path hexasymmetry of zero by construction. Therefore, these trajectory types do not contribute to the neural hexasymmetry. For random walks, the path hexasymmetry depends on the number of time steps $M$ and the tortuosity parameter $\sigma}_{\theta$ (Figure 1—figure supplement 4). To compare the path hexasymmetry to the neural hexasymmetry, the path hexasymmetry needs to be multiplied by the timeaveraged population activity. Since the path hexasymmetry ranges from 0 to 1 (Equation 14), its contribution to the neural hexasymmetry ranges from 0 (for a uniform sampling of movement directions) to the timeaveraged population activity. For M=900,000 time steps and a tortuosity of ${\sigma}_{\theta}=0.5$ rad/s^{1/2}, randomwalk trajectories contribute ∼0.8% of the timeaveraged population activity to the neural hexasymmetry (Figure 1—figure supplement 4). In our simulations with an average population activity in the range of 10^{3} spk/s, the 0.8% noise floor leads to a path hexasymmetry of about 8 spk/s (Figure 5).
Another factor that contributes noise is the finite sampling of grid phase offsets. For all our simulations of the conjunctive grid by HD cell hypothesis and the repetition suppression hypothesis (and for control simulations of the structurefunction mapping hypothesis), we sample a finite number of grid phase offsets from a twodimensional von Mises distribution with a clustering concentration parameter ${\kappa}_{s}=0$ (equivalent to a uniform distribution on the unit rhombus). In this scenario, random fluctuations give rise to realizations with a nonzero sample mean clustering concentration parameter $\hat{\kappa}}_{s$, which led to a corresponding mean neural hexasymmetry of 0.7 spk/s (labeled ‘standard’ and ‘star’ in Figure 4—figure supplement 1), even in the absence of conjunctive tuning and repetition suppression. For 1024 grid cells, this contribution is 10 times smaller than the neural hexasymmetry resulting from the ‘realistic’ parameters for both the conjunctive grid by HD cell hypothesis and the repetition suppression hypothesis (Figure 5). In the case of the structurefunction mapping hypothesis, these fluctuations would not be considered as ‘noise’ as they directly contribute to the hexasymmetry from the clustering of grid phase offsets.
On the relation of firing rates and fMRI/iEEG signals
A topic that this study did not investigate is the question of how the sum signal of single neurons translates into fMRI and iEEG signals. In neocortical regions such as the auditory cortex, a clear linear relationship between singleneuron activity and fMRI activity has been observed (Mukamel et al., 2005), but it remains elusive whether this linear relationship also applies to the EC in general and to entorhinal grid cells in particular. In the neighboring hippocampus, for example, the relationship between singleneuron activity and the fMRI signal is highly complex (Ekstrom, 2010; Kunz et al., 2019). Future studies are needed to detail the relationship between singleneuron firing and fMRI/iEEG signals in the EC. This would allow us to clarify whether a hexadirectional modulation of sum gridcell activity directly results in a hexadirectional modulation of fMRI/iEEG activity or whether currently unknown factors modulate the expression of hexadirectional fMRI/iEEG signals.
Conclusion
Using numerical simulations and analytical derivations we showed that a hexadirectional neural population signal can emerge from the activity of grid cells given the ideal conditions of three different hypotheses. Whether a given hypothesis leads to a hexadirectional population signal is significantly influenced by the subjects’ type of navigation through the spatial environment and by the exact biological properties of human grid cells.
Methods
Trajectory modeling
To describe gridcell activity as a function of time $t$ during which a subject (animal or human) is exploring an environment, we model three distinct trajectory types. We first describe trajectory types in environments without bounds, which are quasiinfinite, and then add rules that account for boundaries.
Environments without boundaries
The first trajectory type is a random walk ${X}_{t}=[{x}_{t},{y}_{t}]$, which is defined by
with $\theta}_{t}={\sigma}_{\theta}\cdot {W}_{t$, where $\sigma}_{\theta$ controls the tortuosity of the trajectory and $W}_{t$ is a standard Wiener process. In a numerical simulation with a time step $\mathrm{\Delta}t$, the angle is updated in each time step by ${\theta}_{t+\mathrm{\Delta}t}={\theta}_{t}+{\sigma}_{\theta}\cdot \mathrm{\Delta}W$, where $\mathrm{\Delta}W$ is a normally distributed random variable with variance $\mathrm{\Delta}t$. The variable $v$ depicts the (constant) speed.
The second type of navigation is a starlike walk, where the subject moves radially outward from a predefined origin in space at a certain angle $\theta$ on a straight line to a maximum distance $r}_{max$ at a constant speed $v$. In simulations, this movement is repeated (with the same predefined origin) for $N}_{\theta$ angles that are equally spaced on the interval $[0,2\pi )$. Within each individual radial path, the subject does not turn around and move back to the origin, i.e., the entire trajectory of $N}_{\theta$ radial paths is not continuous.
Finally, we introduce a piecewise linear walk, which is constructed by placing all the radial paths of the starlike walk endtoend such that they form one single continuous trajectory of length $N}_{\theta}\phantom{\rule{thinmathspace}{0ex}}{r}_{\mathrm{m}\mathrm{a}\mathrm{x}$. The trajectory thus consists of successive straight runs for the simulated subject, which can be interpreted as a random walk with a time step $\mathrm{\Delta}t={r}_{\mathrm{m}\mathrm{a}\mathrm{x}}/v$ and directions that are sampled uniformly without replacement from a predetermined set of angles. In comparison to the random walk and the starlike walk, this procedure presumably reflects the situation in human virtualreality setups most closely, as participants often move along straight trajectories with intermittent turns (Doeller et al., 2010; Kunz et al., 2015; Horner et al., 2016).
Environments with boundaries
Most virtualreality studies in humans use finite instead of infinite spatial environments to examine gridlike representations. We wondered whether the size and the shape of these finite environments might modulate the strength of macroscopic gridlike representations obtained through one or more of the three hypotheses. Hence, we performed simulations not only in infinite but also in finite environments with a given size (between one and six times the grid spacing) and shape (circle and square).
For randomwalk trajectories, we enforce that the navigating subject stays within the circular or square environment by performing an ‘outofbounds check’ at each time point. This means that, after every time step $\mathrm{\Delta}t$, we measure the distance that the subject has moved outside of the boundary. This is done differently in square and circular environments, both of which are centered at the origin $(x=0,y=0)$. In the square environment, we define the variables $\mathrm{\Delta}x$ and $\mathrm{\Delta}y$ as the distance the subject has moved out of bounds in the $x$ and $y$ coordinates, respectively. Let $L$ be the half of the length of the square’s sides. $\mathrm{\Delta}x$ and $\mathrm{\Delta}y$ are then defined as $\mathrm{\Delta}x=\text{max}[xL,0]$ and $\mathrm{\Delta}y=\text{max}[yL,0]$. For circular environments, let $R$ be the radius of the circle, and let $\overrightarrow{r}=(x,y)$ be the position vector of the subject. We then introduce the measure $\mathrm{\Delta}r:=\text{max}[\overrightarrow{r}R,0]$, such that $\mathrm{\Delta}r$ is nonzero only when the subject has moved outside of the circular boundary. If at any time point either $\mathrm{\Delta}x$, $\mathrm{\Delta}y$, or $\mathrm{\Delta}r$ are nonzero, the outofbounds check fails. In this case, we reject the movement in this time step and keep resampling a new angle $\theta}_{t$ (Kropff and Treves, 2008; Si et al., 2012) until the check succeeds, meaning that the subject has made a move that is within the boundaries of the environment. If $\mathrm{\Delta}x$, $\mathrm{\Delta}y$, or $\mathrm{\Delta}r$ remain nonzero for 50 consecutive samples of $\theta}_{t$, we temporarily increase the tortuosity $\sigma}_{\theta$ by a factor 1.1. Without this increase in tortuosity, the subject tends to get stuck when approaching the boundary at angles close to the perpendicular or at the corners of the square boundary. The tortuosity is reset to its initial value once a valid move is made. We visually checked the randomwalk trajectories, which show some oversampling along the boundaries, and found that they were comparable to the navigation trajectories in rodent studies (e.g. Hafting et al., 2005).
Implementation of grid cell activity
Following Burgess et al., 2007, the activity profile $G}_{i$ of grid cell i (for $i=1,...,N$ in a population of $N$ grid cells) is modeled as the product of three cosine waves rotated by 60° ($=\pi /3$) from each other:
where $A}_{\text{max}$ is the grid cell’s maximal firing rate, $s}_{i$ depicts the cell’s grid spatial scale (‘grid spacing’), $x}_{\text{off},i$ and $y}_{\text{off},i$ are the phase offsets of the grid (‘grid phase’) in the two spatial dimensions (called $x$ and $y$ here), and $\gamma}_{i$ is the orientation of the grid (‘grid orientation’); see Table 1 for numerical values of parameters and Figure 1A for an illustration of the three grid characteristics.
To describe the activity of many grid cells (e.g. in a voxel for an MRI scan), we sum up the firing rates of $N$ grid cells in Equation 2. For a given trajectory ${X}_{t}=[{x}_{t},{y}_{t}]$, the macroscopic activity as a function of time $t$ is then simply described by the sum $\sum _{i=1}^{N}{G}_{i}({x}_{t},{y}_{t})$.
Implementation of the three hypotheses to explain macroscopic gridlike representations
Here, we summarize how the activity in a population of $N$ grid cells can be described if they also exhibit (i) HD tuning, (ii) repetition suppression (i.e. firingrate adaptation), or (iii) grid phases that are clustered across grid cells.
In all our models, the activity $G}_{i$ of a grid cell is described by Equation 2. Different grid cells typically have different phase offsets ($x}_{\text{off},i$, $y}_{\text{off},i$) but the same grid spacing ${s}_{i}:=s,\phantom{\rule{thickmathspace}{0ex}}\mathrm{\forall}i$ and grid orientation ${\gamma}_{i}:=\gamma ,\phantom{\rule{thickmathspace}{0ex}}\mathrm{\forall}i$ (Hafting et al., 2005; Boccara et al., 2010; Stensola et al., 2012; Gardner et al., 2022). We assumed that all grid cells have the same grid spacing of 30 cm.
Conjunctive grid by HD cell hypothesis
To include HD tuning in our model, we note that a given trajectory $X}_{t$ has an angle $\theta}_{t$ at time $t$. The summed firing rate, i.e., the population activity, $A}^{c$ from $N$ such conjunctive cells can then be described by
where the upper index ‘c’ indicates ‘conjunctive’ and where we incorporate conjunctive gridHD tuning via the (scaled by factor $2\pi$) von Mises distribution
with concentration parameter $\kappa}_{c$ and preferred angle $\mu}_{c$. The symbol $I$ represents the modified Bessel function of the first kind of order 0. The parameter $\kappa}_{c$ describes the width of the HD tuning: if $\kappa}_{c}\u27f6\mathrm{\infty$, the HD tuning is sharpest; the smaller $\kappa}_{c$, the wider the HD tuning (see Figure 2); for ${\kappa}_{c}=0$, there is no HD tuning, and our scaling leads to $h(\theta )\equiv 1$. We choose the preferred angle as ${\mu}_{c}=\frac{\pi}{3}k+\eta$, where $k$ is randomly drawn from $\{0,1,2,3,4,5\}$ and $\eta$ is randomly drawn from a normal distribution with mean 0 and standard deviation $\sigma}_{c$. For ${\sigma}_{c}=0$, the directional tuning is thus centered around a multiple of 60°. The parameter $\sigma}_{c$ introduces jitter in the alignment of directional tuning to one of the grid axes.
We modeled the cases in which all grid cells show HD tuning (‘ideal’ case, fraction of conjunctive cells: ${\mathrm{p}}_{\mathrm{c}}=1$) as well as a more ‘realistic’ case in which only a third of the cells is conjunctive (${\mathrm{p}}_{\mathrm{c}}=1/3$; Boccara et al., 2010; Sargolini et al., 2006). We note that this is an approximation, since the proportion of conjunctive cells is highly variable across layers of the EC, with up to 90% conjunctive cells in layer V.
Repetition suppression hypothesis
To incorporate repetition suppression in the model, we add an explicit dependence of gridcell activity on time $t$. Specifically, we subject the firing rate $G}_{i$ of a grid cell to an adaptation mechanism:
where $a$ depicts the adaptation variable, and $\tau}_{r$ and $w}_{r$ are the repetition suppression time constant and the weight of the suppression, respectively. The upper index ‘r’ in $G}_{i}^{r$ indicates ‘repetition suppression’.
The adaptation time constant $\tau}_{r$ is on the order of seconds. We restrict the adaptation weight $w}_{r$ to the interval $[0,1]$ (Figure 3B). The maximum operation ‘$max$’ in Equation 5 ensures that the output firing rate $G}_{i}^{r$ is always positive. Together, the summed firing rate $A}^{r$ from $N$ such adapting cells can then be described by
We note that the explicit dependence of the firing rate $G}_{i}^{r$ of grid cell i on time $t$ needs to be considered separately for every cell for repetition suppression, which makes numerical simulations more computationally expensive. In contrast, the functions $G}_{i$ (in Equation 3) and $h}_{i$ (in Equation 4) for the conjunctive hypothesis depend also on time but only implicitly via the location $[{x}_{t},{y}_{t}]$ and the direction $\theta}_{t$ of movement at time $t$—and therefore the explicit time dependence of individual cells can be disregarded, which makes numerical simulations computationally cheaper.
Structurefunction mapping hypothesis
The structurefunction mapping hypothesis relies on a preferred grid phase for neighboring cells. We use two possible choices for the set of grid phases $({x}_{\text{off},i},{y}_{\text{off},i})$: they are either drawn from a uniform distribution on the unit rhombus or clustered. For clustered spatial phases, we draw $x}_{\text{off},i$ and $y}_{\text{off},i$ independently, each from a univariate von Mises distribution (with a defined central phase $\mu}_{s$ and concentration parameter $\kappa}_{s$). For a grid phase drawn from a uniform distribution, we note that random fluctuations lead to a certain degree of clustering of the gridphase sample (Figure 4—figure supplement 1). We can describe the resulting summed activity of $N$ grid cells simply by
with the upper index ‘s’ representing ‘structurefunction’.
Quantification of hexasymmetry of neural activity and trajectories
Combining the mathematical descriptions of grid cell activity for the three hypotheses (‘conjunctive’, ‘repetition suppression’, and ‘structurefunction’), we can denote the resulting population activity $A$ of $N$ grid cells by
where ${A}_{i}(t):={G}_{i}^{r}({x}_{t},{y}_{t},t)\phantom{\rule{thinmathspace}{0ex}}{h}_{i}({\theta}_{t})$ is the firing rate of cell i. To derive from $A$ the activity as a function of movement (or heading) direction $\theta}_{t$, we focus on time steps of length $\mathrm{\Delta}t$ in which the trajectory is linear. In time step $m$, i.e., for time $t$ in the time interval $[{t}_{m},{t}_{m+1})$ where ${t}_{m}=m\mathrm{\Delta}t$ and $m$ is an integer, the trajectory has the fixed angle $\theta}_{m$. The timediscrete mean activity $\overline{A}({t}_{m})$ associated to this interval is the average of $A(t)$ along the linear segment of the trajectory:
The integral in Equation 9 is either calculated analytically, as derived in the following section, or numerically. For a total number of $M$ time steps in a trajectory, the (normalized) mean activity $\stackrel{~}{A}(\varphi )$ as a function of some head direction $\varphi$ is then
where $\delta$ is the Dirac delta distribution. With complex Fourier coefficients $c}_{n$ (with $n\in \mathbb{N}$) defined as
we can quantify the hexasymmetry $H$ of the activity of a population of grid cells as
where $\stackrel{~}{{A}_{i}}}_{6$ is the sixth Fourier coefficient of cell i. The hexasymmetry is nonnegative ($H\ge 0$) and it has the unit spk/s (spikes per second). Furthermore, the average (over time) population activity can be expressed as
The hexasymmetry $H$ could be generated by various properties of the cells, but $H$ as defined above contains also contributions from the hexasymmetry of the underlying trajectory. This is due to the fact that we sum up in Equation 10 the population activities $\overline{A}({t}_{m})$ without taking into account the distribution of movement directions $\theta}_{m$. In this way, a hexasymmetry that is potentially contained in the subject’s navigation trajectory contributes to the hexasymmetry $H$ of the neural activity.
Empirical (fMRI/iEEG) studies (e.g. Doeller et al., 2010; Kunz et al., 2015; Maidenbaum et al., 2018) addressed this problem of trajectories spuriously contributing to hexasymmetry by fitting a generalized linear model (GLM) to the timediscrete fMRI/EEG activity. In contrast, our new approach to hexasymmetry in Equation 12 quantifies the contribution of the path to the neural hexasymmetry explicitly, and has the advantage that it allows an analytical treatment (see next section). Comparing our new method with previous methods for evaluating hexasymmetry led to qualitatively identical statistical effects (Figure 5—figure supplement 4).
To nevertheless be able to quantify how much a specific trajectory contributed to the neural hexasymmetry $H=\left{\stackrel{~}{A}}_{6}\right$, we also explicitly quantified the hexasymmetry of navigation trajectories and interpreted $H$ relative to the path hexasymmetry. To quantify the hexasymmetry of a trajectory, we used the same approach as in Equation 10 and defined the distribution of movement directions of the trajectory by
The path hexasymmetry of the trajectory is then ${\stackrel{~}{T}}_{6}:=\left\frac{1}{M}\sum _{m=0}^{M1}\mathrm{exp}(6j{\theta}_{m})\right$, which is similar to Equation 12. The path hexasymmetry is a unitless quantity, and ranges from 0 (for a uniform sampling of movement directions) to 1 (when only one movement direction is sampled).
To be able to estimate how much of the hexasymmetry $H$ of neuronal activity is due to the hexasymmetry of the trajectory, we compare the relative hexasymmetry of the activity, $\left{\stackrel{~}{A}}_{6}/{\stackrel{~}{A}}_{0}\right$, with the relative hexasymmetry of the trajectory, $\left{\stackrel{~}{T}}_{6}/{\stackrel{~}{T}}_{0}\right$, noting that ${\stackrel{~}{T}}_{0}\equiv 1$. The two terms being similar in magnitude, i.e., $\left{\stackrel{~}{A}}_{6}/{\stackrel{~}{A}}_{0}\right\approx \left{\stackrel{~}{T}}_{6}\right$, indicates that the trajectory is a major source of hexasymmetry whereas $H=\left{\stackrel{~}{A}}_{6}\right\gg {\stackrel{~}{A}}_{0}\left{\stackrel{~}{T}}_{6}\right$ suggests that hexasymmetry has a neural origin.
Analytical derivation of mean activity
In the following, we provide derivations that allow us to analytically integrate Equation 9 for the conjunctive grid by HD cell hypothesis and the structurefunction mapping hypothesis but not for the repetition suppression hypothesis. The respective results are shown in Figure 2B–H and Figure 4A–C, demonstrating that they are very similar to the numerical results.
We analytically calculate the mean activity $\overline{A}$ by averaging $A$ along a linear segment of a trajectory (Equation 9). For convenience, the following abbreviations are used in Equation 2 with the same grid spacing, ${s}_{i}=s$, and the same grid orientation, ${\gamma}_{i}=\gamma$, for all cells i:
A single grid cell
We start with a single grid cell without HD tuning and without repetition suppression. Equation 2 can be described in polar coordinates
In order to integrate along a piece of a straight line through the origin (similarly to the starlike walk), the angle $\psi}_{t}\equiv \overline{\psi$ can now be kept fixed (for that particular straight line) and we only have to consider ${G}_{i}({r}_{t},\overline{\psi})$. If we define $r}_{m$ and $r}_{m+1$ as the distances from zero that the subject is located at times $t}_{m$ and $t}_{m+1$ respectively, integration by substitution gives us
Since the speed of movement ${r}_{t}^{\mathrm{\prime}}\equiv v$ is assumed to be constant along the whole trajectory and $\mathrm{\Delta}r=v\mathrm{\Delta}t$, for $\mathrm{\Delta}r:={r}_{m+1}{r}_{m}$, we obtain
Thus, we have
where ${z}_{k}(r):=\mathrm{cos}[a\phantom{\rule{thinmathspace}{0ex}}{c}_{k,x}(r\mathrm{cos}(\overline{\psi}){b}_{x,i})+a\phantom{\rule{thinmathspace}{0ex}}{c}_{k,y}(r\mathrm{sin}(\overline{\psi}){b}_{y,i})]$. The four parts (A)–(D) are integrated separately. We obtain
where $\begin{array}{l}[[Q]]:=\{\begin{array}{l}1\phantom{\rule{1em}{0ex}}\text{if}Q\text{is true}\\ 0\phantom{\rule{1em}{0ex}}\text{else}\end{array}\end{array}$ is the Iverson bracket, and with the abbreviations
Note that $e}_{k$ actually depends on the cell index i which is omitted in Equation 25 and Equation 28 in order to keep the notation simpler.
HD tuning
For a conjunctive grid by HD cell, the HD tuning depends only on the angle (which is fixed when integrating along a straight line through zero) and not at all on the distance from zero. The mean activity is thus obtained from the mean activity of a grid cell without HD tuning by multiplying it with $h(\theta )$ defined in Equation 4:
Many cells
For more than one cell, the total mean activity (population rate) can simply be calculated as the sum of the mean activities of the single cells
Trajectories
The derived analytical description of the mean activity can be applied only to pieces of linear trajectories through zero. For the starlike walk, we can simply integrate
Piecewise linear trajectories and randomwalk trajectories consist of segments of straight lines that do not necessarily pass through zero. In order to integrate along the $m$th segment of a trajectory from $({x}_{{t}_{m}},{y}_{{t}_{m}})$ to $({x}_{{t}_{m+1}},{y}_{{t}_{m+1}})$ with movement direction $\theta}_{m$, we shift this path segment to the origin by subtracting $({x}_{{t}_{m}},{y}_{{t}_{m}})$ from the grid offset $({x}_{\text{off}},{y}_{\text{off}})$ and then integrating
where $r}_{m+1}=\sqrt{({x}_{{t}_{m+1}}{x}_{{t}_{m}}{)}^{2}+({y}_{{t}_{m+1}}{y}_{{t}_{m}}{)}^{2}$.
Upper bound for hexasymmetry of path trajectories
In the following we derive approximations for the expected value of the hexasymmetry $\left{\stackrel{~}{T}}_{6}\right$ of a path trajectory described by the number of time steps $M$, the movement tortuosity $\sigma}_{\theta$, and the time step size $\mathrm{\Delta}t$. These approximations can be used to assess the contribution of the underlying trajectory to the overall hexasymmetry of neural activity.
From the definition of the Fourier coefficients in Equation 11 and the movement direction distribution in Equation 14, we get the sixth Fourier coefficient of a trajectory:
The hexasymmetry of the trajectory can thus be expressed as
We simplify the sum of squares in Equation 36
with the help of the multinomial theorem and the trigonometric identity $\mathrm{cos}(\alpha \beta )=\mathrm{cos}(\alpha )\mathrm{cos}(\beta )+\mathrm{sin}(\alpha )\mathrm{sin}(\beta )$. If $\mathbb{E}(\mathrm{cos}(6({\theta}_{{m}_{1}}{\theta}_{{m}_{2}})))$ is known, we can compute the expected value of the square of the hexasymmetry of the path trajectory as
As $\mathbb{E}({X}^{2})=(\mathbb{E}(X){)}^{2}+\mathrm{V}\mathrm{a}\mathrm{r}(X)$ for any random variable $X$, we can use the result in Equation 38 to obtain the upper bound
In the following, we focus on the derivation of $\mathbb{E}(\mathrm{cos}(6({\theta}_{{m}_{1}}{\theta}_{{m}_{2}})))$. Using the movement statistics
that were introduced below Equation 1, the distribution of $\theta}_{m$ (after $m$ time steps) can be derived when we start at some angle $\theta}_{0$:
where $\mathcal{W}\mathcal{N}(\mu ,{\sigma}^{2})$ denotes the wrapped normal distribution with parameters μ and $\sigma}^{2$, which correspond to the mean and variance of the corresponding unwrapped distribution (Jammalamadaka and SenGupta, 2001).
In the following, we will derive the probability distribution of $\theta}_{{m}_{1}}{\theta}_{{m}_{2}$. We first define ${X}_{{m}_{1}}\sim \mathcal{N}({\theta}_{0},{m}_{1}{\sigma}_{\theta}^{2}\mathrm{\Delta}t)$ and ${X}_{{m}_{2}}\sim \mathcal{N}({\theta}_{0},{m}_{2}{\sigma}_{\theta}^{2}\mathrm{\Delta}t)$ as the unwrapped versions of $\theta}_{{m}_{1}$ and $\theta}_{{m}_{2}$, respectively, and $G}_{{X}_{{m}_{1}}{X}_{{m}_{2}}$ as the distribution function of their difference $X}_{{m}_{1}}{X}_{{m}_{2}$. Then, the distribution function (Fisher, 1995) of $\theta}_{{m}_{1}}{\theta}_{{m}_{2}$ reads as
where ${f}_{{X}_{{m}_{1}},{X}_{{m}_{2}}}(x,y)$ is the joint probability distribution.
In order to calculate the distribution of $\theta}_{{m}_{1}}{\theta}_{{m}_{2}$, we have to take the dependence between the two angles $\theta}_{{m}_{1}$ and $\theta}_{{m}_{2}$ into account. We will first consider the case $m}_{1}>{m}_{2$. In this case, the conditional distribution of $\theta}_{{m}_{1}$ given ${\theta}_{{m}_{2}}=y$ is wrapped normal with conditional mean $\mathbb{E}({\theta}_{{m}_{1}}{\theta}_{{m}_{2}}=y)=y$ and conditional variance $\mathrm{V}\mathrm{a}\mathrm{r}({\theta}_{{m}_{1}}{\theta}_{{m}_{2}}=y)=({m}_{1}{m}_{2}){\sigma}_{\theta}^{2}\mathrm{\Delta}t$. The same can be said about the unwrapped versions of $\theta}_{{m}_{1}$ and $\theta}_{{m}_{2}$. We will use this result to calculate the joint probability distribution
By applying the Leibniz integral rule in $(\ast )$ we obtain the probability density of $\theta}_{{m}_{1}}{\theta}_{{m}_{2}$:
For $m}_{2}>{m}_{1$, we derive in the same way as above
Finally, for $m}_{2}={m}_{1$, we have $\theta}_{{m}_{2}}={\theta}_{{m}_{1}$. Altogether, we thus get
In order to calculate an upper bound for the average path hexasymmetry, we will now use Equation 55 in Equation 38. Since $\mathbb{E}(\mathrm{\Re}(Z))=\mathrm{\Re}(\mathbb{E}(Z))$ for a complex random variable $Z$ and the $n$ th moment of a wrapped normal distribution with parameters μ and $\sigma}^{2$ is $\mathbb{E}(\mathrm{exp}(j\phantom{\rule{thinmathspace}{0ex}}X{)}^{n})=\mathrm{exp}(j\phantom{\rule{thinmathspace}{0ex}}n\phantom{\rule{thinmathspace}{0ex}}\mu \frac{1}{2}{n}^{2}{\sigma}^{2})$, we can derive
where ${\mu}_{{m}_{1},{m}_{2}}:=0$ and ${\sigma}_{{m}_{1},{m}_{2}}^{2}:={m}_{1}{m}_{2}{\sigma}_{\theta}^{2}\mathrm{\Delta}t$. We thus obtain
The solid lines in Figure 1—figure supplement 4 show the square root of Equation 57 (Equation 39).
From Equation 57 we can derive simplified approximations for two limiting cases. For convenience, we use $\alpha :=\frac{1}{2}\phantom{\rule{thinmathspace}{0ex}}36\phantom{\rule{thinmathspace}{0ex}}{\sigma}_{\theta}^{2}\mathrm{\Delta}t$. First, if $\alpha \gg 1$, i.e., if the new direction after one step is almost uniformly distributed or independent of the previous direction, we can neglect the double sum and we have
The corresponding line is shown in red in Figure 1—figure supplement 4. Note that in Figure 1—figure supplement 4 (with simulation step size $\mathrm{\Delta}t=0.01$ s), $\alpha \gg 1$ corresponds to ${\sigma}_{\theta}\gg 2.36$. A comparison with results from numerical simulations shows that for any ${\sigma}_{\theta}>3.5$ Equation 58 constitutes a viable upper bound of the mean path hexasymmetry.
Second, we assume $M\alpha \gg 1$, i.e., the direction after $M$ steps is almost independent from the original direction. The double sum in Equation 57 can then be approximated:
where the first series in Equation 61 is a geometric series and the second series is the polylogarithm function of order –1. We further approximate
For $\alpha \ll 1$, the error in the last approximation is $\frac{\mathrm{exp}(\alpha )}{{\left(\mathrm{exp}(\alpha )1\right)}^{2}}\approx \frac{1+\alpha}{{\alpha}^{2}}$, which, if $M\alpha \gg 1$, is negligible compared to the remaining term $\frac{M}{\mathrm{exp}(\alpha )1}$. Since $\frac{\mathrm{exp}(\alpha )}{\mathrm{exp}(\alpha )1}$ is a strictly monotonically decreasing function of $\alpha$, this approximation does not only hold for $\alpha \ll 1$ but is good in general.
Inserting Equation 63 into Equation 57 gives
Hence, we get the approximation (for $M\alpha \gg 1$)
This expression is used to compute the dashed lines in Figure 1—figure supplement 4, which all have slope $1/\sqrt{M}$ but a prefactor that depends on $\alpha =18{\sigma}_{\theta}^{2}\mathrm{\Delta}t$.
For $\alpha \ll 1$ (but still $M\alpha \gg 1$), we can use in Equation 64 the firstorder Taylor expansion of the exponential function at 0 to obtain
Hence, the path hexasymmetry for $\alpha \ll 1$ and $M\alpha \gg 1$ can be approximated by
which allows us to see how the key variables $M$, $\sigma}_{\theta$, and $\mathrm{\Delta}t$ interact in this limiting case. For instance, given a certain trajectory $A$ with $M}_{A$ steps and randomwalk parameters $\sigma}_{\theta A$ and $\mathrm{\Delta}{t}_{A}$, we can use Equation 67 to derive how many steps $M}_{B$ are necessary in a second trajectory $B$ with parameters $\sigma}_{\theta B$ and $\mathrm{\Delta}{t}_{B}$ to achieve the same mean path hexasymmetry. From Equation 67, we know that the two path hexasymmetries will have the same upper bound if
Hence, the number of time steps $M}_{A$ has to be multiplied by a factor $\mathrm{\Delta}M:=\frac{{\sigma}_{\theta A}^{2}\mathrm{\Delta}{t}_{A}}{{\sigma}_{\theta B}^{2}\mathrm{\Delta}{t}_{B}}$:
We illustrate the above considerations with an example: Let $\sigma}_{\theta A}=1{\text{rad/s}}^{1/2$, $\sigma}_{\theta B}=0.5{\text{rad/s}}^{1/2$, and $\mathrm{\Delta}{t}_{A}=\mathrm{\Delta}{t}_{B}=0.01$ s. From the given values, we obtain
Thus, the mean hexasymmetry value of a trajectory with ${\sigma}_{\theta A}=1$ after $M}_{A$ time steps is the same as the mean hexasymmetry value of a trajectory with ${\sigma}_{\theta B}=0.5$ after $M}_{B}=4\cdot {M}_{A$ time steps. Results from numerical simulations of path hexasymmetries, shown in Figure 1—figure supplement 4, support the derived theoretical approximations.
Randomfield simulations
To quantitatively evaluate the structurefunction mapping hypothesis, we set out to simulate a set of grid cells in threedimensional anatomical space. The grid phases associated with these grid cells follow the correlation structure suggested by Gu et al., 2018, and Heys et al., 2014. Our aim is to quantify the clustering of grid phases for a realistically sized fMRI voxel given this correlation structure.
We use a threedimensional representation of a voxel with a volume of 3 × 3 × 3 mm^{3}. Within this voxel, we define a grating of 200^{3} potential grid cells that are equally spaced in the voxel, with a distance between neighboring cells of 15 μm along the axes of the grating. To generate a set of random but spatially correlated grid phases on this area, we first define two random unit vectors in the complex plane, $\mathit{Z}}_{1$ and $\mathit{Z}}_{2$, for each of the 200^{3} potential grid cells in the voxel; angles of the unit vectors are thus drawn from a uniform distribution on the interval $[0,2\pi )$. $\mathit{Z}}_{1$ and $\mathit{Z}}_{2$ are further resolved into their real and imaginary components $\text{Re}({\mathit{Z}}_{i})$ and $\text{Im}({\mathit{Z}}_{i})$, respectively, where $i\in \{1,2\}$. To generate correlations between grid phases, we then convolve the two resulting gratings of 200^{3} components separately with either a Gaussian kernel (Figure 4K) or a grid kernel (Figure 4O) to yield the convolved components $\text{Re}({\hat{\mathit{Z}}}_{i})$ and $\text{Im}({\hat{\mathit{Z}}}_{i})$. The grid phases can be obtained by first calculating the angles of the new set of complex numbers and normalizing the result by $2\pi$:
We note that $\hat{x}}_{\text{off}$ and $\hat{y}}_{\text{off}$ are defined on the interval $[0,1)$ and correspond to the grid phases of a single grid cell mapped to the unit square. Transforming the result to the unit rhombus yields the grid phases $x}_{\text{off}$ and $y}_{\text{off}$ in the $x$ and $y$ direction respectively:
To find the average pairwise grid phase distances as a function of the pairwise anatomical distances, 10^{8} pairs of grid cells are sampled randomly from the uniform distribution defined on the discrete space of grating cell positions. The Euclidean distance in anatomical space between the two grid cells in each pair is calculated and sorted into 50 bins of equal width on the interval [10, 500] μm. Then, for each pair of grid cells, $n}_{1$ and $n}_{2$, eight copies of the grid phase $({x}_{\text{off},2},{y}_{\text{off},2})$ of the second cell $n}_{2$ are made, which are offset from the initial position of the grid phase such that they are positioned at the same phase within unit rhombi laid endtoend on a 3 × 3 grid. The minimum distance between the grid phase of the cell $n}_{1$ and the grid phase of each of the copies of the cell $n}_{2$ is taken as the pairwise phase distance. Finally, the pairwise distance between grid phase offsets per distance bin is obtained by taking the mean over all grid cell pairs whose Euclidean distance in anatomical space falls into the corresponding bin (Figure 4H, L, and P).
To estimate the clustering concentration parameter $\kappa}_{s$ in Figure 4, the phases $\hat{x}}_{\text{off}$ and $\hat{y}}_{\text{off}$ are mapped to a circular distribution by multiplying them with $2\pi$. The sets of grid phases $\{2\pi {\hat{x}}_{\text{off}}^{i}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}i\in \{1,2,...,N\}\}$ and $\{2\pi {\hat{y}}_{\text{off}}^{i}\phantom{\rule{thinmathspace}{0ex}}\phantom{\rule{thinmathspace}{0ex}}i\in \{1,2,...,N\}\}$ are then each separately fit to a onedimensional von Mises distribution to obtain a clustering concentration parameter for each axis. The final value of $\kappa}_{s$ is taken as the average of these two values.
Parameter estimation
The hexasymmetry values of the tested hypotheses depend on the respective parameters. We investigate two scenarios for each hypothesis: an ‘ideal’ set of parameters that results in a large (basically maximal) hexasymmetry value and a ‘realistic’ set of parameters that we try to derive from experimental data. All sets of parameters are summarized in Table 1, i.e., the symbols are introduced and defined, and parameter values and units are stated. We comment here on the ‘realistic’ sets of parameters and how we derived them.
Conjunctive grid by HD hypothesis
For the ‘realistic’ value of $\sigma}_{c}={3}^{\circ$, we use the 95% confidence interval of 12.4° from the third row of the table in the Supplementary Figure 5b of Doeller et al., 2010, and translate this value to a Gaussian standard deviation, assuming a 95% confidence interval of $4{\sigma}_{c}$.
For the ‘realistic’ value of ${\kappa}_{c}=4$, we refer to the Supplementary Figure 3 of Doeller et al., 2010. Visually, the tuning curves usually cover between one sixth and one third of the circle. Stronger HD tuning contributes the most to the resulting hexasymmetry, so we choose a value of ${\kappa}_{c}=4$, which corresponds to a tuning width of roughly one sixth of the circle. For further discussion on $\kappa}_{c$, see also Figure 2—figure supplement 1.
The ‘realistic’ value of ${\mathrm{p}}_{\mathrm{c}}=1/3$ was chosen based on the values reported by Sargolini et al., 2006, and Boccara et al., 2010. In our simulations, we observed that the hexasymmetry shows a linear dependence on the fraction of conjunctive cells up to a noise floor (Figure 2—figure supplement 2).
Repetition suppression hypothesis
To obtain ‘realistic’ parameter values, we simply divide the ‘ideal’ ones by two, which results in a ‘realistic’ adaptation time constant ${\tau}_{r}=1.5$ s and a ‘realistic’ adaptation weight ${w}_{r}=0.5$. ‘Ideal’ parameter values attenuate the firing rate by ≈20% (Figure 3A, with 24% when running aligned with the grid axes and 18% when running misaligned to the grid axes). ‘Realistic’ parameter values attenuate firing rates by ≈16% (for a dependence of hexasymmetry on the two parameters, see Figure 3B). Such a decrease has a similar order of magnitude as the attenuation of 5% of firing rates in Reber et al., 2023 (their Figure 2B, right). The value of the ‘realistic’ time constant is also comparable to the time scales found in two electrophysiological studies: first, Giocomo and Hasselmo, 2008, showed that the time constant in the slow component of the hyperpolarizationactivated current ($I}_{h$) in the ventral mEC is 2–3 s (their Figure 2F). Second, Magistretti and Alonso, 1999, found a secondlong inactivation of a persistent sodium current in EC layer II cells.
Structurefunction mapping hypothesis
The ‘realistic’ value of the concentration parameter for clustering, ${\kappa}_{s}=0.1$, is derived from the randomfield simulations in Figure 4G–S and a comparison to results in Gu et al., 2018, and Heys et al., 2014. Figure 4G in Gu et al., 2018, shows for grid cells in the rat mEC the relationship between their pairwise physical distance and their pairwise phase distance. Similarly, for grid cells in the mouse mEC, Heys et al., 2014, in their Figure 4B show the correlation of (onedimensional) grid firing maps as a function of anatomical pairwise distance between grid cells. Together, these data serve as evidence for more similar grid cell firing for anatomically close grid cells (distance <100 μm)—with decreasing similarity for increasing anatomical distance. In our randomfield simulations (Figure 4), we mimic these experimental results by choosing the width of a spatial correlation kernel accordingly (Figure 4K, Gaussian kernel with standard deviation: 30 μm). Consequently, in the simulations, the initial rise of phase distance as a function of anatomical distance (Figure 4L) reflects the experimental curves. In Gu et al.’s Figure 4G there is a slight increase in gridfiring similarity (decrease of phase distance) for anatomical distances of approximately 300 μm. In our simulations, we reflect such a nonmonotonous behavior by introducing a gridlike correlation kernel (Figure 4O) with a grid spacing of 300 μm. Applying this gridlike correlation kernel to the random field in Figure 4G results in a phasedistance curve with a slight dip around 300 μm (Figure 4P), just like in the experimental curves in Gu et al., 2018. For both cases, the Gaussian kernel and the gridlike kernel, we find that the resulting phase clustering for an anatomical voxel of realistic size is quite low, i.e., $\kappa}_{s$ is well below 0.1 (Figure 4M and Q). As it poses an upper bound for the realistic values of $\kappa}_{s$, we use 0.1 as an estimate of the ‘realistic’ value for $\kappa}_{s$ in our further simulations. For the central phase of the cluster, we always use (without loss of generality) the origin: ${\mu}_{s}=(0,0)$.
Implementation of previously used metrics
We applied three previously used metrics to our framework: the GLM method by Doeller et al., 2010; the GLM method with binning by Kunz et al., 2015, the circularlinear correlation method by Maidenbaum et al., 2018; Figure 5—figure supplement 4.
In brief, in the GLM method (e.g. used in Doeller et al., 2010), the hexasymmetry is found in two steps: the orientation of the hexadirectional modulation is first estimated on the first half of the data by using the regressors ${\beta}_{1}\mathrm{cos}(6{\theta}_{t})$ and ${\beta}_{2}\mathrm{sin}(6{\theta}_{t})$ on the timediscrete fMRI activity (Equation 9), with $\theta}_{t$ being the movement direction of the subject in time step $t$. The amplitude of the signal is then estimated on the second half of the data using the single regressor $\beta \mathrm{cos}[6({\theta}_{t}\varphi )]$, where $\varphi =\mathrm{arctan}({\beta}_{2}/{\beta}_{1})/6$. The hexasymmetry is then evaluated as $H=\beta /2$.
The GLM method with binning (e.g. used in Kunz et al., 2015) uses the same procedure as the GLM method for estimating the grid orientation in the first half of the data, but the amplitude is estimated differently on the second half by a regressor that has a value 1 if $\theta}_{t$ is aligned with a peak of the hexadirectional modulation (aligned if ${\theta}_{t}\varphi \phantom{\rule{thinmathspace}{0ex}}\mathrm{\%}\phantom{\rule{thinmathspace}{0ex}}{60}^{\circ}<{15}^{\circ}$, % modulo operator) and a value of −1 if $\theta}_{t$ is misaligned. The hexasymmetry is then calculated from the amplitude in the same way as the GLM method.
The circularlinear correlation method (e.g. used in Maidenbaum et al., 2018) is similar to the GLM method in that it uses the regressors ${\beta}_{1}\mathrm{cos}(6{\theta}_{t})$ and ${\beta}_{2}\mathrm{sin}(6{\theta}_{t})$ on the timediscrete mean activity, but instead of using $\beta}_{1$ and $\beta}_{2$ to estimate the orientation of the hexadirectional modulation, the beta values are directly used to estimate the hexasymmetry using the relation $H=\sqrt{{\beta}_{1}^{2}+{\beta}_{2}^{2}}/2$.
Regarding the statistical evaluation, each method evaluates the size of the neural hexasymmetry differently. Specifically, the new method developed in our manuscript compares the neural hexasymmetry to path hexasymmetry to test whether neural hexasymmetry is significantly above path hexasymmetry. For the two GLM methods, we compare the hexasymmetry to zero (using the MannWhitney U test) to establish significance. Hexasymmetry values can be negative in these approaches, allowing the statistical comparison against 0. Negative values occur when the estimated grid orientation from the first data half does not match the grid orientation from the second data half. Regarding the statistical evaluation of the circularlinear correlation method, we calculated a zscore by comparing each empirical observation of the hexasymmetry to hexasymmetries from a set of surrogate distributions (as in Maidenbaum et al., 2018). We then calculate a pvalue by comparing the distribution of zscores versus zero using a MannWhitney U test. We use the zscores instead of the hexasymmetry for the circularlinear correlation method to match the procedure used in Maidenbaum et al., 2018. We obtained the surrogate distributions by circularly shifting the vector of movement directions relative to the timedependent vector of firing rates. For random walks, the vector is shifted by a random number drawn from a uniform distribution defined with the same length as the number of time points in the vector of movement directions. For the starlike walks and piecewise linear walks, the shift is a random integer multiplied by the number of time points in a linear segment. Circularly shifting the vector of movement directions scrambles the correlations between movement direction and neural activity while preserving their temporal structure.
Tuning widths and the large$\kappa$ approximation
The numerical values of our tuning width of HD modulated cells should not be directly compared to the numbers reported by Sargolini et al., 2006, with values in the range of 55° because Sargolini et al., 2006 reported the tuning widths as the angular standard deviation $\sigma$ of the mean vector of HDdependent firing rates whereas we derive the tuning width from the large$\kappa}_{c$ approximation of a von Mises distribution: $\sigma =1/\sqrt{{\kappa}_{c}}$. This simple approximation is reasonably good only for large enough $\kappa}_{c$ (or small enough $\sigma$), and 55° is already beyond the range where the approximation is reasonable (Figure 2—figure supplement 1B).
To examine the range of values where the large$\kappa}_{c$ approximation holds, we can compare it to the angular standard deviation expressed as $\sigma =\sqrt{2\mathrm{ln}[{I}_{1}({\kappa}_{c})/{I}_{0}({\kappa}_{c})]}$, where ${I}_{n}({\kappa}_{c})$ is the modified Bessel function of the first kind of order ‘n’. This equation is valid for all values of the concentration parameter $\kappa}_{c$. The results of this comparison are summarized in Figure 2—figure supplement 1B. In brief, the large$\kappa}_{c$ approximation closely matches the angular standard deviation for ${\kappa}_{c}>4$, which corresponds to $\sigma <{29}^{\circ}$. The angular standard deviation $\sigma ={55}^{\circ}$ reported by Sargolini et al., 2006 corresponds to a concentration parameter $\kappa}_{c$ of about 1.7 and (in our large$\kappa}_{c$ approximation) to a tuning width $1/\sqrt{{\kappa}_{c}}={44}^{\circ}$ (upward triangle in Figure 2—figure supplement 1B and C).
Numerical simulations
All simulations were implemented in Python 3.10 using the packages NumPy, SciPy, Numba, and JobLib. Matplotlib was used for plotting. Inkscape was used for the final adjustment to the figures.
Data availability
The code used to generate the data is available at https://github.com/ikhwankhalid/grid_bold, copy archived at Khalid, 2024.
References

ConferenceA HighResolution Model of the Human Entorhinal Cortex in the ‘BigBrain’ – Use Case for Machine Learning and 3D AnalyseBraininspired computing: 4th international workshop, BrainComp 2019, Cetraro, taly, July 15–19, 2019, Revised Selected Papers 4. pp. 3–21.https://doi.org/10.1007/9783030824273_1

Grid cells in pre and parasubiculumNature Neuroscience 13:987–994.https://doi.org/10.1038/nn.2602

An oscillatory interference model of grid cell firingHippocampus 17:801–812.https://doi.org/10.1002/hipo.20327

A singlecell spiking model for the origin of gridcell patternsPLOS Computational Biology 13:e1005782.https://doi.org/10.1371/journal.pcbi.1005782

BookStatistical Analysis of Circular DataCambridge University Press.https://doi.org/10.1017/CBO9780511564345

Grid cells are modulated by local head directionNature Communications 11:4228.https://doi.org/10.1038/s41467020175001

Impaired path integration in mice with disrupted grid cell firingNature Neuroscience 21:81–91.https://doi.org/10.1038/s4159301700393

Repetition and the brain: neural models of stimulusspecific effectsTrends in Cognitive Sciences 10:14–23.https://doi.org/10.1016/j.tics.2005.11.006

Gridlike processing of imagined navigationCurrent Biology 26:842–847.https://doi.org/10.1016/j.cub.2016.01.042

Direct recordings of gridlike neuronal activity in human spatial navigationNature Neuroscience 16:1188–1190.https://doi.org/10.1038/nn.3466

Human entorhinal cortex represents visual space using a boundaryanchored gridNature Neuroscience 21:191–194.https://doi.org/10.1038/s4159301700491

Nonhexagonal neural dynamics in vowel spaceAIMS Neuroscience 7:275–298.https://doi.org/10.3934/Neuroscience.2020015

Mesoscopic neural representations in spatial navigationTrends in Cognitive Sciences 23:615–630.https://doi.org/10.1016/j.tics.2019.04.011

Spatial representation in the hippocampal formation: a historyNature Neuroscience 20:1448–1464.https://doi.org/10.1038/nn.4653

Hexadirectional coding of visual space in human entorhinal cortexNature Neuroscience 21:188–190.https://doi.org/10.1038/s4159301700508

Structural modularity and grid activity in the medial entorhinal cortexJournal of Neurophysiology 119:2129–2144.https://doi.org/10.1152/jn.00574.2017

Singleneuron mechanisms of neural adaptation in the human temporal lobeNature Communications 14:2496.https://doi.org/10.1038/s41467023381905

Grid alignment in entorhinal cortexBiological Cybernetics 106:483–506.https://doi.org/10.1007/s0042201205137

Microcircuits for spatial coding in the medial entorhinal cortexPhysiological Reviews 102:653–688.https://doi.org/10.1152/physrev.00042.2020
Article and author information
Author details
Funding
Deutsche Forschungsgemeinschaft (327654276  SFB 1315)
 Richard Kempter
Bundesministerium für Bildung und Forschung (01GQ1705)
 Richard Kempter
Einstein Center for Neurosciences Berlin
 Ikhwan Bin Khalid
Bundesministerium für Bildung und Forschung (01GQ1705A)
 Lukas Kunz
National Institutes of Health (U01 NS11319801)
 Lukas Kunz
German Research Foundation (447634521)
 Lukas Kunz
German Research Foundation (527084865)
 Lukas Kunz
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We would like to thank Tiziano D’Albis for discussions.
This work was funded by the German Research Foundation (DFG; project nos. 327654276  SFB 1315 to RK), the German Federal Ministry of Education and Research (01GQ1705 to RK), and the Einstein Foundation Berlin (to IK). LK received funding from the German Research Foundation (DFG; project nos. 447634521 and 527084865), the return program of the Ministry of Culture and Science of the State of North RhineWestphalia, the Federal Ministry of Education and Research (BMBF; 01GQ1705A), and by NIH/NINDS grant U01 NS11319801.
Copyright
© 2024, Bin Khalid, Reifenstein et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 227
 views

 15
 downloads

 0
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
Female sexual receptivity is essential for reproduction of a species. Neuropeptides play the main role in regulating female receptivity. However, whether neuropeptides regulate female sexual receptivity during the neurodevelopment is unknown. Here, we found the peptide hormone prothoracicotropic hormone (PTTH), which belongs to the insect PG (prothoracic gland) axis, negatively regulated virgin female receptivity through ecdysone during neurodevelopment in Drosophila melanogaster. We identified PTTH neurons as doublesexpositive neurons, they regulated virgin female receptivity before the metamorphosis during the thirdinstar larval stage. PTTH deletion resulted in the increased EcRA expression in the whole newly formed prepupae. Furthermore, the ecdysone receptor EcRA in pC1 neurons positively regulated virgin female receptivity during metamorphosis. The decreased EcRA in pC1 neurons induced abnormal morphological development of pC1 neurons without changing neural activity. Among all subtypes of pC1 neurons, the function of EcRA in pC1b neurons was necessary for virgin female copulation rate. These suggested that the changes of synaptic connections between pC1b and other neurons decreased female copulation rate. Moreover, female receptivity significantly decreased when the expression of PTTH receptor Torso was reduced in pC1 neurons. This suggested that PTTH not only regulates female receptivity through ecdysone but also through affecting female receptivity associated neurons directly. The PG axis has similar functional strategy as the hypothalamic–pituitary–gonadal axis in mammals to trigger the juvenile–adult transition. Our work suggests a general mechanism underlying which the neurodevelopment during maturation regulates female sexual receptivity.

 Neuroscience
Theoretical computational models are widely used to describe latent cognitive processes. However, these models do not equally explain data across participants, with some individuals showing a bigger predictive gap than others. In the current study, we examined the use of theoryindependent models, specifically recurrent neural networks (RNNs), to classify the source of a predictive gap in the observed data of a single individual. This approach aims to identify whether the low predictability of behavioral data is mainly due to noisy decisionmaking or misspecification of the theoretical model. First, we used computer simulation in the context of reinforcement learning to demonstrate that RNNs can be used to identify model misspecification in simulated agents with varying degrees of behavioral noise. Specifically, both prediction performance and the number of RNN training epochs (i.e., the point of early stopping) can be used to estimate the amount of stochasticity in the data. Second, we applied our approach to an empirical dataset where the actions of low IQ participants, compared with high IQ participants, showed lower predictability by a wellknown theoretical model (i.e., Daw’s hybrid model for the twostep task). Both the predictive gap and the point of early stopping of the RNN suggested that model misspecification is similar across individuals. This led us to a provisional conclusion that low IQ subjects are mostly noisier compared to their high IQ peers, rather than being more misspecified by the theoretical model. We discuss the implications and limitations of this approach, considering the growing literature in both theoretical and datadriven computational modeling in decisionmaking science.