To evaluate the different hypotheses on the emergence of grid-like representations, we considered three different types of navigation trajectories: ‘star-like 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 star-like 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 real-world experiments, this type of navigation including teleportation is unusual, but it can be implemented in virtual-reality experiments (Vass et al., 2016; Deuker et al., 2016).

Figure 2 with 2 supplements see all

Download asset Open asset

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 star-like 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 star-like 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 (${\displaystyle {\sigma }_{\theta }=0.5}$ rad/s^1/2) and a time step ${\displaystyle \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 ${\displaystyle 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 (${\displaystyle M=9\cdot {10}^5}$ steps).

As we describe below, the emergence of grid-like 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.

To test how the activity of grid cells could give rise to hexasymmetry of a macroscopic signal, we used a firing-rate model of grid cell activity (Equation 2). Furthermore, we developed a new measure ${\displaystyle 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 ${\displaystyle H}$ has a value ${\displaystyle 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 ${\displaystyle H}$. Using the same approach, we quantified the hexasymmetry of a trajectory and called this path hexasymmetry (see Methods after Equation 14).

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 grid-cell activity. This hypothesis proposes that grid-cell activity is subject to firing-rate adaptation and thus leads to reduced grid-cell 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).

Figure 3

Download asset Open asset

In our model, the repetition suppression hypothesis depends on two adaptation parameters: the adaptation time constant ${\displaystyle {\tau }_r}$ and the adaptation weight ${\displaystyle 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 ${\displaystyle v}$ divided by the grid scale ${\displaystyle s}$ (Figure 3B). We constrained values of the adaptation weight to the full range of reasonable values (${\displaystyle 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 star-like 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 grid-cell activity because the tortuosity of the random walk that we typically used (${\displaystyle {\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 movement-direction 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 ${\displaystyle {\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 grid-cell 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 grid-cell activity) is shifted by ${\displaystyle {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 grid-cell activity more strongly when moving aligned with a grid axis as compared to when moving misaligned with a grid axis (Figure 3A).

Structure-function mapping hypothesis

We next investigated the structure-function 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 ‘structure-function 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 ${\displaystyle {\kappa }_s=10}$) resulted in significant hexadirectional modulations of sum grid-cell activity when the subject performed star-like walks starting at a phase offset of ${\displaystyle (0,0)}$, i.e., the center of the cluster of firing fields of grid cells (Figure 4A). Interestingly, the hexasymmetry values during star-like 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 ${\displaystyle (0,0)}$ or ${\displaystyle (0.3,0.3)}$ led to clear hexasymmetry whereas others, e.g., ${\displaystyle (0.6,0)}$, did not (Figure 4D, left). Additionally, the ‘apparent preferred grid orientation’ (i.e. the movement directions associated with the highest sum grid-cell activity) was shifted by ${\displaystyle {30}^{\circ }}$ for certain offsets in the unit rhombus illustrating the subject’s starting locations relative to the firing-field locations (Figure 4D, right). It was furthermore notable that the summed grid-cell activity as a function of movement direction exhibited relatively sharp peaks at multiples of ${\displaystyle {60}^{\circ }}$ with additional small peaks in between (Figure 4A, right). This pattern is clearly distinct from the more sinusoidal modulation of sum grid-cell activity resulting from the conjunctive grid by HD cell hypothesis (Figure 2) and the more full-wave rectified sinusoidal modulation resulting from the repetition suppression hypothesis (Figure 3).

Figure 4 with 1 supplement see all

Download asset Open asset

When examining the structure-function mapping hypothesis for piecewise linear walks and random walks, hexasymmetry appeared to be considerably lower as compared to simulations with star-like walks (Figure 4B and C). For piecewise linear walks (but not for random walks), population grid-cell activity as a function of movement direction again exhibited sharp peaks at multiples of ${\displaystyle {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 star-like walks (Figure 4D).

In the above-mentioned simulations for the structure-function 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 ${\displaystyle {\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, ${\displaystyle {\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 (${\displaystyle {\kappa }_s<0.1}$; Figure 4G–S) than to strong, ideal clustering ${\displaystyle ({\kappa }_s=10)}$. Specifically, we implemented a ‘simple’ clustering using a Gaussian kernel that led to a clustering concentration parameter ${\displaystyle {\kappa }_s\approx 0.014}$ (Figure 4M). The clustering concentration parameter remained low when we added an empirically guided ‘higher-order’ spatial clustering of grid phase offsets by adding some longer-range spatial autocorrelations to the grid-phase offsets (Figure 4N–Q; Gu et al., 2018), which led to a clustering concentration parameter ${\displaystyle {\kappa }_s\approx 0.052}$ (Figure 4Q). Accordingly, we found that these more ‘realistic’ clustering concentration parameters resulted in substantially reduced hexasymmetries for the structure-function 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.

Figure 5 with 4 supplements see all

Download asset Open asset

Our above-mentioned results for the structure-function mapping hypothesis indicate that the navigation pattern has a strong influence on the neuronal hexasymmetry ${\displaystyle 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 ${\displaystyle 1/\sqrt{M}}$ for a large number ${\displaystyle M}$ of steps (Figure 1—figure supplement 4). For the random walk in Figure 4C with tortuosity ${\displaystyle {\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 (${\displaystyle M=0.9\cdot {10}^6}$ steps for ${\displaystyle \mathrm{\Delta }t=0.01}$ s) is about 0.007, which results in a contribution to the neural hexasymmetry of ${\displaystyle H\approx 7}$ spk/s for a population rate of about 1024 spk/s. This number is close to the obtained neural hexasymmetry (${\displaystyle H=8.0}$ spk/s) in Figure 4C. We thus believe that for random walks the hexasymmetry ${\displaystyle H}$ obtained in the structure-function mapping case is mainly determined by the path hexasymmetry.

We then turned to piecewise linear walks and examined neural hexasymmetries ${\displaystyle H}$ across a broad range of total trajectory distances. Surprisingly, we observed that ${\displaystyle H}$ also decreased as ${\displaystyle 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 ${\displaystyle H}$ of sum grid-cell 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 grid-field locations.

Taken together, the structure-function mapping hypothesis with very strong (‘ideal’) clustering produces neural hexasymmetry values that are larger than expected from path hexasymmetries only with respect to star-like 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).

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 Mann-Whitney U tests between the neural hexasymmetries (${\displaystyle H:=|{\tilde{A}}_6|}$; see also Equation 12) and the product ${\displaystyle |{\tilde{T}}_6|\cdot {\tilde{A}}_0}$ with the multipliers path hexasymmetry ${\displaystyle |{\tilde{T}}_6|}$ and average population activity ${\displaystyle {\tilde{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 (Mann-Whitney U tests, all U=0, all p<0.001) and for more realistic parameters (Mann-Whitney U tests, all U=0, all p<0.001). For star-like 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 star-like walks and piecewise linear walks resulted in significant neural hexasymmetries for both ideal and weaker parameters (Mann-Whitney U tests, all U=0, all p<0.001), whereas random walks did not (Mann-Whitney 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 star-like 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 structure-function 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 (Mann-Whitney U tests, all U<736, all p<0.001) and more realistic parameters (Mann-Whitney U tests for star-like 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 (Mann-Whitney U test, U=1786, p=0.472). However, these results regarding the structure-function mapping hypothesis should be treated with great caution. First, the neural hexasymmetries for star-like 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 structure-function 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 subject-specific path distances with the subject-specific 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 structure-function mapping hypothesis leads to true neural hexasymmetry only in the case of star-like walks.

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 grid-cell 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 grid-cell activity.

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 grid-like representations. This includes the possibility that neural hexasymmetry is completely unrelated to grid-cell 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 grid-like 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°. Attractor-network structures may be involved in such a mechanism, but this remains speculative at the current stage.

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 random-walk 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 star-like 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 star-like trajectories and piecewise linear trajectories (but not for random trajectories); and that for the structure-function mapping hypothesis ‘true’ hexadirectional population signals can only be observed for star-like 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 star-like 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 star-like or piecewise linear trajectories (in favor of the repetition suppression hypothesis); or whether they are mainly visible for star-like 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 structure-function 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 structure-function mapping hypothesis predicts that in fMRI studies involving star-like 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).

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 10-fold 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 structure-function 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 longer-range 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 grid-cell properties (i.e. the adaptation time constant ${\displaystyle {\tau }_r}$ and the adaptation weight ${\displaystyle 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.

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 grid-cell 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 structure-function mapping hypothesis and star-like 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 structure-function mapping hypothesis, random walks and piecewise-linear 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³ 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 structure-function mapping hypothesis, significant neural hexasymmetry can be achieved only for star-like 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 structure-function 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 grid-field 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.

The neural hexasymmetry as defined by Equation 12 also contains contributions from the path hexasymmetry of the underlying trajectory. Star-like 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 ${\displaystyle M}$ and the tortuosity parameter ${\displaystyle {\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 time-averaged 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 time-averaged population activity. For M=900,000 time steps and a tortuosity of ${\displaystyle {\sigma }_{\theta }=0.5}$ rad/s^1/2, random-walk trajectories contribute ∼0.8% of the time-averaged population activity to the neural hexasymmetry (Figure 1—figure supplement 4). In our simulations with an average population activity in the range of 10³ 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 structure-function mapping hypothesis), we sample a finite number of grid phase offsets from a two-dimensional von Mises distribution with a clustering concentration parameter ${\displaystyle {\kappa }_s=0}$ (equivalent to a uniform distribution on the unit rhombus). In this scenario, random fluctuations give rise to realizations with a non-zero sample mean clustering concentration parameter ${\displaystyle {\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 structure-function mapping hypothesis, these fluctuations would not be considered as ‘noise’ as they directly contribute to the hexasymmetry from the clustering of grid phase offsets.

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 single-neuron 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 single-neuron activity and the fMRI signal is highly complex (Ekstrom, 2010; Kunz et al., 2019). Future studies are needed to detail the relationship between single-neuron firing and fMRI/iEEG signals in the EC. This would allow us to clarify whether a hexadirectional modulation of sum grid-cell activity directly results in a hexadirectional modulation of fMRI/iEEG activity or whether currently unknown factors modulate the expression of hexadirectional fMRI/iEEG signals.

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.

To describe grid-cell activity as a function of time ${\displaystyle 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 quasi-infinite, and then add rules that account for boundaries.

Following Burgess et al., 2007, the activity profile ${\displaystyle G_i}$ of grid cell i (for ${\displaystyle i=1,...,N}$ in a population of ${\displaystyle N}$ grid cells) is modeled as the product of three cosine waves rotated by 60° (${\displaystyle =\pi /3}$) from each other:

where ${\displaystyle A_{\text{max}}}$ is the grid cell’s maximal firing rate, ${\displaystyle s_i}$ depicts the cell’s grid spatial scale (‘grid spacing’), ${\displaystyle x_{\text{off},i}}$ and ${\displaystyle y_{\text{off},i}}$ are the phase offsets of the grid (‘grid phase’) in the two spatial dimensions (called ${\displaystyle x}$ and ${\displaystyle y}$ here), and ${\displaystyle {\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 ${\displaystyle N}$ grid cells in Equation 2. For a given trajectory ${\displaystyle X_t=[x_t,y_t]}$, the macroscopic activity as a function of time ${\displaystyle t}$ is then simply described by the sum ${\displaystyle \sum _{i=1}^N{G}_i(x_t,y_t)}$.

Here, we summarize how the activity in a population of ${\displaystyle N}$ grid cells can be described if they also exhibit (i) HD tuning, (ii) repetition suppression (i.e. firing-rate adaptation), or (iii) grid phases that are clustered across grid cells.

In all our models, the activity ${\displaystyle G_i}$ of a grid cell is described by Equation 2. Different grid cells typically have different phase offsets (${\displaystyle x_{\text{off},i}}$, ${\displaystyle y_{\text{off},i}}$) but the same grid spacing ${\displaystyle s_i:=s,\mathrm{\forall }i}$ and grid orientation ${\displaystyle {\gamma }_i:=\gamma ,\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.

Combining the mathematical descriptions of grid cell activity for the three hypotheses (‘conjunctive’, ‘repetition suppression’, and ‘structure-function’), we can denote the resulting population activity ${\displaystyle A}$ of ${\displaystyle N}$ grid cells by

where ${\displaystyle A_i(t):=G_i^r(x_t,y_t,t)h_i({\theta }_t)}$ is the firing rate of cell i. To derive from ${\displaystyle A}$ the activity as a function of movement (or heading) direction ${\displaystyle {\theta }_t}$, we focus on time steps of length ${\displaystyle \mathrm{\Delta }t}$ in which the trajectory is linear. In time step ${\displaystyle m}$, i.e., for time ${\displaystyle t}$ in the time interval ${\displaystyle [t_m,t_{m+1})}$ where ${\displaystyle t_m=m\mathrm{\Delta }t}$ and ${\displaystyle m}$ is an integer, the trajectory has the fixed angle ${\displaystyle {\theta }_m}$. The time-discrete mean activity ${\displaystyle \overline{A}(t_m)}$ associated to this interval is the average of ${\displaystyle 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 ${\displaystyle M}$ time steps in a trajectory, the (normalized) mean activity ${\displaystyle \tilde{A}(\varphi )}$ as a function of some head direction ${\displaystyle \varphi }$ is then

where ${\displaystyle \delta }$ is the Dirac delta distribution. With complex Fourier coefficients ${\displaystyle c_n}$ (with ${\displaystyle n\in \mathbb{N}}$) defined as

we can quantify the hexasymmetry ${\displaystyle H}$ of the activity of a population of grid cells as

where ${\displaystyle {\widetilde{A_i}}_6}$ is the sixth Fourier coefficient of cell i. The hexasymmetry is nonnegative (${\displaystyle 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 ${\displaystyle H}$ could be generated by various properties of the cells, but ${\displaystyle 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 ${\displaystyle \overline{A}(t_m)}$ without taking into account the distribution of movement directions ${\displaystyle {\theta }_m}$. In this way, a hexasymmetry that is potentially contained in the subject’s navigation trajectory contributes to the hexasymmetry ${\displaystyle 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 time-discrete 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 ${\displaystyle H=\left|{\tilde{A}}_6\right|}$, we also explicitly quantified the hexasymmetry of navigation trajectories and interpreted ${\displaystyle 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 ${\displaystyle |{\tilde{T}}_6|:=\left|\frac{1}{M}\sum _{m=0}^{M-1}\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 ${\displaystyle H}$ of neuronal activity is due to the hexasymmetry of the trajectory, we compare the relative hexasymmetry of the activity, ${\displaystyle \left|{\tilde{A}}_6/{\tilde{A}}_0\right|}$, with the relative hexasymmetry of the trajectory, ${\displaystyle \left|{\tilde{T}}_6/{\tilde{T}}_0\right|}$, noting that ${\displaystyle {\tilde{T}}_0\equiv 1}$. The two terms being similar in magnitude, i.e., ${\displaystyle \left|{\tilde{A}}_6/{\tilde{A}}_0\right|\approx \left|{\tilde{T}}_6\right|}$, indicates that the trajectory is a major source of hexasymmetry whereas ${\displaystyle H=\left|{\tilde{A}}_6\right|\gg {\tilde{A}}_0\left|{\tilde{T}}_6\right|}$ suggests that hexasymmetry has a neural origin.

In the following we derive approximations for the expected value of the hexasymmetry ${\displaystyle \left|{\tilde{T}}_6\right|}$ of a path trajectory described by the number of time steps ${\displaystyle M}$, the movement tortuosity ${\displaystyle {\sigma }_{\theta }}$, and the time step size ${\displaystyle \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 ${\displaystyle \cos (\alpha -\beta )=\cos (\alpha )\cos (\beta )+\sin (\alpha )\sin (\beta )}$. If ${\displaystyle \mathbb{E}(\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 ${\displaystyle \mathbb{E}(X^2)=(\mathbb{E}(X){)}^2+\mathrm{V}\mathrm{a}\mathrm{r}(X)}$ for any random variable ${\displaystyle X}$, we can use the result in Equation 38 to obtain the upper bound

In the following, we focus on the derivation of ${\displaystyle \mathbb{E}(\cos (6({\theta }_{m_1}-{\theta }_{m_2})))}$. Using the movement statistics

that were introduced below Equation 1, the distribution of ${\displaystyle {\theta }_m}$ (after ${\displaystyle m}$ time steps) can be derived when we start at some angle ${\displaystyle {\theta }_0}$:

where ${\displaystyle \mathcal{W}\mathcal{N}(\mu ,{\sigma }^2)}$ denotes the wrapped normal distribution with parameters μ and ${\displaystyle {\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 ${\displaystyle {\theta }_{m_1}-{\theta }_{m_2}}$. We first define ${\displaystyle X_{m_1}\sim \mathcal{N}({\theta }_0,m_1{\sigma }_{\theta }^2\mathrm{\Delta }t)}$ and ${\displaystyle X_{m_2}\sim \mathcal{N}({\theta }_0,m_2{\sigma }_{\theta }^2\mathrm{\Delta }t)}$ as the unwrapped versions of ${\displaystyle {\theta }_{m_1}}$ and ${\displaystyle {\theta }_{m_2}}$, respectively, and ${\displaystyle G_{X_{m_1}-X_{m_2}}}$ as the distribution function of their difference ${\displaystyle X_{m_1}-X_{m_2}}$. Then, the distribution function (Fisher, 1995) of ${\displaystyle {\theta }_{m_1}-{\theta }_{m_2}}$ reads as

where ${\displaystyle f_{X_{m_1},X_{m_2}}(x,y)}$ is the joint probability distribution.

In order to calculate the distribution of ${\displaystyle {\theta }_{m_1}-{\theta }_{m_2}}$, we have to take the dependence between the two angles ${\displaystyle {\theta }_{m_1}}$ and ${\displaystyle {\theta }_{m_2}}$ into account. We will first consider the case ${\displaystyle m_1>m_2}$. In this case, the conditional distribution of ${\displaystyle {\theta }_{m_1}}$ given ${\displaystyle {\theta }_{m_2}=y}$ is wrapped normal with conditional mean ${\displaystyle \mathbb{E}({\theta }_{m_1}|{\theta }_{m_2}=y)=y}$ and conditional variance ${\displaystyle \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 ${\displaystyle {\theta }_{m_1}}$ and ${\displaystyle {\theta }_{m_2}}$. We will use this result to calculate the joint probability distribution

By applying the Leibniz integral rule in ${\displaystyle (\ast )}$ we obtain the probability density of ${\displaystyle {\theta }_{m_1}-{\theta }_{m_2}}$: