\chapter{Learning Grid Cells by Predictive Coding}\label{chap:grid}

\section{Introduction}
Our brain contains a rich set of neural representations of space that help us navigate in an ever-changing world. These include hippocampal place cells \citep{o1976place}, which fire when an animal is at a specific spatial position, and grid cells observed in the medial entorhinal cortex (MEC) \citep{hafting2005microstructure}, which fire when an animal occupies multiple positions on a hexagonal or triangular grid. Grid cells have been observed across various species \citep{fyhn2008grid,yartsev2011grid,doeller2010evidence}, and their remarkable regularity has raised extensive interest in the computational mechanism underlying their emergence. Earlier models have focused on how mechanisms, such as membrane potential oscillation \citep{o2005dual,hasselmo2007grid} and specialized recurrent connectivity, can generate grid-like firing patterns \citep{fuhs2006spin,burak2009accurate}. More recently, research has shown that grid cells can emerge in recurrent neural networks (RNNs) trained using backpropagation through time (BPTT) for path integration tasks. The models are trained to predict their current location by integrating velocity inputs \citep{cueva2018emergence,banino2018vector,whittington2020tolman,sorscher2023unified}, providing a normative, task-driven account of the computational problem that the MEC grid cells address. However, the process by which the MEC circuit acquires, or \emph{learns} the grid cells in a biologically plausible way has been largely neglected, despite the fact that grid cells are known to be learned, rather than hardwired at birth \citep{langston2010development,wills2010development}. Although models for synaptic plasticity underlying grid cells do exist (e.g. \citet{weber2018learning}), they are highly specialized for learning grid cells, and do not capture of dynamic nature of the learning process through navigation.

In this chapter, we directly tackle the learning problem underlying the emergence of grid cells using predictive coding. Our approach to modeling grid cell emergence through predictive coding is motivated by three key factors: Firstly, the predictive coding algorithm can be implemented in predictive coding networks (PCNs) with local computations and Hebbian plasticity \citep{bogacz2017tutorial}, making it more biologically plausible than learning rules such as backpropagation. Secondly, PCNs have been successful in replicating representations in other regions of the brain, such as the visual cortex \citep{rao1999predictive,olshausen1996emergence,millidge2024predictive}. Thirdly, PCNs have demonstrated the ability to perform hippocampus-related functions, such as associative and sequential memories (as discussed in previous chapters and in \citet{salvatori2021associative}, \citet{tang2023recurrent} and \citet{tang2024sequential}).

The primary contribution of this chapter is to demonstrate for the first time that grid cells naturally emerge in PCNs trained to represent spatial inputs with biologically plausible plasticity rules. In this chapter we:

\begin{itemize}[leftmargin=*]
    \item show that hexagonal grid cells develop as the latent representations of place cells in classical PCNs \citep{rao1999predictive,olshausen1996emergence} with sparse and non-negative constraints;
    \item train temporal predictive coding network (tPCN) (Chapter \ref{chap:time} and \citet{millidge2024predictive}) in path integration tasks and observe that the latent activities of the tPCN develop hexagonal, grid-like representations, similar to what has been discovered in RNNs;
    \item develop an understanding of grid cell emergence in tPCN, by showing analytically that the Hebbian learning rule of tPCN implicitly approximates truncated BPTT \citep{williams1990efficient};
    \item show that tPCN can robustly develop grid cells under different architectural choices, and even without velocity inputs in path integration.
\end{itemize}

Overall, our results present an effective and plausible learning rule for hexagonal grid cells in the MEC based on predictive coding. We offer a novel extension of predictive coding theory, which has traditionally been used to model visual representations \citep{rao1999predictive, olshausen1996emergence}, to encompass spatial representations in the MEC. Our findings therefore offer a novel understanding of how a single, unified learning algorithm can be employed by different brain regions to represent inputs of various levels of abstraction.

\section{Background: Computational Models of Grid Cells} 
The periodicity of grid cells inspired early models of grid cells based on membrane potential oscillations, where the periodic firing of grid cells results naturally from the interference between somatic and dendritic oscillators in MEC pyramidal neurons \citep{o2005dual, hasselmo2007grid}. These models were subsequently extended to incorporate multiple networks of oscillatory neurons \citep{zilli2010coupled}. However, these models lack biological plausibility as they require an unrealistically large number of networks \citep{giocomo2011computational}. Another major family of models leverages the recurrent attractor networks and obtains grid firing patterns \citep{fuhs2006spin, burak2009accurate, ocko2018emergent} by hand-tuning the recurrent connectivity to form a center-surround structure. These networks perform robust and accurate path integration \citep{burak2009accurate} and can explain experimental observations such as the deformation of grid cells in irregular environments \citep{ocko2018emergent}. However, as pointed out by \citet{sorscher2023unified}, these models lack an explanation for the underlying spatial task that gives rise to the specific recurrent connectivity.

To address this gap, recent studies have explored the question \emph{`If grid cell is the answer, what is the question?'}. \citet{dordek2016extracting} showed that grid cells emerge as the non-negative principal components of place cells, while \citet{stachenfeld2017hippocampus} proposed that grid cells form a basis for predicting future observations. Other studies have focused on the multi-modularity of grid cells by optimizing biologically constrained objective functions \citep{dorrell2022actionable,schaeffer2024self}. Notably, multiple research tracks have found that RNNs trained to perform path integration tasks will develop hexagonal grid representations in their latent states \citep{cueva2018emergence,banino2018vector,whittington2020tolman}, suggesting that grid cells emerge as a result of successful navigation. These findings were further reinforced by \citet{sorscher2023unified}, who analytically demonstrated that path integration with certain implementation choices, such as non-negativity, is a sufficient condition for the emergence of grid cells, clarifying earlier controversies \citep{schaeffer2022no}. However, none of these works have addressed how the MEC/hippocampal network learns the grid cells. The RNN models are trained by BPTT, a learning rule unlikely to be employed by the brain \citep{lillicrap2019backpropagation}. Even though the principal component model by \citet{dordek2016extracting} can be learned with the plausible Sanger's rule \citep{sanger1989optimal}, it has been shown that principal component analysis (PCA) cannot be applied to other brain regions such as the visual areas \citep{olshausen1996emergence}, and Sanger's rule cannot be generalized to dynamical tasks such as path integration. Earlier models of the learning process of grid cells have explored plausible learning rules such as spike time-dependent plasticity \citep{widloski2014model} and variants of Hebbian learning rules \citep{kropff2008emergence} within networks of excitatory and inhibitory neurons \citep{weber2018learning}. However, these learning rules are highly specialized, and have not been shown to reproduce representations from other brain regions with non-spatial tasks. Recent works have also modeled the hippocampal formation using generative models with plausible learning rules similar to predictive coding \citep{george2024generative,bredenberg2021impression}, though these studies did not address 2D spatial learning.

% \subsection{Predictive Coding} Predictive coding has been an influential theory in understanding cortical computations \citep{friston2005theory,rao1999predictive,bogacz2017tutorial} and has been applied to modeling various cortical functions (see \citet{millidge2021predictive} for a review). Specifically, in the visual cortex, PCNs develop realistic visual representations such as Gabor-like receptive fields in response to both static \citep{rao1999predictive,olshausen1996emergence} and moving stimuli \citep{millidge2024predictive}. Recently, theories have been developed to describe hippocampo-neocortical interactions using predictive coding \citep{barron2020prediction}, and PCNs have demonstrated the ability to memorize and retrieve static and dynamic visual patterns, a key function of the hippocampus \citep{salvatori2021associative, tang2023recurrent, tang2024sequential}. Our work explores whether the representational learning capabilities of predictive coding can be extended to the hippocampal formation, which has so far only been functionally modeled by PCNs.

% The computations of PCNs use only local neural dynamics and Hebbian plasticity, making it biologically more plausible than backpropagation \citep{whittington2017approximation}. It has also been shown that predictive coding approximates backpropagation both in theory and practice \citep{whittington2017approximation, song2024inferring, pinchetti2024benchmarking}. Unlike many other Hebbian learning rules, predictive coding can be extended to temporal predictive coding networks (tPCNs), which use recurrent connections to process dynamic stimuli \citep{millidge2024predictive}. However, while \citet{millidge2024predictive} demonstrated that tPCNs approximate Kalman filtering, the relationships between tPCNs and RNNs remain unclear. In this chapter, we train tPCNs for path integration and compare their performance with RNNs both analytically and experimentally in this context.

\section{Models}

\begin{figure}[t]
\begin{center}
\includegraphics[width=\textwidth]{figs/chap6/grid-cell-circuit.pdf}
\end{center}
\vspace{-3mm}
\caption{\textbf{Architecture and circuit implementation of PCNs.} A: Two ways of obtaining grid cells within PCNs. Left: Sparse, non-negative PCN as a generative model. Right: tPCN performing path integration. B: Circuit implementation of path-integrating tPCN with a mapping to MEC and hippocampus.}
\label{fig:networks}
\end{figure}

\subsection{Non-negative Sparse PCN}
We first investigate the classical PCN \citep{rao1999predictive} for its ability to form grid representations. Assuming a place cell input $\bfp \in \mathbb{R}^{N_p}$ that represents a location in $2$D space as an $N_p$-dimensional vector, a simple $2$-layer PCN generates predictions of $\bfp$ using its latent activities $\bfg \in \mathbb{R}^{N_g}$ (which will develop grid-like representations) and a weight matrix $\bfW$ (Fig~\ref{fig:networks}A). The generative model minimizes the following loss function subject to two constraints:
\begin{equation}\label{eq:pcn}
    \mathcal{F}_{\textnormal{PCN}} =  \Vert \bfp - \bfW\bfg \Vert_2^2 + \Vert \bfg \Vert_2^2 + 2\lambda \Vert \bfg \Vert_1
\end{equation}
where $\Vert \bfg \Vert_2^2$ constrains the $l$2 norm of the latent $\bfg$ and $\lambda \Vert \bfg \Vert_1$ enforces sparsity, similar to the sparse coding model \citep{olshausen1996emergence}. This loss function is minimized via an expectation-maximization (EM) algorithm, alternating between the optimization over $\bfg$ (inference) and $\bfW$ (learning) (see Appendix \ref{appendix:alg} for the training algorithm):
\begin{equation}\label{eq:pc-inf}
    \Delta \bfg \propto -\frac{\partial\mathcal{F}_{\textnormal{PCN}}}{\partial\bfg}  = - \bfg - \lambda\textnormal{sgn}(\bfg) + \bfW ^ {\top} \bm{\epsilon}^{\bfp}; \quad \bfg \leftarrow \texttt{ReLU}(\bfg + \Delta\bfg) 
\end{equation}
\begin{equation}\label{eq:pc-learn}
    \Delta \bfW \propto -\frac{\partial\mathcal{F}_{\textnormal{PCN}}}{\partial\bfW} =  \bm{\epsilon}^{\bfp} \bfg ^ {\top}
\end{equation}
where $\bm{\epsilon}^{\bfp}:= \bfp - \bfW\bfg$ and we apply a $\texttt{ReLU}$ to the inference dynamics to constrain the latent activities to be non-negative. The inference and learning dynamics can be implemented in a plausible circuit \citep{bogacz2017tutorial}. After convergence, we examine the firing fields of the latent activities $\bfg$.

\subsection{Path Integrating tPCN}
To account for the learning of spatial representations in moving animals, we also investigate tPCN that extends the classical PCNs to the temporal domain \citep{millidge2024predictive,tang2024sequential} in path integration tasks (Fig.~\ref{fig:networks}B). The model receives a series of place cell activities $\bfp_1, ..., \bfp_T$ and velocity inputs $\bfv_1, ..., \bfv_T$ that represent the trajectory of an agent moving in a $2$D space, and minimizes the following loss function at each time step $t$\footnote{Note that the position of the nonlinear functions $h$ and $f$ is different from previous chapters. In this chapter, they are applied after the activities have been linearly weighted. This is to ensure the non-negativity of the latent, grid activities via $h=\texttt{ReLU}$. The latent and observed variables are also denoted by $\bfg$ and $\bfp$ in this Chapter to reflect grid and place cells.}:
\begin{equation}\label{eq:loss-tpc}
    \mathcal{F}_{\textnormal{tPCN}, t} = \Vert \bfp_t - f(\bfF \bfg_t) \Vert_2^2 + \Vert \bfg_t - h(\bfW \hat{\bfg}_{t-1} + \bfB \bfv_t) \Vert_2^2
\end{equation}
where $f$ and $h$ are both nonlinear activation functions, and $\bfB$, $\bfW$ and $\bfF$ are weight matrices projecting the predictions. We define $\bm{\epsilon}^{\bfp}_t:= \bfp_t - f(\bfF \bfg_t)$, $\bm{\epsilon}^{\bfg}_t:=\bfg_t - h(\bfW \hat{\bfg}_{t-1} + \bfB \bfv_t)$. The model learns by first optimizing the loss function with respect to $\bfg_t$ via gradient descent:
\begin{equation}\label{eq:tpc_inf}
    \Delta \bfg_t \propto -\frac{\partial \mathcal{F}_{\textnormal{tPCN}, t}}{\partial \bfg_t} = -\bm{\epsilon}^{\bfg}_t + \bfF^{\top} f'(\bfF \bfg_t)  \bm{\epsilon}^{\bfp}_t
\end{equation}
and then optimizing weights by:
\begin{equation}\label{eq:tpc_learn}
\begin{split}
    \{\Delta \bfF, \Delta \bfW, \Delta \bfB\} 
    & \propto \frac{\partial \mathcal{F}_{\textnormal{tPCN}, t}}{\partial {\{\bfF, \bfW, \bfB\}}}  
    \\ &= \{
        f'(\bfF \bfg_t) \bm{\epsilon}^{\bfp}_t \bfg_t^{\top}, 
        h'(\tilde{\bfg}_t) \bm{\epsilon}^{\bfg}_t \hat{\bfg}_{t-1}^{\top}, 
        h'(\tilde{\bfg}_t) \bm{\epsilon}^{\bfg}_t \bfv_t^{\top}\}
\end{split}
\end{equation}
where $f'$ and $h'$ are Jacobians of the nonlinear functions $f$ and $h$, and $\tilde{\bfg}_t := \bfW \hat{\bfg}_{t-1} + \bfB \bfv_t$. After the inference (Eq.~\ref{eq:tpc_inf}) converges, we set $\hat{\bfg}_t$ to the converged value of $\bfg_t$, which will be used for optimizing the objective function at the next time step i.e., $\mathcal{F}_{\textnormal{tPCN}, t+1}$. The model is trained on a large number of trajectories $\{\bfv_t, \bfp_t\}$ and after training, a set of velocity inputs from unseen trajectories is presented to the model. The model then performs a forward pass through time and layers to predict the positions encoded by place cells (see Appendix \ref{appendix:alg} for the training and testing algorithms of tPCN):
\begin{equation}\label{eq:tpc-forward}
    \bfg_t = h(\bfW \bfg_{t-1} + \bfB \bfv_t), \quad \hat{\bfp}_t = f(\bfF\bfg_t)
\end{equation}
The model is evaluated on 1) the accuracy of path integration position prediction $\hat{\bfp}_t$ and 2) the firing fields of the latent $\bfg$. When both $f$ and $h$ are linear, these computations can be plausibly implemented in a neural circuit shown in Fig.~\ref{fig:networks}C, with local inference computations (Eq.~\ref{eq:tpc_inf}) and Hebbian learning rules (Eq.~\ref{eq:tpc_learn}) (see Chapters \ref{chap:time} and \ref{chap:seqmem} for detailed discussions of this implementation). When the activation functions involve only local nonlinearity, such as \texttt{tanh} or \texttt{ReLU}, the Jacobians are diagonal and the inference and learning rules remain local and Hebbian \citep{millidge2022predictive}, and additional circuitry components can be included to plausibly implement the nonlinearities \citep{whittington2017approximation}. Within the context of spatial representation learning, this circuit implementation can be naturally mapped to the circuitry of the hippocampal formation, the the latent layer represents the MEC and the observation layer represents the hippocampus.

\subsection{Input of the Model}
In models discussed in this chapter, we assume that grid cells are inferred as latent representations of place cells. Although previous models have followed the opposite direction of the relationship, several strands of experimental evidence have suggested the emergence of grid cells as a result of place cells, including the earlier development of place cells \citep{bush2014grid,langston2010development,wills2010development}. In both PCN and tPCN models, the place cell inputs are constructed as $2$D difference-of-softmaxed-Gaussian (DoS) curves flattened into $1$D vectors, which have been shown to yield hexagonal grid representations in RNNs \citep{schaeffer2022no,sorscher2023unified}. The firing centers of the place cells are uniformly distributed across a $2$D environment. For PCN, the inputs are $N_x$ evenly distributed locations in the environment ($N_x$ large enough to cover the whole environment) represented by the $N_p$ place cells. For tPCN, the trajectories for the path integration task are obtained by simulating an agent performing a smooth random walk in the square environment. At each point in time, the $N_p$ place cells will be uniquely activated, representing the agent's current location. The velocity inputs $\bfv_t$ are $2$D vectors representing the speed of the simulated agent on the $x$ and $y$ coordinates at time step $t$. The effect of boundaries is simulated by slowing down the agent and reverting its moving direction near the borders of the environment. We sample a large number of trajectories to cover the whole simulated environment for training.

\section{Results}
\subsection{Sparse non-negative PCN develops latent grid cells}

\begin{figure}[t]
\begin{center}
\includegraphics[width=\textwidth]{figs/chap6/pcn-pca.pdf}
\end{center}
\caption{\textbf{Grid cells developed in PCN.} A: Latent representations of a sparse, non-negative PCN, resembling hexagonal grid cells in the MEC. Numbers in the title reflect the grid scores. B: Grid cells obtained via the pattern formation theory/non-negative PCA discussed in \citet{sorscher2023unified,dordek2016extracting}. C, D: Latent representations without sparsity or non-negativity, respectively. E: Distribution of grid scores of the representations in A and B.}
\label{fig:pcn-pca}
\end{figure}

Here we examine whether the sparse non-negative PCN can develop hexagonal, grid-like latent representations of the space after training, by plotting each latent neuron's responses to the $N_x=900$ locations in the $2$D space. We use $N_p=512$ and $N_g=256$. The ``gridness" of the 2D latent representations is evaluated using the \emph{grid score} metric, commonly employed in both experimental and computational studies \citep{sargolini2006conjunctive,banino2018vector} (see Appendix \ref{appendix:exp} for grid score calculation).  We found that this simple, $2$-layer PCN can develop hexagonal grid cells similar to those observed in the MEC (Fig.~\ref{fig:pcn-pca}A). For comparison, we reproduce the results from \cite{dordek2016extracting} and \cite{sorscher2023unified}, which show theoretically that performing non-negative PCA on the place cell inputs is guaranteed to produce hexagonal grid representations as the principal components of the $N_x \times N_p$ place cell input matrix. The visual results of the reproduction are shown in Fig.~\ref{fig:pcn-pca}B, and we compare the distribution of grid scores of the PCN's latent neuron firing fields with those of the non-negative principal components in Fig.~\ref{fig:pcn-pca}E. The grid scores between our sparse non-negative PCN and non-negative PCA are similarly distributed.

Why does the sparse, non-negative PCN develop hexagonal grid cells? While a precise analytical explanation is left for future work, we offer an intuitive hypothesis here. When presented with a batch $N_x$ of place cell inputs, the objective of PCN (Eq.~\ref{eq:pcn}) can be written compactly as:
\begin{equation}
    \mathcal{F}_{\textnormal{PCN}} = \Vert \bfP -  \bfG \bfW^{\top}\Vert_F^2 + \sum_{i=1}^{N_x} \Vert \bfg_i \Vert_2^2 + 2\lambda\Vert \bfg_i \Vert_1
\end{equation}
where $\bfP\in\mathbb{R}^{N_x\times N_p}$ is the place cell activities across $N_x$ locations, and $\bfG\in\mathbb{R}^{N_x\times N_g}$ represents grid cell responses. On the other hand, the objective function of PCA is (see Appendix \ref{app:pca} for how this squared error minimization formulation of PCA is related to the more well-known maximum variance formulation):
\begin{equation}\label{eq:PCA}
    \mathcal{F}_{\textnormal{PCA}} = \Vert \bfP -  \bfG \mathbf{M}\Vert_F^2 \quad \textnormal{s.t.} \quad \bfG^{\top} \bfG = \bfI_{N_g}
\end{equation}
where $\mathbf{M}$ is the $N_g \times N_p$ readout matrix\footnote{Note that although in PCA, $\bfG$ is usually framed as the matrix that projects $\bfP$ into its lower ($N_g$)-dimensional \textit{column} representations $\mathbf{M}$, itself can also be understood as lower ($N_g$)-dimensional \textit{row} representations of $\bfP$.}. The constraint $\bfG^{\top} \bfG = \bfI_{N_g}$ in Eq.~\ref{eq:PCA} enforces orthonormality of the grid cell matrix $\bfG$ columns, meaning they are orthogonal and have unit norm. We hypothesize that the constraint $\Vert \bfg_i \Vert_2^2 + 2\lambda\Vert \bfg_i \Vert_1$ for our sparse PCN achieves this orthonormality implicitly: while the constraints are imposed on the rows of $\bfG$, the overall sparsity of entries in $\bfG$ could induce orthogonality among its columns, with the $l_2$ term constraining the norm of the columns to achieve normality implicitly. Indeed, Fig.~\ref{fig:pcn-pca}C shows that if we remove the sparsity constraint, the latent neurons' firing fields will no longer be hexagonal. Similarly, without $\texttt{ReLU}$ i.e., non-negativity applied to the inference dynamics, we also could not obtain hexagonal grid cells (Fig.~\ref{fig:pcn-pca}D). It is worth noting that although (non-negative) PCA can be learned with the biologically plausible Sanger's rule \citep{sanger1989optimal}, it lacks PCN's generalizability to different architectures \citep{salvatori2022learning} and to other brain regions such as the visual cortex \citep{olshausen1996emergence,rao1999predictive}. 
% In fact, it has been precisely shown that PCA \emph{cannot} develop brain-like representations that PCNs are capable of generating \citep{olshausen1996emergence}.
However, it can be noticed that the grid cells by PCN lack the multi-modularity of the grid cells by non-negative PCA i.e., grid cells with different firing periods. 
% We suspect that although sparse PCNs can approximate the orthonormality of latent variables, they cannot approximate the process by which PCA extracts the importance of latent variables, which naturally corresponds to the multi-modularity of grid cells.
We suspect that although sparse PCNs can approximate the orthonormality of latent variables, they lack PCA's ability to extract latent variables ordered by the amount of explained variance in data, with higher variance naturally corresponding to larger spatial scales and vice versa. 

\subsection{tPCN develops grid cells by path integration}

\begin{figure}[t]
\begin{center}
\includegraphics[width=\textwidth]{figs/chap6/tpc_rnn.pdf}
\end{center}

\caption{\textbf{tPCN in path integration}. A: Visual demonstration of the performance of tPCN and RNN in path integration. B: RMSEs between the decoded and ground-truth 2D positions by tPCN and RNN with different agent moving speed. C: Grid score distributions of tPCN and RNN with different agent moving speed. D, E: Firing fields of latent neurons in a tPCN and an RNN respectively, when $dt=0.02$. F, G: Firing fields of latent neurons in a tPCN and an RNN respectively, when $dt=0.1$.}

\label{fig:tpc-rnn}
\end{figure}

Although training a static PCN with a large number of place cell activations can already give rise to brain-like hexagonal grid cells,
% it is unlikely that animals develop grid cells by receiving place cell representations of the entire environment altogether. 
the emergence of grid cells is known to rely on dynamic motion of animals \citep{mcnaughton2006path,winter2015passive}. Therefore, we investigate tPCN in a path integration task, where the simulated agent uses dynamic velocity inputs to determine its current position. As a reference, we compare tPCN with RNNs trained in path integration, which have been shown to develop hexagonal grid cells \citep{cueva2018emergence,banino2018vector,sorscher2023unified} and share the same graphical structure as tPCN (Fig.~\ref{fig:networks}B). However, it is important to note that RNNs are trained with the biologically implausible backpropagation-though-time (BPTT) algorithm, which requires ``unrolling'' of the network through time, a process unlikely to occur in the brain \citep{lillicrap2019backpropagation}.

We first evaluate whether tPCN can learn to perform the path integration task using local and Hebbian learning rules. We trained a tPCN model with $N_g=2048$ latent neurons on trajectories within a $1.4\textnormal{m}\times 1.4$m environment represented by $N_p=512$ place cells. After training, we tested the model on a set of unseen trajectories with velocity input $\bfv_t$, and assessed whether the tPCN and RNN models could predict the correct positions using Eq.~\ref{eq:tpc-forward}. As the output of the networks is the $N_p$-dimensional population activity of the place cells, we calculate the predicted $2$D positions by averaging the center positions of the $3$ most active place cells in the output $\hat{\bfp}_t$, and calculate the root mean square error (RMSE) between the decoded and ground-truth $2$D positions. The visual and numerical results are shown in Fig.~\ref{fig:tpc-rnn}A and B, where we also varied a scaling factor $dt$ of the simulated agent's speed, sampled from a Rayleigh distribution with mean $1$, to test the robustness of the results. Note that we do not intend to model physiologically realistic speed of animals with these values. The performance of tPCN is comparable to that of the RNN, though it slightly deteriorates when the agent moves at higher speeds.

Next, we examine whether the tPCN model develops grid-like representations in its latent layer during path integration. We plot the firing fields of the $2048$ latent neurons given an unseen set of trajectories covering the entire space. The neurons with the highest grid scores are shown in  Fig.~\ref{fig:tpc-rnn}C, which reveals a grid-like, hexagonal firing pattern with high grid scores. Visually, these grid cells are similar to those in a trained RNN with the same architecture shown in Fig.~\ref{fig:tpc-rnn}E, replicating the results from \citep{sorscher2023unified}. To systematically compare the grid cells in tPCN and RNN, we plot the distribution of grid scores in both models as a function of the movement speed of the agent in the environment in Fig.~\ref{fig:tpc-rnn}C. When the movement is slow, the grid score distributions are similar between tPCN and RNN. However, as the $dt$ increases to $0.05$ and $0.1$, tPCN tends to have higher grid scores than RNN. This is visually reflected in Fig.~\ref{fig:tpc-rnn}F (tPCN) and G (RNN), which shows the latent representations developed by tPCN largely retain the grid-like pattern whereas firing centers of many of the RNN neurons no longer form a grid when $dt=0.1$. Interestingly, the band-like representations present in both models in this case are observed in MEC \citep{krupic2012neural}, although their existence is controversial \citep{navratilova2016grids}.


\subsection{tPCN approximates truncated BPTT}

\begin{figure}[t]
\begin{center}
\includegraphics[width=\textwidth]{figs/chap6/tbptt.pdf}
\end{center}
\vspace{-3mm}
\caption{\textbf{Comparing tPCN and tBPTT.} A: Dependencies of latent grid cells in tPCN and RNN trained with $1$-step tBPTT. Black arrows indicate the flow of computations during a forward pass and red arrows indicate the dependency of latent variables. B: Firing fields of the latent neurons of an RNN trained by $1$-step tBPTT. C, D: Path integration RMSE and grid score distributions of $1$-step tBPTT, BPTT and tPCNs with different inference iterations. ``tPCNk" indicates tPCN trained with k inference iterations.}
\label{fig:tbptt}
\end{figure}

Next, we asked why hexagonal grid representations emerge both when training a tPCN using a BPTT-free Hebbian learning rule and when training an RNN using BPTT. We provide an analytical comparison between the learning rules of tPCN and RNN. Assuming a vanilla, sequence-to-sequence RNN with exactly the same graphical structure as in Fig.~\ref{fig:networks}A, its dynamics can be recursively described as:
\begin{equation}\label{eq:rnn_forward}
    \bfg_t = h(\bfW \bfg_{t-1} + \bfB \bfv_t); \quad \hat{\bfp}_t =  f(\bfF \bfg_t)
\end{equation}
The loss that this RNN is trained to minimize is the cumulative prediction error:
\begin{equation}
    \mathcal{F}_{\textnormal{RNN}} = \textstyle{\sum_{t=1}^T \mathcal{F}_{\textnormal{RNN}, t} = \sum_{t=1}^T \Vert \bfp_t - \hat{\bfp}_t \Vert_2^2}
\end{equation}

Suppose BPTT is performed at every step $t$ to update weights in this RNN, the learning rule for $\bfW$ at step $t$ can be expressed as (see Appendix \ref{appendix:bptt} for derivations):
\begin{equation}\label{eq:bptt-learning-rule}
    \Delta \bfW^{\textnormal{RNN}} = \textstyle{\sum_{k=1}^t \tfrac{\partial \bfg_t}{\partial \bfg_k} h'(\tilde{\bfg}_t) \bfF^{\top} f'(\bfF\bfg_t) \bm{\epsilon}_t^{\bfp} \underline{\bfg_{k-1}^{\top}}}
\end{equation}
where $\bm{\epsilon}_t^{\bfp}$ denotes the prediction error $\bfp_t - \hat{\bfp}_t$ and the $\frac{\partial \bfg_t}{\partial \bfg_k}$ terms correspond to the unrolling in BPTT, which can be factorized into a chain of partial derivatives \citep{bellec2020solution}. On the other hand, for tPCN, if we assume that the inference dynamics in Eq.~\ref{eq:tpc_inf} has fully converged ($\Delta \bfg_t=0$) at the time of weight update, the learning rule of tPCN can be written as (see Appendix \ref{appendix:bptt} for derivations):
\begin{equation}\label{eq:tpc-learning-rule}
    \Delta \bfW^{\textnormal{tPCN}} = h'(\tilde{\bfg}_t) \bfF^{\top} f'(\bfF\bfg_t) \bm{\epsilon}_t^{\bfp} \underline{\hat{\bfg}_{t-1}^{\top}}
\end{equation}
Two key differences between these learning rules stand out. First, tPCN does not involve the recursive unrolling term, thereby avoiding the need to maintain a perfect memory of all preceding hidden states. Second, instead of using the forward-propagated $\bfg_{t-1}$ as in Eq.~\ref{eq:rnn_forward}, tPCN employs the inferred $\hat{\bfg}_{t-1}$ from Eq.~\ref{eq:tpc_inf} (underlined). The first difference suggests an equivalence between tPCN and RNN trained with truncated BPTT (tBPTT) with a truncation window of size 1 (1-step tBPTT) \citep{williams1990efficient}, where the RNN does not backpropagate \emph{any} hidden states through time when updating the weights. This characteristic could potentially harm the RNN's performance as it cannot effectively perform temporal credit assignment. However, the second difference partially solves this problem, as $\hat{\bfg}_{t-1}$ is inferred following Eq.~\ref{eq:tpc_inf}, which includes the term $\bm{\epsilon}^{\bfp}_{t-1}$ that communicates the place cell prediction error at step $t-1$. Therefore, when $\bfW$ is updated at step $t$, the $\hat{\bfg}_{t-1}^{\top}$ term in $\Delta \bfW^{\textnormal{tPCN}}$ will effectively form an eligibility trace \citep{bellec2020solution} that allows the model to access historical prediction errors on the place cell level. Fig.~\ref{fig:tbptt}A illustrates this difference between tPCN and RNN trained by $1$-step tBPTT, highlighting the dependency of tPCN hidden states on past place cell activations. In Appendix \ref{appendix:bptt} we also discuss the relationship between the update rules for $\bfB$ and $\bfF$ in these two models.

To verify this theoretical difference, we compare tPCN with RNNs trained by tBPTT in the path integration task. Since $\hat{\bfg}_t$ in tPCN is initialized by a forward pass $f(\bfW\hat{\bfg}_{t-1})$ and then updated by the iterative inference (Appendix \ref{appendix:alg}), the behavior of 1-step tBPTT, which computes its latent states via a forward pass at each time step, should be closer to tPCN with fewer inference iterations. Therefore, we evaluate tPCN with various inference iterations. Fig.~\ref{fig:tbptt}B shows the grid cells learned by an RNN trained with $1$-step tBPTT, which still exhibit hexagonal grid firing fields, though with lower grid scores than those from full BPTT. This suggests that backpropagating the error through all time steps is not entirely necessary for RNNs to generate grid cell-like representations. In Fig.~\ref{fig:tbptt}C we show the path integration performance of RNN by 1-step tBPTT and BPTT, as well as tPCNs with different inference iterations from $1$ to $50$. As can be seen, tPCN with a single inference iteration has identical performance to RNN trained by tBPTT, and its performance will improve as we increase the number of inference iterations but will saturate around $20$ iterations. Overall, this graph suggests that tPCN with $5$ or more inference iterations can effectively perform temporal credit assignment that improves upon tPCN$1$ or $1$-step tBPTT, potentially due to the eligibility trace. However, this eligibility trace arises from local inference dynamics (Eq.~\ref{eq:tpc_inf}) rather than from unrolling the RNN graph as in \cite{bellec2020solution}. This improvement is also reflected in the grid scores (Fig.~\ref{fig:tbptt}D), although increasing the inference iterations does not necessarily result in better grid score representations. We suspect that although the gridness of latent representations is somewhat related to path integration performance, their relationship is not linear. It is also worth noting that to fully evaluate the similarities and differences between BPTT and tPCN, an in-depth comparison is needed across different tasks and versions of tBPTT. We aim to investigate this question in future works as it is beyond the scope of this chapter.

\subsection{Robustness of grid cell representations in tPCN}\label{sec:robust}

\begin{figure}[t]
\begin{center}
\includegraphics[width=\textwidth]{figs/chap6/tpc_robustness.pdf}
\end{center}
\vspace{-0.5cm}
\caption{\textbf{Robust emergence of grid cells in tPCN.} A, B: Path integration RMSE and grid scores of tPCN in different setups. ``Stationary baseline'' refers to a model that always predicts the initial position regardless of movement. C-I: Firing fields of latent neurons in tPCNs with C: Gaussian place cells; D: $f=$\texttt{tanh}; E: $h=$\texttt{tanh}; F: $1.8\tn{m}\times 1.8\tn{m}$ environment; G: $1.2\tn{m}\times 1.2\tn{m}$ environment; H: $256$ latent neurons; I: tPCN without velocity input.}
\label{fig:tpc-robust}
\vspace{-0.2cm}
\end{figure}

Inspired by \citet{schaeffer2022no}, we examine the robustness of our results against different architectural choices of the tPCN model, to understand what contributes to the emergence of grid cells within tPCN. Specifically, we vary the following components of the model: 1) Encoding of the place cell activities; 2) Output nonlinearity $f$; 3) Recurrent nonlinearity $h$; 4) Environment sizes; 5) Latent sizes and 6) Velocity input to the model. The baseline model has DoS place cell encodings, $h=\texttt{ReLU}$, $f=\texttt{softmax}$, $1.4\tn{m}\times 1.4\tn{m}$ environment and latent size $2048$ with velocity inputs.

We first examine whether replacing the place encoding with Gaussian curves affects the model's performance. As shown in Fig.~\ref{fig:tpc-robust}A, B and C, the Gaussian place cells do not affect the path integration performance, but the latent representations are no longer hexagonal. This is consistent with earlier findings that the DoS place cell encoding is necessary for hexagonal grid cells \citep{dordek2016extracting,sorscher2023unified,schaeffer2022no}.

The choices of $f$ and $h$ are particularly interesting: as discovered by earlier works \citep{dordek2016extracting, sorscher2023unified}, a choice of $h$ that imposes non-negativity constraint on the latent activities, such as \texttt{ReLU}, is necessary for the emergence of hexagonal grid cells. In our tPCN model, the activation functions are also important for biological plausibility: in both Eq.~\ref{eq:tpc_inf} and Eq.~\ref{eq:tpc_learn}, the multiplication with the Jacobians $h'$ and $f'$ can be reduced to local, element-wise multiplications if $h$ and $f$ are element-wise nonlinearities such as $\texttt{ReLU}$ and $\texttt{tanh}$. Although it is possible to design a circuit to perform the computations in $\texttt{softmax}$ \citep{snow2022biological}, it is unclear how the Jacobian matrix of $\texttt{softmax}$ can be computed in a biological circuit. Therefore, we first replace $f$ with a $\texttt{tanh}$ function in our tPCN model and evaluate the model's performance in both path integration and its latent representations. As shown in Fig.~\ref{fig:tpc-robust}A, replacing $f$ with $\texttt{tanh}$ results in slightly worse path integration performance and lower grid scores than the $\texttt{softmax}$ baseline. However, visually, the latent representations are hexagonal and grid-like (Fig.~\ref{fig:tpc-robust}D), suggesting that using a biologically more plausible $f$ would not significantly affect the emergence of grid cells within tPCN. On the other hand, replacing the non-negative constraint (\texttt{ReLU}) on the latent activities with $h=\texttt{tanh}$ results in the amorphous latent representations (Fig.~\ref{fig:tpc-robust}E), which is consistent with \cite{sorscher2023unified}.

We next investigate the impact of the size of the environment, by training tPCN within a square environment of size $1.8\tn{m} \times 1.8\tn{m}$ (big) and an environment of size $1.2\tn{m} \times 1.2\tn{m}$ (small). Changing environment sizes does not affect the path integration performance, and does not affect tPCN's capability of developing grid cells either (Fig.~\ref{fig:tpc-robust}F for big environment and G for small environment). We also vary the number of latent neurons in the model from $256$ to $512$ and $1024$, which does not affect the grid cell representations (Fig.~\ref{fig:tpc-robust}H shows the latent representations learned by a tPCN with $256$ latent neurons). However, with fewer latent neurons, the performance in path integration becomes worse as the model has fewer number of parameters to perform the task (Fig.~\ref{fig:tpc-robust}A).

Earlier studies using PCNs to model visual representations have mostly used unsupervised PCNs \citep{rao1999predictive,olshausen1996emergence,millidge2024predictive}, which corresponds to blocking the velocity input $\bfv_t$ into tPCN in Fig.~\ref{fig:networks}B. Here we asked how removing velocity input would affect the path integration performance and grid cell emergence of tPCN. Mathematically, this is achieved simply by re-defining $\tilde{\bfg}_t := \bfW_{\tn{r}} \hat{\bfg}_{t-1}$ without changing any inference or learning dynamics.
% , and this velocity-free model has been shown to approximate Kalman filtering \citep{millidge2024predictive}. 
It can be seen from Fig.~\ref{fig:tpc-robust}A that the path integration performance is significantly affected by the absence of velocity input, with an RMSE even higher than the stationary baseline, where the model does not predict any movement at all. Intriguingly, the latent representations developed by this unsupervised tPCN are still grid cell-like (Fig.~\ref{fig:tpc-robust}I) with a similar grid score distribution to the baseline model. This result demonstrates that grid cells can still emerge even in a model unable to perform path integration at all. Therefore, our model predicts that path integration is not a sufficient condition for the emergence of grid cells, which resonates with \citet{schaeffer2022no}. In other words, it predicts that animals unable to navigate due to impaired velocity encoding may still develop grid cells as a result of self-motion.
% It also strengthens the argument in \citet{schaeffer2022no} that path integration should not serve as a sufficient condition for grid cells, and the exact relationship between path integration and grid cells should be more cautiously evaluated.

\section{Discussion} 
\subsection{Summary}
In this chapter, we have demonstrated a biologically plausible learning rule for grid cells based on predictive coding. We have shown that with sparsity and non-negative constraints, classical PCNs can develop grid cell-like representations of batched place cell inputs. With inputs representing trajectories of moving agents, tPCN can also develop grid cell activations while performing path integration. We have developed a theoretical understanding of this property of tPCN by deriving and comparing its learning dynamics with that of BPTT, showing that unrolling a recurrent network is unnecessary for it to learn grid cells, and a more plausible approach with recursive inference dynamics should suffice. Furthermore, we have examined the robustness of our results by varying hyper-parameters of the model, and found that grid cells can be learned even without velocity inputs. Overall, our work demonstrates that predictive coding can serve as an effective and biologically plausible plasticity rule for neural networks to learn grid cells observed in the MEC. Importantly, compared with earlier learning rules specialized for grid cells, predictive coding is a general learning rule able to reproduce many other cortical functions and representations. Thus, our findings suggest that a single, unified plausible learning rule can be employed by the brain to find the most appropriate representation of cortical inputs in different regions.

\subsection{Relationship to RNNs}
In this chapter we have demonstrated the relationship between tPCN and single-step truncated BPTT, which partly resolves the question left in Chapter \ref{chap:time} on whether tPCN can effectively learn long-term dependencies. Our results here suggest that in highly Markovian tasks, such as path integration, where long-range dependency does not play an important role, tPCN can match BPTT. However, it may under-perform in tasks that requires long-term credit assignment, such as those involving language \citep{graves2013speech}. An interesting future line of research can thus explore tPCN's performance in these more challenging machine learning tasks. A possible solution to the long-term learning issue with tPCN can be multiple layers of recurrently connected neurons, where higher layers may account for long-term dependencies, similar to what was discussed in another temporal version of predictive coding \citep{jiang2022dynamic}. Another route is to apply tPCN to neuroscience-related tasks that have been studied using BPTT-trained RNNs, such as those studied in \citet{yang2019task}, which could provide a more plausible account of how the brain performs these tasks.

\subsection{Relationship to biological data}
This chapter is mainly focused on demonstrating that grid cells can emerge in predictive coding models. The abundant experimental findings, on the other hand, provide many future directions for this research to follow. For example, it is known that the development of grid cells is abrupt \citep{wills2010development}. It would be interesting to investigate if predictive coding can train networks to abruptly learn grid cells, and if the model could provide a mechanistic account for the abrupt learning. Additionally, grid cells emerged in predictive coding models does not present the multi-scale property observed in rats \citep{hafting2005microstructure}. In particular, it is known that grid cell scales have a multi-peak, discretized distribution \citep{stensola2012entorhinal}. One possible approach to reproduce this observation is to use place cell inputs with distributed place field sizes, which has a known continuous distribution \citep{dabaghian2012topological}, instead of the fixed field size used in experiments in this chapter. The hexagonal firing pattern of grid cells is also known to distort in polarized environments such as trapezoidal ones \citep{krupic2015grid}. This property has been replicated in reinforcement learning models also employing the concept of predictive coding \citep{stachenfeld2017hippocampus}, and it would be intriguing to explore whether predictive coding networks can also reproduce this phenomenon.