Grid cells encode local head direction

Grid and head direction codes represent cognitive spaces for navigation and memory. Pure grid cells generate grid codes that have been assumed to be independent of head direction, whereas conjunctive cells generate grid representations that are tuned to a single head direction. Here, we demonstrate that pure grid cells also encode head direction, but through distinct mechanisms. We show that individual firing fields of pure grid cells are tuned to multiple head directions, with the preferred sets of directions differing between fields. This local directionality is not predicted by previous continuous attractor or oscillatory interference models of grid firing but is accounted for by models in which pure grid cells integrate inputs from co-aligned conjunctive cells with firing rates that differ between their fields. Local directional signals from grid cells may contribute to downstream computations by decorrelating different points of view from the same location.

These representations include grid codes, in which neurons are active at the vertices of a grid of tessellating triangles that cover the environment an animal explores 5 , and head direction codes, in which neuronal activity is tuned to head direction 6 . Effective navigation and memory require observed data with shuffled datasets generated by allocating spikes to the behavioural trajectory according to the location-dependent average firing rate ( Fig. 2b and Figs. S1, S2). We refer to the resulting plots as 'distributive plots' as they enable assessment of a neuron's directional tuning in comparison to predictions made under the distributive hypothesis. When we tested whether the experimental firing rate for each directional bin in the distributive plots differed significantly from the shuffled data (threshold p < 0.05 after correcting for multiple comparisons made across bins), we found 7.3 ± 3.9 and 4.1 ± 4.5 significant bins out of 20 bins / cell for mice and rats respectively (Fig. 2c). For all mice (n = 34 / 34 cells) and most rats (n = 56 / 68 cells) there were more significant bins than expected based on the shuffled data (for both datasets p < 10 -16 , U > 10, Mann-Whitney U test vs the shuffled data) (Fig. 2c). We obtained similar results when we analysed firing as a function of movement direction rather than head direction, although the effects of movement direction were smaller (Fig. S3a). In contrast to the unimodal directional tuning of conjunctive cells 8 , the directionally binned firing of pure grid cells had multiple peaks and troughs. The orientation of the peaks differed substantially between pure grid cells indicating that they were not driven by common external cues (Fig. S4).
Variation in running speed between different parts of the environment is also unlikely to account for directional tuning as, in agreement with previous studies 19 , firing of most pure grid cells had speed scores below the threshold previously used to identify speed cells (cf. [19][20][21] )(median speed score for mouse grid cells = 0.068 ± 0.18, n = 34, n = 26 / 34 with speed score < 0.1; median for rat grid cells from the rat dataset = 0.038 ± 0.048, n = 68, n = 60 / 68 with speed score < 0.1)(see Fig S5a and Supplemental Data) and directional tuning was independent of a neuron's speed score (Fig S5b).
Together, these analyses indicate that firing of pure grid cells has a multimodal directional structure that is qualitatively distinct from the unidirectional tuning of conjunctive cells .

Pure grid cell firing fields are locally tuned to head direction
If firing by pure grid cells encodes head direction then we expect this to also manifest at the level of individual firing fields. To test this we isolated spikes from each field using a watershed algorithm (44 fields isolated from 13 pure grid cells in 4 mice and 83 fields from 25 pure grid cells in 5 rats) (Fig. 3a) and analysed directional firing separately for each field (Fig. 3b-c). We used the watershed algorithm to avoid potential bias from manual selection of fields and only selected cells for further analysis when the algorithm identified at least two fields.
Individual fields demonstrated clearer directional peaks than for the arena as a whole suggesting a greater degree of directional modulation at smaller spatial scales (cf. Fig. 1b and   3c). To assess whether these peaks could have arisen by chance we generated distributive plots for each field (Fig. 3b, c, Fig. S6). For almost all pure grid cells, the firing rate in at least one directional bin of the plot differed significantly from the corresponding shuffled data after correcting for multiple comparisons across bins (mice: 12 / 13 cells; rats 24 / 25 cells). The number of significant bins per field was substantially greater than predicted by the shuffled data (mice: 4.3 ± 3.2 bins / field; rats: 2.1 ± 2.8 bins / field, for mice and rats p < 10 -16 , U > 8, vs shuffled data, Mann-Whitney U-test) (Fig. 3d) and did not correlate with bias in the behavioural head direction within the field (Fig. S7). The proportion of fields per grid cell with directional bins that remained significant after correcting for multiple comparisons was 84.6 ± 27.1 % for mice (38 / 44 fields) and 61.7 ± 28.2 % for rats (47 / 83 fields). This head direction dependence at the level of individual fields could not be explained by speed-dependence of neuronal firing (Fig.   S5c). In contrast to unimodal tuning of conjunctive cells and head direction cells, when fields from pure grid cells had significantly modulated bins within their distributive plot, their relative orientation was often consistent with multidirectional tuning (cf. Fig. 3c and Supplemental data).
Although the presence of multiple preferred directions in the fields of pure grid cells limits the ability to directly compare the strength of directional modulation with unidirectional conjunctive cells, each cell type nevertheless demonstrated a similar proportion of directionally modulated fields ( Fig S8). Together, these analyses demonstrate directional modulation of pure grid cell firing that is not explained by the distributive hypothesis (cf. Fig. 2a) or by potential confounds introduced by speed-dependent neuronal firing.
Inspection of individual fields from the same cell suggests that their directional modulation differs from one another (Fig. 3c). If this is the case, then directional tuning of fields from the same cell would on average show little or no correlation. In contrast, if directional modulation is non-local, then head direction tuning of fields from the same cell should be correlated. We found that fields from the same pure grid cell on average showed weak or no correlation in their directional firing (median correlation for mice: 0.016 ± 0.36; for rats: 0.12 ± 0.28) (Fig. 4a-b), consistent with local directional modulation. In contrast, correlations between fields from the smaller number of spikes available for the analysis. Directional firing remained uncorrelated between fields when considering only fields adjacent to the walls of the arena, or only fields in the centre of the arena, indicating that the location dependence of directional tuning is also not related to the proximity of fields to the borders of the arena ( Fig. 4d and Fig S9). Examination of correlations between fields as a function of their distance from one another also did not reveal any local patterns in their directional tuning (Fig. S10).
Together, these analyses indicate that individual firing fields from pure grid cells are tuned to head direction, with many fields having multiple preferred directions. As directional modulation within each grid field is independent of direction modulation of other fields from the same grid cell, these data indicate that directional modulation of grid cell firing is location specific.

Stability of directional modulation
To assess the stability of directional tuning of pure grid cells we first quantified the correlation between global classic head direction plots generated from the first and second half of each recording session (Fig. 5a-c). Across the population of pure grid cells from mice the correlation coefficients were skewed to positive values (median Pearson correlation = 0.47 ± 0.39 ) and differed significantly from zero (Fig 5c). To address the possible contribution of biases in the running direction between different parts of the firing field (cf. 18 ), we compared correlations for shuffled and observed data from distributive plots for the two halves of each recording session ( Fig. 5d). We found that 16 / 34 mouse grid cells had percentile scores in the top 95 % of correlation scores generated from the shuffled data. For individual fields, their classic head direction plots were also positively correlated (median Pearson correlation = 0.29 ± 0.31)( Fig.   5e-g) and 8 / 44 fields had correlation scores in the 95th percentile of the shuffled data (Fig. 5h).
These results are similar to predictions from 'ground truth' data generated by simulations of spiking models of directionally modulated pure grid cells (see below and Fig. S14). These simulations predicted correlation scores between the first and second half of each session for global (median = 0.38 ± 0.27 ) and local (median = 0.58 ± 0.30 ) directional tuning that were comparable to the experimental data ( Fig S14). They also predicted a similar proportion of correlation scores in the 95th percentile of the corresponding scores generated from the shuffled data (all cells: 12 / 19; fields: 29 / 106)( Fig S14). Thus, global directional tuning of pure grid cell firing appears stable at the time scale of single recording sessions.

Integration of inputs from conjunctive cells accounts for directional modulation of grid cells
Our analyses imply that spatially specific multi-directional firing is a core feature of the activity of pure grid cells that does not appear to be predicted by existing models for grid firing [9][10][11][12][13]22 (Supplemental Table 1). For example, a requirement of many continuous attractor models is that each grid cell has a single preferred direction of input 9,11 . Consistent with this intuition, simulations using experimentally recorded trajectories as input to each of two previously described continuous attractor models 9,13 and two previous oscillatory interference models 12,23 did not generate directionally modulated grid firing fields ( Fig. 6 and Fig. S11). Because these simulations used real trajectories, their results also provide further strong evidence against the hypothesis that directional tuning identified in our analyses of experimental data is a result of movement-related variables.
What mechanisms might then explain the local directional firing of pure grid cells? We reasoned that three ingredients may be important (Fig. 7a). First, upstream directional grid signals, which may originate from conjunctive cells 8 . Second, differences between fields of the upstream cells in their mean firing rate. This is consistent with grid fields encoding positional information in their local firing rate 24 . Third, convergence of multiple upstream signals that have co-aligned spatial firing fields and different head direction preferences. This is consistent with the prominent projections from deep layers of MEC, where conjunctive cells are most frequently found, to superficial layers where pure grid cells are most frequently found 8,[25][26][27] . Models have been proposed previously that combine the first and third ingredients to account for firing of pure grid cells (cf. 15,25 ). We refer to these as spatially uniform conjunctive cell input models as each of the upstream conjunctive cells' fields have similar mean firing rates. We consider here the possibility that introducing the second ingredient would account for the experimentally observed local directional tuning of pure grid cells. We refer to this new model as a spatially non-uniform conjunctive cell input model (Fig. 7a). We simulated both classes of model using an anatomically reconstructed stellate cell configured to account for experimentally recorded excitable properties and with synaptic inputs from conjunctive cells distributed across its dendritic tree ( Fig. 7a and S12). Because the strength of in vivo synaptic inputs to stellate cells is unknown, we simulated versions of the uniform and non-uniform input model with relatively low and relatively high synaptic conductances.
We found that the spatially non-uniform conjunctive cell input model accounted well for the global and local directional modulation of pure grid cell firing fields. When we simulated this architecture with experimentally recorded trajectories as inputs, the firing of the model grid cell across an environment was directionally modulated (Fig. 7b) and each of its fields had a distinct directional modulation (Fig. 7c), resembling the experimentally recorded directional tuning of pure grid cells ( Fig. S6 and cf. Fig. 7c with 3c). The proportion of directionally modulated fields and the number of significant bins per field in their distributive plots generated by the low and high conductance versions of the non-uniform model spanned the range of the experimental observations from rats and mice (Fig. 7d). Similar to our experimental observations (cf. Fig. 4b), correlations between fields from the same simulated pure grid cell were weak in this model (median correlation = -0.0061 ± 0.11) (Fig. 8b) and between field correlations differed from the within field correlations (Fig 8c), also resembling the experimental data (cf. Fig. 4c)(p = 2.2 x 10 -16 , D = 0.62). In contrast to the non-uniform input models, versions of the models with spatially uniform conjunctive cell input generated fewer directional fields with fewer significant directional bins per fields (Fig. S13). Despite their weaker directional modulation, the fields generated by the spatially uniform conjunctive cell input versions of the model were nevertheless positively correlated (median correlation = 0.13 ± 0.26), which is consistent with a lack of local directional modulation (Fig. 8b), and the between field correlations did not differ from the corresponding within field correlations ( Fig. 8b)(p = 0.14, D = 0.13). Thus, integration of input from co-aligned conjunctive cells can account for the local direction selectivity of grid cell firing, but only when the conjunctive cells have fields that differ from one another in their mean firing rates. Because the models with uniform and non-uniform conjunctive cell inputs generate spikes through biophysically plausible integrative mechanisms, we used their outputs to further validate the analyses of our experimental data. Systematic variation in the length of the simulated recording indicated that the typical duration of the recording sessions we used experimentally is sufficient to reliably detect directional tuning at a population and local level (Fig. S14). Detection of correlations between the first and second half of each recording session was less reliable, with large variability in the correlation scores obtained. This is also consistent with our experimental data (cf. Fig. 5c, g and Fig. S14g, h).

Discussion
Our analyses show that the firing fields of pure grid cells are locally tuned to head direction. This property of grid cells is not predicted by existing models but can be accounted for by integration of signals from conjunctive cells with spatially non-uniform firing rates. Spatial signals that distinguish points of view within an environment may be critical to downstream computations, for example pattern separation functions of the dentate gyrus and encoding of head direction by place cells.
The directional modulation of grid cell firing that we describe here is qualitatively distinct from previously described properties of conjunctive cells 8 . Whereas conjunctive cells are tuned to a single preferred head direction, we find that pure grid cells are tuned to multiple directions (Figs. 1 -4). Whereas for conjunctive cells directional tuning of individual firing fields is strongly correlated and therefore global, for pure grid cells directional tuning of individual fields differs and is therefore local (Fig. 4). These differences suggest that directional tuning of conjunctive cells and grid cells result from distinct circuit mechanisms. They also suggest reasons why directional tuning of pure grid cells has not previously been described. For grid cells summed directional tuning from several fields generates global multi-directional modulation (cf. Fig. 1a(iii) that is not typically considered in analyses of directional firing and is challenging to detect with shorter duration recordings (cf. Fig. S14).
Embedding of location-dependent directional information within grid cell firing fields has implications for the organisation of spatial computations in entorhinal circuits and downstream structures. While more complex models that account for location-dependent directional firing of pure grid cells could be envisaged (Table S1), integration of input from conjunctive cells with spatially non-uniform firing rates provides a parsimonious explanation for the location-dependent directional firing of grid cells (Fig. 7). This explanation is consistent with models in which grid codes originate in conjunctive cells 15,25,28 , but differs in that in the model we propose here each of a conjunctive cell's firing fields has a different mean firing rate. It is this difference in firing rate that leads to the location-dependent directional tuning of pure grid cell firing. According to this scenario, investigation of the mechanisms underlying the grid code would most productively focus on conjunctive rather than pure grid cells. This scheme may also resolve discrepancies between functionally and anatomically defined cell types in MEC 29 . For example, entorhinal stellate cells, which are a major input to the hippocampal dentate gyrus and are required for spatial memory 30,31 , appear to be the most numerous grid cell type, but less than half of stellate cells are grid cells 32 . This functional divergence is consistent with all stellate cells implementing similar cellular computations 33 , but with the emergence of grid firing patterns depending on the identity of their dominant synaptic inputs.
A key future question is how directional information encoded by pure grid cells is used by downstream neurons. One possibility is that by disambiguating different views at the same location, directional firing by pure grid cells may facilitate pattern separation functions of the immediately downstream dentate gyrus [34][35][36] or CA1 place cells 37 . For example, recently described directional selectivity of place cells is also multi-modal and location-dependent 37 , and could be driven by the location-dependent directional tuning of pure grid cells that we describe here. Thus, direction-dependent grid firing may be a powerful computational feature of representations of visual and conceptual scenes 4,38,39 .  determined by the spatial firing rate map. Circles represent recorded spikes and squares represent shuffled spikes. 2. Observed and shuffled spike rates binned by head direction (bin width 18°, distributive plot) and visualized either with cartesian (left) or polar coordinates (right). Error bars / shaded areas indicate the 90 % confidence intervals of the shuffled data. Asterisks indicate bins in which the observed data differs significantly from the shuffled data (p < 0.05, two-tailed p value calculated from the shuffled distribution and corrected for multiple comparisons with the Benjamini-Hochberg procedure). 3. As for step 2 but comparing individual shuffles to the overall shuffled distribution. 4. The number of significant bins in the observed data (determined as in step 2) and the shuffled data (determined by repeating step 3 for all shuffles). The example data is from the cell in Fig. 1b. (c) The distributions of number of significant bins per cell differed significantly between observed and shuffled data (n = 34 grid cells from mice, n = 68 grid cells from rats, for both comparisons p < 10 -16 , U > 10, Mann-Whitney U test). , example directional firing rate histogram for observed and shuffled spikes (middle) for the highlighted field from (c)(yellow box) and distribution of number of significant bins from the shuffled data (grey) and the observed data (blue)(lower).The error bars represent the 90 % confidence interval of the shuffled distribution. (c) Distributive plots for each firing field (coloured according to (a), field in yellow box also shown in (b)). The maximum firing rates are shown above the head direction plots. Significantly directional bins are marked with an asterix (*). (d) The number of directional bins (out of 20 bins) that differ significantly (threshold p < 0.05 after correction for multiple comparisons) from the shuffled data. The numbers of significant bins differed between observed and shuffled data (n = 44 fields from 13 grid cells, n = 83 fields from 25 grid cells, p < 10 -16 , U > 8, for all comparisons, Mann-Whitney U test).  head direction plots for spikes from the first and second half of the session for the neuron in (a) and comparison of the correlation coefficient for the observed data (blue) with the distribution of coefficients generated from the shuffled data (grey). (f) Classic head direction plot for the first (green) and second half (blue) of the recording session for the pink grid field shown from (e). (g) Distributions of Pearson correlation coefficients for all fields calculated as in (f). The median (0.29 ± 0.31 for mice) differed from zero (Wilcoxon test p = 3.37 x 10 -9 , T = 321, n = 44). (h) Distributive head direction plots for the first and second half of the session for the field in (f) and comparison of the correlation coefficient for the observed data (blue) with the distribution of coefficients generated from the shuffled data (grey). , and binned firing rate as a function of head direction for the colour coded fields for two types of grid models ran on our trajectory data. (c) Number of significant bins per field compared across previous grid cell models and experimental data. Each of the previous models differed significantly from the rat and mouse experimental data (p < 0.05, one-way ANOVA with Games-Howell test). (d) Proportion of directionally modulated fields compared across existing grid cell models and experimental data. are the corresponding shuffled distributions. (d) Comparison of conjunctive cell input models with experimental data, showing significant bins per field (left), and proportion of directional fields (right). The number of significant bins per field in the high conductance (high g ex ) non-uniform model was greater than for the rat data (p = 5.7 x 10 -9 , one-way ANOVA with Games-Howell test) but did not differ significantly from the mouse data (p = 0.299), whereas the low conductance (low g ex ) non-uniform model differed had fewer significant fields that in the mouse data (p = 6.1 x 10 -6 ) but not did not differ significantly from the rat data (p = 0.172). Both high and low conductance uniform models had fewer significant bin per field than the mouse and rat data (p < 0.001, one-way ANOVA with Games-Howell test).  (Fig. 4b). The distribution of correlation coefficients differed between uniform and non-uniform conjunctive input models (p = 3.99 x 10 -6 , D = 0.21, Kolmogorov-Smirnov test). (c) Correlations between firing rate histograms for the same field (grey) and between different fields of the same cell (blue) generated as in Fig. 4. The between and within field distributions differed from one another for the non-uniform model (p < 2.2 x 10 -16 , D = 0.62, Kolmogorov-Smirnov test) but not for the uniform model (p = 0.14, D = 0.13). The half-session between field distributions differed between the non-uniform and uniform models (p = 7.1 x 10 -8 , D = 0.17, KS test).

Animals
All animal procedures were performed under a UK Home Office project license (PC198F2A0) in accordance with The University of Edinburgh Animal Welfare committee's guidelines. All procedures complied with the Animals (Scientific Procedures) Act, 1986, and were approved by the Named Veterinary Surgeon and local ethical review committee.
We report data from 8 mice (4 males and 4 females) from a total of 16 that were implanted with recording electrodes. For 7 of the 8 mice the location of the tetrodes in the MEC was confirmed anatomically after experiments were complete. For the other mouse the tetrode location could not be determined. The age of the mice when tetrodes were implanted was 7 -13 weeks (mean = 10.6 ± 1.7 weeks). Animals were excluded either because the recordings did not identify grid cells (n = 6) or because they had to be terminated before the end of the experiment, because their implant came off (n = 2). The mice used were from the p038 line 40 . Before surgery animals were group housed (3 -5 mice per cage) in a standard holding room on a standard 12 hour on / off light cycle (dark from 7 pm to 7 am). After surgery, mice were single housed in a different holding room in otherwise similar conditions (average temperature 20°C, relative humidity 50 %). Mice were kept in individually ventilated cages containing sawdust, tissues, chewing sticks and cardboard tubes, which after surgery were replaced by a larger cardboard igloo. Two days after the surgery, a running wheel was placed in the cages. Standard laboratory chow and water were given ad libitum .

Microdrive design
We modified previous designs for 16-channel microdrives consisting of 4 tetrodes and an optic fiber 41,42 . A 21 gauge 9 mm long inner cannula (Stainless Tube & Needle Co. LTD) was attached to an EIB-16 board (Neuralynx) using epoxy (RS components 132-605). Tetrodes were made with 18 µm HML-coated 90 % platinum 10 % iridium wire (Neuralynx). We connected grounding wires (1.5 cm long insulated part) to the reference and ground pins of the EIB. The 4 tetrodes were threaded through the inner cannula, connected to the pinholes of the EIB-16 board and fixed with gold pins (Neuralynx, EIB Pins). A 13 mm long optic fibre stub (Plexon, PX.OPT-FS-Flat-200/230-13L) was threaded through the inner cannula between the tetrodes. The wires on the board and the optic ferrule were covered to about 2/5 th of the length of the optic ferrule with epoxy and glued (RS components, 473-455) the tetrodes to four sides of the optic fibre. The tetrodes were cut using ceramic scissors (Fine Science Tools, Germany) to a length of 0.5 mm from the tip of the optic fibre. The tips of the tetrodes were plated in a non-cyanide gold plating solution (Neuralynx). Tetrodes were cleaned, by applying three 1 second 4 µA pulses with the tetrodes as an anode and then plated by passing 2 µA 1 second pulses with the tetrodes as a cathode until their impedance was between 150 and 200 KΩ. General surgical procedures and stereotaxic viral injections were carried out as described previously 43 . We induced inhalation anaesthesia using 5 % isoflurane / 95 % oxygen, and sustained at 1 -2 % isoflurane / 98 -99 % oxygen throughout the procedure (1 L / minute). Before implanting the drive, we injected AAV9-tre-ChR2-mCherry (Gene Therapy Center, University of Massachusetts Medical School)(800 -2000 nl total injection volume, 3 -5 injections sites, 200 -400 nl / site) for additional opto-tagging experiments (data not shown). All animals were injected 3.4 mm lateral relative to Bregma (Table S2). For electrical grounding, we drilled two small craniotomies, and implanted M1 x 4 mm screws (AccuGroup SFE-M1-4-A2) on both sides about 3.4 mm lateral, and 1 mm rostral relative to Bregma.

Microdrive implantation
We stereotaxically lowered the tetrodes 1.5 mm into the brain, beginning 3.4 mm lateral from Bregma (right hemisphere of two mice, and left hemisphere of 14 mice) and targeting the dorsal third of the medial entorhinal cortex. We sealed the outer cannula with sterile Vaseline and fixed the implant by putting dental acrylic (Simplex Rapid powder) around the drive frame and the outer cannula, wrapped the grounding wires around the grounding screws, fixed the wires with silver paint (RS components 101-5621) and applied another layer of dental acrylic to cover the skull and the grounding screws. Mice recovered on a heat mat for approximately 20 minutes and were given Vetergesic jelly (0.5 mg / kg of body weight buprenorphine in raspberry jelly) for 12 hours after surgery.

Open field exploration task
All recordings from mice were performed in an open field arena consisting of a box made from a black metal frame (parts from Kanya UK, C01-1, C20-10, A33-12, B49-75, B48-75, A39-31, ALU3), with removable black metal walls, a polarizing cue (white A4 paper) on one wall, and a floor area of 1 m 2 . A camera (Logitech B525, 1280 x 720 pixels Webcam, RS components 795-0876) was mounted on the top of the frame for motion tracking. To record head direction we used a custom Bonsai script 44 to track red and green polystyrene balls attached to either side of the mouse's head on the left and right sides. The distance between the polystyrene balls was approximately 3.5 cm.
Mice were handled three times a week for 5 -10 minutes for four weeks following surgery. For 3 consecutive days before recording we habituated the mice by allowing them to explore the open field arena for 5 -10 minutes. For recording sessions, mice explored the open field arena unrewarded until they covered the whole area, or for a maximum of 90 minutes. An opto-tagging experiment was performed at the end of each session (data not shown). After each recording session we lowered the tetrodes by 50 µm using the drive mechanism on the implant.
For electrophysiological recording, the 16 channel optetrode was connected to an Open Ephys acquisition board 45 and computer (HP Z440 Tower Workstation i7, 16GB, 512GB SSD, Cat.: J9CO7EA#ABU) using an SPI cable (Intan Technologies, RHD2000 6-ft (1.8 m) Ultra Thin SPI interface cable C3216) and via a commutator (SPI cable adapter board, Intan Technologies C3430 and custom 3D printed holder). Signals were filtered between 2.5 Hz -7603.8 Hz using a second order Butterworth filter implemented in Open Ephys. We aligned position and electrophysiology data using light pulses generated at random intervals (20 s to 60 s) by an LED (light-emitting diode) attached to the side of the open field arena hidden from the mouse but in the field of view of the camera.

Post recording assessment of tetrode locations
To enable determination of tetrode locations, after the last recording day we anesthetized mice using an isoflurane chamber and pentobarbital (100 -150 μl), and applied a 2 second ~20 µA current to burn the tissue at the tip of the electrodes. We then intracardially perfused PBS (phosphate buffered saline, Gibco, 70011044, 10 times diluted with distilled water) for 2 minutes, then 4 % PFA (paraformaldehyde, Sigma Aldrich, 30525-89-4) in 0.1 M PB (phosphate buffer, Sigma Aldrich, P7994) for 4 minutes at a 10 mL / minute flow rate. We left the brains in 4 % PFA in 0.1 M PB for 16 hours, then transferred them to 30 % sucrose (Sigma Aldrich, S0389) in PBS until they sank.
We cut 50 µm sagittal sections of the fixed brains using a freezing microtome. Sections were processed to label them with primary antibody rat anti-mCherry (Invitrogen M11217, 1:1000) followed by secondary antibody goat anti-rat alexa 555 (Invitrogen A-21434,1:1000) and stained with either NeuroTrace 640/660 (Invitrogen N21483, 1:500) or NeuroTrace 435/455 (Invitrogen N21479, 1:500) following procedures described previously 43 . Images were taken on a Zeiss Axio Scan Z1 using a 10x objective and visually inspected to determine the final position of the recording electrodes (see Supplemental Data).

Data analyses
Analyses were carried out using Python (version 3.5.1 in Anaconda environment 4.0) and R version: 3.3.1 (2016-06-21). All code will be available at https://github.com/MattNolanLab . Summary values are reported as mean or median ± standard deviation.

Spike sorting
To isolate spikes from electrophysiological data we used an automated analysis and clustering pipeline based around MountainSort 3 (v 0.11.5 and dependencies) 46 . Python scripts pre-processed the data by converting Open Ephys files to mda format and organized these files together with spike sorting input parameter files. We defined the four channels of each tetrode to be in the same 'sorting neighbourhood'. We excluded broken channels identified during data acquisition.
MountainSort filtered the data from 600 Hz -6000 Hz using a bandpass filter and then performed spatial whitening over all 16 channels to remove correlated noise. Events with peaks three standard deviations above average and at least 0.33 ms away from other events on the same channel were detected. The first 10 principal components of the detected waveforms were calculated. A spike sorting algorithm, ISO SPLIT, was applied to the resulting feature space 46 .
Cluster quality was evaluated using metrics for isolation, noise-overlap, and peak signal to noise 586 ratio 46 . Units that had a firing rate higher than 0.5 Hz, isolation more than 0.9, noise overlap less than 0.05, and peak signal to noise ratio more than 1 were accepted for further analysis. Any units that did not have a refractory period or hyperpolarization component of their spike waveform were discarded. These exclusions were based on visually assessing output figures generated for sorted clusters. No additional manual curation was used.

Classification of functional cell types
To classify recorded neurons we used established grid and head direction scores 8 . Grid scores were defined as the difference between the minimum correlation coefficient for rate map autocorrelogram rotations of 60 and 120 degrees and the maximum correlation coefficient for autocorrelogram rotations of 30, 90 and 150 degrees 47 . The firing rate map was calculated by summing the number of spikes in each location, dividing that by the time the animal spent there and then smoothing the surface with a Gaussian centred on each location bin ( 34 ). Autocorrelograms were calculated by shifting the binned firing rate map 34 into every possible binned position along both horizontal and vertical axes and calculating correlation scores for each of these positions. This autocorrelogram was converted into a binary array using a 20 % threshold on normalized data. If the binary array had more than 7 local maxima, a grid score was calculated. Subsequent parts of the analysis, where correlations between the rotated autocorrelograms were calculated, only included the ring containing 6 local maxima closest to the centre of the binary array, excluding the maximum at the centre. The ring was detected based on the average distance of the 6 fields near the centre of the autocorrelogram (middle border = 1.25 * average distance, outer border = 0.25 * average distance).
To calculate head direction scores, the head direction angles corresponding to the firing events of each neuron were first binned into 360 bins between 0 and 2π. The obtained polar histogram was smoothed by calculating a rolling sum over a 23 degree window. For angles between -179 and 180 degrees in steps of 1 degree, dx and dy were calculated in a unit circle (radius = 1), as and . To obtain the x and y components of the head direction vector, y d = radius sin (angle) x d = radius cos(angle) the head direction polar histogram was multiplied by the dx and dy values, respectively, and normalized to the number of observations in the polar head direction histogram, so that and . The head direction score was then calculated histogram using the Pythagorean theorem as head direction score= .
We defined grid cells as cells with a grid score ≥ 0.4, which was chosen as a conservative threshold (cf. 24 ). We defined head direction cells as cells with a head direction score ≥ 0.5. We defined conjunctive grid cells as cells that passed both head direction and grid cell criteria.

Identification and analysis of individual fields
We identified individual firing fields using methods similar to those used previously to detect place fields 48 . The open field arena was divided into 42 x 42 bins, where each bin contained a smoothed firing rate value calculated by summing the number of spikes at the locations corresponding to each bin, dividing this by the time the animal spent in the bin and then smoothing the surface with a Gaussian ( ) centred on each bin 34 . We next identified the e 2 −x2 bin of the rate map with the highest firing rate. If the rate was higher than the average firing rate plus the standard deviation of the rest of the rate map, we then added adjacent bins if they had a firing rate higher than 35 % of the maximum firing rate within the field. We recursively added to the field further bins that satisfied these criteria with respect to the newly added bins. We accepted a detected field if it had more than 45 bins, but it was smaller than half of the arena. After successfully detecting a field, it was removed from the rate map by replacing the values with zeros and the analysis was repeated until we found no more fields. All detected fields were visually assessed and if a detected field appeared to be a combination of two fields or only part of a field, it was tagged as a false positive to be excluded from the analyses.

Analysis of head direction
Classic head direction polar plots (Figs. 1b, 5b,f) were generated by dividing the histogram of head directions when the cell fired by the histogram of head directions from the trajectory to obtain a firing rate (Hz) for each angle. Each histogram has 1 degree bins with 23 degree smoothing. (Figs. 2, 3, 5d,h, 6, 7) were generated by first dividing the histogram of head directions when the cell fired by the histogram of head directions from the trajectory using 20 bins of width 18 degrees to generate a directional firing rate plot. The shuffling procedure was designed to rule out the possibility that directional firing is a result of different sampling of directions in different bins of the rate map 18 . Before generating shuffled data, the cell's spatial firing rate map was used to allocate firing probabilities to each point on the trajectory such that the total probability was 1 and the local probability was proportionate to that of the firing rate in the given bin of the firing rate map. Shuffled spike locations were generated by random selection with replacement of locations from the trajectory of the animal using the allocated probabilities. Each shuffled dataset was sorted into directional bins in the same way as the experimental data. After generating 1000 shuffled datasets, the median and 90 % confidence intervals were calculated for each bin by ranking the shuffled data.

Two-sample Watson test.
To evaluate whether head direction when the cell fired differed from the head direction of the animal during the time spent in a field, or arena, we performed Watson's two sample test 49,50 for homogeneity on the two distributions using the R package circular ( https://r-forge.r-project.org/projects/circular/ ). To illustrate the distributions of head directions from the trajectory and the normalised firing rate, we made smooth head directions plots (1 degree bins and 23 degree smoothing). On plots where the firing rate and trajectory histograms are shown, the trajectory histogram was re-scaled to fit the plot.
Shuffle tests. Analyses compared the observed data with 1000 shuffled data sets generated as described above. For each bin the percentile position of the observed data relative to the randomized data was used to calculate a p value. The twenty p values within a field were corrected for multiple comparisons using the Benjamini-Hochberg procedure and the number of intervals where the corrected p values were < 0.05 were counted. To obtain null distributions the same analyses were performed for each of the 1000 randomized data sets, where each shuffle was treated like the observed data and was compared to the distributions from the 1000 shuffles. Fields where the percentile score of the observed number of significant bins relative to the null distribution from shuffles was > 97.5 were considered directional.
Correlation between directional firing histograms. To evaluate stability of firing within a session, Pearson correlation coefficients were calculated using the SciPy function scipy.stats.pearsonr for pairs of classic head direction polar plots. Comparison of pure grid cell between and within field correlations used fields with > 500 spikes. To evaluate whether correlations are influenced by the distributive hypothesis, we generated shuffled head direction polar plots for each half of the session as described above for the whole session. We calculated Pearson correlation coefficients between the distributive plots for the first and second halves of the session for the experimental data and separately for all possible combinations of the shuffled data (1000 x 1000 comparisons). We then determined the percentile score of the Pearson correlation coefficient corresponding to the observed data relative to the population of coefficients derived from the shuffled data.

Rat data
Data from rats 8 was downloaded from the Kavli Institute's online database ( https://www.ntnu.edu/kavli/research/grid-cell-data ). The data was available in a format that contained the trajectory of the animal and firing times of sorted cells in MATLAB files. The MATLAB files were converted into a spatial data frame similar to the mouse data so the same analysis scripts could be used to perform all analyses. Our analyses use 32 conjunctive cells and 68 pure grid cells identified in this dataset.

Simulation of previous grid cell models
We simulated four previous grid cell models with experimentally recorded behavioural trajectories as inputs. Code for the Guanella, Giocomo, and Burgess models was acquired from (Zilli 2012) (hosted at ModelDB, accession number: 144006), and the Pastoll model was hosted on ModelDB (accession number: 150031). We extracted spike times as outputs from the simulations and matched these to corresponding head directions and x, y positions from the trajectory. These model results were converted into a format used by our analysis pipeline and underwent the same standardised analysis as experimental data. Code for these models was adapted from the authors and from implementations provided by (Zilli 2012) that were available on ModelDB.

Stellate cell model
All simulations were performed in the NEURON simulation environment 51 . Simulation code will be available at https://github.com/MattNolanLab .

712
The model stellate cell used a previous morphological reconstruction of a mouse MEC layer 2 stellate cell 52 . Voltage-gated sodium and potassium channels were inserted into the soma and axons (channel models from 53 ), HCN channels were inserted into the dendrites and soma (channel models from 54 ) and leak channels were inserted into all compartments. Maximum channel conductances were adjusted so that electrophysiological properties of the neuron were similar to the experimentally determined properties of stellate cells 55 , with the resting membrane potential (RMP), input resistance, sag, rheobase spike peak and half-width fit within the range of experimental values 33 . The best fit for voltage threshold that could be achieved was within 30% of the mean experimentally determined values. RMP was defined as the average membrane potential over 4 s with no current input. Spike peak and half-width were determined from a single suprathreshold response to a 20 ms depolarizing current, where peak potential was the maximum voltage and half-width was the width of the spike at a voltage halfway between RMP and the peak. Sag ratio was determined in response to a 40 pA hyperpolarizing current step and defined as the ratio between the peak decrease and steady state decrease in voltage. Voltage threshold was defined as the highest voltage reached without spiking and was determined from responses to a series of 3 s duration current steps with progressively increasing amplitude. Rheobase was defined as the current required to initiate a spike and was determined using a current ramp increasing linearly from 0 -150 pA over 2 s. Input resistance was determined as the slope of the line fitted to the voltage increase resulting from current injections between 0.0016 nA to 0.048 nA.
The modeled cell received simulated synaptic input from conjunctive cells. Synapses were randomly localized to dendritic locations with 9 synapses per input cell. The probability of synapse placement on a dendrite was given by the ratio of the dendrite length to the total basal or apical dendritic length. Synapses generated fast conductance changes with an instantaneous rise and exponential decay of 2 ms (cf. 52 ). All synapses had the same maximal conductance. We evaluate a 'low conductance' version of the model in which maximal synaptic conductances was 0.603 nS, and a 'high conductance' version in which the maximal synaptic conductances was 250 nS.

Simulation of conjunctive cell firing
To simulate firing of conjunctive cells we first generated for each cell a grid pattern and a head direction tuning curve. For the grid pattern, the centre of each field was specified from the vertices of equilateral triangles with a length of 50 cm that were aligned to tessellate the simulated environment. Each cell had the same spatial phase. Each firing field was described by a circular Gaussian distribution with a full width at half maximum (FWHM) of 20.0 cm. To generate non-uniform field maxima (cf. 24 ), the peak for each field was scaled by a random value from a uniform distribution between 0 and 1. The grids were then scaled to have a peak height of 1. The head direction tuning curve consisted of a Gaussian centered around a preferred direction with a FWHM of 141 degrees and a height of 1. The preferred head directions were evenly distributed between 0 and 360 degrees. For each simulated conjunctive cell, firing probability distributions were calculated from the grid pattern and a head direction tuning curve at a spatial resolution of 1 cm. The probability of firing was determined as 0.28 x the probability 755 of grid x 0.12 x the probability of head direction firing . These values resulted in peak spatial firing rates of around 12 Hz in the low conductance version and 28 Hz in the high conductance version. Peak head direction firing rates were around 5 Hz in the low conductance version, and 10 Hz in the high conductance version, which are comparable to experimental data 8 .
Conjunctive cell spike times were generated for a behavioural trajectory within a 1 x 1 m open field. The x and y position was determined every ms and if the probability of each cell firing was greater than a random number between 0 and 1 the cell spiked. For 100 simulated conjunctive cells, the peak spatial firing rate was 12.3 ± 2.71 Hz and the normalized peak head direction firing rate was 4.94 ± 1.04 Hz.
Each conjunctive cell was connected at 9 synapses to the downstream stellate cell. To evaluate representations generated in the stellate neuron model, we simulated 20 trials for each model configuration. Each trial differed in the randomly determined synapse placement of each conjunctive input and in the randomly determined peaks of the conjunctive cell fields.

Data Availability
All data will be made available from the Nolan Lab repository on the University of Edinburgh's DataShare site ( https://datashare.is.ed.ac.uk/handle/10283/777 ).

Code Availability
All analysis code will be made available from the Nolan Lab GitHub site ( https://github.com/orgs/MattNolanLab/ ).   25 cells for rats)( mice: 4.3 ± 3.2 bins / field; rats: 2.1 ± 2.8 bins / field ). The number of significant bins per field differed from the shuffled data (Mann-Whitney U test comparing the number of significant bins, p < 10 -20 for mice and p < 10 -14 for rats).  of significant directional bins / cell plotted as a function of speed score for grid cells from mice and rats. (c) The number of significant directional bins / field plotted as a function of speed score for grid cells from mice and rats. Note that data from (c) are from the subset of the grid cells in (b) from which fields could be automatically detected. with the corresponding shuffled data (grey) for each pure grid cell. The lighter grey region indicates the range of the 5 th and 95 th percentiles. The maximum firing rate of the observed data is shown above each plot. Significantly directional bars are marked with an asterix (*). Fields classified as 'directional' had at least one bin that differed statistically from the shuffled distribution (p < 0.05 after correcting for multiple comparisons), whereas field classified as 'not directional' did not.  The directional dependence of spiking globally and at the level of single field was determined as in Fig. 2 and Fig. 3 respectively. The plots show the cumulative proportion of bins that differ significantly from their corresponding shuffled distribution globally (upper) and per field (lower) for pure grid cells (left) and conjunctive cells (right). (b) Summary data for proportion of fields with significant bins and average number of significant bins for pure and conjunctive grid cells.      Relationship between field size and field spacing and the proportion of directional bins detected in non-uniform conjunctive input models. Simulations of the non-uniform conjunctive cell input models were run with different field spacing or field sizes. Each point is a cell, because field spacing and field size are metrics calculated per-cell, rather than per field. (a) Field size had a weak but statistically significant effect on the proportion of significant bins detected (Pearson coefficient = -0.12, p = 0.0091). (b) Field spacing had no effect on the proportion of significant bins detected (Pearson coefficient = 0.079, p = 0.34).  Table S1. Comparison of predictions for directional modulation of grid cell firing made by grid models.
References are to example publications and are not exhaustive. Notes for each model: 1. Continuous attractor network models assume that grids result from activity bumps driven around a recurrent network by inputs that encode speed and head direction. Depending on their configuration the models predict either unidirectional head direction modulation (e.g. 59 ) or spatially uniform omnidirectional modulation (e.g. 9 ). Models with unidirectional head direction firing can account for the properties of conjunctive cells. Because spatially uniform directional modulation is a requirement for these models to generate grid fields they are unlikely to explain the local head direction modulation reported here. 2. Oscillatory interference models generate grid fields through summation (or multiplication) of multiple oscillatory inputs that are tuned to a particular direction and have phase modulated by running speed. Because the effect of direction on firing rate is very weak in these models and as uniform directional modulation of input phase is a requirement for these models to generate grid fields they are unlikely to explain the local head direction modulation reported here. 3. In plasticity models grid fields emerge through synaptic plasticity mechanisms in conjunction with adaptation rules (e.g. 57 ) or training signals (e.g. 60 ). To date grid cells generated through these mechanisms appear to have omnidirectional firing fields with no local modulation of directional tuning (cf. 60 ).
4. It is conceivable that grid cell firing patterns are generated through a purely sub-threshold mechanism. In other words, grid patterns are generated by membrane potential changes that alone are insufficient to trigger action potentials. In this scenario, additional input from head direction cells with local spatial firing fields, would convert the silent cell into a grid cell with local-directionally modulated firing fields. A challenge for this scenario is to establish biophysical mechanisms that would ensure that alone neither the grid pattern generator, or the head direction input will drive action potential firing. 5. The scenario we proposed here assumes that co-aligned conjunctive cells make convergent synaptic input onto the common postsynaptic neurons. We show here that these postsynaptic neurons will have grid firing fields that are locally modulated by head direction. This model contrasts with the previous scenario (4) which requires multiple relatively complex and untested assumptions.  Table S2. Estimated position of tips of recording electrodes at the beginning and end of experiments and estimated recording sites in deep and superficial layers of the MEC. Estimated angles are relative to the straightened skull and are based on histology images. Two animals were terminated for health reasons and their brains were not processed. Animals with grid cells used in the analyses are highlighted in bold.

Mouse grid cells with more than two successfully detected firing fields.
Each series of plots summarises data from a single grid cell. Top row left to right: action potential waveforms overlaid for the four channels of the tetrode, autocorrelograms of spike times, histogram of spike counts as a function of recording time, plot of spike rate as a function of movement speed, and coverage heat map for the position of the animal in the arena. Second row left to right: trajectory of the animal (black line) and firing events (red dots), firing rate map, autocorrelation matrix for the rate map, smoothed polar histogram of head direction when the cell fired (red, Hz) and from the whole session (black), scatter plot of firing events colour-coded for head direction on trajectory. Third row left to right: Detected firing fields on rate map and polar histograms of head direction in detected fields.

Grid cells
Animal ID: F