Angular and linear speed cells in the parahippocampal circuits

An essential role of the hippocampal region is to integrate information to compute and update representations. How this transpires is highly debated. Many theories hinge on the integration of self-motion signals and the existence of continuous attractor networks (CAN). CAN models hypothesise that neurons coding for navigational correlates – such as position and direction – receive inputs from cells conjunctively coding for position, direction, and self-motion. As yet, very little data exist on such conjunctive coding in the hippocampal region. Here, we report neurons coding for angular and linear velocity, uniformly distributed across the medial entorhinal cortex (MEC), the presubiculum and the parasubiculum, except for MEC layer II. Self-motion neurons often conjunctively encoded position and/or direction, yet lacked a structured organisation. These results offer insights as to how linear/angular speed – derivative in time of position/direction – may allow the updating of spatial representations, possibly uncovering a generalised algorithm to update any representation.


Introduction
The brain is constantly bombarded with all types of information that need to be related to each other in order to build a coherent picture of one's surroundings (or of the situation one is experiencing at a given time). A broad range of studies have led to assert that one of the main roles of the hippocampal region is to do that: integrate multimodal information -coming from a variety of sensory and associative cortices, as well as deeper structures -to build a dynamic representation of an environment or an event [1][2][3][4] . The accurate updating of these representations -and their comparison with previously stored representations as well as projections of likely near futures -could allow one to evaluate the outcome of a range of decisions in order to react adequately. In the context of spatial cognition, information integration is implemented in interconnected subareas of the hippocampal region through neurons coding for specific instantaneous navigational features such as position (place cells) 5 , direction (head direction cells) 6 , local metrics (grid cells) 7 and boundaries (border cells) 8,9 . Successful navigation is also thought to depend on the accurate updating of spatial representations, which themselves would hinge on self-motion signals and their integration with both positional and directional information [10][11][12][13] . Despite their crucial role, where and how self-motion signals are integrated remains largely elusive.
Linear speed modulation has so far mainly been reported in the CA1 region of the hippocampus and the medial entorhinal cortex (MEC) in conjunction with positional information or as an (apparently) self-standing code (speed cells) 14,15 . In addition, speed has been reported to influence oscillatory activity recorded in the hippocampal field potential where the theta power seems correlated to locomotory activity [16][17][18] . In contrast, angular velocity coding has not been yet established in principal neurons of the hippocampal region. Most reports come from recordings of subcortical structures (e.g., lateral mammillary nuclei, dorsal tegmental nucleus), linked to the processing of vestibular information [19][20][21][22] .
Most crucially, it has remained unclear whether angular velocity coding reflects an intelligent design. In fact, some theoretical models within the class of continuous attractor networks (CANs) require that neurons coding for instantaneous navigational correlates -such as position and direction -receive inputs from cells conjunctively coding for position, direction and selfmotion 11,13,[23][24][25][26][27] . These conjunctive cells -sometimes refers to as the 'hidden layer' -are hypothesised to mediate the shift of activity from position (or direction) at time t to the next position (or direction) at time t+1 (Figs. 1c and 2b). Exciting evidence compatible with such CAN models was recently provided by investigations of the Drosophila melanogaster central complex, where head direction cells whose activity was modulated by angular velocity were shown to be organised in a ring according to their preferred direction 28,29 . These reports and the underlying mechanisms are likely restricted, in insects, to directional coding.
To understand the circuit mechanism by which spatial representations can be updated in mammal, we carefully mapped the activity of individual neurons recorded in three main areas of the rat parahippocampal region: the MEC, the parasubiculum and the presubiculum. We specifically investigated whether these neurons could respond to both linear and angular selfmotions signals. Our study revealed the existence of parahippocampal neurons coding conjunctively for direction, position and self-motion, possible evidence for the elusive "hidden layer", pillar of many CAN models. However, the thorough examination of self-motion neurons exposed several inconsistencies with such models, highlighting a need for their revision. In particular, it appears that direction, position and speed selectivity are randomly admixed with each other, perhaps pointing to a capacity for self-organization of these circuits.
Given that linear and angular speed are the derivative in time of position and direction respectively, this may reflect a general self-organizing derivative algorithm for the updating of any type of information.  30 .

Angular velocity coding in the parahippocampal region
First, we sought to determine whether parahippocampal neurons responded to the angular head velocity (AHV) signal, which is the derivative of head direction in time. To that end, we computed, for each cell, its angular velocity score as the Pearson product moment correlation between the instantaneous value of angular velocity and the firing rate of the cell across the recording session (see Methods and extended data Fig. 4a). We defined cells as AHV modulated when their score was greater than the 99th percentile of the shuffled distribution. As per our definition, AHV cells are neurons whose firing rate is positively modulated by angular velocity, meaning that these cells are more active when the animal is turning its head.
About half of the AHV cells had their activity modulated solely when the animal had its head turning only in one direction, either clockwise (CW) or counterclockwise (CCW). The other half of the AHV cells were bidirectional (BiDir) and modulated by angular motions in both directions (MEC: 26 Fig. 1h).
To reduce the dependence of our results on a unique scoring method, we further analysed all recorded cells with two additional methods. Given that the Pearson correlation would present higher scores for linear modulation, we selected methods that would allow us, for example, to select cells responding to a specific speed band. The first additional score was inspired by a method commonly used to characterized spatial modulation 31 . The "Skaggs score" of each cell was obtained by computing the information per spike that each neuron conveyed about the running speed or AHV (see Methods). The second additional method relied on a generalized linear model (GLM) approach 14 to calculate a neuron firing profile as a function of velocity values. For this approach, velocity values were treated as categorical variables; therefore, no linear dependence was implied (see Methods and extended data Fig. 4b). While these additional methods revealed a substantial fraction of AHV modulated cells not captured by linear scoring method (and vice-versa), the cell populations yielded by the three methods were significantly overlapping (Binomial test, p < 0.001, extended data Fig. 4c). Out of 182 AHV cells solely picked up by the GLM method, about half of them (54%) were anticorrelated with the absolute value of angular velocity -meaning that their activity was maximal when the animal was not turning its head. Only a small number of them (11%) were responsive to a particular value of angular head velocity (measured with a gaussian fit), while the rest did not show a modulation profile with a simple shape. Given the dominance of linear coding and to facilitate the comparison with linear speed analyses, we decided to use the more conservative Pearson scoring methods for all further AHV analyses.

Parahippocampal neurons upstream of the entorhinal cortex code for linear speed
Once established that angular velocity coding was widespread across several parahippocampal areas, we tested whether linear speed coding could also extend beyond the medial entorhinal cortex (MEC). To that end, we determined the speed score of each cell as a Pearson product moment correlation between the instantaneous value of rectilinear speed and the firing rate of the cells across the recording session (see Methods and extended data Fig. 4a As for AHV cells, speed cells were stable across time (Fig. 2f).
We further analysed speed modulation with two "non-linear biased" methods (i.e., "Skaggs information per spike" and GLM) similar to those we used for AHV cells. The population of speed cells detected with these methods significantly overlapped with those obtained with our linear scoring method (binomial test, pvalue < 0.001; extended data for GLM in Fig. 4c).
Similar to AHV scoring, the Pearson method yielded the lowest number of speed-modulated neurons (Skaggs and GLM approaches respectively resulted in 334 and 551 selected neurons).
Out of 296 speed-modulated cells solely picked up by the GLM method, about one fifth were anticorrelated with the absolute value of speed and about one sixth showed a very low linear correlation. Only a very small number of them (8%) were responsive to a particular speed band (measured with a gaussian fit), while the rest did not show a modulation profile with a simple shape. In order to allow for comparison with the latest reports of speed cells in the MEC, we decided to use the more conservative Pearson scoring methods for all further analyses.

Conjunctive coding of primary (place, direction) and derivative (velocity) signals
Based on the fact that angular head velocity (AHV) and speed are the derivative in time of head direction and position respectively, we defined positional and directional signals as primary and self-motion signals as derivative.
Given the key role of a "hidden layer" of cells presenting conjunctive coding for continuous attractor network (CAN) theories 11, 32 , we next sought to determine to which degree selfmotion signals are co-existing in conjunction with other types of coding at the unit level. To that end, we computed the grid and head direction (HD) scores of each recorded unit. We labelled as 'significantly modulated', cells whose score exceeded the 99th percentile of the score distribution calculated on shuffled data (see Methods). According to these parameters, the majority (80.4%) of AHV cells coded for at least one other feature (HD: 55.2%; grid: 13.4%; grid-by-HD: 5.6%; rectilinear speed: 33.7%; Figs. 3a-b, 3e and 4a-c). These percentages were similar to the percentages observed in the general population (HD: 53.3%; grid: 17.8%; grid-by-HD: 7.7%, rectilinear speed: 19,2%; Fig. 3a). A similar distribution of conjunctive coding was observed among speed cells (HD: 35.4%; grid: 18.7%; grid-by-HD: 3.2%; AHV: 29.9%; Figs. 3c-d, 3f and 4a-c). We observed all possible types of conjunction of code including AHV-by-HD (the hidden layer of directional CAN models) and grid-by-HDby-speed (the hidden layer of positional CAN models). This result contrasted significantly with previously published studies of a predominantly self-standing code in speed cells recorded in the MEC superficial layers 15 .
To test whether these differences could be explained by regional variations and an over-

Derivative signals seemed encoded differently than primary signals
In order to grasp whether derivative signals (i.e., speed and AHV modulation) were encoded in a similar fashion as primary signals (i.e., position and direction modulation), we compared the firing properties of each class of neurons. We observed that cells coding for derivative signals exhibited higher average firing rates than cells coding for primary signals (t-test: pvalue < 0.001; extended data Fig. 6a). They also showed a shorter average inter-spike interval (t-test: pvalue <0.001; extended data Fig. 6c) and a larger peak firing (defined as the fifth quintile of the rate distribution, t-test: pvalue <0.001; extended data Fig. 6b). It is important to notice that, contrary to what would be expected from CAN models, unidirectional AHV cells were not silent when the rats were not turning their head. Instead, they decrease their rate in response to movement in one direction and increase it in the other. The differences in firing properties between primary and derivative neurons could be explained by the fact that the monotonic firing profiles used to encode motion signals is less sparse than the receptive-field coding of grid and HD cells which are largely silent outside their firing field. To test this hypothesis, we calculated the percentage of the correlate values at which each cell fired more than its average firing rate (see Methods and extended data Fig 6d). Derivative cells showed a significantly larger percentage (mean: 61.3%, std: 15%) than primary cells (mean: 29.8%, std: 9.4%; t-test: To further characterize whether derivative signals are monotonically encoded in the parahippocampal region, we fitted the rate-response tuning curve of both AHV-and speed modulated neurons with either a linear or a sigmoid function (see Methods). The majority of AHV cells (68 %) were better described by a linear fit, compared to a sigmoidal fit (extended data Figs. 3a-c). As the steepness parameter of the sigmoidal fits was generally low (mean: 0.47 (rad/s) -1 , std: 0.16 (rad/s) -1 ), we concluded that most AHV cells followed a quasi-linear rate function. In contrast with the AHV population, the sigmoidal fit with low steepness was slightly more predominant among speed cells (56%). This result seemed to be due to a saturation observed at high speeds and was probably related to low sampling in that speed band (extended data Figs. 3a-c). Yet, also in this case, the fitted sigmoids generally exhibited low steepness (mean: 0.05 (cm/s) -1 , std: 0.025 (cm/s) -1 ). We thus concluded that most speed cells also followed a quasi-linear rate function.

Velocity coding is independent of theta modulation
Because the theta rhythm of the local field potential has historically been strongly associated to running speed, we investigated whether both AHV and speed cells had their activity modulated by theta. Following previous work 30 , we defined that a cell is theta modulated when its mean spectral power around the peak in the 5-11 Hz range was at least fivefold greater than the average spectral power in the 0-125 Hz range (see Methods). We observed that only around 40% of the AHV and the speed cells passed these criteria for theta modulation (Fig. 5a), while many of the remaining cells did not show any modulation by theta (Fig. 5b-c). The theta modulation of AHV/speed cells was observed in comparable proportions to those observed in the general population (Fig. 5e). There was no significant correlation between AHV score and theta score, or between speed score and theta score (Pearson correlation: pvalue > 0.05).
Conjunctive coding for grid or head direction did not influence the proportions of velocity cells that were theta modulated (Fig. 5f, t-test: pvalue < 0.05). Theta modulation was uniformly distributed across all layers of each area except for MEC LII which showed more theta modulation and MEC LVI which showed less (extended data Fig. 7a, t-test: pvalue <0.01). The proportion of velocity cells theta modulated in each layer did not differ from what was expected based on the theta modulation observed in the general population, except for MEC LVI speed cells, which showed less theta modulation than expected (extended data Fig. 7b-c, t-test: pvalue <0.05). Together, these results suggest that the code for self-motion in the parahippocampal region seems largely independent of theta modulation at the single cell level.

Discussion
Here we reveal the existence of a network of parahippocampal principal neurons whose activity is modulated by both angular and linear self-motion signals. Extensive mapping of the rat parahippocampal region showed that this network was spread homogenously across all layers of several interconnected areas upstream of the hippocampal formation: the medial entorhinal cortex (MEC), the presubiculum (PrS) and the parasubiculum (PaS). We observed that some neurons modulated by self-motion seemed to only respond to either angular head velocity or linear speed (i.e., an apparent self-standing code). Yet, the majority of those self-motion neurons were also concomitantly responding to spatial or directional information (i.e., a conjunctive code). Such integration at the unit level may be a crucial mechanism underlying the generation and the updating of the representation of position (place and grid cells) and direction (head direction cells) in the hippocampal/parahippocampal circuits. Moreover, given that linear and angular speed are the derivative in time of position and direction respectively, our results may uncover a general algorithm for the updating of any type of information.
In support of this idea of a general parahippocampal mechanism that could compute the derivative in time of other correlates, we demonstrated that both angular and linear self-motion (derivative) signals were encoded in a different manner with respect to positional and directional (primary) signals. It is well documented that head direction, grid, and place cells tend to be active only when an individual is either in a given position or with its head in a given direction. In contrast, we showed that only a negligible proportion of self-motion neurons responded preferentially for a given speed. The vast majority presented a much higher baseline activity and had their firing rate linearly (or quasi-linearly) ramped up in response to increasing speed. This could be the hallmark of a general strategy in the parahippocampal region for the neural coding of scalar quantities whose magnitude has a well-defined meaning -speed and angular velocity in this case -as opposed to neural activity manifolds used to encode position and direction, where coordinates are always relative to a reference frame.
Previously, angular head velocity (AHV) cells had been mainly characterized upstream of the hippocampal circuits, in subcortical structures linked to the processing of vestibular information (e.g., lateral mammillary nuclei, dorsal tegmental nucleus, thalamic nuclei and striatum) 19-22, 33, 34 . In addition, some example of AHV modulation was reported in the retrosplenial cortex 35,36 and has been linked to the accurate processing of visual inputs in the primary visual cortex 37 . Besides, a few examples of hippocampal neurons modulated by whole body motion have been reported in the primate 38 . Furthermore, angular head velocity has been shown to influence the preferred orientation or the selectivity for pitch and azimuth of some presubicular head direction cells 36,37,39 , as well as to modulate the firing rate of a small fraction of fast spiking presubicular interneurons 40 . We observed both bilateral and unilateral (i.e., CW and CCW) AHV modulation in the presubiculum, the parasubiculum and the MEC. We hypothesise that such signals were not reported in previous studies either because of a restricted scope in analyses or because recordings were often clustered in the most dorsolateral part of the presubiculum, at the border of the retrosplenial cortex. To avoid such bias, we recorded from a much larger population spread out across the medio-lateral and antero-posterior axis.
Nevertheless, it is worth noting that our recordings did not show any topographic organisation of self-motion modulation. The lack of topography, in areas near the top of the cortical hierarchy is consistent with AHV being one of potentially several high-level signals, which could be extracted by derivation with respect to time 41 .
Since the discovery of grid cells, many have attempted to understand how such a strikingly regular signal could be generated by individual neurons 11,23,[42][43][44][45] . A speed code is central to much of this theoretical work, and a break-through was the characterization of speed cells and speed modulation in the MEC 14,15,[46][47][48][49] . Interestingly, recent experimental work has demonstrated that grid cell activity is dependent on the integrity of the speed signal 50 , which itself seems driven by the brainstem locomotor circuit 51 . Likewise, the stability in head direction coding seems dependent on angular head velocity (AHV) signal 52 and vestibular inputs 20,53 . Therefore, self-motion signals could be similarly involved in the generation and the maintenance of both position and direction signals. It is interesting to note that MEC LII contains neither HD nor AHV signals, a major exception to the apparent random assignation of correlates, arguing for the importance of a local network in some cases.
Our findings offer robust experimental evidence particularly relevant for the reappraisal of grid and head direction cell generation theories based on continuous attractor neural networks (CAN). CAN are popular theoretical models of how reciprocally interconnected neurons may extract, refine and sustain stable representations of continuous behavioural variables. This function is particularly crucial when afferent signals are weak, noisy and/or intermittent. How can such representations be updated? A common view is that neurons coding for "primary" navigational variables -such as position and direction -are connected by a so-called "hidden layer" of cells conjunctively coding for position, direction and self-motion. This conjunctive layer is postulated to drive the activity of the output layer to remain congruent, at any time, with the external world. Although several variants of these mechanisms have been proposed, in general they posit a categorical division of labour between the output units representing the instantaneous variables and the hidden units updating them 23,54 . Recent studies in drosophila melanogaster have indeed shown that a ring attractor network with local excitation and global inhibition underlies the representation of head direction 28,29 . It is possible that similar mechanisms operate in mammals. Yet, it is difficult to imagine how a similar neural structure, neatly laid out and possibly hard-wired in the fly brain, could be engineered in the absence of topography in the mammalian brain, and be retained in areas with such prominent plasticity.
Here we report the existence of cells conjunctively coding for position, direction and selfmotion in both the pre-and the parasubiculum as well as in the medial entorhinal cortex, providing the first recording in the hippocampal region of the conjunctive coding postulated to be found in the "hidden layer" of CAN models. Yet, our findings should not be perceived as a blanket endorsement of such models. Indeed, our results challenge theories assuming very specific architectures, such as those requiring that all speed and AHV cells must conjunctively code for position and head direction. While we observed more conjunctive motion-by-primary signal coding than what was reported in previous studies of speed cells in MEC layer II 15,47,49 , we still observed a significant number of self-standing AHV and speed cells. Furthermore, the distinct types of selectivity not only appeared randomly admixed, but also did not present evidence of clear-cut categorical distinctions between conjunctive and non-conjunctive cells.
We observed that scores (i.e., grid, HD, AHV and linear speed score) showed no correlation, and the observed percentages of conjunctive cells were compatible with a scenario of independent assignment of coding properties. Such results are thus more consistent with selforganizing models than with precisely engineered ones 45,55 . Ultimately, they challenge the very concept of a well-defined representation of the instantaneous variables. In terms of attractor structures, they support extending CAN models to incorporate dynamical variables in what their output "represents", for example including units firing at a baseline rate for zero AHV, and quiescent for high either CW or CCW AHV. However, theoretical work on such dynamical attractors or "structured flow manifolds" is only just starting 56,57 . It will have to take into account the remarkable lack of segregation of the parahippocampal regions in the representation of different correlates. One exception to this observation was the absence of AHV and HD cells in MEC LII, suggesting that the AHV signal is needed locally, among the same cells coding for the primary signal it serves to update. In the framework of attractor network models, this leads us to hypothesize that the neural circuitry representing and updating dynamical variables is one and the same, not only at the level of the brain area, but also locally within a layer and even within individual cells. Interventional studies will be necessary to test how such randomly assorted conjunctive coding of primary and self-motion signals can be necessary for accurate spatial cognition.
Historically, running speed has been reported to show a strong correlation with the amplitude of theta oscillations recorded in the local field potential of freely behaving rodents [16][17][18] .
Likewise, many place and grid cells exhibit a strong modulation of their firing rate following those theta oscillations, either in a phase-locked or in a phase-precessing manner [58][59][60] . In line with these observations, many models point to theta oscillations as inherent to the generation of the grid signal. Among them, the oscillatory interference-based models of grid cells assume a velocity input to the grid network composed of translational speed and movement direction 42,44 . Here we show that only 40% of our speed cells and AHV cells show a strong modulation by theta, a finding that should be taken into account in refined versions of such models. This apparent decoupling between speed/velocity signals and theta may be surprising for some.
However, the role of theta in spatial coding had already been recently challenged by two lines of evidence. One is that modulation of the septal oscillatory activity had no apparent consequence on grid signal maintenance 61 , therefore suggesting that non-theta septum correlates -such as attention -were involved in the grid cell signal disruption observed after septal inactivation 62,63 . Second is that no stable theta oscillation had been recorded in the hippocampal local field potentials of bats, who do exhibit both place and grid coding 64 . Yet, recent results show that bats seem to exhibit theta-band modulation of grid firing and matching phase precession when considering a non-stable theta frequency 65 . Further targeted interventional studies would be essential to determine the ambiguous role of theta in the generation and maintenance of the grid signal. Discriminating between the theta modulated and the non-theta modulated self-motion neurons, revealed by our work, would be a chief target for these studies.
In conclusion, we provide clear evidence of a widespread parahippocampal network involved in linear and angular speed coding that could have a crucial role in the updating of the cognitive map, or perhaps be part of the map itself. The existence in the hippocampal region of neurons conjunctively coding for self-motion, position and direction would prima facie appear to fill a gap in the framework of continuous attractor network models. Yet, the analysis of these neurons reveals features divergent from those expected of units serving solely to update the instantaneous representation of static variables. Our work urges for the revision of such models so that they can (i) express dynamical continuous attractors and (ii) account for the apparent random nature of the spatial code, as well as its peculiar lack of a clear organization. We hypothesize that derivative algorithms may have a generalized role in the updating of continuously varying information, not just of a spatial nature. Further studies, with either targeted inactivation of neurons or testing of non-spatial correlates will be necessary to establish whether one of the main roles of the parahippocampal region is to ensure the accurate updating of the hippocampal representation.

Acknowledgement
We thank Edvard and May-Britt Moser for consent to use of these data 30     Comment on scoring method: The correlation method and the GLM approach have rather complementary advantages and weaknesses. The former assumes that the modulation of the firing rate is linear or close to linear but does not require any further arbitrariness in the fixing of additional parameters. The latter is able to detect more general forms of modulation but is subject to choices (e.g., binning size, regularization procedure, form of the link function) that are, at least to some degree, arbitrary. Since the majority of the modulated cells are found with both methods, we proceeded in the further analysis with the simpler one, which also allows a more direct comparison with previous literature.

Subject and surgeries
All the data presented here have been previously published 30 but was re-analysed for this manuscript. The neuronal activity was recorded from 28 male Long-Evans rats (3-5 months old, 350-450 g at implantation, housed and food deprived as in described previously). All experiments were approved by the National Animal Research Authority of Norway. Tetrodes configuration and surgical implantations are described in previously published work 30 .

Data acquisition and training procedures
General data acquisition procedures have been described previously 30 . In brief, rats were trained to collect food crumbs thrown randomly into a 50-cm-high square or circular box with black floor and black walls surrounded by black curtains. Each trial lasted 10 mins.

Spike sorting
The spike detection in the local field potential and sorting was performed as previously described 30 . Spike sorting was performed offline using graphical cluster-cutting software.

Position
The position of the animal was estimated from the coordinates of two light-emitting diodes (LEDs) on the head of the animal. The X and Y coordinates of both the LEDs where smoothed with a gaussian filter with a 250 ms standard deviation, chosen to match the smoothing performed on the firing rate (see below), and the average between the two LED positions was used as the position of the animal.

Head Direction
The head direction (HD) was calculated as the angle between the line connecting the small LED to the big LED and the x axis. HD is expressed in radians, 0 meaning that the rat head is aligned with the x-axis, facing right.

Linear speed and angular head velocity
Speed was calculated as the modulus of the vector difference between the position at time t and the position at time t+1. Angular head velocity was calculated as the signed difference between the head direction at time t and the head direction at time t+1.
The absolute value of the angular head velocity was used for the scoring of bidirectional angular head velocity cells (see below). No further smoothing was applied.

Firing rate calculation
Instantaneous firing rate was obtained dividing the whole session in bins of 20 ms, coinciding with the frames of the tracking cameras. The spike count in each time bin was then calculated and divided by the temporal width to obtain the rate. The rate profile was smoothed with a 250 ms wide Gaussian filter.

Speed filtering
The analysis on speed and angular velocity was performed on movement periods, defined as the ones in which the animal speed was >2 cm/s. A speed filter was applied on the timeseries of each correlate, discarding the time points for which the instantaneous speed was below 2 cm/s, that were excluded in the subsequent analysis.

Spatial rate maps
The histograms for spike count and time spent in each location were constructed using equally spaced bins of 2-cm linear size. Each bin of the rate map was obtained as the ratio between spike count and time spent, smoothed with a gaussian filter with standard deviation of 4 cm.

Directional rate maps
The histograms for spike count and time spent facing each direction were constructed using equally spaced bins of size 6 degrees. Each bin of the rate map was obtained as the ratio between spike count and time spent, smoothed with a gaussian filter with standard deviation of 6 degrees.

Shuffling
Chance-level statistics was calculated for a given variable W through a shuffling procedure.
For each repetition, the firing rate time series was time shifted of a random interval of at least 30 seconds, with the end of the trial wrapped to the beginning.
This procedure was repeated 100 times for each cell, and the shuffled score for variable W was calculated for each instance to compose the chance level statistics.
For cell classification, all shuffled data from the same region were pooled together and the 99th percentile of the distribution was used as a classification threshold.

Speed Score
The speed score was defined as the Pearson product-moment correlation between the cell's instantaneous firing rate and the instantaneous speed of the animal, across the whole recording session. This yields a score ranging from -1 to +1.

Unidirectional angular velocity score
The unidirectional angular velocity score was defined as the Pearson product-moment correlation between the cell's instantaneous firing rate and the instantaneous angular velocity of the animal. Positive values of angular velocity correspond to clockwise head movement.
Cells that had a score greater than the 99th percentile of the shuffled distribution were classified as clockwise modulated (CW), while cells whose score was lower than the 1st percentile were classified as counterclockwise modulated (CCW): they significantly code for head movement in the counterclockwise direction. CW and CCW populations are mutually exclusive by construction.

Bidirectional angular velocity score
The bidirectional angular velocity score was defined as the Pearson product-moment correlation between the cell's instantaneous firing rate and the absolute value of the instantaneous angular velocity of the animal. Cells in this population increase their firing rate in response to head movement regardless of the direction. The unidirectional and bidirectional angular velocity scores are not mutually exclusive by construction. A cell with a strong ramping up of the activity for positive value of AHV, and a mild sensitivity to negative values, for example, could be selected as both a bidirectional and a unidirectional CW cell. A strong modulation from both negative and positive angular velocities would be picked up only by the bidirectional score, while a strictly monotonic increase (or decrease) along the whole range of velocities would give a high unidirectional score.

Mean vector length (head-direction score)
The mean vector length score is calculated from the head-direction tuning map of a given cell as the sum: , where ! is the orientation in radians associated with bin i and ! is the firing rate in the bin.
The sums run over all N directional bins, and the modulus of the resulting complex number is taken. Head direction was binned in bins of 6 degrees and smoothed with a gaussian filter with a standard deviation of 6 degrees.

Grid score
The grid score was calculated from the spatial autocorrelogram of a given cell with a procedure similar to 66 . After exclusion of the centre of the autocorrelogram, the Pearson correlation of the autocorrelogram rotated by 30,60,90, 120 and 15 degrees (+-3 degrees offsets) was considered. Only bins closer to the centre than an outer radius s were included in the calculation of the correlation. Given s, the grid score was defined as the difference between the average of the maximum correlations around 60 and 120 degrees (+-3 degrees offsets) and the average of the minimum correlations around 30,90 and 150 degrees (+-3 degrees offsets). The final grid score of the cell was then defined as the maximum grid score over values of s ranging from twenty to forty bins, computed at intervals of one bin.

Theta index
Theta modulation of individual cells was estimated from the frequency power spectrum of the spike-train autocorrelation histogram of the cell.
A cell was defined to be theta modulated if the mean power in a 2 Hz window centred in the peak in the 5-to 11-Hz frequency range was at least fivefold greater than the mean spectral power in the 0-to 125-Hz range.

Estimation of the significance of overlaps between cell populations
The observed overlaps between cell populations were compared to the ones that would result from the statistical null hypothesis of independent random assignment with a two-sided binomial test. The probability of observing an overlap of size k between two populations of sizes # and $ , independently drawn from a total number of cell N is given by: Where #$ = # $ and # = # / , $ = $ / .

Information analysis
The information per spike conveyed by each cell about the correlate of interest (speed or angular head velocity) was calculated using the formula: Where i is the index of the correlate bin, ! is the probability of observing the correlate in bin i (i.e. the normalized occupancy), ! is the average firing rate of the cell in bin i, and is the average firing rate of the cell.
Speed was divided in 2 cm/s bins in the range 2-50 cm/s (as in all analysis, stillness periods were excluded), while angular head velocity was divided in 0.5 rad/s bins, in the range -5-5 rad/s. Cells were considered to carry significant information about the correlate if the observed information rate exceeded the 99th percentile threshold of the null distribution obtained by shuffling the cell firing rate values (1000 shuffles per cell).

Generalised linear model (GLM) analysis
We analysed the effect of each correlate (speed and angular head velocity) with a linearnonlinear Poisson spiking GLM model.
This model assumes that the firing rate of the cell depends on the value of the correlate as: Where ! ( ) is a one-hot vector (i.e., a vector with only one non-null element) indexing which value the correlate is taking at time t, and ! are the coefficients of a linear filter quantifying the contribution of each value of the correlate to the firing rate of the cells.
The model is fitted using the python module statsmodel.api, which finds the set of parameters ! maximising the log-likelihood of the observed spikes, subject to an elastic-net regularization constraint.
To perform the fitting procedure, the speed values have been binned in 10 bins in the range 2-50 cm/s, and the angular velocity values divided in 10 bins in the range -3-3 rad/s. Cells were considered significantly modulated by a correlate if the log-likelihood of the best fit was significantly larger than the value obtained with only the average firing rate as a predictor. Significance was estimated with a 10-fold bootstrapping procedure to extract the confidence interval of the observed log-likelihood.

Tuning curve fitting
Two different functional forms were fitted and compared to the tuning curve of modulated cells.
A linear model:

= +
And a sigmoid model: Here x is the value of the correlate (speed or angular head velocity) and r is the average firing rate of the cell at that value of x. Tuning curves were rescaled by their maximum value, in order to match the two model by number of parameters.
The R2 fitting scores were then compared for each cell. Cells with a linear R2 greater than the sigmoid R2 were classified as linear, and vice versa.
To quantify the number of cells that showed modular coding (i.e., an increased rate in a particular speed or angular velocity band), we fitted their rescaled tuning curves with a gaussian profile: Cells were labelled as gaussian if the fit yielded a R2 score > 0.5 and the average a of the gaussian laid within the tuning curve interval (2-50 cm/s for speed, -3-3 rad/s for AHV).

Sparsity calculation
To quantify the sparsity of the firing of primary and derivative cells, we calculated their tuning curves (as described above) as a function of speed (for speed cells), of AHV (for AHV cells), head direction (for HD cells) and position (for grid cells).
We then calculated as the sparsity the percentage of the bins of the tuning curve that had a firing rate value larger than the average overall firing rate of the cell.