Grid pattern development, distortions and topological defects may depend on distributed anchoring

The firing pattern of grid cells in rats has been shown to exhibit elastic distortions that compresses and shears the pattern and suggests that the grid is locally anchored. Anchoring points may need to be learned to account for different environments. We recorded grid cells in animals encountering a novel environment. The grid pattern was not stable but moved between the first few sessions predicted by the animals running behavior. Using a learning continuous attractor network model, we show that learning distributed anchoring points may lead to such grid field movement as well as previously observed shearing and compression distortions. The model further predicted topological defects comprising a pentagonal/heptagonal break in the pattern. Grids recorded in large environments were shown to exhibit such topological defects. Taken together, the final pattern may be a compromise between local network attractor states driven by self-motion signals and distributed anchoring inputs from place cells.


Introduction 19
Cortex may represent variables in discrete or continuous low dimensional attractor states. Although 20 point attractors may be stably learned, for instance by Hopfield networks (Hopfield 1982), continuous 21 attractors networks (CANs) are notably difficult to balance as small perturbations will lead them to 22 wander away from the proper state. Thus, CANs need to monitor and reset their state either by lower 23 level sensory inputs or by learned associations with representations found in higher cortical levels. Crystalline materials may undergo similar elastic distortions to a certain degree, after which increased 33 shear strain is released by slippage of a grain boundary, or nucleation of topological defects such as 34 dislocations. Whether strong local elastic distortions may lead to similar topological defects in the grid 35 pattern is unknown, but they would provide strong evidence of anchoring mechanics. 36 A concern for any system that relies solely on internal locomotor cues to update a position estimate is 37 sensitivity to noise. Incomplete or imperfect sensory input may impart drift to the representation and 38 undermine the possibility to derive stable spatial information from it over time. This caveat could be 39 mitigated by a resetting mechanism that would ensure the systems long-term stability. Thus, to be 40 useful as a spatial reference the grid needs to anchor to its environment. Anchoring has been 41 suggested as a process where stable features of the environment provides a reference frame that will 42 reset the position estimate at these various locations (Campbell et  anchoring in an open field may be most conspicuous during initial encounters with a novel 49 environment where the pattern may depend more strongly on plastic interactions with the 50 environmental layout. The grid has been shown to get increasingly hexagonal over the first few 51 sessions after being subjected to a novel environment (Barry et al. 2012), possibly as a process of 52 increasingly acquired anchoring points. 53 Here we demonstrate how the grid pattern develops in relation to behavior when subjected to a novel 54 environment, how anchoring may explain details of the development as well as previously described 55 distortions of the grid, and how anchoring may affect the topology of the grid pattern. 56

GRID FIELDS MOVE OPPOSITE TO THE RUNNING DIRECTION DURING NOVELTY 58
To investigate the impact of learning on anchoring and the form of the grid we implanted rats with 59 multi-tetrode devices and recorded grid cells in the medial entorhinal cortex. Animals foraged for 60 rewards in a square 1.5 m box, a 1.8 m equilateral triangle or a 2.4 m equilateral triangle for a novel 61 session followed by several subsequent sessions in the same box between 5 h to 3 days later (5 animals 62 in 9 novelty exposures, 40 grid cells, 2-6 grid cells recorded simultaneously in each experiment). Cells 63 that had a gridness score (Langston et al. 2010) above that of the 95 th percentile of a shuffled 64 distribution was chosen, if this were true for at least one of the sessions. Thus, some cells in the 65 analyses did not qualify as a grid cell in the initial session. Only cells with grid spacing smaller than 75 66 cm was chosen to improve the resolution of the local analyses below. 67 Local changes in the grid were studied by subdividing firing rate maps into 9 equal squares for the 68 square environments ( Figure 1A respectively. For each cell, each subdivision from session 1 was cross-correlated with that from session 71 2 ( Figure 1A, blue vs red respectively). Spatial differences in the grid between session 1 and 2 was 72 measured from the offset of the peak closest to the middle of the cross-correlogram in each 73 subdivision, if the peak was higher than 0.4 r and the peak had not moved more than 25% of the grid 74 spacing of the cell. Grid fields were not stable between the 1 st and 2 nd session but had shifted a mean 75 of 5.1 cm (sd = 3.2 cm, Figure 1B). There were no differences in the magnitude of the shift in the 76 corners, along the walls nor in the middle subdivisions of the experiments in the square environment 77 (mean magnitude 5.0 cm, sd = 2.9; 5.1 cm, sd = 3.6; and 5.4 cm, sd 2.7, in the mid, corner and wall 78 compartments respectively, Wilcoxon's rank sum test of the magnitude in the mid vs corner 79 compartments, p = 0.96, mid vs wall, p = 0.35, and corner vs wall compartments, p = 0.21). 80 Across animals, the grid did not shift in any coherent pattern (N = 303 subdivision measurements, 81 mean resultant vector length (MVL) = 0.052, weighted circular mean direction = 21°, Rayleigh's test for 82 circular non-uniformity, p = 0.18). However, when correcting for the behavior within each subdivision 83 by rotating the grid shift direction according to the prevailing running direction in that subdivision from 84 the first session, the shift of the grid was significantly distributed oppositely to the average local 85 running direction ( Figure 1B right shows example and Figure 1C, left, shows the entire data set; MVL 86 = 0.17, weighted circular mean direction = 197°, Rayleigh's test for circular non-uniformity, p = 0.031, 87 weighted circular V-test for non-uniformity around the specific direction 180° from the average 88 running direction (V-test180), v = 179, p = 9.4 × 10 -15 ). Using a subset of the same cells recorded at later 89 3 90 Figure 1. Individual grid fields move oppositely to the average running velocity during novelty. 91 (A) Two example cells recording during novelty (blue) and a subsequent session (red) are shown 92 together with the average running velocity (right).

93
(B) The cross-correlation between sub compartments in session 1 and 2 shows how the grid moves 94 locally (left, black arrows, 6 cells) displayed in relation to the average running velocity (blue arrows).

95
Normalizing the grid movement to the running direction in this example shows that most shifts occur 96 in the opposite direction of the average running direction (right).

97
(C) Pooling the shifts normalized to average running direction shows that the shifts from 40 grid cells 98 from 9 experiments run in the opposite direction of the behavior (left). Black arrows denote field 99 shifts, the red arrow displays the mean vector and the blue running direction. The circular histogram is 100 divided into 10° bins and displays counts weighted by shift magnitude (bottom). The 2 nd to last vs last 101 session grid movement was smaller in magnitude and did not shift in any significant direction (C).

102
(D) The mean magnitude of movement dropped during subsequent exposure to the environment 103 (mean and SEM). 104 sessions demonstrated that the effect was due to novelty (Figure 2C, right; 26 grid cells comparing  105  sessions 4 and 5, 7 and 8, 3 and 4, 5 and 6, 2 and 3, 5 and 6 and finally 5 and 6 for the respective  106 experiments, MVL = 0.089, weighted circular mean direction after rotating according to running 107 direction = 284°, V-test180, v = -13, p = 0.78). Although the direction of the shift in relation to the 108 running direction was significant, the resultant mean vector length was short, suggesting that the grid 109 may move locally also in other directions, possibly to consolidate a less regular initial pattern (Barry et 110 al. 2012). The movement of the grid fields opposite to behavior was significant between session 1 and 111 2 and between session 2 and 3 (MVL = 0.10, V-test180, V = 58, p = 0.0036) after which any movement 112 was statistically undetectable by these means ( Figure 1D, session 3 to 4, MVL 0.07, V-test180, V = -8.3, 113 p = 0.67). 114 Although the shift was undetectable after a few sessions, there is a possibility that movement 115 continues albeit at a slower pace. We therefore analyzed a set of previously recorded data (Stensola et  116 al. 2012) that spanned over weeks. Since cells were only recorded once, we here examine whether the 117 population representation of the shape of the grid changed by measuring the offset of the grid axis 118 that was most closely aligned to a wall of the environment. Out of 39 modules, 15 modules were 119 recorded at more than 10 occasions. The correlation between the absolute offset of the axis and 120 recording day showed that in 6 of these cases (40%), the orientation changed over time. In 4 out of the 121 6 cases, the orientation rotated towards 0° offset (r = -0.44, p = 0.025; r = -0.71, p = 0.0010; r = -0.45, 122 p = 2.3 × 10 -5 ; r = -0.30, p = 0.029) and in the remaining two cases it rotated away from 0° (r = 0.45, p 123 = 0.01; r = 0.27, p = 0.0055). 124 Thus, these data show that the grid is initially volatile and moves between the first and subsequent 125 encounter with an unfamiliar environment and that this movement occurs preferentially in the 126 direction opposite to the animals typical running direction. This movement slows down within a few 127 sessions so that it is undetectable by these means but on a population level the grid may continue 128 moving over weeks. 129

AN ANCHORING CAN MODEL RECAPITULATES FIELD MOVEMENT AND DISTORTIONS 130
The movement of the pattern suggests that the grid may change its points of anchoring. To investigate 131 whether plastic anchoring inputs to a CAN would permit these dynamics we modelled a small CAN of 132 672 grid cells connected to plastic inputs from 289 place cells with place fields in a square lattice 133 covering the entire environment. Place cells closer than 20 cm to a wall had an 20% lower firing rate 134 than the rest to mimic the real distribution of firing rates of hippocampal CA1 place cells (Muessig et  135 al. 2015). Each grid cell received inputs from every place cell through Hebbian synapses modelled by 136 the BCM-rule (Bienenstock, Cooper, and Munro 1982) where the sliding threshold was modulated by 137 the firing rate during the last 300 ms. The patterning on the cortical sheet within the CAN aligned to 138 an edge because of the square architecture of the modeling environment. To allow for arbitrary 139 orientations, an offset variable µ was added to the head direction. Every 100 seconds, the grid 140 orientation was measured from an autocorrelogram of the last 300 seconds of activity and µ was 141 updated according to the prevailing orientation. This rather artificial operation allows the CAN to 142 slowly realign the orientation of the attractor to the pattern expressed during the last few minutes. In 143 the brain, this operation could be performed by plastic interactions between the head direction system 144 and the grid cells. 145 The balance between the collateral connections within the CAN and the inputs from the place cells 146 was set by α so that α = 1 meant that grid cells only received inputs from other grid cells, and at α = 0, 147 they only received place cell inputs. As α was increased, the lower inputs from the place cells made 148 the grid smaller and less stable and at α = 1, the grid pattern was completely abolished, the firing rate 149 was approximately 80% lower and the cells activity strongly reflected the head direction inputs 150 5 (Supplementary Figure 1), similar to experimental data (Bonnevie et al. 2013). The simulated animal 151 moved with anisotropic behavior so that whenever it reached a wall there was 0.02 higher probability 152 per iteration of turning to the left than the right, which also leads to a higher average velocity and 153 higher occupancy along the walls (Figure 2A). This stereotypical behavior of rats is correlated to 154 compression distortions (Hagglund et al. 2019). 155 The model was run 100 times for 500 000 iterations (83 minutes) in a 1.5 m box with 289 place cells, α 156 = 0.92, an initial position in the corner of the box and the initial orientation of the CAN at 0°. We chose 157 a subset of 10 random cells from each simulation to produce an average that we compared between 158 the simulations. After rotating and reflecting the grid pattern from different simulations so that their 159 primary axis was between 0° and 15° the final pattern displayed a median primary axis offset of 5.3°, 160 sd = 3.7 ° with few simulations with offsets close to 0° and 15° ( Figure 2B The development of the pattern was measured by performing a local cross-correlation in each of 9 167 sub-compartments between an initial epoch ( Figure 2C, blue) and a subsequent epoch (red), with 168 epochs being 125000 iterations (21 minutes) long. The grid had moved in directions that were highly 169 distributed overall ( Figure 2C right shows an example simulation, red arrows denote the shift of 170 individual cells, 20 cells plotted), but when adjusting the shifts for the average running velocity in each 171 sub-compartment, the shift of the grid was directed in the opposing direction ( Figure 2D shows the 172 example in C, for the whole population the circular mean direction was 167.5°, MVL = 0.63, V-test180, 173 p < 10 × 10 -20 , mean shift magnitude = 3.3 cm). The deviation from 180° was smaller if half of the 174 simulations were reflected so that the running anisotropy direction was unbiased (circular mean 175 direction 179.5°, MVL 0.35, V-test180, p < 10 × 10 -20 ). Although larger than in the data, the modest MVL 176 suggest that the grid pattern also here moved in other directions. The shifts were most conspicuous in 177 the first epochs and later subsided ( Figure 2E, shows epoch 3 and 4, and Figure 2F shows the average 178 magnitude of the shifts for 7 different simulations with measurements starting every 5 min of the 179 duration of the simulation) but had not disappeared at epoch 3 vs 4. 180 To further elucidate the effects of anchoring on the finer spatial details of the grid we ran the model 181 in a larger 2.2 m environment with a smaller spacing and a lower α to uncover whether anchoring may 182 also reproduce local spatial compression distortions found in such an environment (Hagglund et al. 183 2019). 100 simulations were run for 750 000 iterations (125 min) with 625 place cells and α = 0.85. To 184 relate the grid to both environmental axes we defined a secondary axis of the grid as the line that 185 passes through the 2 nd and 6 th field of an autocorrelogram when counting counter clock wise from 0° 186 and runs orthogonal to the primary axis ( Figure 2A, right). Local grid spacing was measured from the 4 187 out of 9 equal square sub-compartments that abutted each mid-wall section. The average 188 autocorrelogram of 10 randomly chosen cells from each simulation showed that the primary axis was 189 compressed by 8% along the walls that aligned to the primary axis ( Figure 2G; spacing difference 190 between the E-W and the N-S mid wall bin 2.3 cm, Wilcoxon's test, p = 8.2 × 10 -12 ). The spacing along 191 the secondary axis was similarly compressed by 12% along the walls that aligned to the secondary axis 192 ( Figure  animal showed that the grid shift was distributed in the opposite direction from the running direction 206 (example from C).

207
(E) During later epochs the movement had subsided. 208 (F) Measuring the shift every 5 minutes displays how the magnitude of the shift goes down with time 209 (7 simulations shown).

210
(G) The model reproduced local spatial distortions. The spacing of the primary and secondary axes was 211 smaller along the walls that ran alongside those axes (arrows denote the axial direction). The boxplot 212 displays the spacing at the mid wall of each axis (notches show the 95% confidence intervals). 213 7 Disabling the differential place cell firing rates did not remove the compression distortion (primary and 214 secondary axis compressed by 9% and 7% respectively, Wilcoxon's test p = 1.2 × 10 -20 and p = 7.3 × 10 -215 17 respectively) but it did abolish the diagonally symmetric form of the grid (Supplementary Figure 2B). 216 However, as the BCM-rule may introduce intrinsic grid formation on the single cell level independently 217 from the CAN(Stepanyuk 2015) we speculated that the compression distortion may stem from border 218 fields not being constrained by fields lying outside the box (Supplementary Figure 2C). The 219 compression stemming from this effect should thus be reduced by lowering the impact of the plasticity 220 by increasing α (and lowering the speed gain to compensate for the smaller spacing stemming from 221 the smaller anchoring input), which is exactly what we saw. With α = 0.96 (as high as possible without 222 losing stability of the grid) and vgain = 0.018, primary and secondary axes were less compressed (both 223 by 4%; smaller than simulations run at α = 0.85, Wilcoxon's test p = 0.0030 and p = 4.5 × 10 -5 for the 224 primary and secondary axes respectively). Shearing depended on both the anisotropic behavior and 225 the compression distortion and was abolished if either one was not present (Supplementary Figure 2D; 226 without anisotropic behavior, median primary axis offset increased from 3.5°, sd = 3.4° to 5.5° sd = 5. approach. The de-elliptified gridness score is a metric that fits an ellipse on the 6 innermost peaks of 241 the autocorrelogram (excluding the center peak) and then transforms the autocorrelogram so that the 242 peaks lie on a circle. The gridness score (Langston et al. 2010) is subsequently measured throughout 243 the environment to provide an estimate for hexagonality that is not biased by ellipticity. Henceforth 244 we will only use this de-elliptified gridness score. 245 Many maps exhibited a homogenously distributed high gridness score throughout the environment 246 (Supplementary Figure 3, top) In some simulations there was a local region where the score was either 247 below zero or undefined because the 6 innermost peaks of the autocorrelogram did not form an ellipse 248 ( Figure 3A). In the same region as these low gridness scores, the local spacing and orientation of the 249 grid underwent a spatial phase transition around a nexus, where one or two of the axes of the grid 250 suddenly jumped to a different configuration ( Figure 3B and 3C). This jump signified that the two peaks 251 (autocorrelograms are 180° rotationally symmetric so peaks come in pairs) closest to the center of the 252 autocorrelogram became farther away than another pair of peaks. Thus, we concluded that in the 253 location of the anomaly, one row of fields locally exhibited a larger spacing than a neighboring row, 254 but farther from the anomaly the pattern was again coherent. Taken together, these transitions are 255 consistent with a phenomenon known as a dislocation. 256 Dislocations have been studied in materials science for over 80 years as they have a large impact on 257 the properties of metals and other crystalline materials (Burgers 1939). They are defined by the 258 8 259 Figure 3. The model produces grids with topological defects. 260 (A) Some simulations contained local dips in local gridness score (blue) and areas where the gridness 261 score could not be defined (white pixels).

262
(B and C) The regions of low gridness score were accompanied by a local phase transition in spacing 263 and orientation respectively showing that these grids had a dislocation defect.

264
(D) A Burgers circuit can be drawn that encapsulates the defect and shows the discrepancy in 265 positional order.

266
(E) Voronoi tesselation analyses was used to identify pentagonal and heptagonal fields.

267
(F) There was a local pentagon/heptagon hot-spot in the same area as the low gridness score (same 268 simulation as in A, B, C and D). Transparent red and blue fields denote those fields that had one 269 Voronoi edge closer than 20 cm to a wall.

270
(G) The gridness scores at the areas with pentagons, hexagons and heptagons shows that non-271 hexagonal Voronoi cells display a lower gridness score (outliers not shown). heptagons in the pattern using Voronoi tessellation by Delaunay triangulation, with grid field locations 292 as seeds ( Figure 3E). All maps had high deviation from hexagonality along the borders ( Figure 3G), 293 which is to be expected from this method, as there are no fields outside the borders to constrain the 294 Delaunay triangulation. Therefor we disregarded any Voronoi cells that had an edge closer than 20 cm 295 to the border ( Figure 3F, transparent colored fields). The number of bins with a zero or lower gridness 296 score, or an undefined gridness score, correlated with the mean number of pentagons (r = 0.57, p = 297 6.1 × 10 -10 ) and heptagons (r = 0.52, p = 2.1 × 10 -8 ) in the population. Moreover, the location of the 298 pentagonal/heptagonal fields coincided with the areas of lower gridness scores ( Figure 3G; median 299 local gridness scores where there were pentagonal fields, 0.99; hexagonal fields, 1.33; heptagonal 300 fields 1.27; 1000 cells from 100 simulations, ANOVA, F(4, 11955) = 1114, p = 2.9 × 10 -182 , Wilcoxon's 301 test of local gridness scores at pentagons and heptagons vs hexagons, p = 6.8 × 10 -99 ). 302 We next asked whether the local phase structure of the dislocated simulations displayed topological 303 breaks in the translational symmetry between grid cell pairs, such as a reflection of the phase-to-phase 304 coupling between cells, or a disclination defect where the rotational order is altered, in contrast to the 305 positional order which is changed in dislocations. To visualize the variations in translational phase 306 coupling, a sliding window cross-correlation was performed between 50 cell pairs that had an average 307 offset distance of 25% +/-10% of the grid spacing (see examples in Figure 3H). The offset between cell 308 pairs were normalized by the mean so that a map could be constructed displaying the local average 309 directional and magnitude deviation as well as cross-correlation amplitude between cells in the 310 population ( Figure 3H, 3 right-most panels) showing that close to areas with low gridness score the 311 map often underwent a local change in the phase structure, but it was not much different from the 312 variation seen in non-dislocated simulations. In each map with a gridness score of zero or below, or 313 where the gridness score was undefined, we measured the amplitude and deviation of magnitude and 314 direction at the minimum and maximum level of gridness score ( Figure 3I). None of these variables 315 were different across the population (difference in deviation in direction 0.3°, magnitude 4.7 cm, 316 amplitude r = 0.0043, Wilcoxon's test p = 0.19, p = 0.55, p = 0.07 respectively). To further quantify the 317 phase structure, we subdivided rate maps into 9 sub-compartments and measured the phase offset 318 between 50 cells if their average distance were 25% +/-10%. The results showed that in locations 319 where the gridness score was zero or smaller or undefined the offset magnitude and the offset 320 direction displayed mean deviations of 0.4 cm and 0.4° respectively, compared to where the gridness 321 score was higher than 0 (Wilcoxon's test, p = 0.35 and p = 0.018, respectively). Although significant, 322 the phase offset deviation within dislocated regions were small and the overall translational symmetry 323 between cells was thus mostly contained, which is not surprising since we impose the phase structure 324 with the synaptic connectivity architecture of the CAN. 325

DISLOCATIONS ARE FOUND IN REAL GRID CELLS 326
To investigate whether the dislocations described above could occur in real grid cells we analyzed a 327 data set of 131 grid cells from 5 modules recorded from three animals foraging in a 2. gridness scores were measured locally throughout the environment using a sliding window with a size 332 of 73 cm ( Figure 4A). Three of the modules (here referred to as module 1-3) described a consistently 333 high gridness score throughout the environment. However, two modules (module 4 and 5) contained 334 a region of negative or undefined gridness score denoting low hexagonality. At these locations, local 335 spacing and orientation was found to exhibit a similar transition as that found in the modelled grid 336 ( Figure 4C and D). 337 To investigate the hexagonality and deviations from it directly, fields were defined as regional maxima 338 in locally normalized oversmoothed rate maps ( Figure 4B, bottom) with a smoothing kernel 339 proportional to the grid spacing (see Methods). The Voronoi method was used as described above, 340 disregarding any Voronoi cells that had an edge closer than 20 cm to the border. Non-hexagonal 341 Voronoi cells were found scattered throughout the environment in each module, likely false positives 342 stemming from the noise of low sampling, individual differences in firing rate and disproportionate 343 coverage. However, in module 4 and 5 there was a local hot-spot of non-hexagonal Voronoi cells 344 (example in Figure 4E). For each rate map the numbers of pentagonal and hexagonal Voronoi cells 345 were counted. Module 4 and 5 contained a higher ratio of grid cells with at least one heptagonal or 346 pentagonal Voronoi cell ( Figure 4F left; chi-square test, chi = 18.9, p = 1.4 × 10 -5 ), higher ratio of 347 heptagons ( Figure 4F mid; chi-square test, chi = 11.5, p = 6.7 × 10 -4 ) and a higher ratio of pentagons 348 ( Figure 4F right; chi-square test, chi = 6.4, p = 0.011). The local gridness score was lower at pentagonal 349 and heptagonal fields compared to hexagonal fields ( Figure  Wilcoxon's test, p = 2.0 × 10 -5 ). Thus, within each rate map, the grid in the quadrant containing a 369 dislocation had a different shape compared to the three other quadrants. 370 As in the model, we asked whether the local phase structure of the dislocated modules displayed any 371 topological breaks in the translational symmetry of the grid, such as a reflection of the phase-to-phase 372 coupling between cells or a disclination defect. Local phase deviation maps showed that similar to the 373 model, the offset displayed a local deviation at the region of the dislocations which was not abrupt but 374 continuous (see example in Figure  To quantify the local translational symmetry, each quadrant of the rate maps was cross-correlated 399 between cells within each of the dislocated modules. The distance and direction offset from the middle 400 to the closest peak in the cross-correlogram was measured. Local population deviation of direction and 401 distance between cells were evaluated by subtracting the mean magnitude and direction of the offset 402 from each measurement ( Figure 4H, bottom right) that the variability in phase relationships between cells is small in a dislocated region. The translational 413 symmetry between cells is thus not topologically disturbed but may stay roughly consistent also across 414 dislocations. 415

Discussion 416
We have shown that the layout of the grid pattern is experience dependent and may be in a state of 417 change depending on the behavior during initial encounters of an environment as well as over longer 418 time spans. The novelty induced re-anchoring could be reproduced by a CAN that learns to anchor to 419 place cell inputs. Moreover, the learning CAN explained previously described spatial distortions and 420 predicted topological defects with retained translational symmetry between cells. These topological 421 defects were found in a data set of animals running in a large environment. 422 Grid cells have been shown to fire in anticipation of grid fields (Kropff et al. 2015;De Almeida et al. 423 2012). This effect may lead to associative re-anchoring during times of high plasticity that likely occurs 424 during an encounter with a novel environment. In our model, upon entering a grid field, a grid cells 425 firing rate increases which leads to potentiation of place cell synapses from place cells with place fields 426 at this grid field entry. However, potentiation according to adapting Hebbian learning, such as with the 427 BCM-rule used here, saturates with time, which means that when the animal exits the grid field place 428 cell synapses from place cells with fields at the exit point will be less potentiated or may even be 429 depressed. Thus, grid cells slowly shifted their representation by changing which place cells they 430 anchored to. Within session movement was not analyzed because of a lack of data. Although possible, 431 we find it less likely that the field shift occurred solely during the inter-session rest period and more 432 likely to be a gradual process unfolding during a session. lack of traversal of fields lying outside the environment which leads to a lesser force separating the 457 fields. The lack of constraining neighboring fields along the secondary axis would impose a stronger 458 force than along the primary axis which is why the secondary axis is relatively more compressed than 459 the primary axis. This effect may further explain why the pattern is not only rotated by the anisotropic 460 behavior but shears because of unequal compression between the primary and the secondary axis. although not addressed directly in the current study, may be biased by several factors. Large 472 environments, small grid spacing, and strong anchoring could all lead to a decrease in long range 473 interactions of the pattern and subsequent breaks in the topology of the grid. Possibly they are more 474 prone to form from two or more disparate anchoring points between which the pattern could not 475 coherently consolidate, which may predispose them in geometrically altered environments 476 (Supplementary Figure 4). If and how dislocations may generalize to other attractor systems of the 477 brain is an intriguing question for future studies. determined by repeated shuffling of the experimental data as previously described (Langston 2010). If 580 the gridness score was larger than that of the 95th percentile of the grid scores from the shuffled data, 581 the cell was defined as a grid cell. 582 Rate maps of well separated neural signals were produced as described earlier (Sargolini et al. 2006). 583 Position estimates were smoothed with a 15-point median filter. The position data were sorted into 584 1.5 cm bins and a smoothing algorithm using a Gaussian kernel with a sigma of 6 cm was applied to 585 the spike position data, which was normalized by the spatial occupancy. Locations further than 3 cm 586 away from the animal's path were considered unvisited. The local shift of the grid from one session to 587 the next was measured by subdividing each map into 9 equal bins and measuring the offset by cross-588 correlation. Average running velocity was also collected from each sub-compartment and the grid shift 589 was normalized to the average running direction by subtraction. 590 Parameters of the grid were measured in autocorrelograms where the location of the six innermost 591 local maxima, excluding the middle peak, were extracted. The mean orientation of three axes of the 592 grid in each module was measured and a primary axis was determined as the axis that was closest to 593 either the N-S or E-W axis of the recording enclosure. All cells within a module/simulation were rotated 594 and/or reflected consistently depending on the mean orientation of the primary axis of each 595 module/simulation so that it ended up between 0° and 15° offset from the E-W axis. 596 To evaluate the ellipticity of the grid, an ellipse was fitted to the six innermost peaks of 597 autocorrelograms (excluding the middle peak) using a least-square criterion and the orientation of the 598 long axis of the ellipse was defined within the interval 0° to 180° (Fitzgibbon 1999). Ellipticity was 599 defined as the length of the semi-major axis divided by the semi-minor axis. 600 In the shearing analyses, the position of the six innermost peaks was incrementally sheared by step-601 wise increasing the shearing factor q in 602 ( ) = ( + ) 603 along the y axis and the value of the shearing factor q and primary axis offset was measured where 604 shearing produced the lowest the ellipticity. 605 To measure the de-elliptified gridness score, autocorrelograms were rotated to align their elliptic 606 orientation with the environmental orientation after which they were compressed along this axis to 607 retain a circular layout of the innermost 6 peaks. After this transformation the gridness score was 608 measured (Langston 2010). 609 Local grid parameter maps were constructed by sliding a 74 × 74 cm window across 10 randomly 610 chosen cells from each simulation, calculating an autocorrelogram at every 7.5 cm. Autocorrelograms 611 were averaged between cells to get a mean representation of the grid at each location of the 612 environment. Local spacing ( Figure 2G, 3B, 4C), orientation ( Figure 2B, 4D, Supplementary Figure 2) 613 and gridness scores (Figure 3 and 4) was measured in each simulation and was then averaged to 614 produce the final maps. However, in the data analysis all cells from respective modules were used and 615 the step-size of the sliding window was 1.5 cm. 616 Local phase deviation maps were constructed by a sliding window cross-correlation in steps of 4.5 cm 617 or 1.5 cm per measurement in the model and data respectively. In the model, 50 randomly chosen cell 618 pair combinations from each simulation was cross-correlated locally in a 74 × 74 bin window if their 619 grid offset was between 15% -35%. At each location the peak closest to the middle was interpreted as 620 the translational phase offset between the cells. The mean offset was measured from the full map 621 cross-correlation which was subsequently subtracted from the local maps to produce maps describing 622 the local deviation from the mean. Magnitude, direction and amplitude of the cross-correlation offset 623 peak was then extracted and plotted as individual parameters in color coded maps. In the data the 624 process was the same except for the number of cell pairs used (293), the size of the window (30 × 30) 625 and the step size = 1.5 cm. 626

GRID POLYGONS 627
Fields were defined as local maxima in rate maps produced by two filtering steps. Rate maps were first 628 locally normalized by being divided by a gaussian over-smoothed version of themselves with a σ = 15 629 cm. Next, in the model data the resulting map was further smoothed with a σ = 9 cm. In the real data 630 the σ for the second smoothing was 631 σ = ̅̅̅̅̅̅̅̅ 2 / 200 632 derived from the average spacing ̅̅̅̅̅̅̅̅ of the entire module. Field locations were extracted as local 633 maxima and used as seeds for Delaunay triangulation. Grid field polygonal structure was analyzed by 634 Voronoi diagrams derived from the Delaunay triangulation. 635 MODEL 636 The model consists of a periodic CAN with 24 x 28 grid cells situated on a twisted torus (Guanella,Kiper,637 and Verschure 2007), connected to 289 or 625 gaussian place cells in the 1.5 m and 2.2 m boxes 638 respectively, distributed evenly every 9 cm with a sigma of 15cm. Each grid cell integrates grid cell 639 inputs, place cell inputs and velocity input, and after this the synaptic weights from the place cells are 640 updated according to the BCM rule with a temporal sliding threshold β based on the average firing rate 641 for each neuron for the last 300 ms. The balance between the place cell and grid cell inputs to a grid 642 cell is set by α. 643 The firing rate g of neuron i at time t is 644     C. Fields in the middle of the environment (red) are constrained from all directions. Fields along the borders lack neighboring elds outside the box (dashed circles) and are thus less constrained (green arrows denote lacking forces) meaning the axis along the wall may compress (blue arrows). The lack of constraint has a stronger e ect along the secondary axis since the lacking border eld (green dashed circle) would have had a larger constraining in uence on the spacing. From 13 animals, three grid modules had more than 3 cells with a grid spacing less than 60 cm. Local gridness score maps (top) revealed that two of the modules contained areas of low delity, and in one of the (Module 7) seemed to contain 2 such areas. The ratio of pentagonal or heptagonal Voronoi cells is signi cantly higher in the modules that had a zero or below or unde ned gridness score (chi = 22.4, p = 2.2e-6).