1. Neuroscience
Download icon

Extracting grid cell characteristics from place cell inputs using non-negative principal component analysis

  1. Yedidyah Dordek
  2. Daniel Soudry Is a corresponding author
  3. Ron Meir
  4. Dori Derdikman Is a corresponding author
  1. Technion – Israel Institute of Technology, Israel
  2. Columbia University, United States
Research Article
Cited
3
Views
2,269
Comments
0
Cite as: eLife 2016;5:e10094 doi: 10.7554/eLife.10094

Abstract

Many recent models study the downstream projection from grid cells to place cells, while recent data have pointed out the importance of the feedback projection. We thus asked how grid cells are affected by the nature of the input from the place cells. We propose a single-layer neural network with feedforward weights connecting place-like input cells to grid cell outputs. Place-to-grid weights are learned via a generalized Hebbian rule. The architecture of this network highly resembles neural networks used to perform Principal Component Analysis (PCA). Both numerical results and analytic considerations indicate that if the components of the feedforward neural network are non-negative, the output converges to a hexagonal lattice. Without the non-negativity constraint, the output converges to a square lattice. Consistent with experiments, grid spacing ratio between the first two consecutive modules is −1.4. Our results express a possible linkage between place cell to grid cell interactions and PCA.

https://doi.org/10.7554/eLife.10094.001

eLife digest

Long before the invention of GPS systems, ships used a technique called dead reckoning to navigate at sea. By tracking the ship’s speed and direction of movement away from a starting point, the crew could estimate their position at any given time. Many believe that some animals, including rats and humans, can use a similar process to navigate in the absence of external landmarks. This process is referred to as “path integration”.

It is commonly believed that the brain’s navigation system is based on such path integration in two key regions: the entorhinal cortex and the hippocampus. Most models of navigation assume that a network of grid cells in the entorhinal cortex processes information about an animal’s speed and direction of movement. The grid cell network estimates the animal’s future position and relays this information to cells in the hippocampus called place cells. Individual place cells then fire whenever the animal reaches a specific location.

However, recent work has shown that information also flows from place cells back to grid cells. Further experiments have suggested that place cells develop before grid cells. Also, inactivating place cells eliminates the hexagonal patterns that normally appear in the activity of the grid cells.

Using a computational model, Dordek, Soudry et al. now show that place cell activity could in principle trigger the formation of the grid cell network, rather than vice versa. This is achieved using a process that resembles a common statistical algorithm called principal component analysis (PCA). However, this only works if place cells only excite grid cells and never inhibit their activity, similar to what is known from the anatomy of these brain regions. Under these circumstances, the model shows hexagonal patterns emerging in the activity of the grid cells, with similar properties to those patterns observed experimentally.

These results suggest that navigation may not depend solely on grid cells processing information about speed and direction of movement, as assumed by path integration models. Instead grid cells may rely on position-based input from place cells. The next step is to create a single model that combines the flow of information from place cells to grid cells and vice versa.

https://doi.org/10.7554/eLife.10094.002

Introduction

The system of spatial navigation in the brain has recently received much attention (Burgess, 2014; Morris, 2015; Eichenbaum, 2015). This system involves many regions, which seem to divide into two major classes: regions such as CA1 and CA3 of the hippocampus, which contain place cells (O'Keefe and Dostrovsky, 1971; O'Keefe and Nadel, 1978), vs. regions, such as the medial-entorhinal cortex (MEC), the presubiculum and the parasubiculum, which contain grid cells, head-direction cells and border cells (Hafting et al., 2005; Boccara et al., 2010; Sargolini et al., 2006; Solstad et al., 2008; Savelli et al., 2008). While the phenomenology of those cells is described in many studies (Derdikman and Knierim, 2014; Tocker et al., 2015), the manner in which grid cells are formed is quite enigmatic. Many mechanisms have been proposed. The details of these mechanisms differ, however, they mostly share in common the assumption that the animal’s velocity is the main input to the system (Derdikman and Knierim, 2014; Zilli, 2012; Giocomo et al., 2011), such that positional information is generated by the integration of this input in time. This process is termed 'path integration' (PI) (Mittelstaedt and Mittelstaedt, 1980). A notable exception to this class of models was suggested in a previous paper by Kropff and Treves (2008); and in a sequel to that paper (Si and Treves, 2013), in which they demonstrated the emergence of grid cells from place cell inputs without using the rat's velocity as an input signal.

We note here that generating grid cells from place cells may seem at odds with the architecture of the network, since it is known that place cells reside at least one synapse downstream of grid cells (Witter and Amaral, 2004). Nonetheless, there is current evidence that the feedback from place cells to grid cells is of great functional importance. Specifically, there is evidence that inactivation of place cells causes grid cells to disappear (Bonnevie et al., 2013), and furthermore, it seems that, in development, place cells emerge before grid cells do (Langston et al., 2010; Wills et al., 2010). Thus, there is good motivation for trying to understand how the feedback from hippocampal place cells may contribute to grid cell formation.

In the present paper, we thus investigated a model of grid cell development from place cell inputs. We showed the resemblance between a feedforward network from place cells to grid cells to a neural network architecture previously used to implement the PCA algorithm (Oja, 1982). We demonstrated, both analytically and through simulations, that the formation of grid cells from place cells using such a neural network could occur given specific assumptions on the input (i.e. zero mean) and on the nature of the feedforward connections (specifically, non-negative, or excitatory).

Results

Comparing neural-network results to PCA

We initially considered the output of a single-layer neural network and of the PCA algorithm in response to the same inputs. These consisted of the temporal activity of a simulated agent moving around in a two-dimensional (2D) space (Figure 1A; see Materials and methods for details). In order to mimic place cell activity, the simulated virtual space was covered by multiple 2D Gaussian functions uniformly distributed at random (Figure 1B), which constituted the input. In order to calculate the principal components, we used a [Neuron x Time] matrix (Figure 1C) after subtracting the temporal mean, generated from the trajectory of the agent as it moved through the place fields. Thus, we displayed a one-dimensional mapping of the two-dimensional activity, transforming the 2D activity into a 1D vector per input neuron. This resulted in the [Neuron X Neuron] covariance matrix (Figure 1D), on which PCA was performed by evaluating the appropriate eigenvalues and eigenvectors.

Construction of the correlation matrix from behavior.

(A) Diagram of the environment. Black dots indicate places the virtual agent has visited. (B) Centers of place cells uniformly distributed in the environment. (C) The [Neuron X Time] matrix of the input-place cells. (D) Correlation matrix of (C) used for the PCA process.

https://doi.org/10.7554/eLife.10094.003

To learn the grid cells, based on the place cell inputs, we implemented a single-layer neural network with a single output (Figure 2). Input to output weights were governed by a Hebbian-like learning rule. As described in the Introduction (see also analytical treatment in the Methods section), this type of architecture induces the output’s weights to converge to the leading principal component of the input data.

Neural network architecture with feedforward connectivity.

The input layer corresponds to place cells and the output to a single cell.

https://doi.org/10.7554/eLife.10094.004

The agent explored the environment for a sufficiently long time allowing the weights to converge to the first principal component of the temporal input data. In order to establish a spatial interpretation of the eigenvectors (from PCA) or the weights (from the converged network) we projected both the PCA eigenvectors and the network weights onto the place cells space, producing corresponding spatial activity maps. The leading eigenvectors of the PCA and the network’s weights converged to square-like periodic spatial solutions (Figure 3A–B).

Results of PCA and of the networks' output (in different simulations).

(A) 1st 16 PCA eigenvectors projected on the place cells' input space. (B) Converged weights of the network (each result from different simulation, initial conditions and trajectory) projected onto place cells' space. Note that the 8 outputs shown here resemble linear combinations of components #1 to #4 in panel A.

https://doi.org/10.7554/eLife.10094.005

Being a PCA algorithm, the spatial projections of the weights were periodic in space due to the covariance matrix of the input having a Toeplitz structure (Dai et al., 2009) (a Toeplitz matrix has constant elements along each diagonal). Intuitively, the Toeplitz structure arises due to the spatial stationarity of the input. In fact, since we used periodic boundary conditions for the agent’s motion, the covariance matrix was a circulant matrix, and the eigenvectors were sinusoidal functions, with length constants determined by the scale of the box (Gray, 2006) [a circulant matrix is defined by a single row (or column), and the remaining rows (or columns) are obtained by cyclic permutations. It is a special case of a Toeplitz matrix - see for example Figure 1D]. The covariance matrix was heavily degenerate, with approximately 90% of the variance accounted for by the first 15% of the eigenvectors (Figure 4B). The solution demonstrated a fourfold redundancy. This was apparent in the plotted eigenvalues (from the largest to the smallest eigenvalue, Figure 4A and C), which demonstrated a fourfold grouping-pattern. The fourfold redundancy can be explained analytically by the symmetries of the system – see analytical treatment of PCA in Methods section (specifically Figure 15C).

Eigenvalues and eigenvectors of the input's correlation matrix.

(A) Eigenvalue size (normalized by the largest, from large to small (B) Cumulative explained variance by the eigenvalues, with 90% of variance accounted for by the first 35 eigenvectors (out of 625). (C) Amplitude of leading 32 eigenvalues, demonstrating that they cluster in groups of 4 or 8. Specifically, the first four clustered groups correspond respectively (from high to low) to groups A,B,C & D In Figure 15C, which have the same redundancy (4,8,4 & 4).

https://doi.org/10.7554/eLife.10094.006

In summary, both the direct PCA algorithm and the neural network solutions developed periodic structure. However, this periodic structure was not hexagonal but rather had a square-like form.

Adding a non-negativity constraint to the PCA

It is known that most synapses from the hippocampus to the MEC are excitatory (Witter and Amaral, 2004). We thus investigated how a non-negativity constraint, applied to the projections from place cells to grid cells, affected our simulations. As demonstrated in the analytical treatment in the Methods section, we could expect to find hexagons when imposing the non-negativity constraint. Indeed, when adding this constraint, the outputs behaved in a different manner and converged to a hexagonal grid, similar to real grid cells. While it was straightforward to constrain the neural network, calculating non-negative PCA directly was a more complicated task due to the non-convex nature of the problem (Montanari and Richard, 2014; Kushner and Clark, 1978).

In the network domain, we used a simple rectification rule for the learned feedforward weights, which constrained their values to be non-negative. For the direct non-negative PCA calculation, we used the raw place cells activity (after spatial or temporal mean normalization), as inputs to three different iterative numerical methods: NSPCA (Nonnegative Sparse PCA), AMP (Approximate Message Passing) and FISTA (Fast Iterative Threshold and Shrinkage) based algorithms (see Materials and methods section).

In both cases, we found that hexagonal grid cells emerged in the output layer (plotted as spatial projection of weights and eigenvectors: Figure 5A–B, Figure 6A–B, Video 1, Video 2). When we repeated the process over many simulations (i.e. new trajectories and random initializations of weights) we found that the population as a whole consistently converged to hexagonal grid-like responses, while similar simulations with the unconstrained version did not (compare Figure 3 to Figure 5Figure 6).

Output of the neural network when weights are constrained to be non-negative.

(A) Converged weights (from different simulations) of the network projected onto place cells space. See an example of a simulation in Video 1. (B) Spatial autocorrelations of (A). See an example of the evolution of autorcorrelation in simulation in Video 2.

https://doi.org/10.7554/eLife.10094.007
Video 1
Evolution in time of the network's weights.

625 Place-cells used as input. Video frame shown every 3000 time steps up to t=1,000,000. Video converges to results similar to those of Figure 5.

https://doi.org/10.7554/eLife.10094.008
Video 2
Evolution of autocorrelation pattern of network's weights shown in Video 1.
https://doi.org/10.7554/eLife.10094.009
Results from the non-negative PCA algorithm.

(A) Spatial projection of the leading eigenvector on input space. (B) Corresponding spatial autocorrelations. The different solutions are outcomes of multiple simulations with identical settings in a new environment and new random initial conditions.

https://doi.org/10.7554/eLife.10094.010

In order to further assess the hexagonal grid emerging in the output, we calculated the mean (hexagonal) Gridness scores ([Sargolini et al., 2006], which measure the degree to which the solution resembles a hexagonal grid [see Materials and methods]). We ran about 1500 simulations of the network (in each simulation, the network consisted of 625 place cell-like inputs and a single grid cell-like output), and found noticeable differences between the constrained and unconstrained cases. Namely, the Gridness score in the non-negatively constrained-weight simulations was significantly higher than in the unconstrained-weight case (Gridness = 1.07 ± 0.003 in the constrained case vs. 0.302 ± 0.003 in the unconstrained case. see Figure 7). A similar difference was observed with the direct non-negative PCA methods (1500 simulations, each with different trajectories, Gridness = 1.13 ± 0.0022 in the constrained case vs. 0.27 ± 0.0023 in the unconstrained case).

Histograms of Gridness values from network and PCA.

First row (A) + (C) corresponds to network results, and second row (B) + (D) to PCA. The left column histograms contain the 60° Gridness scores and the right one the 90° Gridness scores.

https://doi.org/10.7554/eLife.10094.011

Another score we tested was a 'Square Gridness' score (see Materials and methods) where we measured the 'Squareness' of the solutions (as opposed to 'Hexagonality'). We found that the unconstrained network had a higher square-Gridness score while the constrained network had a lower square-Gridness score (Figure 7); for both the direct-PCA calculation (square-Gridness = 0.89 ± 0.0074 in the unconstrained case vs. 0.1 ± 0.006 in the constrained case) and the neural-network (square-Gridness = 0.073 ± 0.006 in the constrained case vs. 0.73 ± 0.008 in the unconstrained case).

All in all, these results suggest that when direct PCA eigenvectors and neural network weights were unconstrained they converged to periodic square solutions. However, when constrained to be non-negative, the direct PCA, and the corresponding neural network weights, both converged to a hexagonal solution.

Dependence of the result on the structure of the input

We investigated the effect of different inputs on the emergence of the grid structure in the networks' output. We found that some manipulation of the input was necessary in orderto enable the implementation of PCA in the neural network. Specifically, PCA requires a zero-mean input, while simple Gaussian-like place cells do not possess this property. In order to obtain input with zero-mean, we either performed differentiation of the place cells’ activity in time, or used a Mexican-hat like (Laplacian) shape (See Materials and methods for more details on the different types of inputs). Another option we explored was the usage of positive-negative disks with a total sum of zero activity in space (Figure 8). The motivation for the use of Mexican-hat like transformations is their abundance in the nervous system (Wiesel and Hubel, 1963; Enroth-Cugell and Robson, 1966; Derdikman et al., 2003).

Different types of spatial input used in our network.

(A) 2D Gaussian function, acting as a simple place cell. (B) Laplacian function or Mexican hat. (C) A positive (inner circle) - negative (outer ring) disk. While inputs as in panel A do not converge to hexagonal grids, inputs as in panels B or C do converge.

https://doi.org/10.7554/eLife.10094.012

We found that usage of simple 2-D Gaussian-functions as inputs did not generate hexagonal grid cells as outputs (Figure 9). On the other hand, time-differentiated inputs, positive-negative disks or Laplacian inputs did generate grid-like output cells, both when running the non-negative PCA directly (Figure 6), or by simulating the non-negatively constrained Neural Network (Figure 5). Another approach we used for obtaining zero-mean was to subtract the mean dynamically from every output individually (see Materials and methods). The latter approach, related to adaptation of the firing rate, was adopted from Kropff & Treves (Kropff and Treves, 2008), who used it to control various aspects of the grid cell's activity. In addition to controlling the firing rate of the grid cells, if applied correctly, the adaptation could be exploited to keep the output's activity stable, with zero-mean rates. We applied this method in our system and in this case the outputs converged to hexagonal grid cells as well, similarly to the previous cases (e.g. derivative in time, or Mexican hats as inputs; data not shown).

Spatial projection of outputs’ weights in the neural network when inputs did not have zero mean (such as in Figure 8A).

(A) Various weights plotted spatially as projection onto place cells space. (B) Autocorrelation of (A).

https://doi.org/10.7554/eLife.10094.013

In summary, two conditions were required for the neural network to converge to spatial solutions resembling hexagonal grid cells: (1) non-negativity of the feedforward weights and (2) an effective zero-mean of the inputs (in time or space).

Stability analysis

Convergence to hexagons from various initial spatial conditions

In order to numerically test the stability of the hexagonal solution, we initialized the network in different ways, randomly, using linear stripes, squares, rhomboids (squares on hexagonal lattice) and noisy hexagons. In all cases, the network converged to a hexagonal pattern (Figure 10; for squares and stripes, other shapes not shown here).

Evolution in time of the networks’ solutions (upper rows) and their autocorrelations (lower rows).

The network was initialized in shapes of (A) Squares and of (B) stripes (linear).

https://doi.org/10.7554/eLife.10094.014

We also ran the converged weights in a new simulation with novel trajectories and tested the Gridness scores, and the inter-trial stability in comparison to previous simulations. We found that the hexagonal solutions of the network remained stable although the trajectories varied drastically (data not shown).

Asymptotic stability of the equilibria

Under certain conditions (e.g., decaying learning rates and independent and identically distributed (i.i.d.) inputs), it was previously proved (Hornik and Kuan, 1992), using techniques from the theory of stochastic approximation, that the system described here can be asymptotically analyzed in terms of (deterministic) Ordinary Differential Equations (ODE), rather than in terms of the stochastic recurrence equations. Since the ODE defining the converged weights is non-linear, we solved the ODEs numerically (see Materials and methods), by randomly initializing the weight vector. The asymptotic equilibria were reached much faster, compared to the outcome of the recurrence equations. Similarly to the recurrence equations, constraining the weights to be non-negative induced them to converge into a hexagonal shape while a non-constrained system produced square-like outcomes (Figure 11).

Numerical convergence of the ODE to hexagonal results when weights are constrained.

(A) + (B): 60° and 90° Gridness score histograms. Each score represents a different weight vector of the solution J. (C) + (D): Spatial results for constrained and unconstrained scenarios, respectively. (E) + (F) Spatial autocorrelations of (C) + (D).

https://doi.org/10.7554/eLife.10094.015

Simulation was run 60 times, with 400 outputs per run. 60° Gridness score mean was 1.1 ± 0.0006 when weights were constrained and 0.29 ± 0.0005 when weights were unconstrained. 90° Gridness score mean was 0.006 ± 0.002 when weights were constrained and 0.8 ± 0.0017 when weights were unconstrained.

Effect of place cell parameters on grid structure

A more detailed view of the resulting grid spacing showed that it was heavily dependent on the field widths of the place cells inputs. When the environment size was fixed and the output calculated per input size, the grid-spacing (distance between neighboring peaks) increased for larger place cell field widths.

To enable a fast parameter sweep over many place cell field widths (and large environment sizes), we took the steady state limit, and the limit of a high density of place cell locations, and used the fast FISTA algorithm to solve the non-negative PCA problem (see Materials and methods section).

We performed multiple simulations, and found that there was a simple linear dependency between the place field size and the output grid scale. For the case of periodic boundary conditions, we found that grid scale was S = 7.5sigma+0.85, where sigma was the width of the place cell field (Figure 12A). For a different set of simulations with zero boundary conditions, we achieved a similar relation: S=7.54sigma+0.62 (figure not shown). Grid scale was more dependent on place field size and less on box size (Figure 12H). We note that for very large environments, the effects of boundary conditions diminishes. At this limit, this linear relation between place field size and grid scale can be explained from analytical considerations (see Materials and methods section). Intuitively, this follows from dimensional analysis: given an infinite environment, at steady state the length scale of the place cell field width is the only length scale in the model, so any other length scale must be proportional to this scale. More precisely, we can provide a lower bound for the linear fit (Figure 12A), which depends only on the tuning curve of the place cells (see Materials and methods section). This lower bound was derived for periodic boundary conditions, but works well even with zero boundary conditions (not shown).

Effect of changing the place-field size in fixed arena (FISTA algorithm used; periodic boundary conditions and Arena size 500);

(A) Grid scale as a function of place field size (sigma); Linear fit is: Scale = 7.4 Sigma+0.62; the lower bound, equal to 2π/k, were k is defined in Equation 32 in the Materials and methods section; (B) Grid orientation as a function of gridness; (C) Grid orientation as a function of sigma – scatter plot (blue stars) and mean (green line); (D) Histogram of grid orientations; (E) Mean gridness as a function of sigma; and (F) Histogram of mean gridness. (G) Gridness as a function of sigma and (arena-size/sigma) (zero boundary conditions). (H) Grid scale for the same parameters as in G.

https://doi.org/10.7554/eLife.10094.016

Furthermore, we found that the grid orientation varied substantially for different place cell field widths, in the possible range of 0–15 degrees (Figure 12C,D). For small environments, the orientation strongly depended on the boundary conditions. However, as described in the Methods section, analytical considerations suggest that as the environment grows, the distribution of grid orientations becomes uniform in the range of 0–15 degrees, with a mean at 7.5°. Intuitively, this can be explained by rotational symmetry – when the environment size is infinite, all directions in the model are equivalent, and so we should get all orientations with equal probability, if we start the model from a uniformly random initialization. In addition, grid orientation was not a clear function of the gridness of the obtained grid cells (Figure 12B). For large enough place cells, gridness was larger than 1 (Figure 12E–G).

Modules of grid cells

It is known that in reality grid cells form in modules of multiple spacings (Barry et al., 2007; Stensola et al., 2012). We tried to address this question of modules in several ways. First, we used different widths for the Gaussian/Laplacian input functions: Initially, we placed a heterogeneous population of widths in a given environment (i.e., uniformly random widths) and ran the single-output network 100 times. The distribution of grid spacings was almost comparable to the results of the largest width if applied alone, and did not exhibit module like behavior. This result is not surprising when thinking about a small place cell overlapping in space with a large place cell. Whenever the agent passes next to the small one, it activates both weights via synaptic learning. This causes the large firing field to overshadow the smaller one. Additionally, when using populations of only two widths of place fields, the grid spacings were dictated by the size of the larger place field (data not shown).

The second option we considered was to use a multi-output neural network, capable of computing all 'eigenvectors' rather than only the principal 'eigenvector' (where by 'eigenvector' we mean here the vectors achieved under the positivity constraint, and not the exact eigenvectors themselves). We used a hierarchical network implementation introduced by Sanger, 1989 (see Materials and methods). Since the 1st output’s weights converged to the 1st 'eigenvector', the network (Figure 13A–B) provided to the subsequent outputs (2nd, 3rd, and so forth) a reduced-version of the data from which the projection of the 1st 'eigenvector' has been subtracted out. This process, reminiscent of Gram-Schmidt orthogonalization, was capable of computing all 'eigenvectors' (in the modified sense) of the input's covariance matrix. It is important to note though that, due to the non-negativity constraint, the vectors achieved in this way were not orthogonal, and thus it cannot be considered a real orthogonalization process, although, as explained in the Methods section, the process does aim for maximum difference between the vectors.

Hierarchial network capable of computing all 'principal components'.

(A) Each output is a linear sum of all inputs weighted by the corresponding learned weights. (B) Over time, the data the following outputs 'see' is the original data after subtration of the 1st 'eigenvector's' projection onto it. This is an iterative process causing all outputs' weights to converge to the 'prinipcal components' of the data.

https://doi.org/10.7554/eLife.10094.017

When constrained to be non-negative, and using the same homogeneous 'place cells' as in the previous network, the networks' weights converged to hexagonal shapes. Here, however, we found that the smaller the 'eigenvalue' was (or the higher the principal component number) the denser the grid became. We were able to identify two main populations of grid-distance 'modules' among the hexagonal spatial solutions with high Gridness scores (>0.7, Figure 14A–B). In addition, we found that the ratio between the distances of the modules was −1.4, close to the value of 1.42 found by Stensola et al. (Stensola et al., 2012). Although we searched for additional such jumps, we could only identify this single jump, suggesting that our model can yield up to two 'modules' and not more. The same process was repeated using the direct PCA method, utilizing the covariance matrix of the data after simulation as input for the non-negative PCA algorithms, and considering their ability to calculate only the 1st 'eigenvector'. By iteratively projecting the 1st 'eigenvector' on the simulation data and subtracting the outcome from the original data, we applied the non-negative PCA algorithm to the residual data obtaining the 2nd 'eigenvector' of the original data. This 'eigenvector' now constituted the 1st eigenvector' of the new residual data (see Materials and methods). Applying this process to as many 'outputs' as needed, we obtained very similar results to the ones presented above using the neural network (data not shown).

Modules of grid cells.

(A) In a network with 50 outputs, the grid spacing per output is plotted with respect to the hierarchical place of the output. (B) The grid spacing of outputs with high Gridness score (>0.7). The centroids have a ratio of close to 2. (C) + (D) Example of rate maps of outputs and their spatial autocorrelations for both of the modules.

https://doi.org/10.7554/eLife.10094.018

Discussion

In our work, we explored the nature and behavior of the feedback projections from place cells to grid cells. We shed light on the importance of this relation and showed, both analytically and in simulation, how a simple single-layer neural network could produce hexagonal grid cells when subjected to place cell-like temporal input from a randomly-roaming moving agent. We found that the network resembled a neural network performing PCA (Oja, 1982), with the constraint that the weights were non-negative. Under these conditions, and also under the requirements that place cells have a zero mean in time or space, the first principal component in the 2D arena had a firing pattern resembling a hexagonal grid cell. Furthermore, we found that in the limit of very large arenas, grid orientation converged to a uniform distribution in range of 0–15°. When looking at additional components, grid scale tends to be discretely clustered, such that two modules emerge. This is partially consistent with current experimental findings (Stensola et al., 2012; 2015). Furthermore, the inhibitory connectivity between multiple grid cells is consistent with the known functional anatomy in this network (Couey et al., 2013).

Place-to-Grid as a PCA network

As a consequence of the requirements for PCA to hold, we found that the place cell input needed to have a zero-mean, otherwise the output was not periodic. Due to the lack of the zero-mean property in 2D Gaussians, we used various approaches to impose zero-mean on the input data. The first, in the time domain, was to differentiate the input and use the derivatives (a random walk produces zero-mean derivatives) as inputs. Another approach was to dynamically subtract the mean in all iterations of the simulation. This approach was reminiscent of the adaptation procedure suggested in the Kropff & Treves paper (Kropff and Treves, 2008). A third approach, applied in the spatial domain was to use inputs with a zero-spatial mean such as Laplacians of Gaussians (Mexican hats in 2D, or differences-of-Gaussians) or negative – positive disks. Such Mexican-hat inputs are quite typical in the nervous system (Wiesel and Hubel, 1963; Enroth-Cugell and Robson, 1966; Derdikman et al., 2003), although in the case of place cells it is not completely known how they are formed. They could be a result of interaction between place cells and the vast number of inhibitory interneurons in the local hippocampal network (Freund and Buzsáki, 1996).

Another condition we found crucial, which was not part of the original PCA network, was a non-negativity constraint on the place-to-grid learned weights. While rather easy to implement in the network, adding this constraint to the non-convex PCA problem was harder to implement. Since the problem is NP-hard (Montanari and Richard, 2014), we turned to numerical methods. We used three different algorithms (Montanari and Richard, 2014; Zass and Shashua, 2006; Beck and Teboulle, 2009) to find the leading 'eigenvector' of every given temporal based input. As shown in the results section, both processes (i.e. direct PCA and the neural network) resulted in hexagonal outcomes when the non-negativity and zero-mean criteria were met. Note that the ease of use of the neural network for solving the positive PCA problem is a nice feature of the neural network implementation, and should be investigated further.

We also note that while our network focused on the projection from place cells to grid cells, we cannot preclude the importance of the reciprocal projection from grid cells to place cells. Further study will be needed to ‘close the loop’ and simultaneously consider both of these projections at once.

Similar studies

We note that similar work has noticed the relation between place-cell-to-grid-cell transformation and PCA. Notably, Stachenfeld et al., (2014) have demonstrated, from considerations related to reinforcement learning, that grid cells could be related to place cells through a PCA transformation. However, due to the unconstrained nature of their transformation, the resulting grid cells were square-like. Furthermore, there has been an endeavor to model the transformation from place cells to grid cells using independent-component-analysis (Franzius et al., 2007).

We also note that there is now a surge of interest in the feedback projection from place cells to grid-cells, which is inverse to the anatomical downstream direction from grid cells to place cells (Witter and Amaral, 2004) that has guided most of the models to-date (Zilli, 2012; Giocomo et al., 2011). In addition to several papers from the Treves group, in which the projection from place cells to grid cells is studied (Kropff and Treves, 2008; Si and Treves, 2013), there has been also recent work from other groups as well exploring this direction (Castro and Aguiar, 2014; Stepanyuk, 2015). As far as we are aware, none of the previous studies noted the importance of the non-negativity constraint and the requirement of zero mean input. Additionally, to the best of our knowledge, the analytic results and insights provided in this work (see Materials and methods) are novel, and provide a mathematically consistent explanation for the emergence of hexagonally-spaced grid cells.

Predictions of our model

Based on the findings of this work, it is possible to make several predictions. First, the grid cells must receive zero-mean input over time to produce hexagonally shaped firing patterns. With all feedback projections from place cells being excitatory, the lateral inhibition from other neighboring grid cells might be the balancing parameter to achieve the temporal zero-mean (Couey et al., 2013). Alternatively, an adaptation method, such as the one suggested in Kropff and Treves, (2008) may be applied. Second, if indeed the grid cells are a lower dimensional representation of the place cells in a PCA form, the place-to-grid neural weights distribution should be similar across identically spaced grid cell populations. This is because all grid cells with similar spacing would have maximized the variance over the same input, resulting in similar spatial solutions. As an aside, we note that such a projection may be a source of phase-related correlations in grid cells (Tocker et al., 2015). Third, we found a linear relation between the size of the place cells and the spacing between grid cells. Furthermore, the spacing of the grid cells is mostly determined by the size of the largest place cell – predicting that the feedback from large place cells is not connected to grid cells with small spacing. Fourth, we found modules of different grid spacings in a hierarchical network with the ratio of distances between successive units close to √2. This result is in accordance with the ratio reported in Stensola et al., (2012). However, we note that there is a difference between our results and experimental results because the analysis predicts that there should only be two modules, while the data show at least 5 modules, with a range of scales, the smallest and most numerous having approximately the scale of the smaller place fields found in the dorsal hippocampus (25–30 cm). Fifth, for large enough environments our model suggests that, from mathematical considerations, the grid orientation should approach a uniform orientation in the possible range of 0–15°. This is in discrepancy with experimental results which measure a peak at 7.5°, and not a uniform distribution (Stensola et al., 2015). As noted, the discrepancies between our results and reality may relate to the fact that a more advanced model will have to take into account both the downstream projection from grid cells to place cells together with the upstream projection from place cells to grid cells discussed in this paper. Furthermore, such a model will have to take into account the non-uniform distribution of place-cell widths (Kjelstrup et al., 2008).

Why hexagons?

In light of our results, we further asked what is special about the hexagonal shape which renders it a stable solution. Past works have demonstrated that hexagonality is optimal in terms of efficient coding. Two recent papers have addressed the potential benefit of encoding by grid cells. Mathis et al., (2015) considered the decoding of spatial information based on a grid-like periodic representation. Using lower bounds on the reconstruction error based on a Fisher information criterion, they demonstrated that hexagonal grids lead to the highest spatial resolution in two dimensions (extensions to higher dimensions were also provided). The solution is obtained by mapping the problem onto a circle packing problem. The work of Wei et al., (2013) also took a decoding perspective, and showed that hexagonal grids minimize the number of neurons required to encode location with a given resolution. Both papers offer insights into the possible information theoretic benefits of the hexagonal grid solution. In the present paper, we were mainly concerned with a specific biologically motivated learning (development) mechanism that may yield such a solution. Our analysis suggests that the hexagonal patterns can arise as a solution that maximizes the grid cell output variance, under non-negativity constraints. In Fourier space, the solution is a hexagonal lattice with lattice constant near the peak of the Fourier transform of the place cell tuning curve (Figures 15 and 16; see Materials and methods).

To conclude, this work demonstrates how grid cells could be formed from a simple Hebbian neural network with place cells as inputs, without needing to rely on path-integration mechanisms.

Materials and methods

All code was written in MATLAB, and can be obtained on https://github.com/derdikman/Dordek-et-al.-Matlab-code.git or on request from authors.

Neural network architecture

We implemented a single-layer neural network with feedforward connections that was capable of producing a hexagonal-like output (Figure 2). The feedforward connections were updated according to a self-normalizing version of a Hebbian learning rule referred to as the Oja rule (Oja, 1982),

(1) ΔJit=εt(ψtrit(ψt)2Jit),

where εt denotes the learning rate, Jit is the ith weight and ψt,rit are the output and the ith input of the network, respectively (all at time t). The weights were initialized randomly according to a uniform distribution and then normalized to have norm 1. The output ψt was calculated every iteration by summing up all pre-synaptic activity from the entire input neuron population. The activity of each output was processed through a sigmoidal function (e.g., tanh) or a simple linear function. Formally,

(2) ψt=fi=1nJitrit,

where n is the number of input place cells. Since we were initially only concerned with the eigenvector associated with the largest eigenvalue, we did not implement a multiple-output architecture. In this formulation, in which no lateral weights were used, multiple outputs were equivalent to running the same setting with one output several times.

As discussed in the introduction, this kind of simple feedforward neural network with linear activation and a local weight update in the form of Oja’s rule (1) is known to perform Principal Components Analysis (PCA) (Oja, 1982; Sanger, 1989; Weingessel and Hornik, 2000). In the case of a single output the feedforward weights converge to the principal eigenvector of the input's covariance matrix. With several outputs, and lateral weights, as described in the section on modules, the weights converge to the leading principal eigenvectors of the covariance matrix, or, in certain cases (Weingessel and Hornik, 2000), to the subspace spanned by the principal eigenvectors. We can thus compare the results of the neural network to those of the mathematical procedure of PCA. Hence, in our simulation, we (1) let the neural networks' weights develop in real time based on the current place cell inputs. In addition, we (2) saved the input activity for every time step to calculate the input covariance matrix and perform (batch) PCA directly.

It is worth mentioning that the PCA solution described in this section can be interpreted differently based on the Singular Value Decomposition (SVD). Denoting by R the T×d spatio-temporal pattern of place cell activities (after setting the mean to zero), where T is the time duration and d is the number of place cells, the SVD decomposition (see Jolliffe, 2002; sec. 3.5) for R is R=ULA'. For a matrix R of rank r, L is a r×r diagonal matrix whose kth element is equal to lk1/2, the square root of the kth eigenvalue of the covariance matrix RR' (computed in the PCA analysis), A is the d×r matrix with kth column equal to the kth eigenvector of RR', and U is the T×r matrix whose kth column is lk1/2Rak. Note that U is a T×r dimensional matrix whose kth column represents the temporal dynamics of the kth grid cell. In other words, the SVD provides a decomposition of the place cell activity in terms of the grid cell activity, as opposed to the grid cell representation in terms of place cell activity we discussed so far. The network learns the spatial weights over place cells (the eigenvectors) as the connections weights from the place cells, and 'projection onto place cell space' (lk1/2Rak) is simply the firing rates of the output neuron plotted against the location of the agent.

The question we therefore asked was under what conditions, when using place cell-like inputs, a solution resembling hexagonal grid cells emerges. To answer this we used both the neural-network implementation and the direct calculation of the PCA coefficients.

Simulation

We simulated an agent moving in a 2D virtual environment consisting of a square arena covered by n uniformly distributed 2D Gaussian-shaped place cells, organized on a grid, given by

(3) rit(X(t))=exp((X(t)Ci)22σi2),i=1,2,...,n

where X(t) represents the location of the agent. The variables rit constitute the temporal input from place cell i at time t, and Ci,σi are the ith place cell’s field center and width, respectively (see variations on this input structure below). In order to eliminate boundary effects, periodic boundary conditions were assumed. The virtual agent moved about in a random walk scheme (see Appendix) and explored the environment (Figure 1A). The place cell centers were assumed to be uniformly distributed (Figure 1B) and shared the same standard deviation σ. The activity of all place cells as a function of time (r(t)1,r(t)2r(t)n) was dependent on the stochastic movement of the agent, and formed a [Neuron x Time] matrix (rRnxT, with T- being the time dimension, see Figure 1C).

The simulation was run several times with different input arguments (see Table 1). The agent was simulated for T time steps, allowing the neural network's weights to develop and reach a steady state by using the learning rule (Equations 1,2) and the input (Equation 3) data. The simulation parameters are listed below and include parameters related to the environment, simulation, agent and network variables.

Table 1

List of variables used in simulation.

https://doi.org/10.7554/eLife.10094.019
Environment:Size of arenaPlace cells field widthPlace cells distribution
Agent:Velocity (angular & linear)Initial position-------------------
Network:# Place cells/ #Grid cellsLearning rateAdaptation variable (if used)
Simulation:Duration (time)Time step-------------------

To calculate the PCA directly, we used the MATLAB function Princomp in order to evaluate the n principal eigenvectors {qk}k=1n and corresponding eigenvalues of the input covariance matrix. As mentioned in the Results section, there exists a near fourfold redundancy in the eigenvectors (X-Y axis and in phase). Figure 3 demonstrates this redundancy by plotting the eigenvalues of the covariance matrix. The output response of each eigenvector qk corresponding to a 2D input location (x,y) is

(4) Φ(x,y)k=j=1nqkj exp((xcxj)22σx2(ycyj)22σy2),k=1,2,...,n

where cxj and cyj are the x,y components of the centers of the individual place cell fields. Unless otherwise mentioned, we used place cells in a rectangular grid, such that a place cell is centered at each pixel of the image (that is – number of place cells equals the number of image pixels).

Non-negativity constraint

Projections between place cells and grid cells are known to be primarily excitatory (Witter and Amaral, 2004), thus if we aim to mimic the biological circuit, a non-negativity constraint should be added to the feedforward weights in the neural network. While implementing a non-negativity constraint in the neural network is rather easy (a simple rectification rule in the weight dynamics, such that weights which are smaller than 0 are set to 0), the equivalent condition for calculating non-negative Principal Components is more intricate. Since this problem is non-convex and, in general, NP-hard (Montanari and Richard, 2014), a numerical procedure was imperative. We used three different algorithms for this purpose.

The first (Zass and Shashua, 2006) named NSPCA (Nonnegative Sparse PCA) is based on coordinate-descent. The algorithm computes a non-negative version of the covariance matrix's eigenvectors and relies on solving a numerical optimization problem, converging to a local maximum starting from a random initial point. The local nature of the algorithm did not guarantee a convergence to a global optimum (recall that the problem is non-convex). The algorithm's inputs consisted of the place cell activities’ covariance matrix, α - a balancing parameter between reconstruction and orthonormality, β – a variable which controls the amount of sparseness required, and an initial solution vector. For the sake of generality, we set the initial vector to be uniformly random (and normalized), α was set to a relatively high value – 104 and since no sparseness was needed, β was set to zero.

The second algorithm (Montanari and Richard, 2014) does not require any simulation parameters except an arbitrary initialization. It works directly on the inputs and uses a message passing algorithm to define an iterative algorithm to approximately solve the optimization problem. Under specific assumptions it can be shown that the algorithm asymptotically solves the problem (for large input dimensions).

The third algorithm we use is the parameter free Fast Iterative Threshold and Shrinkage algorithm FISTA (Beck and Teboulle, 2009). As described later in this section, this algorithm is the fastest of the three, and allowed us rapid screening of parameter space.

Different variants of input structure

Performing PCA on raw data requires the subtraction of the data mean. Some thought was required in order to determine how to perform this subtraction in the case of the neural network.

One way to perform the subtraction in the time domain was to dynamically subtract the mean during simulation by using the discrete 1st or 2nd derivatives of the inputs in time [i.e. from Equation 3r(t+1)=r(t+1)rt]. Under conditions of an isotropic random walk (namely, given any starting position, motion in all directions is equally likely) it is clear that E[r(t)]=0. Another option for subtracting the mean in the time domain was the use of an adaptation variable, as was initially introduced by Kropff and Treves, (2008). Although originally exploited for control over the firing rate, it can be viewed as a variable that represents subtraction of a weighted sum of the firing rate history. Instead of using the inputs rit directly in Equation 2 to compute the activation ψt, an intermediate adaptation variable ψadpt(δ) was used (δ being the relative significance of the present temporal sample) as

(5) ψadpt=ψtψ¯t,
(6) ψ¯t=(1δ)ψ¯t1+δψt.

It is not hard to see that for i.i.d. variables ψadpt, the sequence ψ¯t converges for large t to the mean of ψt. Thus, when t we find that E[ψadpt]0, specifically, the adaptation variable is of zero asymptotic mean.

The second method we used to enforce a zero mean input was simply to create it in advance. Rather than using 2D Gaussian functions (i.e. [Equation 3]) as inputs we used 2D difference-of-Gaussians (all σ are equal in x and y axis):

(7) rit(X(t))=c1,i exp((X(t)Ci)22σ1,i2)c2,i exp((X(t)Ci)22σ2,i2),i=1,2,...,n

where the constants c1 and c2 are set so the integral of the given Laplacian function is zero (if the environment size is not too small, then c1,i/c2,i  σ2,i/σ1,i). Therefore, if we assume a random walk that covers the entire environment uniformly, the temporal mean of the input would be zero as well. Such input data can be inspired by similar behavior of neurons in the retina and the lateral-geniculate nucleus (Wiesel and Hubel, 1963; Enroth-Cugell and Robson, 1966). Finally, we implemented another input data type; positive-negative disks (see Appendix). Analogously to the difference-of-Gaussians function, the integral over input is zero so the same goal (zero-mean) was achieved. It is worthwhile noting that subtracting a constant from a simple Gaussian function is not sufficient since at infinity it does not reach zero.

Quality of solution and Gridness

In order to test the hexagonality of the results we used a hexagonal Gridness score (Sargolini et al., 2006). The Gridness score of the spatial fields was calculated from a cropped ring of their autocorrelogram including the six maxima closest to the center. The ring was rotated six times, 30 per rotation, reaching in total angles of 30,60,90,120,150. Furthermore, for every rotated angle the Pearson correlation with the original un-rotated map was obtained. Denoting by Cγ the correlation for a specific rotation angle γ, the final Gridness score was (Kropff and Treves, 2008):

(8) Gridness=12(C60+C120)13(C30+C90+C150).

In addition to this 'traditional' score we used a Squareness Gridness score in order to examine how square-like the results are spatially. The special reference to the square shape was driven by the tendency of the spatial solution to converge to a rectangular shape when no constrains were applied. The Squareness Gridness score is similar to the hexagonal one, but now the cropped ring of the autocorrelogram is rotated 45 every iteration to reach angles of 45,90,135. As before, denoting Cγ as the correlation for a specific rotation angle γ the new Gridness score was calculated as:

(9) Square Gridness=C9012(C45+C135).

All errors calculated in gridness measures are SEM (Standard Error of the Mean).

Hierarchical networks and modules

As described in the Results section, we were interested to check whether a hierarchy of outputs could explain the module phenomenon described for real grid cells. We replaced the single-output network with a hierarchical, multiple outputs network, which is capable of computing all 'principal components' of the input data while maintaining the non-negativity constraint as before. The network, introduced by Sanger, 1989, computes each output as a linear summation of the weighted inputs similar to Equation 2. However, the weights are now calculated according to:

(10) ΔJijt=εt(rjtψitψitk=1iJkjtψkt).

The first term in the parenthesis when k=1 was the regular Hebb-Oja derived rule. In other words, the first output calculated the first non-negative 'principal component' (in inverted commas due to the non-negativity) of the data. Following the first one, the weights of each output received a back projection from the previous outputs. This learning rule applied to the data in a similar manner to the Gram-Schmidt process, subtracting the 'influence' of the previous 'principal components' on the data and recalculating the appropriate 'principal components' of the updated input data.

In a comparable manner, we applied this technique to the input data X in order to obtain non-negative 'eigenvectors' from the direct nonnegative-PCA algorithms. We found V2 by subtracting from the data the projection of V1 on it,

(11) X~=XV1T(V1X).

Next, we computed V2, the first non-negative 'principal component' of X~, and similarly the subsequent ones.

Stability of hexagonal solutions

In order to test the stability of the solutions we obtained under all types of conditions, we applied the ODE method (Kushner and Clark, 1978; Hornik and Kuan, 1992; Weingessel and Hornik, 2000) to the PCA feature extraction algorithm introduced in pervious sections. This method allows one to asymptotically replace the stochastic update equations describing the neural dynamics by smooth differential equations describing the average asymptotic behavior. Under appropriate conditions, the stochastic dynamics converge with probability one to the solution of the ODEs. Although originally this approach was designed for a more general architecture (including lateral connections and asymmetric updating rules), we used a restricted version for our system. In addition, the following analysis is accurate solely for linear output functions. However, since our architecture works well with either linear or non-linear output functions, the conclusions are valid.

We can rewrite the relevant updating equations of the linear neural network (in matrix form), (see [Weingessel and Hornik, 2000] Equations 15–19):

(12) ψt+1=QJt(rt)T,
(13) ΔJt=εt(ψt(rt)TΦ(ψt(ψt)T)Jt).

In our case we set

Q=I,Φ=diag.

Consider the following assumptions

  1. The input sequence rt consists of independent identically distributed, bounded random variables with zero-mean.

  2. {εt} is a positive number sequence satisfying: t εt=,t (εt)2<.

A typical suitable sequence is εt=1t,t=1,2.

For long times, we denote

(14) E[ψt(rt)T]E[JrrT]=E[J]E[rrT]=JΣ,
(15) limtE[ψψT]=E[J]E[rrT]E[JT]=JΣJT.

The penultimate equalities in these equations used the fact that the weights converge with probability one to their average value, resulting from the solution of the ODEs. Following Weingessel and Hornik, (2000), we can analyze Equations 12,13 under the above assumptions, via their asymptotically equivalent associated ODEs

(16) dJdt=JΣdiag(JΣJT) J,

with equilibria at

(17) JΣ=diag(JΣJT) J.

We solved it numerically by exploiting the same covariance matrix and initializing with random weights J. In line with our previous findings, we found that constraining J to be non-negative (by a simple cut-off rule) resulted in a hexagonal shape (in the projection of J onto the place cells space; Figure 11). In contrast, when the weights were not constrained they converged to square-like results.

Steady state analysis

From this point onwards, we focus on the case of a single output, in which J is a row vector, unless stated otherwise. In the unconstrained case, from Equation 17 any J which is a normalized eigenvector of Σ would be a fixed point. However, from Equation 16, only the principal eigenvector, which is the solution to the following optimization problem

(18) maxJ:JTJ=1 JΣJT

would correspond to a stable fixed point. This is the standard PCA problem. By adding the constraint J0 we get the non-negative PCA problem.

To speed up simulation and simplify analysis we make further simplifications.

First, we assume that the agent’s random movement is ergodic (e.g., an isotropic random walk in a finite box as we used in our simulation), uniform and covering the entire environment, so that

(19) JΣJT=Eψ2(X(t))=1|S|Sψ2(x)dx,

where x denotes location vector (in contrast to X(t), which is the random process corresponding to the location of the agent), S is the entire environment, and |S| is the size of the environment.

Second, we assume that the environment S is uniformly and densely covered by identical place cells, each of which has the same a tuning curve r(x) (which integrates to zero). In this case, the activity of the linear grid cell becomes a convolution operation

(20) ψ(x)=SJ(x')r(xx')dx',

where J(x) is the synaptic weight connecting to the place cell at location x.

Thus, we can write our objective as

(21) 1|S|Sψ2(x)dx=1|S|S(SJ(x')r(xx')dx')2dx

under the constraint that the weights are normalized

(22) 1SSJ2 (x) dx=1,

where either J(x) (PCA) or J(x)  0 ('non-negative PCA').

Since we expressed the objective using a convolution operation (different boundary conditions can be assumed), it can be solved numerically considerably faster. In the non-negative case, we used the parameter free Fast Iterative Threshold and Shrinkage algorithm [FISTA (Beck and Teboulle, 2009); in which we do not use shrinkage, since we only have hard constraints], where the gradient was calculated efficiently using convolutions.

Moreover, as we show in the following sections, if we assume periodic boundary conditions and use Fourier analysis, we can analytically find the PCA solutions, and obtain important insight on the non-negative PCA solutions.

Fourier notation

Any continuously differentiable function f(x), defined over S[0,L]D, a 'box' region in D dimensions, with periodic boundary conditions, can be written using a Fourier series

(23) f^(k)1|S|Sf(x)eikxdx,f(x)kS^f^(k)eikx,

where |S|=LD is the volume of the box and

S^{(2m1πL,,2mdπL)}(m1,,md)D

is the reciprocal lattice of S in k-space (frequency space).

PCA solution

Assuming periodic boundary conditions, we use Parseval’s identity, and the properties of the convolution, to transform the steady state objective (Equation 21) to its simpler form in the Fourier domain,

(24) 1|S|Sψ2(x)dx=kS^|J^(k)r^(k)|2.

Similarly, the normalization constraint can also be written in the Fourier domain,

(25) 1|S|SJ2(x)dx=kS^|J^(k)|2=1.

Maximizing the objective Equation 24 under this constraint in the Fourier domain, we immediately get that any solution is a linear combination of the Fourier components,

(26) J^(k) ={1,     k=k*0,     kk*.

where

(27) k*argmaxkS^r^ (k),

and J^(k) satisfies the normalization constraint. In the original space, the Fourier components are

(28) J(x)=eik·x+iϕ,

where ϕ[0,2π) is a free parameter that determines the phase. Also, since J(x) should assume real values, it is composed of real Fourier components

(29) J(x)=12(eik*x+iϕ+eik*xiϕ)=2cos(k*x+ϕ).

This is a valid solution, since r(x) is a real-valued function, r^(k)=r^(k) and therefore -k* argmaxkS^r^(k).

PCA solution for a difference of Gaussians tuning curve

In this paper we focused on the case where r(x) has the shape of a difference of Gaussians (Equation 7),

(30) r(x)c1 exp(x22σ12)c2 exp(x22σ22)

where c1 and c2 are some positive normalization constants, set so that Sr(x)dx=0 (see appendix). The Fourier transform of r(x) is also a difference of Gaussians

(31) r^(k)exp(12σ12k2)exp(12σ22k2)

kS^, as we show in the appendix. Therefore the value of the Fourier domain objective only depends on the radius k, and all solutions k* have the same radius k*. If L, then the k-lattice S^ becomes dense (S^D) and this radius is equal to

(32) k=argmaxk0 exp(12σ12k2)exp(12σ22k2)

which is a unique maximizer, that can be easily obtained numerically.

Notice that if we multiply the place cell field width by some positive constant c, then the solution k will be divided by c. The grid spacing, proportional to 1k, would therefore also be multiplied by c. This entails a linear dependency between the place cell field width and the grid cell spacing, in the limit of a large box size (L). When the box has a finite size, k-lattice discretization also has a (usually small) effect on the grid spacing.

In that case, all solutions k* are restricted to be on the finite lattice S^. Therefore, the solutions k* are the points on the lattice S^ for which the radius k* is closest to k (see Figure 15B,C).

PCA k-space analysis for a difference of Gaussians tuning curve.

(A) The 1D tuning curve r(x). (B) The 1D tuning curve Fourier transform r^(k). The black circles indicate k-lattice points. The PCA solution, k*, is given by the circles closest to k, the peak of r^(k) (red cross). (C) A contour plot of the 2D tuning curve Fourier transform r^(k). In 2D k-space the peak of r^(k) becomes a circle (red), and the k-lattice S^ is a square lattice (black circles). The lattice point can be partitioned into equivalent groups. Several such groups are marked in blue on the lattice. For example, the PCA solution Fourier components lie on the four lattice points closest to the circle, denoted A1-4. Note the grouping of A,B,C & D (4,8,4 and 4, respectively) corresponds to the grouping of the 20 highest principal components in Figure 4. Parameters: 2σ1=σ2=7.5,L=100.

https://doi.org/10.7554/eLife.10094.020

The degeneracy of the PCA solution

The number of real-valued PCA solutions (degeneracy) in 1D is two, as there are exactly two maxima, k* and k*. The phase ϕ, determines how the components at k* and k* are linearly combined.

However, there are more maxima in the 2D case. Specifically, given a maximum k*, we can write (m,n)=L2πk*, where (m,n)2. Usually there are 7 other different points with the same radius: (m,n),(m,n),(m,n),(n,m),(n,m),(n,m) and (n,m), so we will have a degeneracy of eight (corresponding to the symmetries of a square box). This is case of points in group B, shown in Figure 15C.

However, we can also get a different degeneracy. First, if either m=±n, n=0 or m=0 we will have a degeneracy of 4, since then some of the original eight points will coincide (groups A,C and D in Figure 15C). Second, additional points (k,r) can exist such that k2+r2=m2+n2, (Pythagorean triplets with the same hypotenuse) – for example, 152+202=252=72+242. These points will also appear in groups of four or eight.

Therefore, we will always have a degeneracy which is some multiple of 4. Note that in the full network simulation, the degeneracy is not exact. This is due to the perturbation noise from the agent’s random walk as well as the non-uniform sampling of the place cells.

The PCA solution with a non-negative constraint

Next, we add the non-negativity constraint J(x)0. As mentioned earlier, this constraint renders the optimization problem NP-hard, and prevents us from a complete analytical solution. We therefore combine numerical and mathematical analysis, in order to gain intuition as to why

  1. Locally optimal 2D solutions are hexagonal.

  2. These solutions have a grid spacing near (4π/(3k) (k is the peak of r^k).

  3. The average grid alignment is approximately 7.5°, for large environments.

  4. Why grid cells have modules, and what is their spacing.

1D Solutions

Our numerical results indicate that the Fourier components of any locally optimal 1D solution of non-negative PCA have the following structure:

  1. There is a non-negative 'DC component' (k = 0).

  2. The maximal non-DC component, (k ≠ 0) is k*, where k* is 'close' (more details below) to k, the peak of r^(k).

  3. All other non-zero Fourier components are mk*n=1, weaker harmonies of k*.

This structure suggests that the component at k* aims to maximize the objective, while the other components guarantee the non-negativity of the solution J(x). In order to gain some analytical intuition as to why this is the case, we first examine the limit case that L and r^(k) is highly peaked at k. In that case the Fourier objective (Equation 24) simply becomes 2|r^(k)|2|J^(k)|2. For simplicity, we will rescale our units so that |r^(k)|2=1/2, and the objective becomes |J^(k)|2. Therefore, the solution must include a Fourier component at k or the objective would be zero. The other components exist only to maintain the non-negativity constraint, since if they increase in magnitude, then the objective, which is proportional to |J^(k)|2, must decrease to compensate (due to the normalization constraint – Equation 25). Note that these components must include a positive 'DC component' at k=0, or else SJ(x)dxJ^(0)0, which contradicts the constraints. To find all the Fourier components, we examine a solution composed of only a few (M) components

J(x)=J^(0)+2m=1MJ^mcos(kmx+ϕm).

Clearly, we can set k1=k, or otherwise, the objective would be zero. Also, we must have

J^(0)=minx(2m=1MJ^mcos(kmx+ϕm))0.

Otherwise, the solution would be either (1) negative or (2) non-optimal, since we can decrease J^(0) and increase |J1|.

For M=1, we immediately get that, in the optimal solution, 2J^1=J^(0)=2/3 (ϕm does not matter). For M=2,3 and 4 a solution is harder to find directly, so we performed a parameter grid search over all the free parameters (km,J^m and ϕm) in those components. We found that the optimal solution (which maximizes the objective |J^(k)|2), had the following form

(33) J(x)=m=MMJ^(mk)cos(mk(xx0)),

where x0 is a free parameter. This form results from a parameter grid search for M=1,2,3 and 4, under the assumption that L and r^(k) is highly peaked. However, our numerical results in the general case (Figure 16A), using the FISTA algorithm, indicate that the locally optimal solution does not change much even if L is finite, and r^(k) is not highly peaked. Specifically, it has a similar form

(34) J(x)=m=J^(mk*)cos(mk*(xx0)).

Since J^(mk*) is rapidly decaying (Figure 16A), effectively only the first few components are non-negligible, as in Equation 33. This can also be seen in the value of the objective obtained in the parameter scan

(35) M1234|J^(k)|2160.23670.24570.2457,

where the contribution of additional high frequency components to the objective quickly becomes negligible. In fact, the value of the objective cannot increase above 0.25, as we explain in the next section.

And so, the main difference between Equations 33 and 34 is the base frequency, k*, which is slightly different from k. As explained in the appendix, the relation between k* and k depends on the k-lattice discretization, as well as on the properties of r^(k).

Fourier components of Non-negative PCA on the k-lattice.

(A) 1D solution (blue) includes: a DC component (k=0), a maximal component with magnitude near k (red line), and weaker harmonics of the maximal component. (B) 2D solution includes: a DC component (k = (0,0)), a hexgaon of strong components with radius near k (red circle), and weaker components on the lattice of the strong components. White dots show underlying k-lattice. We used a difference of Gaussians tuning curve, with parameters 2σ1=σ2=7.5,L=100, and the FISTA algorithm.

https://doi.org/10.7554/eLife.10094.021

2D Solutions

The 1D properties, described in the previous section, generalize to the 2D case in the following manner:

  1. There is a non-negative DC component (k=(0,0)).

  2. A small 'basis set' of components k*(i)i=1Bwith similar amplitudes, and with similar radii k*(i) which are all 'close' to k (details below).

  3. All other non-zero Fourier components are weaker, and restricted to the lattice

kS^ | k=i=1Bnik*(i)(n1,...,nB)B

Interestingly, given these properties of the solution we already get hexagonal patterns, as we explain next.

Similarly to the 1D case, the difference between k*(i) and k is affected by lattice discretization, and the curvature of r^(k) near k. To simplify matters, we focus first on the simple case that L and r^(k) is sharply peaked around k. Therefore, the Fourier objective becomes i=1B|J^(k*(i))|2, so the only Fourier components that appear in the objective are {k*(i)}i=1B, which have radius k. We examine the values this objective can have.

All the base components have the same radius. This implies, according to the Crystallographic restriction theorem in 2D, that the only allowed lattice angles (in the range between 0 and 90 degrees) are 0, 60 and 90 degrees. Therefore, there are only three possible lattice types in 2D. Next, we examine the value of the objective for each of these lattice types:

1) Square lattice, in which k*(1)=k(1,0), k*(2)=k(0,1), up to a rotation. In this case,

J(x,y)=mx=my=J^mx,mycos(k(mxx+myy)+ϕmx,my)

and the value of the objective is bounded above by 0.25 (see proof in appendix).

2) 1D lattice, in which k*(1)=k(1,0), up to a rotation. This is a special case of the square lattice, with a subset of J^mx,myequal to zero, so we can write, as we did in the 1D case

J(x)=m=J^mcos(kmx+ϕm)

Therefore, the same objective upper bound, 0.25, holds. Note that some of the solutions we found numerically are close to this bound (Equation 35).

3) Hexagonal lattice, in which the base components are

k*(1)=k(1,0),k*(2)=k(12,32),k*(3)=k(12,32)

up to a rotation by some angle α. Our parameter scans indicate that the objective value cannot surpass 0.2 in any solution composed of only the base hexgonal components {k*(m)}m=13 and a DC component. However, taking into account also some higher order lattice components, we can find a better solution, with an objective value of 0.2558. Though this is not necessarily the optimal solution, it surpasses any possible solutions on the other lattice types (bounded below 0.25, as we proved in the appendix). Specifically, this solution is composed of the base vectors {k*(m)}m=13 and their harmonics

J(x)=J^0+2m=18J^mcos(k*(m)x)

with k*(4)=2k*(1), k*(5)=2k*(2), k*(6)=2k*(3)k*(7)=k*(1)+k*(2), k*(8)=k*(1)+k*(3). Also, J^0=0.6449, J^1=J^2=J^3=0.292, J^4=J^5=J^6=0.0101 and J^7=J^8=0.134.

Thus, any optimal solution must be on the hexagonal lattice, given our approximations. In practice, the lattice hexagonal basis vectors do not have exactly the same radius, and, as in the 1D case, this radius is somewhat smaller then k, due to the lattice discretization, and due to that r^(k) is not sharply peaked. However, the resulting solution lattice is still approximately hexagonal in k-space. For example, this can be seen in the numerically obtained solution in Figure 16B – where the strongest non-DC Fourier components form an approximate hexagon near k, from the Fourier components A, defined in Figure 17.

The modules in Fourier space.

As in Figure 15C, we see a contour plot of the 2D tuning curve Fourier transform r^(k) and the k-space the peak of r^(k) (red circle), and the k-lattice S^ (black circles). The lattice points can be divided into approximately hexgonal shaped groups. Several such groups are marked in blue on the lattice. For example, group A and B are optimal since they are nearest to the red circle. The next best (with the highest-valued contours) group of points, which have an approximate hexgonal shape, is C. Note that group C has a k-radius of approximately the optimal radius times 2 (cyan circle). Parameters: 2σ1=σ2=7.5,L=100.

https://doi.org/10.7554/eLife.10094.022

Grid spacing

In general, we get a hexagonal grid pattern in x-space. If all base Fourier components have a radius of k, then the grid spacing in x-space would be 4π/(3k). Since the radius of the basis vectors can be smaller than k, the value of 4π/(3k) is a lower bound to the actual grid spacing (as demonstrated in Figure 12A), up to lattice discretization effects.

Grid alignment

The angle of the hexagonal grid, α, is determined by the directions of the hexagonal vectors. An angle α is possible, if there exists a k-lattice point k=2πL(m,n) (with m,n integers), for which πL>2πL(m2+n2)k*2, and then α=arctannm. Since the hexagonal lattice has rotational symmetry of 60, we can restrict α to be in the range 30α30. The grid alignment, which is the minimal angle of the grid with the box boundaries is given by

(36) Gridaligment=min(|α|,90(|α|+60))=min(|α|,30|α|)

which is limited to the range [0,15], since 30α30. There are usually several possible grid alignments which are (approximately) rotated versions of each other (i.e., different α). Note that, due to the k-lattice discretization, different alignments can result in slightly different objective values. However, the numerical algorithms we used to solve the optimization problem reached many possible grid alignments with a positive probability (Figure 12C), since we started from a random initialization and converged to a local minimum.

In the limit L, the grid alignment will become uniform in the range [0,15], and the average grid alignment is 7.5.

Hierarchical networks and modules

There are multiple routes to generalize non-negative PCA with multiple vectors. In this paper we chose to do so using a 'Gramm-Schmidt' like process, which can be written in the following way. First we define,

(37) r1(x,y)=r(xy)=c1 exp(xy22σ12)c2 exp(xy22σ22)J0(x)=0,,

and then, recursively, this process recovers non-negative 'eigenvectors' by subtracting out the previous components, similarly to Sanger’s multiple PCA algorithm (Sanger, 1989), and enforcing the non-negativity constraint.

(38) rn+1(x,y)=rn(x,y)-srn(x,z)Jn(z)Jn(y)dzJn(y)=arg maxJ(x)0,1SsJ2(x)dx=1sdxsrn(x,y)J(y)dy2.

To analyze this, we write the objectives we maximize in the Fourier domain, using Parseval’s Theorem.

For n =1, we recover the old objective (Equation 24):

(39) kS^r^(k)J^(k)2.

For n =2, we get

(40) kS^|r^(k)(J^(k)J^1(k)(qS^J^1*(q)J^(q)))|2,

where J^* is the complex conjugate of J^. This objective is similar to the original one, except that it penalizes J^(k) if its components are similar to those of J^1(k). As n increases the objective becomes more and more complicated, but as before, it contains terms which penalize J^n(k) if its components are similar to any of the previous solutions (i.e., J^m(k) for m<n). This form suggests that each new 'eigenvector' tends to occupy new points in the Fourier lattice (similarly to unconstrained PCA solutions).

For example, the numerical solution shown in Figure 16B is composed of the Fourier lattice components in group A, defined in Figure 17. A completely equivalent solution would be in group B (it is just a 90 degrees rotation of the first). The next 'eigenvectors' should then include other Fourier-lattice components outside groups A and B. Note that components with smaller k-radius cannot be arranged to be hexagonal (not even approximately), so they will have a low gridness score. In contrast, the next components with higher k-radius (e.g., group C) can form an approximately hexagonal shape together, and would appear as an additional grid cell 'module'. The grid spacing of this new module will decrease by 2, since the new k-radius is about 2 times larger than the k-radius of groups A and B.

Appendix

Movement schema of agent and environment data

The relevant simulation data used:

Size of arena 10X10Place cells field width: 0.75Place cells distribution:
uniform
Velocity: 0.25 (linear),
0.1-6.3 (angular)
# Place cells: 625Learning rate: 1/(t+1e5)

The agent was moved around the virtual environment according to:

Dt+1=modulo(Dt+ωZ,2π)xt+1=xt+νcos(Dt+1)yt+1=yt+νsin(Dt+1)

where Dt is the current direction angle, ω is the angular velocity, Z-N(0,1) where N is the standard normal Gaussian distribution, ν is the linear velocity, and (xt,yt) is the current position of the agent. Edges were treated as periodic – when agent arrived to one side of box it was teleported to the other side.

Positive-negative disks

Positive-negative disks are used with the following activity rules:

(rjt(xt,yt))={1ρ12ρ22ρ120,         (xtcx,j)2+(ytcy,j)2<ρ1,   ρ1<(xtcx,j)2+(ytcy,j)2<ρ2,   ρ2<(xtcx,j)2+(ytcy,j)2

Where cx, cy are the centers of the disks, and ρ1, ρ2 are the radii of the inner circle and the outer ring, respectively. The constant value in the negative ring was chosen to yield zero integral over the disk.

Fourier transform of the difference of Gaussians function

Here we prove that if we are given a difference of Gaussians function,

r(x)=i=12(1)i+1ci exp(x·x2σi2),

in which the appropriate normalization constants for a bounded box are

ci=[LLexp(z22σi2)dz]1=[2πσi2[2F(L/σi)1]]D,

with F(x)=12πσi2xexp(z22σi2)dz being the cumulative normal distribution, then kS^=m12πL,,mD2πLm1,,mDD, we have

r^(k)1SS r(x) eik·x dx = 1Si=12(-1)i+1exp(-12σi2 k·k).

Note that the normalization constants ci vanish after the Fourier transform, and that in the limit Lσ1,σ2 we obtain the standard result for the Fourier transform of an unbounded Gaussian distribution.

Proof: For simplicity of notation, we assume D=2. However, the calculation is identical for any D.

r^(k)=1|S|Sr(x)eikxdx=1|S|i=12LLdxLLdy(1)i+1ci exp(x2+y22σi2)exp(i(kxx+kyy))=1|S|i=12LLdx(1)i+1ci exp(x22σi2)exp(ikxx) · LLdy(1)i+1exp(y22σi2)exp(ikyy)=2D|S|i=12(1)i+1ciz{x,y}0Ldz exp(z22σi2)exp(ikzz)+exp(ikzz)=2D|S|i=12(1)i+1ciz{x,y}0Ldz exp(z22σi2)cos(kzz)

To solve this integral, we define

Ii (a)=0L dz exp-z22σi2 cos (az).

Its derivative is

I'i(a)=0Ldzz exp(z22σi2)sin(az)=σi2exp(z22σi2)sin(az)|z=0Laσi20Ldz exp(z22σi2)cos(az)=aσi2Ii(a)

where in the last equality we used the fact that aL=2πn (where n is an integer) if aS^. Solving this differential equation, we obtain

Ii(a)=Ii(0)exp(12σi2a2)

where

Ii(0)=0Ldz exp(-z22σi2)=2πσi2 F(Lσi)-F(0)=2πσi2 F(Lσi)-12.

substituting this into our last expression for r(k) we obtain

r(k)=1|S|i=12(1)i+1exp(12σi2z{x,y}kz2),

which is what we wanted to prove.

Why the 'unconstrained' k is a lower bound on 'constrained' k*

In section 'The PCA solution with a non-negative constraint – 1D Solutions', we mention that the main difference between the Equations (Equations 33 and 34) is the base frequency, k*, which is slightly different from k.

Here we explain how the relation between k* and k depends on the k-lattice discretization, as well as the properties of r^(k). The discretization effect is similar to the unconstrained case, and can cause a difference of at most π/L between k* to k. However, even if L, we can expect k* to be slightly smaller than k. To see that, suppose we have a solution as described above, with k*=k+δk and δk is a small perturbation (which does not affect the non-negativity constraint). We write the perturbed k-space objective of this solution

m=|J^(mk)|2|r^(m(k+δk))|2

Since k is the peak of |r^(k)|2, if |r^(k)|2 is monotonically decreasing for k>k (as is the case for the difference of Gaussians function), then any positive perturbation δk>0 would decrease the objective.

However, a sufficiently small negative perturbation δk<0 would improve the objective. We can see this from the Taylor expansion of the objective,

m=|J^(mk)|2(|r^(mk)|2+d|r^(k)|2dk|k=mkmδk+O(δk)2)

In which the derivative d|r^(k)|2dk|k=mk is zero for m=1 (since k is the peak of |r^(k)|2) and negative for m>1, if |r^(k)|2 is monotonically decreasing for k>k. If we gradually increase the magnitude of this negative perturbation δk<0, at some point the objective will stop increasing and start decreasing, because, for the difference of Gaussians function, |r^(k)|2 increases more sharply for k<k than it decreases for k>k. We can thus treat the 'unconstrained' k as an upper bound for the 'constrained' k*, if we ignore discretization effects (i.e., the limit L).

Upper bound on the constrained objective – 2D square lattice

In the Methods section 'The PCA solution with a non-negative constraint – 2D Solutions' we examine different 2D lattice-based solutions, including a square lattice, which is a solution of the following constrained minimization problem

max{Jm}m=0J^1,02+J^0,12s.t.:(1)mx= my=J^mx,my2=1    (2)mx= my=J^mx,mycos(k(mxx+myy)+ϕmx,my)0.

Here we prove that the objective is bounded above by 14, i.e., J^1,02+J^0,1214.

Without loss of generality we assume ϕ1,0=ϕ0,1=0,k=1 (we can always shift and scale x,y). We denote

P(x,y)=2J^0,1cos(x)+2J^1,0cos(y)

and

A(x,y)=mx=my=J^mx,mycos((mxx+myy)+ϕmx,my)P(x,y)

We examine the domain [0,2π]2. We denote by C± the regions in which P(x,y) is positive or negative, respectively. Note that

P((x+π)mod 2π,(y+π)mod 2π)=P(x,y)

so both regions have the same area, and

(A1) CP2(x,y)dxdy=C+P2(x,y)dxdy.

From the non-negativity constraint (2), we must have

(A2) (x,y)C:A(x,y)P(x,y).

Also, note that P(x,y) and A(x,y) do not share any common Fourier (cosine) components. Therefore, they are orthogonal, with the inner product being the integral of their product on the region [0,2π]2

0=[0,2π]2P(x,y)A(x,y)dxdy=C+P(x,y)A(x,y)dxdy+CP(x,y)A(x,y)dxdyC+P(x,y)A(x,y)dxdyC+P2(x,y)dxdy,

where in the last line we used the bound from Equation (A2), and then Equation (A1). Using this result, together with Cauchy-Schwartz inequality, we have

C+P2(x,y)dxdyC+P(x,y)A(x,y)dxdyC+P2(x,y)dxdyC+A2(x,y)dxdy,

So,

C+P2(x,y)dxdyC+A2(x,y)dxdy,

Summing this equation with Equation (A2), squared and integrated over the region C, and dividing by (2π)2, we obtain

1(2π)2[0,2π]2P2(x,y)dxdy1(2π)2[0,2π]2A2(x,y)dxdy,

using the orthogonality of the Fourier components (i.e., Parseval’s theorem) to perform the integrals over P2(x,y) and A2(x,y), we get

2J^1,02+2J^0,12mx=my=J^mx,my22J^1,022J^0,12.

Lastly, plugging in the normalization constraint mx=my=J^mx,my2=1, we find

14J^1,02+J^0,12

as required.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
    Imaging spatiotemporal dynamics of surround inhibition in the barrels somatosensory cortex
    1. D Derdikman
    2. R Hildesheim
    3. E Ahissar
    4. A Arieli
    5. A Grinvald
    (2003)
    Journal of Neuroscience 23:3100–3105.
  10. 10
  11. 11
  12. 12
    The contrast sensitivity of retinal ganglion cells of the cat
    1. C Enroth-Cugell
    2. JG Robson
    (1966)
    The Journal of Physiology 187:517–552.
  13. 13
  14. 14
  15. 15
  16. 16
    Toeplitz and circulant matrices: A review
    1. RM Gray
    (2006)
    Now Publishers Inc 2:155–239.
  17. 17
  18. 18
  19. 19
    Principal component analysis
    1. I Jolliffe
    (2002)
    Wiley Online Library, 10.1002/0470013192.bsa501.
  20. 20
  21. 21
  22. 22
    Applied Mathematical Sciences, Vol. 26
    1. HJ Kushner
    2. DS Clark
    (1978)
    Stochastic Approximation Methods for Constrained and Unconstrained Systems, Applied Mathematical Sciences, Vol. 26, New York, NY, Springer New York, 10.1007/978-1-4684-9352-8.
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
    The mantle of the heavens: Reflections on the 2014 Nobel Prize for medicine or physiology
    1. RG Morris
    (2015)
    Hippocampus, 25, 10.1002/hipo.22455.
  28. 28
    The hippocampus as a spatial map. Preliminary evidence from unit activity in the freely-moving rat
    1. J O'Keefe
    2. J Dostrovsky
    (1971)
    Brain Research 34:171–175.
  29. 29
    The hippocampus as a cognitive map, Oxford, Oxford University Press
    1. J O'Keefe
    2. L Nadel
    (1978)
    The hippocampus as a cognitive map, Oxford, Oxford University Press.
  30. 30
    A simplified neuron model as a principal component analyzer
    1. E Oja
    (1982)
    Journal of Mathematical Biology 15:267–273.
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
    Design Principles of the Hippocampal Cognitive Map
    1. KL Stachenfeld
    2. M Botvinick
    3. SJ Gershman
    (2014)
    NIPS Proceedings pp. 2528–2536.
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
    Local PCA algorithms
    1. A Weingessel
    2. K Hornik
    (2000)
    Neural Networks, IEEE Transactions 11:1242–1250.
  43. 43
    Effects of visual deprivation on morphology and physiology of cells in the cat's lateral geniculate body
    1. TN Wiesel
    2. DH Hubel
    (1963)
    Journal of Neurophysiology 26:978–993.
  44. 44
  45. 45
    Hippocampal Formation
    1. MP Witter
    2. DG Amaral
    (2004)
    In: G Paxinos, editors. The Rat Nervous System (3rd ed). San Diego, CA: Elsevier. pp. 635–704.
    https://doi.org/10.1016/B978-012547638-6/50022-5
  46. 46
    Nonnegative Sparse PCA
    1. R Zass
    2. A Shashua
    (2006)
    NIPS Proceedings pp. 1561–1568.
  47. 47
    Models of grid cell spatial firing published 2005-2011
    1. EA Zilli
    (2012)
    Frontiers in Neural Circuits, 6, 10.3389/fncir.2012.00016.

Decision letter

  1. Michael J Frank
    Reviewing Editor; Brown University, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your work entitled "Extracting grid characteristics from spatially distributed place cell inputs using non-negative PCA" for peer review at eLife. Your submission has been favorably evaluated by Timothy Behrens (Senior editor), Reviewing editor Michael Frank, and two reviewers.

The reviewers have discussed the reviews with one another and the Reviewing editor has drafted this decision to help you prepare a revised submission.

Summary:

The authors present analytical and neural network simulations which demonstrate that, under certain constraints, the first principal component of place cell firing patterns in a 2D arena resemble hexagonal grid cell firing patterns. Specifically, these constraints require that the principal components be constrained to non-negativity, and that the place cell inputs have a zero mean in time or space. In addition, they suggest that interesting features of grid cell firing, namely the tendency for grid scale to be discretely clustered, and for the grid to align at approximately 7 degrees to the walls of a square arena. This is an important result which adds to a growing body of recent studies that have emphasised the importance of the projection from place cells to grid cells, as opposed to the grid cell to place cell projection that has been the subject of several previous theoretical studies.

As you will see however, the Reviewers felt that it was important to provide some intuition and cogent explanation about how the results were obtained, and several questions arose along these lines. It is not clear how the dual constraints of (eigenvector) non-negativity and orthogonality are actually satisfied, nor is it clear why they get 7degree alignment or a ratio of 1.4 for the jump in grid scale.

In addition to the revisions articulated below, the following suggestion arose from discussion amongst the Reviewers/editors.

Figure 13 of the manuscript tracks the evolution of the network from different initial conditions.

One might go some way in clarifying the questions of how and why the observed phenomena occur by tracking individual objective measures, in a manner similar to Figure 13.

To wit, one of the algorithms you use trades off two objectives:

(i) maximize the variance of the projection of the data onto the [non-negative] principal component

(ii) enforce approximate orthogonormality by minimizing || I – eT e||, where I is the identity matrix, and e is the normalized "eigenvector," subject to the non-negativity constraint.

In addition, one ought to quantify the degree of non-negativity. There are several ways one might go about this, but one approach would be to normalize the eigenvector to unit length, then take the norm of the vector composed of only the positive entries and subtract the norm of the vector of negative entries.

Essential revisions:

1) This manuscript is of sufficient general interest for publication in eLife, although there is a general lack of detail when describing the methodology. There is also significant overlap with a recent paper at NIPS 2014 (Stachenfeld et al., 2014). I find both papers (Dordek et al., submitted, and Stachenfeld et al.) extremely interesting, and they both have different focusses. But I do think that there needs to be some link to/discussion of Stachenfeld et al.

2) Perhaps the most important issue is that they provide little intuition into the interesting subsidiary results. There is a section on "Why Hexagons?" but none on why other aspects of the results occur. This becomes more important given the fact that Stachenfeld also see the basic result. Why do grid scales have a ratio of 1.4 (indeed is this figure reliable)? Why do they see a non-zero alignment of grids to borders (and is this figure reliable)? Alignment angle appears to depend on place field size, but how does overall place cell firing coverage (and covariance) vary with place field size? How does the alignment angle covary with gridness (the very nice grid-like patterns they show tend to look aligned)?

3) "Due to the isotropic nature of the data (generated by the agent's motion), there was a 2D redundancy in the X-Y plane in conjunction with a dual phase redundancy in 1D, resulting in a fourfold total redundancy of every solution" – this statement requires unpacking. Presumably the dual phase redundancy is caused by the periodic boundary conditions? What are the implications for non-periodic boundary conditions? This sentence needs to be explained to the non-expert reader.

4) Similarly, Figure 4A is cryptic – how does the figure show fourfold redundancy? The authors also use a reduced number of place cells to obtain this result (225 vs. 625 in other simulations) – what are the reasons for this change? This should be explained and justified in the text.

5) More generally, when changing place field size, what happens to the number of place cells – presumably the density of firing overlap (covariance) needs to be maintained as this will strongly affect the results?

6) It is not clear how the non-negative weights are implemented alongside zero mean inputs, and we could not gain any additional insight from the Methods section. This issue needs to be elaborated upon in the text, otherwise it seems implicitly contradictory. Does it imply that the output GC can have negative firing rates? i.e. are there positive and negative place cell firing rates (with a mean of zero), coupled with positive place cell to grid cell synaptic weights, to give both negative and positive inputs to the output grid cell? Presumably the time step to time step changes in grid cell firing rate between negative and positive values will have a large impact on the direction of synaptic weight change – is this the case?

7) Figure 7C, D; also 14B – the analytical PCA results for constrained weights appear to exhibit a bimodal distribution of 90 degree gridness scores – what in the data accounts for this? Can the authors provide any comment or insight?

"The dependency was almost linear" – presumably there is a mathematical reason for this, based on the relationship between the covariance matrix and eigenvectors? Can the authors provide any comment or insight into this relationship?

8) The methods are generally lacking in detail that would allow these simulations to be replicated, based on details provided in the manuscript (see specific details below). The authors should ensure that sufficient detail to allow replication is incorporated. They may also like to consider uploading their code to ModelDB, or a similar online repository, for the sake of transparency and to allow independent replication.

9) The representation of space lies on a torus in this model, which means that the eigenvectors should be periodic and the spatial frequencies quantized. We fail to grasp why, in Figure 3, the first spatial frequency should be 2 and not 1 (i.e., we have two peaks in each direction), and why the eigenvectors 3B is periodic, but in 3A they are not.

10) The eigenvectors of a correlation matrix are orthogonal. If we normalize these eigenvectors and furthermore enforce an additional constraint that the entries of each eigenvector must be positive semidefinite, then the "synaptic weights" of different eigenvectors should not overlap. In other words, if an entry in the first eigenvector is non-zero, it ought to be zero in all other eigenvectors.

11) While the algorithms the authors use relax this last constraint, it is not clear how skewing or shearing a square pattern that results from PCA helps satisfy the dual constraints of non-negativity and orthogonality.

This is a key point in the paper, and begs for an intuitive explanation, or at least more analysis or quantification of the degree to which the constraints are satisfied.

12) The finding that the spacings corresponding to successive eigenvectors obey ~ 2 also calls for an explanation. Why should this ratio of spacings best satisfy non-negativity, orthogonality, and periodicity?

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Extracting grid characteristics from spatially distributed place cell inputs using non-negative PCA" for further consideration at eLife. Your revised article has been favorably evaluated by Timothy Behrens (Senior editor), Reviewing editor Michael Frank, and two reviewers, one of whom is Neil Burgess. The manuscript has been improved but both reviewers have identified some remaining issues that need to be addressed before acceptance, as articulated in the comments below:

Reviewer #1:

The authors have worked hard to provide further intuition and interpretation to their results, with much success, although there are still some points that should be clarified, see comments below.

There is also a plain error that should be removed. The authors state "the distribution of grid orientations becomes uniform in the range of 0-15 degrees, with a mean (and median) at 7.5 degrees, similar to experimental data [Stensola et al., 2015]". On the contrary the experimental data shows a very non-uniform distribution clearly peaked at about 7.5 degrees, which is quite different. Also, as a circular variable limited to [0,15) – albeit one for which 16 degrees corresponds to 14 degrees rather than 1 degree – the 'mean' and 'median' they refer to are undefined if the distribution is uniform. Thus the authors should acknowledge that they do not provide an explanation for this aspect of the experimental data, and some of the text/ results on this point could be removed, most importantly in the Abstract.

At the start of describing the main result (Results, first paragraph) is not clear that the "1d mapping of the 2d activity" is caused by the trajectory of the rat moving through the place fields.

The next part of the description of their PCA analysis as Eigenvectors and their spatial interpretation could also be improved in terms of intuition. It would be worth mentioning that the PCA, in addition to being calculated as the Eigenvectors of the matrix of temporal covariance between inputs, can also be thought of in terms of singular value decomposition of the (appropriately normalised) input activity. I.e. the input to the network (temporal sequences of place cell activity) can be decomposed into an outer product of (spatial) weights over place cell inputs (their Eigenvectors) and temporally varying weights over these Eigenvectors. The network learns the spatial weights over place cells (the Eigenvectors) as the connections weights from the place cells, and "projection onto place cell space" is simply the firing rates of the output neuron plotted against the location of the rat.

In the second paragraph of the subsection “Adding a non-Negativity constraint to the PCA “it is important to remind the reader that the "raw place cells activity" is not the final input to the PCA numerical methods – spatial or temporal mean normalization has to happen first.

It would be useful to have a surface plot showing Gridness vs. PC Density vs. PC Width. I have found it hard to replicate the Gridness distribution and wonder how general it is to the precise place cell inputs used? The reference to Figure 10A in the third paragraph of the subsection “Effect of place cell parameters on grid structure “is incorrect (I think) – not showing the dependence on the width of the place fields.

Regarding the presence of "modules" of grid cells (subsection “Modules of grid cells”), it is not clear if the interpretation is that different modules correspond to different grid scales, or different Eigenvectors (even if having the same scale). Assuming the latter, it seems that the analysis predicts that there should only be two modules, with the scale of the first module approximately 7.5x the largest place field inputs, and the second module (~2 smaller) existing only to enable non-negative PCA solutions. This deserves discussion, given that the data show at least 5 modules, with a range of scales, the smallest and most numerous having approximately the scale of the smaller place fields found in the dorsal hippocampus (25-30 cm).

Reviewer #2:

The new material, starting with the subsection “Steady state analysis “is great, and just the thing I was asking for, namely an understanding as to why hexagonal patterns result under a non-negative PCA learning rule.

There are a few things I'd like to understand better, though:

The last equation in the section “2D Solutions”: should read ϕ1 = ϕ2 = ϕ3n 2 π/3, n ϵ Z. e.g. ϕ1 = ϕ2 = ϕ3= 0 would be a solution.

As written, ϕ1 = ϕ2 = ϕ3, is not the maximum. Besides, this more general case would allow solutions that are unlike the ones observed in experiments.

This brings me to

In the subsection “2D solutions”:

"To do this [maximize the Fourier components inside the basis set,] we must maximize the minimal value of the cosine components for the basis set"

I am afraid this argument escapes me, even though its role is crucial.

While I think I understand the Fourier domain argument up to this point in the text, how exactly does maximizing the minimum value in real space of the function guarantee that the Fourier components in Eq. 24 are maximized?

Eq. 34: why is the upper limit here ∞, when it was M in Eq. 33?

I'd like some clarification from the authors before I can make a judgment call as to whether the argument is likely to be correct.

https://doi.org/10.7554/eLife.10094.027

Author response

As you will see however, the Reviewers felt that it was important to provide some intuition and cogent explanation about how the results were obtained, and several questions arose along these lines. It is not clear how the dual constraints of (eigenvector) non-negativity and orthogonality are actually satisfied, nor is it clear why they get 7degree alignment or a ratio of 1.4 for the jump in grid scale. Re-analysis demonstrates that we tend to get values which are on average a little lower than 7 degrees (Author response images 3 and 4). However, the range of possible values is from 0 to 15 degrees, and we demonstrate analytically (in the new extended Methods section) that in the case of a very large box, the distribution of orientation angles becomes uniform in the interval 0-15, with a mean at 7.5 degrees, consistent with experimental data.

The jump of 1.4 results from the change in the circular components of the solution in Fourier space (See especially new Figure 17 in Methods section of paper). From running more simulations and also from analytical considerations, we believe that within our theoretical framework there is only one jump, and thus this model can really only explain 2 modules and not more (see analytical derivation for more intuition on this point).

In addition to the new extended Methods section, we have updated the Abstract and Results (throughout) to account for these changes.

In addition to the revisions articulated below, the following suggestion arose from discussion amongst the Reviewers/editors. […]

In addition, one ought to quantify the degree of non-negativity. There are several ways one might go about this, but one approach would be to normalize the eigenvector to unit length, then take the norm of the vector composed of only the positive entries and subtract the norm of the vector of negative entries.

We thank the editors and reviewers for raising this essential issue. In our understanding, the main issue is to comprehend to what extent we can combine maximal -variance, orthogonality, and positivity together, when looking at multiple outputs (such as in old Figures 11 and 12, or new Figures 13 and 14). As we demonstrate now, the outputs from the method are not orthogonal (Author response image 1), and this we believe solves the conundrum. We now added a short comment about this into the Results section: “It is important to note though that, due to the non-negativity constraint, the vectors achieved in this way were not orthogonal, and thus it cannot be considered a real orthogonalization process, although as explained in the Methods section, the process does aim for maximum difference between the vectors”

Author response image 1
Demonstrates that output vectors are not orthogonal.

Matrix of cosine values between every two vectors.

https://doi.org/10.7554/eLife.10094.023

Essential revisions: 1) This manuscript is of sufficient general interest for publication in eLife, although there is a general lack of detail when describing the methodology.

Additional detail has been added to the Methods. Furthermore, extensive sections have been added to the Methods in order to provide an analytical motivation to the whole paper (from the subsection “Steady state analysis”). See detailed points below in the section related to Methods. Furthermore, we have uploaded our MATLAB code to a public server, as described below.

There is also significant overlap with a recent paper at NIPS 2014 (Stachenfeld et al., 2014). I find both papers (Dordek et al., submitted, and Stachenfeld et al) extremely interesting, and they both have different focusses. But I do think that there needs to be some link to/discussion of Stachenfeld et al.

Stachnfeld et al. arrive to similar ideas in a very nice paper based on a reinforcement learning model of place cells. Unlike our paper, their paper does not contain the positivity constraint, and thus the “grid cells” that results from the PCA procedure they use are square rather than hexagonal. Furthermore, they did not provide a theoretical analysis of the nature provided in the revision. We now discuss their paper in the Discussion, as follows: “We note that similar work has noticed the relation between place-cell-to-grid-cell transformation and PCA.[…] However, due to the unconstrained nature of their transformation, the outputted grid cells were square-like.”

2) Perhaps the most important issue is that they provide little intuition into the interesting subsidiary results. There is a section on "Why Hexagons?" but none on why other aspects of the results occur. This becomes more important given the fact that Stachenfeld also see the basic result. Why do grid scales have a ratio of 1.4 (indeed is this figure reliable)?

Following the new analytical derivation in the Methods section (from the subsection “Steady state analysis”), we can now explain the jump of scale as a hop on the Fourier-space lattice (new Figure 17 in paper).

We have now re-run the multiple solutions in order to demonstrate a jump between two consecutive scales. However, we find that it is harder to see the jump to the next smaller third scale, and thus we conclude that the method only partially accounts for the effect of modules. This is now noted in the Results as follows: “[…]we found that the ratio between the distances of the modules was ~1.4, close to the value of 1.42 found by Stensola et al. [Stensola et al., 2012]. Although we searched for additional such jumps, we could only identify this single jump, suggesting that our model can yield up to two “modules” and not more.”

Why do they see a non-zero alignment of grids to borders (and is this figure reliable)?

As we now show in our analytical derivation, the non-zero alignment can arise from the relation between the lattice points in the Fourier domain and the optimal-solution circle. This is now explainedthoroughly in the new sections added to the Methods from the subsection “Steady state analysis” and in Figure 1516 in paper.

Alignment angle appears to depend on place field size, but how does overall place cell firing coverage (and covariance) vary with place field size?

In our previous submission we manipulated the place field size while keeping the size of the arena and number of place cells constant. Thus the coverage has increased as a function of the increase in size. Here (Author response image 2) we show an example of increasing both the place field size and the size of the arena, while keeping the place-cell-size/arena size constant, such that the total coverage increases (a place cell on every pixel).

How does the alignment angle covary with gridness (the very nice grid-like patterns they show tend to look aligned)?

In Author response image 2b we plot grid orientation as a function of gridness. As can be seen the two quantities are not strongly related:

Author response image 2
Example of changing width of place cells (σ) and correspondingly changing arena size, such that there is a fixed ratio: arena_size/σ = 40.

(FISTA algorithm was used for positive PCA, with 0 boundary conditions.) (a) Grid scale as a function of σ (place cell width); (b) Grid orientation as a function of σ. (c) Histogram of grid orientations (d) Mean gridness as a function of σ. (e) Histogram of mean gridness.

https://doi.org/10.7554/eLife.10094.024
Author response image 3
Effect of changing the place-field size in a constant boundary (FISTA algorithm used; zero boundary conditions and Arena size 500);

(A) Grid scale as a function of place field size (σ); Linear fit is: Scale = 7.4. Σ+0.62; the lower bound is equal to , where k is defined in Eq. 32 in the Methods section; (B) Grid orientation as a function of gridness (green line is the mean); (C) Grid orientation as a function of σ – scatter plot (blue stars) and mean (green line); (D) Histogram of grid orientations; (E) Mean gridness as a function of σ; and (F) Histogram of mean gridness.

https://doi.org/10.7554/eLife.10094.025

3) "Due to the isotropic nature of the data (generated by the agent's motion), there was a 2D redundancy in the X-Y plane in conjunction with a dual phase redundancy in 1D, resulting in a fourfold total redundancy of every solution" – this statement requires unpacking. Presumably the dual phase redundancy is caused by the periodic boundary conditions?

We have now devoted a large section in the Methods to analytical considerations, which can explain the fourfold redundancy. Thus we changed the wording of the above paragraph to the following: “The solution demonstrated a fourfold redundancy (Figure 4C). This was apparent in the plotted eigenvalues (from largest to the smallest eigenvalue, Figure 4A and 4C), which demonstrated a fourfold grouping-pattern. The fourfold redundancy can be explained analytically by the symmetries of the system – see PCA analysis in Methods section.”

What are the implications for non-periodic boundary conditions? This sentence needs to be explained to the non-expert reader.

In the case of very large arenas, it does not matter what the boundary conditions are. In smaller arenas, the nature of the boundary conditions (i.e. how the place cells behave at the boundary) has an effect. It seems from our simulations that the boundary conditions matter less for the scale of the grid and more for the grid orientation (Compare Author response images 3 and 4). We now note in the Results section: “[…]for very large environments the effects of boundary conditions diminish.”.

Author response image 4
Effect of changing the place-field size in a constant boundary (FISTA algorithm used; periodic boundary conditions and Arena size 500);

(A) Grid scale as a function of place field size (σ); Linear fit is: Scale = 7.4 Σ+0.62; the lower bound, equal to 2π/k, were k is defined in Eq. 32 in the Methods section; (B) Grid orientation as a function of gridness; (C) Grid orientation as a function of σ – scatter plot (blue stars) and mean (green line); (D) Histogram of grid orientations; (E) Mean gridness as a function of σ; and (F) Histogram of mean gridness.

https://doi.org/10.7554/eLife.10094.026

Author response image 4 title/legend>

4) Similarly, Figure 4A is cryptic – how does the figure show fourfold redundancy? The authors also use a reduced number of place cells to obtain this result (225 vs. 625 in other simulations) – what are the reasons for this change? This should be explained and justified in the text. We have improved the figure by increasing the resolution and using the same number of points for both sub-figures. Now the fourfold redundancy is quite clear (Figure 4C). We updated Figure 4 to contain all 625 values, and have also added an additional panel to emphasize the fourfold redundancy.

5) More generally, when changing place field size, what happens to the number of place cells – presumably the density of firing overlap (covariance) needs to be maintained as this will strongly affect the results?

Author response image 2 above demonstrates that results are not strongly dependent on the place cell overlap.

6) It is not clear how the non-negative weights are implemented alongside zero mean inputs, and we could not gain any additional insight from the Methods section. This issue needs to be elaborated upon in the text, otherwise it seems implicitly contradictory. Does it imply that the output GC can have negative firing rates? i.e. are there positive and negative place cell firing rates (with a mean of zero), coupled with positive place cell to grid cell synaptic weights, to give both negative and positive inputs to the output grid cell?

The grid cells were not clipped to zero so they did obtain negative values. In order to realize this in real neurons we will need to assume a general baseline firing which will account for the fact that real neurons do not have negative firing rates.

Presumably the time step to time step changes in grid cell firing rate between negative and positive values will have a large impact on the direction of synaptic weight change – is this the case?

Yes, this is indeed the case, as far as we understand the comment.

7) Figure 7C, D; also 14B – the analytical PCA results for constrained weights appear to exhibit a bimodal distribution of 90 degree gridness scores – what in the data accounts for this? Can the authors provide any comment or insight?

The bi-modal distribution was a result of the formation of two populations, one with small orientations to one of the walls and another with small orientation to the other wall. We believe this is not an essential result, but cannot provide a full explanation for this phenomenon.

"The dependency was almost linear" – presumably there is a mathematical reason for this, based on the relationship between the covariance matrix and eigenvectors? Can the authors provide any comment or insight into this relationship?

This relation is now explained from the added analytical model. The predicted value provides a tight lower bound on the actual relation, as can be seen from the results of simulations in Author response image 3 and 4 above, where the actual values were that the grid scale is about 7.5 times the width of the place field.

Intuitively, the linear dependency between the place cell width and the grid cell spacing must be true in the limit of infinite box size, from dimensional analysis: the length units of the grid cell spacing must be the same length units of the place cell width, since these are the only length units in the model. This is explained mathematically in the unconstrained case by Equation 32, as we detail below that equation:

"Notice that if we multiply the place cell width by some positive constant c, then the solution 2pi/k† will be divided by c. […] When the box has finite size, k-lattice discretization also has a (usually small) effect on the grid spacing."

A similar reasoning can be applied to the constrained case.

8) The methods are generally lacking in detail that would allow these simulations to be replicated, based on details provided in the manuscript (see specific details below). The authors should ensure that sufficient detail to allow replication is incorporated. They may also like to consider uploading their code to ModelDB, or a similar online repository, for the sake of transparency and to allow independent replication.

We have followed the advice given and have now uploaded our code to t

This is also noted at the beginning of the Methods section.

9) The representation of space lies on a torus in this model, which means that the eigenvectors should be periodic and the spatial frequencies quantized. We fail to grasp why, in Figure 3, the first spatial frequency should be 2 and not 1 (i.e., we have two peaks in each direction),

Indeed the spatial frequencies are quantized by the box dimensions. However, the eigenvectors periodicity is not determined by the box, as we now explain mathematically in extended Methods section. Even in the case that the box size is infinite, we get a finite periodic solution, with spatial frequency equal to the peak of the Fourier transform of the place cell tuning curve – given in Equation 32. When the box size is finite this spatial frequency is also affected by the quantization due to the box dimensions (Figure 15).

and why the eigenvectors 3B is periodic, but in 3A they are not.

In Figure 3B we show multiple examples of the 1st principal component as output from the network, while in Figure 3A we show the first 16 components of the PCA algorithm. Thus both figures do not show the same thing. Generally, the PCA components need not appear periodic (i.e., repeat more than once in some directions), though they do obey periodic boundary conditions. For example, PCA components 5-12 in Figure 3A do not appear periodic since they are composed of the Fourier components B1-8 in Figure 15C.

10) The eigenvectors of a correlation matrix are orthogonal. If we normalize these eigenvectors and furthermore enforce an additional constraint that the entries of each eigenvector must be positive semidefinite, then the "synaptic weights" of different eigenvectors should not overlap. In other words, if an entry in the first eigenvector is non-zero, it ought to be zero in all other eigenvectors.

In the constrained case, the solutions maximizing the variance subject to the non-negativity constraint are no longer eigenvectors of any matrix. Thus, the second such solution described in the section “Hierarchical networks and modules” need not be orthogonal to the first, and the issues raised by the reviewers does not arise. We have added a note on this issue in the Results section: “It is important to note though that due to the non-negativity constraint the vectors achieved this way were not orthogonal, and thus it cannot be considered a real orthogonalization process.”. See Author response image 1, demonstrating that the different vectors outputted from the simulation are indeed not orthogonal.

11) While the algorithms the authors use relax this last constraint, it is not clear how skewing or shearing a square pattern that results from PCA helps satisfy the dual constraints of non-negativity and orthogonality. This is a key point in the paper, and begs for an intuitive explanation, or at least more analysis or quantification of the degree to which the constraints are satisfied.

We believe the above reply addresses this issue. See also Author response image 1.

12) The finding that the spacings corresponding to successive eigenvectors obey ~ 2 also calls for an explanation. Why should this ratio of spacings best satisfy non-negativity, orthogonality, and periodicity?

We now provide an analytical motivation for this in the new revised Methods section of the paper. In short, the change of scale occurs when “jumping” from the first optimal circle in the Fourier space to next available Fourier lattice points which can form a hexagonal shape (Figure 17 in paper).

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Reviewer #1: The authors have worked hard to provide further intuition and interpretation to their results, with much success, although there are still some points that should be clarified, see comments below. There is also a plain error that should be removed. The authors state "the distribution of grid orientations becomes uniform in the range of 0-15 degrees, with a mean (and median) at 7.5 degrees, similar to experimental data [Stensola et al., 2015]". On the contrary the experimental data shows a very non-uniform distribution clearly peaked at about 7.5 degrees, which is quite different. Also, as a circular variable limited to [0,15) – albeit one for which 16 degrees corresponds to 14 degrees rather than 1 degree – the 'mean' and 'median' they refer to are undefined if the distribution is uniform. Thus the authors should acknowledge that they do not provide an explanation for this aspect of the experimental data, and some of the text/ results on this point could be removed, most importantly in the Abstract.

We agree that the experimental data is non-uniform, and clearly peaked around 7.5o. We have thus erased this point from the abstract, and removed it throughout the paper.

At the start of describing the main result (Results, first paragraph) is not clear that the "1d mapping of the 2d activity" is caused by the trajectory of the rat moving through the place fields.

This point is now clarified in the second paragraph of subsection “Comparing Neural-Network results to PCA).

The next part of the description of their PCA analysis as Eigenvectors and their spatial interpretation could also be improved in terms of intuition. It would be worth mentioning that the PCA, in addition to being calculated as the Eigenvectors of the matrix of temporal covariance between inputs, can also be thought of in terms of singular value decomposition of the (appropriately normalised) input activity. I.e. the input to the network (temporal sequences of place cell activity) can be decomposed into an outer product of (spatial) weights over place cell inputs (their Eigenvectors) and temporally varying weights over these Eigenvectors. The network learns the spatial weights over place cells (the Eigenvectors) as the connections weights from the place cells, and "projection onto place cell space" is simply the firing rates of the output neuron plotted against the location of the rat.

Thanks for this good point. This alternative interpretation using SVD has now been added to the in the Methods section.

In the second paragraph of the subsection “Adding a non-Negativity constraint to the PCA “it is important to remind the reader that the "raw place cells activity" is not the final input to the PCA numerical methods – spatial or temporal mean normalization has to happen first.

This reminder has now been added in the third paragraph of the section “Adding a non-Negativity constraint to the PCA”.

It would be useful to have a surface plot showing Gridness vs. PC Density vs. PC Width. I have found it hard to replicate the Gridness distribution and wonder how general it is to the precise place cell inputs used?

We have now added the following subpanels to Figure 12, in which we vary the ratio between the box size and the size of the input place fields. As can be seen, as long as the box and place field sizes are large enough, the Gridness is > 1 (left panel). Furthermore, the scale of the grid is not highly dependent on the size of the box (right panel).

The reference to Figure 10A in the third paragraph of the subsection “Effect of place cell parameters on grid structure “is incorrect (I think) – not showing the dependence on the width of the place fields.

Typo has been corrected to “Figure 12A”.

Regarding the presence of "modules" of grid cells (subsection “Modules of grid cells”), it is not clear if the interpretation is that different modules correspond to different grid scales, or different Eigenvectors (even if having the same scale). Assuming the latter, it seems that the analysis predicts that there should only be two modules, with the scale of the first module approximately 7.5x the largest place field inputs, and the second module (~2 smaller)

existing only to enable non-negative PCA solutions. This deserves discussion, given that the data show at least 5 modules, with a range of scales, the smallest and most numerous having approximately the scale of the smaller place fields found in the dorsal hippocampus (25-30cm).

We agree. We believe that a thorough understanding of the module phenomenon will require further study. Mainly two points need to be added to a future model: (a) Looking at a non-uniform distribution of place-cell sizes, similar to reality and (b) adding the feedforward projection from grid cells to place cells, thus “closing the loop”. This point has now been added to the Discussion.

Reviewer #2: The last equation in the section “2D Solutions”: should read ϕ1 = ϕ2 = ϕ3n 2 π/3, n ϵ Z. e.g. ϕ1 = ϕ2 = ϕ3

= 0 would be a solution. As written, ϕ1 = ϕ2 = ϕ3, is not the maximum. Besides, this more general case would allow solutions that are unlike the ones observed in experiments.

This part was removed due to our revision, in response to the next comment. In any case, this was a typo (which happened due to an inaccurate change of notation from cos(k i ·(x-xo)) to cos(k i ·x+ϕi). It should have been “ϕi = ki ·x0, for some x0”.

This brings me to In the subsection “2D solutions”: "To do this [maximize the Fourier components inside the basis set,] we must maximize the minimal value of the cosine components for the basis set" I am afraid this argument escapes me, even though its role is crucial.

While I think I understand the Fourier domain argument up to this point in the text, how exactly does maximizing the minimum value in real space of the function guarantee that the Fourier components in Eq. 24 are maximized?

This argument was indeed incorrect (it is true if can we neglect the higher order harmonics in each lattice, unfortunately this cannot be done, as we verified numerically). We thank the reviewer for noting this. We changed the sub-section “2D solutions” (in section “The PCA solution with a non-negative constraint”) to correct this. First, we prove (in the appendix section “Upper bound on the constrained objective – 2D square lattice”) that the objective is bounded from above by 0.25 in any non-hexagonal 2D lattice with a single grid length (i.e., a square or 1D lattice). Note this remains true even if we consider higher order harmonics. Second, we give an example for a hexagonal lattice (found via a numerical scan) which achieves an objective value higher than 0.25. Thus, given our approximations, this shows that any optimal solution must be a hexagonal lattice

Eq. 34: why is the upper limit here ∞, when it was M in Eq. 33?

The two equations result from two different numerical procedures and assumptions. We revised the text near the two equations to explain this point:

“This form (Eq. 33) results from a parameter grid search for M =1, 2,3 and 4, under the assumption that L ⟶ ∞ and rˆ (k) is highly peaked. […] Specifically, it has similar form (Eq. 34).”

The main point is to show that both give very similar results. In other words, that our approximate description, using the simplifying assumptions and finite M, is closely related to the solution without approximations. Additional details were added to the text to clarify this important point. Also, to improve focus, we relegated to the appendix some less important details from that section

(on the origin of the difference between k and k*).

https://doi.org/10.7554/eLife.10094.028

Article and author information

Author details

  1. Yedidyah Dordek

    1. Faculty of Electrical Engineering, Technion – Israel Institute of Technology, Haifa, Israel
    2. Rappaport Faculty of Medicine and Research Institute, Technion – Israel Institute of Technology, Haifa, Israel
    Contribution
    YD, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Daniel Soudry
    Competing interests
    The authors declare that no competing interests exist.
  2. Daniel Soudry

    1. Department of Statistics, Columbia University, New York, United States
    2. Center for Theoretical Neuroscience, Columbia University, New York, United States
    Contribution
    DS, Conception and design, Analysis and interpretation of data, Drafting or revising the article, Analytical derivations in Materials and methods and in Appendix
    Contributed equally with
    Yedidyah Dordek
    For correspondence
    1. daniel.soudry@gmail.com
    Competing interests
    The authors declare that no competing interests exist.
  3. Ron Meir

    1. Faculty of Electrical Engineering, Technion – Israel Institute of Technology, Haifa, Israel
    Contribution
    RM, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  4. Dori Derdikman

    1. Rappaport Faculty of Medicine and Research Institute, Technion – Israel Institute of Technology, Haifa, Israel
    Contribution
    DD, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    1. derdik@technion.ac.il
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0003-3677-6321

Funding

Ollendroff center of the Department of Electrical Engineering, Technion (Research fund)

  • Yedidyah Dordek
  • Ron Meir

Gruss Lipper Charitable Foundation

  • Daniel Soudry

Intelligence Advanced Research Projects Activity

  • Daniel Soudry

Israel Science Foundation (Personal Research Grant, 955/13)

  • Dori Derdikman

Israel Science Foundation (New Faculty Equipment Grant, 1882/13)

  • Dori Derdikman

Rappaport Institute (Personal Research Grant)

  • Dori Derdikman

Allen and Jewel Prince Center for Neurodegenrative Disorders (Research Grant)

  • Dori Derdikman

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to thank Alexander Mathis and Omri Barak for comments on the manuscript, and Gilad Tocker and Matan Sela for helpful discussions and advice. The research was supported by the Israel Science Foundation grants 955/13 and 1882/13, by a Rappaport Institute grant, and by the Allen and Jewel Prince Center for Neurodegenerative Disorders of the Brain. The work of RM and Y D was partially supported by the Ollendorff Center of the Department of Electrical Engineering, Technion. DD is a David and Inez Myers Career Advancement Chair in Life Sciences fellow. The work of DS was partially supported by the Gruss Lipper Charitable Foundation, and by the Intelligence Advanced Research Projects Activity (IARPA) via Department of Interior/Interior Business Center (DoI/IBC) contract number D16PC00003. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. Disclaimer: The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of IARPA, DoI/IBC, or the U.S. Government.

Reviewing Editor

  1. Michael J Frank, Reviewing Editor, Brown University, United States

Publication history

  1. Received: July 15, 2015
  2. Accepted: March 8, 2016
  3. Accepted Manuscript published: March 8, 2016 (version 1)
  4. Version of Record published: April 13, 2016 (version 2)

Copyright

© 2016, Dordek et al

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,269
    Page views
  • 605
    Downloads
  • 3
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Comments

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Microbiology and Infectious Disease
    Kenneth J Ertel et al.
    Research Article