A local measure of symmetry and orientation for individual spikes of grid cells

Grid cells have attracted broad attention because of their highly symmetric hexagonal firing patterns. Recently, research has shifted its focus from the global symmetry of grid cell activity to local distortions both in space and time, such as drifts in orientation, local defects of the hexagonal symmetry, and the decay and reappearance of grid patterns after changes in lighting condition. Here, we introduce a method that allows to visualize and quantify such local distortions, by assigning both a local grid score and a local orientation to each individual spike of a neuronal recording. The score is inspired by a standard measure from crystallography, which has been introduced to quantify local order in crystals. By averaging over spikes recorded within arbitrary regions or time periods, we can quantify local variations in symmetry and orientation of firing patterns in both space and time.


Generating artificial spike locations
To generate artificial grid patterns, we draw spikes from a probability distribution that is comprised of two-dimensional Gaussians whose centers lie on a-possibly perturbedhexagonal grid. The center location of each Gaussian is the location of a grid field. For all generated data, we use an arena size of 1m × 1m. To comply with the experimental observation that the size of grid fields scales with the grid spacing (Brun et al., 2008), we use a constant ratio of the width of the Gaussians and the grid spacing before adding perturbations: grid spacing / standard deviation of Gaussian = 7 .
To assess the robustness of the grid score, we study several perturbations of hexagonal firing patterns. To this end, we first generate grid field locations on a perfect hexagonal grid. To avoid edge effects, this grid is generated on a much larger area that contains a perfect hexagonal arrangement of 271 (i.e., side length 10) grid fields, and the actual arena size (1m × 1m) is cut out in the center after the perturbations are applied.
In the scenario of random perturbations (Fig 3a), we add a different random vector to each grid field location. Each such two dimensional vector is drawn from a two dimensional normal distribution. The standard deviation of this distribution dictates the noise on the field locations (Fig 3c).
In the scenario of shearing, we apply a horizontal shearing transformation to each field location of the perfect lattice, by multiplying the two dimensional location of each firing field with the matrix where s is the strength of the shearing (Fig 4c).
In the scenario of different levels of background noise, we create grid patterns with 2000 spikes each. The locations of the 2000 spikes are drawn from two distributions at different fractions. One fraction of spikes is drawn from the above-mentioned perfect hexagonal arrangement of Gaussians. The other fraction is drawn from a spatially homogeneous distribution over the accessible arena. The larger this fraction, the stronger the background noise.
In the scenario of a drift in orientation along the south-north axis, we rotate the location of the i-th grid field (x i , y i ) T via where (x i , y i ) T denotes the location of the grid field after rotation and φ i increases along the south-north axis in the following way: where y min and y max are the y-locations of the southernmost and northernmost grid fields, respectively. The factor d determines the strength of the drift. We use d = π/3, in the example shown in Fig 7c. Finally, in the scenario of a sudden change in orientation (Fig 7b), we extracted the field locations from simulations by Rosay et al. (2018).

The correlogram-based grid score
To compute the correlogram-based grid score, we use the method suggested by Langston et al. (2010). To this end, we first calculate a "rate map" by summing over Gaussian kernels with a width of 5% of the arena side length, centered at the locations of the spikes.
We then determine the spatial autocorrelogram, i.e., the Pearson correlation coefficients for all spatial shifts of the rate map against itself. From this correlogram, we crop out an annulus that contains the six fields that are arranged around the central peak.
To get the inner radius of this annulus, we clip all values in the correlogram that are smaller than 0.1 to 0. We obtain the resulting clusters of values that are larger than 0.1 using scipy.ndimage.measurements.label from the SciPy package for Python with a quadratic filter structure, ((1, 1, 1), (1, 1, 1), (1, 1, 1)), for a correlogram with 201 × 201 pixels. We use the distance from the center to the outermost pixel of the innermost cluster as the inner radius of the annulus. We obtain the outer radius of the annulus by trying 50 values, linearly increasing from the inner radius to a corner of the arena. With each of the resulting 50 annuli, we crop the correlogram, rotate it and correlate it with an unrotated copy of the correlogram. We determine the correlation coefficient for 30, 60, 90, 120 and 150 degrees. We define the grid score as the minimum of the correlation values at 60 and 120 degrees minus the maximum of the correlation values at 30, 90 and 150 degrees. After trying all 50 annuli, we take the highest resulting grid score as the grid score of the cell. A hexagonal symmetry thus leads to positive values whereas a quadratic symmetry leads to negative values. Since the correlogram-based score uses Pearson correlation values that are typically expressed with the Greek letter ρ, we refer to it as the ρ score. Note that the match of the grid score symbols with the first names of the authors that introduced them (ρ-samund and Ψ-mon) is coincidental.

Extracting trials with specific light conditions
In the experimental protocol of Pérez-Escobar et al. (2016), different lights on different sides of the arena walls are followed by epochs of darkness. The authors measured different grid patterns when different sides of the wall were illuminated (Pérez-Escobar et al., 2016). In our analysis, we were only interested in the effect of darkness on an individual pattern. We thus concatenated the spikes from trials with one lighting condition with subsequent trials of darkness and did not mix spikes from trials with two different lighting conditions.

Box plot settings
In all shown box plots, each box extends from the first to the third quartile, with a dark blue line at the median. The lower whisker reaches from the lowest data point still within 1.5 IQR of the lower quartile, and the upper whisker reaches to the highest data point still within 1.5 IQR of the upper quartile, where IQR is the inter quartile range between the third and first quartile.