On Path Integration of Grid Cells: Group Representation and Isotropic Scaling
Year: 2021

Abstract:
Understanding how grid cells perform path integration calculations remains a fundamental problem. In this paper, we conduct theoretical analysis of a general representation model of path integration by grid cells, where the 2D self-position is encoded as a higher dimensional vector, and the 2D self-motion is represented by a general transformation of the vector. We identify two conditions on the transformation. One is a group representation condition that is necessary for path integration. The other is an isotropic scaling condition that ensures locally conformal embedding, so that the error in the vector representation translates conformally to the error in the 2D self-position. Then we investigate the simplest transformation, i.e., the linear transformation, uncover its explicit algebraic and geometric structure as matrix Lie group of rotation, and explore the connection between the isotropic scaling condition and a special class of hexagon grid patterns. Finally, with our optimization-based approach, we manage to learn hexagon grid patterns that share similar properties of the grid cells in the rodent brain. The learned model is capable of accurate long distance path integration. Code is available at https://github.com/ruiqigao/grid-cell-path.

heading: 1 Introduction
text: Imagine walking in the darkness. Purely based on the sense of self-motion, one can gain a sense of self-position by integrating the self motion - a process often referred to as path integration [10, 14, 21, 15, 27]. While the exact neural underpinning of path integration remains unclear, it has been hypothesized that the grid cells [21, 17, 40, 24, 23, 12] in the mammalian medial entorhinal cortex (mEC) may be involved in this process [20, 30, 22]. The grid cells are so named because individual neurons exhibit striking firing patterns that form hexagonal grids when the agent (such as a rat) navigates in a 2D open field [18, 21, 16, 6, 34, 5, 7, 11, 29, 1]. The grid cells also interact with the place cells in the hippocampus [28]. Unlike a grid cell that fires at the vertices of a lattice, a place cell often fires at a single (or a few) locations.
The purpose of this paper is to understand how the grid cells may perform path integration calculations. We study a general optimization-based representational model in which the 2D self-position is
∗The author is now a Research Scientist at Google Brain team.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
represented by a higher dimensional vector and the 2D self-motion is represented by a transformation of the vector. The vector representation can be considered position encoding or position embedding, where the elements of the vector may be interpreted as activities of a population of grid cells. The transformation can be realized by a recurrent network that acts on the vector. Our focus is to study the properties of the transformation.
Specifically, we identify two conditions for the transformation: a group representation condition and an isotropic scaling condition, under which we demonstrate that the local neighborhood around each self-position in the 2D physical space is embedded conformally as a 2D neighborhood around the vector representation of the self-position in the neural space.
We then investigate the simplest special case of the transformation, i.e., linear transformation, that forms a matrix Lie group of rotation, under which case we show that the isotropic scaling condition is connected to a special class of hexagonal grid patterns. Our numerical experiments demonstrate that our model learns clear hexagon grid patterns of multiple scales which share observed properties of the grid cells in the rodent brain, by optimizing a simple loss function. The learned model is also capable of accurate long distance path integration.
Contributions. Our work contributes to understanding the grid cells from the perspective of representation learning. We conduct theoretical analysis of (1) general transformation for path integration by identifying two key conditions and a local conformal embedding property, (2) linear transformation by revealing the algebraic and geometric structure and connecting the isotropic scaling condition and a special class of hexagon grid patterns, and (3) integration of linear transformation model and linear basis expansion model via unitary group representation theory. Experimentally we learn clear hexagon grid patterns that are consistent with biological observations, and the learned model is capable of accurate path integration.
2 General transformation
2.1 Position embedding
Consider an agent (e.g., a rat) navigating within a 2D open field. Let x = (x1,x2) be the selfposition of the agent. We assume that the selfposition x in the 2D physical space is represented by the response activities of a population of d neurons (e.g., d = 200), which form a vector v(x) = (vi(x), i = 1, ...,d)> in the ddimensional “neural space”, with each element vi(x) representing the firing rate of one neuron when the animal is at location x.
v(x) can be called position encoding or position embedding. Collectively, (v(x),∀x) forms a codebook of x ∈ R2, and (v(x),∀x) is a 2D manifold in the d-dimensional neural space, i.e., globally we embed R2 as a 2D manifold in the neural space. Locally, we identify two condi-
tions under which the 2D local neighborhood around each x is embedded conformally as a 2D neighborhood around v(x) with a scaling factor. See Fig. 1. As shown in Section 3.3, the conformal embedding is connected to the hexagon grid patterns.

heading: 2.2 Transformation and path integration
text: At self-position x, if the agent makes a self-motion ∆x = (∆x1,∆x2), then it moves to x+∆x. Correspondingly, the vector representation v(x) is transformed to v(x+∆x). The general form of the transformation can be formulated as:
v(x+∆x) = F(v(x),∆x). (1) The transformation F(·,∆x) can be considered a representation of ∆x, which forms a 2D additive group. We call Eq. (1) the transformation model. It can be implemented by a recurrent network to
derive a path integration model: if the agent starts from x0, and makes a sequence of moves (∆xt , t = 1, ...,T ), then the vector is updated by vt = F(vt−1,∆xt), where v0 = v(x0), and t = 1, ...,T .

heading: 2.3 Group representation condition
text: The solution to the transformation model (Eq. (1)) should satisfy the following condition. Condition 1. (Group representation condition) (v(x),∀x) and (F(·,∆x),∀∆x) form a representation of the 2D additive Euclidean group R2 in the sense that
F(v(x),0) = v(x), ∀x; (2) F(v(x),∆x1 +∆x2) = F(F(v(x),∆x1),∆x2), ∀x,∆x1,∆x2. (3)
(F(·,∆x),∀∆x) is a Lie group of transformations acting on the codebook manifold (v(x),∀x). The reason for (2) is that if ∆x = 0, then F(·,0) should be the identity transformation. Thus the codebook manifold (v(x),∀x) consists of fixed points of the transformation F(·,0). If F(·,0) is furthermore a contraction around (v(x),∀x), then (v(x),∀x) are the attractor points. The reason for (3) is that the agent can move in one step by ∆x1 +∆x2, or first move by ∆x1, and then move by ∆x2. Both paths would end up at the same x+∆x1 +∆x2, which is represented by the same v(x+∆x1 +∆x2).
The group representation condition is a necessary self-consistent condition for the transformation model (Eq. (1)).

heading: 2.4 Egocentric self-motion
text: Self-motion ∆x can also be parametrized egocentrically as (∆r,θ), where ∆r is the displacement along the direction θ ∈ [0,2π], so that ∆x= (∆x1 = ∆r cosθ ,∆x2 = ∆r sinθ). The egocentric self-motion may be more biologically plausible where θ is encoded by head direction, and ∆r can be interpreted as the speed along direction θ . The transformation model then becomes
v(x+∆x) = F(v(x),∆r,θ), (4) where we continue to use F(·) for the transformation (with slight abuse of notation). (∆r,θ) form a polar coordinate system around x.

heading: 2.5 Infinitesimal self-motion and directional derivative
text: In this subsection, we derive the transformation model for infinitesimal self-motion. While we use ∆x or ∆r to denote finite (non-infinitesimal) self-motion, we use δx or δ r to denote infinitesimal self-motion. At self-position x, for an infinitesimal displacement δ r along direction θ , δx= (δx1 = δ r cosθ ,δx2 = δ r sinθ). See Fig. 1 (a) for an illustration. Given that δ r is infinitesimal, for any fixed θ , a first order Taylor expansion of F(v(x),δ r,θ) with respect to δ r gives us
v(x+δx) = F(v(x),δ r,θ) = F(v(x),0,θ)+F ′(v(x),0,θ)δ r+o(δ r) = v(x)+ fθ (v(x))δ r+o(δ r), (5)
where F(v(x),0,θ) = v(x) according to Condition 1, and fθ (v(x)) := F ′(v(x),0,θ) is the first derivative of F(v(x),∆r,θ) with respect to ∆r at ∆r = 0. fθ (v(x)) is the directional derivative of F(·) at self-position x and direction θ . For a fixed θ , (F(·,∆r,θ),∀∆r) forms a one-parameter Lie group of transformations, and fθ (·) is the generator of its Lie algebra.

heading: 2.6 Isotropic scaling condition
text: With the directional derivative, we define the second condition as follows, which leads to locally conformal embedding and is connected to hexagon grid pattern. Condition 2. (Isotropic scaling condition) For any fixed x, ‖ fθ (v(x))‖ is constant over θ .
Let f0(v(x)) denote fθ (v(x)) for θ = 0, and fπ/2(v(x)) denote fθ (v(x)) for θ = π/2. Then we have the following theorem:
Theorem 1. Assume group representation condition 1 and isotropic scaling condition 2. At any fixed x, for the local motion δx= (δ r cosθ ,δ r sinθ) around x, let δv = v(x+δx)−v(x) be the change of vector and s = ‖ fθ (v(x))‖, then we have ‖δv‖= s‖δx‖. Moreover,
δv = fθ (v(x))δ r+o(δ r) = f0(v(x))δ r cosθ + fπ/2(v(x))δ r sinθ +o(δ r), (6) where f0(v(x)) and fπ/2(v(x)) are two orthogonal basis vectors of equal norm s.
See Supplementary for a proof and Fig. 1(b) for an illustration. Theorem 1 indicates that the local 2D polar system around self-position x in the 2D physical space is embedded conformally as a 2D polar system around vector v(x) in the d-dimensional neural space, with a scaling factor s (our analysis is local for any fixed x, and s may depend on x). Conformal embedding is a generalization of isometric embedding, where the metric can be changed by a scaling factor s. If s is globally constant for all x, then the intrinsic geometry of the codebook manifold (v(x),∀x) remains Euclidean, i.e., flat. Why isotropic scaling and conformal embedding? The neurons are intrinsically noisy. During path integration, the errors may accumulate in v. Moreover, when inferring self-position from visual image, it is possible that v is inferred first with error, and then x is decoded from the inferred v. Due to isotropic scaling and conformal embedding, locally we have ‖δv‖= s‖δx‖, which guarantees that the `2 error in v translates proportionally to the `2 error in x, so that there will not be adversarial perturbations in v(x) that cause excessively big errors in x. Specifically, we have the following theorem. Theorem 2. Assume the general transformation model (Eq. (4)) and the isotropic scaling condition. For any fixed x, let s = ‖ fθ (v(x))‖, which is independent of θ . Suppose the neurons are noisy: v = v(x)+ ε , where ε ∼N (0,τ2Id) and d is the dimensionality of v. Suppose the agent infers its 2D position x̂ from v by x̂ = argminx′ ‖v−v(x′)‖2, i.e., v(x̂) is the projection of v onto the 2D manifold formed by (v(x′),∀x′). Then we have
E‖x̂−x‖2 = 2τ2/s2. (7)
See Supplementary for a proof.
Connection to continuous attractor neural network (CANN) defined on 2D torus. The group representation condition and the isotropic scaling condition appear to be satisfied by the CANN models [2, 6, 7, 29, 1] that are typically hand-designed on a 2D torus. See Supplementary for details.

heading: 3 Linear transformation
text: After studying the general transformation, we now investigate the linear transformation of v(x), for the following reasons. (1) It is the simplest transformation for which we can derive explicit algebraic and geometric results. (2) It enables us to connect the isotropic scaling condition to a special class of hexagon grid patterns. (3) In Section 4, we integrate it with the basis expansion model, which is also linear in v(x), via unitary group representation theory.
For finite (non-infinitesimal) self-motion, the linear transformation model is: v(x+∆x) = F(v(x),∆x) =M(∆x)v(x), (8)
where M(∆x) is a matrix. The group representation condition becomes M(∆x1 +∆x2)v(x) = M(∆x2)M(∆x1)v(x), i.e., M(∆x) is a matrix representation of self-motion ∆x, and M(∆x) acts on the coding manifold (v(x),∀x)). For egocentric parametrization of self-motion (∆r,θ), we can further write M(∆x) =Mθ (∆r) for ∆x = (∆r cosθ ,∆r sinθ), and the linear model becomes v(x+∆x) = F(v(x),∆r,θ) =Mθ (∆r)v(x).

heading: 3.1 Algebraic structure: matrix Lie algebra and Lie group
text: For the linear model (Eq. (8)), the directional derivative is: fθ (v(x)) = F ′(v(x),0,θ) = M ′θ (0)v(x) =B(θ)v(x), where B(θ) =M ′ θ (0), which is the derivative of Mθ (∆r) with respect to ∆r at 0. For infinitesimal self-motion, the transformation model in Eq. (5) becomes v(x+δx) = (I+B(θ)δ r)v(x)+o(δ r), (9)
where I is the identity matrix. It can be considered a linear recurrent network where B(θ) is the learnable weight matrix. We have the following theorem for the algebraic structure of the linear transformation.
Theorem 3. Assume the linear transformation model so that for infinitesimal self-motion (δ r,θ), the model is in the form of Eq. (9), then for finite displacement ∆r,
v(x+∆x) =Mθ (∆r)v(x) = exp(B(θ)∆r)v(x). (10)
Proof. We can divide ∆r into N steps, so that δ r = ∆r/N→ 0 as N→ ∞, and
v(x+∆x) = (I+B(θ)(∆r/N)+o(1/N))Nv(x)→ exp(B(θ)∆r)v(x) (11)
as N→ ∞. The matrix exponential map is defined by exp(A) = ∑∞n=0 An/n!.
The above math underlies the relationship between matrix Lie algebra and matrix Lie group in general [38]. For a fixed θ , the set of Mθ (∆r) = exp(B(θ)∆r) for ∆r ∈ R forms a matrix Lie group, which is both a group and a manifold. The tangent space of Mθ (∆r) at identity I is called matrix Lie algebra. B(θ) is the basis of this tangent space, and is often referred to as the generator matrix.
Path integration. If the agent starts from x0, and make a sequence of moves ((∆rt ,θt), t = 1, ...,T ), then the vector representation of self-position is updated by
vt = exp(B(θt)∆rt)vt−1, (12)
where v0 = v(x0), and t = 1, ...,T .
Approximation to exponential map. For a finite but small ∆r, exp(B(θ)∆r) can be approximated by a second-order (or higher-order) Taylor expansion
exp(B(θ)∆r) = I+B(θ)∆r+B(θ)2∆r2/2+o(∆r2). (13)

heading: 3.2 Geometric structure: rotation, periodicity, metic and error correction
text: If we assume B(θ) = −B(θ)>, i.e., skew-symmetric, then I +B(θ)δ r in Eq. (9) is a rotation matrix operating on v(x), due to the fact that (I+B(θ)δ r)(I+B(θ)δ r)> = I+O(δ r2). For finite ∆r, exp(B(θ)∆r) is also a rotation matrix, as it equals to the product of N matrices I+B(θ)(∆r/N) (Eq. (11)). The geometric interpretation is that, if the agent moves along the direction θ in the physical space, the vector v(x) is rotated by the matrix B(θ) in the neural space, while the `2 norm ‖v(x)‖2 remains fixed. We may interpret ‖v(x)‖2 = ∑di=1 vi(x)2 as the total energy of grid cells. See Fig. 1(b).
The angle of rotation is given by ‖B(θ)v(x)‖δ r/‖v(x)‖, because ‖B(θ)v(x)‖δ r is the arc length and ‖v(x)‖ is the radius. If we further assume the isotropic scaling condition, which becomes that ‖ fθ (v(x))‖= ‖B(θ)v(x)‖ is constant over θ for the linear model, then the angle of rotation can be written as µδ r, where µ = ‖B(θ)v(x)‖/‖v(x)‖ is independent of θ . Geometrically, µ tells us how fast the vector rotates in the neural space as the agent moves in the physical space. In practice, µ can be much bigger than 1 for the learned model, thus the vector can rotate back to itself in a short distance, causing the periodic patterns in the elements of v(x). µ captures the notion of metric.
For µ 1, the conformal embedding in Fig. 1 (b) magnifies the local motion in Fig. 1 (a), and this enables error correction [34]. More specifically, we have the following result, which is based on Theorem 2.
Proposition 1. Assume the linear transformation model (Eq. (9)) and the isotropic scaling condition 2. For any fixed x, let µ = ‖B(θ)v(x)‖/‖v(x)‖. Suppose v= v(x)+ε , where ε ∼N (0,τ2Id) and τ2 = α2(‖v(x)‖2/d), so that α2 measures the variance of noise relative to the average magnitude of (vi(x)2, i= 1, ...,d). Suppose the agent infers its 2D position x̂ from v by x̂= argminx′ ‖v−v(x′)‖2. Then we have
E‖x̂−x‖2 = 2α2/(µ2d). (14)
See Supplementary for a proof. By the above proposition, error correction of grid cells is due to two factors: (1) higher dimensionality d of v(x) for encoding 2D positions x, and (2) a magnifying µ 1 (our analysis is local for any fixed x, and µ may depend on x).

heading: 3.3 Hexagon grid patterns formed by mixing Fourier waves
text: In this subsection, we make connection between the isotropic scaling condition 2 and a special class of hexagon grid patterns created by linearly mixing three Fourier plane waves whose directions are 2π/3 apart. We show such linear mixing satisfies the linear transformation model and the isotropic scaling condition. Theorem 4. Let e(x) = (exp(i〈a j,x〉), j = 1,2,3)>, where (a j, j = 1,2,3) are three 2D vectors of equal norm, and the angle between every pair of them is 2π/3. Let v(x) =Ue(x), where U is an arbitrary unitary matrix. Let B(θ) =U ∗D(θ)U , where D(θ) = diag(i〈a j,q(θ)〉, j = 1,2,3), with q(θ) = (cosθ ,sinθ)>. Then (v(x),B(θ)) satisfies the linear transformation model (Eq. (9)) and the isotropic scaling condition 2. Moreover, B(θ) is skew-symmetric.
See Supplementary for a proof. We would like to emphasize that the above theorem analyzes a special case solution to our linear transformation model, but our optimization-based learning method does not assume any superposition of Fourier basis functions as in the theorem. Our experimental results are learned purely by optimizing a loss function based on the simple assumptions of our model with generic vectors and matrices.
We leave it to future work to theoretically prove that the isotropic scaling condition leads to hexagon grid patterns in either the general transformation model or the linear transformation model. The hexagon grid patterns are not limited to superpositions of three plane waves as in the above theorem.

heading: 3.4 Modules
text: Biologically, it is well established that grid cells are organized in discrete modules [4, 37] or blocks. We thus partition the vector v(x) into K blocks, v(x) = (vk(x),k = 1, ...,K). Correspondingly the generator matrices B(θ) = diag(Bk(θ),k = 1, ...,K) are block diagonal, so that each sub-vector vk(x) is rotated by a sub-matrix Bk(θ). For the general transformation model, each sub-vector is transformed by a separate sub-network. By the same argument as in Section 3.2, let µk = ‖Bkvk(x)‖/‖vk(x)‖, then µk is the metric of module k.

heading: 4 Interaction with place cells

heading: 4.1 Place cells
text: For each v(x), we need to uniquely decode x globally. This can be accomplished via interaction with place cells. Specifically, each place cell fires when the agent is at a specific position. Let A(x,x′) be the response map of the place cell associated with position x′. It measures the adjacency between x and x′. A commonly used form of A(x,x′) is the Gaussian adjacency kernel A(x,x′) = exp(−‖x−x′‖2/(2σ2)). The set of Gaussian adjacency kernels serve as inputs to our optimizationbased method to learn grid cells.
4.2 Basis expansion
A popular model that connects place cells and grid cells is the following basis expansion model (or PCA-based model) [13]:
A(x,x′) = 〈v(x),u(x′)〉= d
∑ i=1 ui,x′vi(x), (15)
where v(x) = (vi(x), i = 1, ...,d)>, and u(x′) = (ui,x′ , i = 1, ...,d)>. Here (vi(x), i = 1, ...,d) forms a set of d basis functions (which are functions of x) for expanding A(x,x′) (which is a function of x for each place x′), while u(x′) is the read-out weight vector for place cell at x′ and needs to be learned. See Fig. 2 for an illustration. Experimental results on biological brains have shown that the connections from grid cells to place cells are excitatory [42, 31]. We thus assume that ui,x′ ≥ 0 for all i and x′.

heading: 4.3 From group representation to basis functions
text: The vector representation v(x) generated (or constrained) by the linear transformation model (Eq. (8)) can serve as basis functions of the PCA-based basis expansion model (Eq. (15)), due to the fundamental theorems of Schur [41] and Peter-Weyl [38], which reveal the deep root of Fourier analysis and generalize it to general Lie groups. Specifically, if M(∆x) is an irreducible unitary representation of ∆x that forms a compact Lie group, then the elements {Mi j(∆x)} form a set of orthogonal basis functions of ∆x. Let v(x) = M(x)v(0) (where we choose the origin 0 as the reference point). The elements of v(x), i.e., (vi(x), i = 1, ...,d), are linear mixings of the basis functions {Mi j(x)}, so that they themselves form a new set of basis functions that serve to expand (A(x,x′),∀x′) that parametrizes the place cells. Thus group representation in our path integration model is a perfect match to the basis expansion model, in the sense that the basis functions are results of group representation.
The basis expansion model (or PCA-based model) (Eq. (15)) assumes that the basis functions are orthogonal, whereas in our work, we do not make the orthogonality assumption. Interestingly, the learned transformation model generates basis functions that are close to being orthogonal automatically. See Supplementary for more detailed explanation and experimental results.

heading: 4.4 Decoding and re-encoding
text: For a neural response vector v, such as vt in Eq. (12), the response of the place cell associated with location x′ is 〈v,u(x′)〉. We can decode the position x̂ by examining which place cell has the maximal response, i.e.,
x̂= argmax x′ 〈v,u(x′)〉. (16)
After decoding x̂, we can re-encode v← v(x̂) for error correction. Decoding and re-encoding can also be done by directly projecting v onto the manifold (v(x),∀x), which gives similar results. See Supplementary for more analysis and experimental results.

heading: 5 Learning
text: We learn the model by optimizing a loss function defined based on three model assumptions discussed above: (1) the basis expansion model (Eq. (15)), (2) the linear transformation model (Eq. (10)) and (3) the isotropic scaling condition 2. The input is the set of adjacency kernels A(x,x′),∀x,x′. The unknown parameters to be learned are (1) (v(x) = (vk(x),k = 1, ...,K),∀x), (2) (u(x′),∀x′) and (3) (B(θ),∀θ). We assume that there are K modules or blocks and B(θ) is skew-symmetric, so that B(θ) are parametrized as block-diagonal matrices (Bk(θ),k = 1, ...,K),∀θ) and only the lower triangle parts of the matrices need to be learned. The loss function is defined as a weighted sum of simple `2 loss terms constraining the three model assumptions: L = L0 +λ1L1 +λ2L2, where
L0 = Ex,x′ [A(x,x′)−〈v(x),u(x′)〉]2, (basis expansion) (17)
L1 = K
∑ k=1 Ex,∆x‖vk(x+∆x)− exp(Bk(θ)∆r)vk(x)‖2, (transformation) (18)
L2 = K
∑ k=1 Ex,θ ,∆θ [‖Bk(θ +∆θ)vk(x)‖−‖Bk(θ)vk(x)‖]2. (isotropic scaling) (19)
In L1, ∆x = (∆r cosθ ,∆r sinθ). λ1 and λ2 are chosen so that the three loss terms are of similar magnitudes. A(x,x′) are given as Gaussian adjacency kernels. For regularization, we add a penalty on ‖u(x′)‖2, and further assume u(x′)≥ 0 so that the connections from grid cells to place cells are excitatory [42, 31]. However, note that u(x′)≥ 0 is not necessary for the emergence of hexagon grid patterns as shown in the ablation studies.
Expectations in L0, L1 and L2 are approximated by Monte Carlo samples. L is minimized by Adam [25] optimizer. See Supplementary for implementation details.
It is worth noting that, consistent with the experimental observations, we assume individual place field A(x,x′) to exhibit a Gaussian shape, rather than a Mexican-hat pattern (with balanced excitatory center and inhibitory surround) as assumed in previous basis expansion models [13, 33] of grid cells.
ReLU non-linearity. We also experiment with a non-linear transformation model where a ReLU activation is added. See Supplementary for details.

heading: 6 Experiments
text: We conduct numerical experiments to learn the representations as described in Section 5. Specifically, we use a square environment with size 1m × 1m, which is discretized into a 40× 40 lattice. For direction, we discretize the circle [0,2π] into 144 directions and use nearest neighbor linear interpolations for values in between. We use the second-order Taylor expansion (Eq. (13)) to approximate the exponential map exp(B(θ)∆r). The displacement ∆r are sampled within a small range, i.e., ∆r is smaller than 3 grids on the lattice. For A(x,x′), we use a Gaussian adjacency kernel with σ = 0.07. v(x) is of d = 192 dimensions, which is partitioned into K = 16 modules, each of which has 12 cells.

heading: 6.1 Hexagon grid patterns
text: Fig. 3 shows the learned firing patterns of v(x) = (vi(x), i = 1, ...,d) over the 40×40 lattice of x. Every row shows the learned units belonging to the same block or module. Regular hexagon grid patterns emerge. Within each block or module, the scales and orientations are roughly the same, but with different phases or spatial shifts. For the learned B(θ), each element shows regular sine/cosine tuning over θ . See Supplementary for more learned patterns.
We further investigate the characteristics of the learned firing patterns of v(x) using measures adopted from the literature of grid cells. Specifically, the hexagonal regularity, scale and orientation of grid-like patterns are quantified using the gridness score, grid scale and grid orientation [26, 32], which are determined by taking a circular sample of the autocorrelogram of the response map. Table 1 summarizes the results of gridness scores and comparisons with other optimization-based approaches [3, 33]. We apply the same threshold to determine whether a learned neuron can be considered a grid cell as in [3] (i.e., gridness score > 0.37). For our model, 73.10% of the learned neurons exhibit significant hexagonal periodicity in terms of the gridness score. Fig. 4 shows the
histogram of grid scales of the learned grid cell neurons (mean 0.33, range 0.21 to 0.49), which follows a multi-modal distribution. The ratio between neighboring modes are roughly 1.52 and 1.51, which closely matches the theoretical predictions [39, 36] and also the empirical results from rodent grid cells [37]. Collectively, these results reveal striking, quantitative correspondence between the properties of our model neurons and those of the grid cells in the brain.
Connection to continuous attractor neural network (CANN) defined on 2D torus. The fact that the learned response maps of each module are shifted versions of a common hexagon periodic pattern implies that the learned codebook manifold forms a 2D torus, and as the agent moves, the responses of the grid cells undergo a cyclic permutation. This is consistent with the CANN models hand-crafted on 2D torus. See Supplementary for a detailed discussion.
Ablation studies. We conduct ablation studies to examine whether certain model assumptions are empirically important for the emergence of hexagon grid patterns. The conclusions are highlighted as follows: (1) The loss term L2 (Eq. (19)) constraining the isotropic scaling condition is necessary for learning hexagon grid patterns. (2) The constraint u(x′)≥ 0 is not necessary for learning hexagon patterns, but the activations can be either excitatory or inhibitory without the constraint. (3) The skew-symmetric assumption on B(θ) is not important for learning hexagon grid pattern. (4) Hexagon patterns always emerge regardless of the choice of block size and number of blocks. (5) Multiple blocks or modules are necessary for the emergence of hexagon grid patterns of multiple scales. See Fig. 5 for several learned patterns and Supplementary for the full studies.

heading: 6.2 Path integration
text: We then examine the ability of the learned model on performing multi-step path integration, which can be accomplished by recurrently updating vt (Eq. (12)) and decoding vt to xt for t = 1, ...,T (Eq. (16)). Re-encoding vt ← v(xt) after decoding is adopted. Fig. 6(a) shows an example trajectory of accurate path integration for number of time steps T = 30. As shown in Fig. 6(b), with re-encoding, the path integration error remains close to zero over a duration of 500 time steps (< 0.01 cm, averaged over 1,000 episodes), even if the model is trained with the single-time-step transformation model (Eq. (18)). Without re-encoding, the error goes slight higher but still remains small (ranging from 0.0 to 4.2 cm, mean 1.9 cm in the 1m × 1m environment). Fig. 6(c) summarizes the path integration performance by fixing the number of blocks and altering the block size. The performance of path integration would be improved as the block size becomes larger, i.e., with more neurons in each module. When block size is larger than 16, path integration is very accurate for the time steps tested.
Error correction. See Supplementary for numerical experiments on error correction, which show that the learn model is still capable of path integration when we apply Gaussian white noise errors or Bernoulli drop-out errors to vt .

heading: 6.3 Additional experiments on path planning and egocentric vision
text: We also conduct additional experiments on path planning and egocentric vision with our model. Path planning can be accomplished by steepest ascent on the adjacency to the target position. For egocentric vision, we learn an extra generator network that generates the visual image given the position encoding formed by the grid cells. See Supplementary for details.

heading: 7 Related work
text: Our work is related to several lines of previous research on modeling grid cells. First, RNN models have been used to model grid cells and path integration. The traditional approach uses simulation-based models with hand-crafted connectivity, known as continuous attractor neural network (CANN) [2, 6, 7, 29, 1]. On the other hand, more recently two pioneering papers [9, 3] developed optimization-based RNN approaches to learn the path integration model and discovered that grid-like response patterns can emerge in the optimized networks. These results are further substantiated in [33, 8]. Our work analyzes the properties of the general recurrent model for path integration, and these properties seem to be satisfied by the hand-crafted CANN models. Our method belongs to the scheme of optimization-based approaches, and the learned response maps share similar properties as assumed by the CANN models.
Second, our work differs from the PCA-based basis expansion models [13, 33, 35] in that, unlike PCA, we make no assumption about the orthogonality between the basis functions, and the basis functions are generated by the transformation model. Furthermore, in previous basis expansion models [13, 33], place fields with Mexican-hat patterns (with balanced excitatory center and inhibitory surround) had to be assumed in order to obtain hexagonal grid firing patterns. However, experimentally measured place fields in biological brains were instead well characterized by Gaussian functions. Crucially, in our model, hexagonal grids emerge from learning with Gaussian place fields, and there is no need to assume any additional surround mechanisms or difference of Gaussians kernels.
In another related paper, [19] proposed matrix representation of 2D self-motion, while our work analyzes general transformations. Our investigation of the special case of linear transformation model reveals the matrix Lie group and the matrix Lie algebra of rotation group. Our work also connects the linear transformation model to the basis expansion model via unitary group representation theory.

heading: 8 Conclusion
text: This paper analyzes the recurrent model for path integration calculations by grid cells. We identify a group representation condition and an isotropic scaling condition that give rise to locally conformal embedding of the self-motion. We study a linear prototype model that reveals the matrix Lie group of rotation, and explore the connection between the isotropic scaling condition and hexagon grid patterns. In addition to these theoretical investigations, our numerical experiments demonstrate that our model can learn hexagon grid patterns for the response maps of grid cells, and the learned model is capable of accurate long distance path integration.
In this work, the numerical experiments are mostly limited to the linear transformation model, with the exception of an experiment with ReLU non-linearity. We will conduct experiments on the other non-linear transformation models, especially the forms assumed by the hand-crafted continuous attractor neural networks. Moreover, we assume that the agent navigates within a square open-field environment without obstacles or rewards. It is worthwhile to explore more complicated environments, including 3D environment.
Acknowledgments and Disclosure of Funding
The work was supported by NSF DMS-2015577, ONR MURI project N00014-16-1-2007, DARPA XAI project N66001-17-2-4029, and XSEDE grant ASC170063. We thank Yaxuan Zhu from UCLA Department of Statistics for his help with experiments on egocentric vision. We thank Dr. Wenhao Zhang for sharing his knowledge and insights on continuous attractor neural networks. We thank Sirui Xie for discussions. We thank the three reviewers for their constructive comments.

Review summaries:
This paper focuses on the find that how the grid cells perform path integration. Speciafically, this paper introduces the isotropic scaling condition to produce the hexagon firing pattern for grid cells. The optimization experiments show that the proposed formulation of place and grid cells demonstrated this.

The paper derives a model of grid cells from first principles, starting with the basic assumption that position is represented as a high-D vector and that it is updated as a function of direction and speed.  The authors stipulate an isotropic scaling condition in order to make the representation robust to neural noise.  The elements of the resulting high-D embedding resemble grid cells, and importantly, the ability for error correction during path integration. 

The paper introduces a general framework for modeling 2D navigation via grid cells. The main idea is to embed a 2D position x into a higher-dimensional positional encoding v(x) and model 2D translations via group transformations in this embedding space. A specific case are linear transformations via rotation matrices that act on the embedding v(x). Under this assumption, it can be shown that hexagonal patterns which are characteristic for neural grid cells naturally emerge as an optimal encoding.