Hippocampal Spike-Timing Correlations Lead to Hexagonal Grid Fields

TL;DR

This paper presents a theory linking hippocampal spike-timing correlations to the formation of hexagonal grid fields through spike-timing-dependent plasticity and pattern formation dynamics.

Contribution

It introduces a novel model explaining how temporal spike correlations transform into spatial grid patterns via synaptic plasticity and pattern formation mechanisms.

Findings

Hexagonal grid patterns emerge under specific parameter regimes.

Spike-timing-dependent plasticity can produce spatial firing fields.

The model aligns with observed grid cell patterns in medial entorhinal cortex.

Abstract

Space is represented in the mammalian brain by the activity of hippocampal place cells as well as in their spike-timing correlations. Here we propose a theory how this temporal code is transformed to spatial firing rate patterns via spike-timing-dependent synaptic plasticity. The resulting dynamics of synaptic weights resembles well-known pattern formation models in which a lateral inhibition mechanism gives rise to a Turing instability. We identify parameter regimes in which hexagonal firing patterns develop as they have been found in medial entorhinal cortex.

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}{J}_n={\int }_{-\infty }^{\infty }\mathrm{d}sW(s)C_n(s)+F(J_n).$$

$$C_{n{n}^{\prime }}(s):={\int }_{-\infty }^{\infty }\mathrm{d}t{H}_n(t)H_{n^{\prime }}(t-s).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{J}_n=\sum _{n^{\prime }}{J}_{n^{\prime }}{G}_{n{n}^{\prime }}+F(J_n),G_{n{n}^{\prime }}:={\int }_{-\infty }^{\infty }\mathrm{d}sW(s)C_{n{n}^{\prime }}(s).$$

$$H_n(t;\mathcal{P})=a{\mathrm{e}}^{\frac{-(x_{\mathcal{P}}(t)-x_n{)}^2}{2{\sigma }^2}}[\cos (\omega t+{\varphi }_n)+1]/2.$$

$$C_{n{n}^{\prime }}(s)$$

$$\frac{\mathrm{d}}{\mathrm{d}t}J(x)$$

$$\Gamma (r)=c\frac{\sqrt{\pi }{\mathrm{e}}^{-\frac{r^2}{4{\sigma }^2}}}{4v/(a\sigma )}\left(1+\alpha J_0\left(\frac{\pi r}{R}\right)+\frac{\beta r}{\sigma }{J}_1\left(\frac{\pi r}{R}\right)\right)$$

$$c[W]:$$

$$\alpha [W]:$$

$$\beta [W]:$$

$$W(s)=W_0(2\pi {\mu }^2{\rho }^2{)}^{-\frac{1}{2}}\left(1-(s/\rho {)}^2\right){\mathrm{e}}^{-\frac{s^2}{2{\rho }^2{\mu }^2}},$$

$$\tilde{\Gamma }(k):={\int }_{\Omega }\mathrm{d}x\Gamma (x)\cos (kx)\Theta [\cos (kx)-\cos (|k|R_g)]$$

Full text

The spatial position of an animal can be reliably decoded from the neuronal activity of several cell populations in the hippocampal formation Harris et al. (2003); Mathis et al. (2012); Kammerer and Leibold (2014). For example, place cells in the hippocampus fire at only few locations in a spatial environment O’Keefe (1976); Muller et al. (1987) and the position of the animal can be readily read out from single active neurons. Grid cells of the medial entorhinal cortex (MEC) fire at multiple distinct places that are arranged on a hexagonal lattice Fyhn et al. (2004); Hafting et al. (2005). Although hexagonal patterns are abundant in nature and there exist well-studied physical theories for their emergence, the mechanistic origin of this neuronal grid pattern is still unclear. Initially it was suggested that they result from continuous attractor dynamics Fuhs and Touretzky (2006); Burak and Fiete (2009) or superposition of plane wave inputs Burgess (2008) and, based on circuit anatomy, place cells would then result from a superposition of many grid cells McNaughton et al. (2006); Moser and Moser (2008). More recent experiments, however, reported place cell activity without intact grid cells, such that grid cells are not the unique determinants of place field firing Koenig et al. (2011); Brandon et al. (2011); Wills et al. (2010); Hales et al. (2014); Schlesiger et al. (2015). Conversely, it would thus be possible that grid fields may arise from place field input as suggested in Kropff and Treves (2008); Dordek et al. (2016); Stepanyuk (2015). The biological mechanisms proposed by these latter theories, however, remain hypothetical. In the present Letter, we propose a learning rule for grid cells based on the individual spike timings of place cells using spike-timing dependent synaptic plasticity (STDP) Markram et al. (1997); Kempter et al. (1999); Bi and Poo (1998). The theory thereby predicts that the observed temporal hippocampal firing patterns (phase precession and theta-scale correlations; see below) O’Keefe and Recce (1993); Dragoi and Buzsaki (2006); Feng et al. (2015) translate the temporal proximity of sequential place field spikes into spatial neighborhood-relations observed in grid-field activity. For our model to work, we only have to assume that the synaptic plasticity rule averages over a sufficiently long time interval.

Model.

We use the classical formulation of pairwise additive STDP Gerstner et al. (1996); Kempter et al. (1999), where the update of a synaptic weight $J_{n}\ ,\ n=1,\ldots,N$ at time $t$ is computed as Kempter et al. (1999)

[TABLE]

$C_{n}(s)$ denotes the time averaged correlation function between the spike train of presynaptic neuron $n$ and the postsynaptic neuron, the learning window $W(s)$ describes the update of the synaptic weight as a function of the time difference $s$ between a pair of pre- and postsynaptic action potentials, and the function $F$ implements soft bounds for the weight increase. The dynamics is further constrained such that weights cannot become negative.

To be able to treat eq. (1) analytically, we use a linear Poisson neuron model, i.e., the mean firing rate of the postsynaptic neuron $E(t)=\mathbf{J}\cdot\mathbf{H}(t)$ results from a weighted sum of hippocampal firing rates $\mathbf{H}=(H_{1}(t),\ldots,H_{N}(t))^{\rm T}$. Under these assumptions $C_{n}(s)$ can be approximated for large $N$ Kempter et al. (1999) as $C_{n}(s)=\sum_{n^{\prime}}J_{n^{\prime}}C_{nn^{\prime}}(s)$ with

[TABLE]

Inserting the correlation functions from eq. (2) into the weight dynamics from eq. (1) yields

[TABLE]

Following van Rossum et al. (2000); Kistler and van Hemmen (2000) we introduce the quadratic stabilization term $F(J)=F_{0}\,J\,(K-J)\ ,F_{0}>0$ that implements a soft upper bound.

As an input to the postsynaptic neuron, we consider a population of $N$ hippocampal place cells. The firing of these neurons is characterized by a bell shaped envelope modulating the spatial path $x_{\cal P}(t)$ and oscillations in time $t$ (Fig. 1A),

[TABLE]

The oscillation frequency $\omega$ of a neuron is slightly higher than the frequency $\omega_{\theta}$ of the theta oscillation ($\sim 8$ Hz) in the local field potential giving rise to a phenomenon called theta phase precession (Fig. 1B): spikes early in the field come at later phases than spikes late in the field O’Keefe and Recce (1993). During traversal of a place field, phase precession spans a whole theta cycle Maurer et al. (2006). Thus the two frequencies have to relate to each other like $\omega=\omega_{\theta}+\frac{\pi}{R}\,v$, with $R$ denoting the distance from the place field center at which the firing rate has decreased to $10\%$ (i.e., $R=\sigma\,\sqrt{2\,\ln(10)}$), and $v$ denoting the running speed, which we fix at 25 cm/s. At each individual entry into the place field, the phase of the cellular oscillation is reset to phase zero with respect to the theta oscillation phase $\phi_{\theta}$ by fixing $\phi_{n}=(\phi_{\theta}-2\pi)\,\omega/\omega_{\theta}$.

To obtain a closed expression for the correlation function $C_{nn^{\prime}}(s)$ (Fig. 1C), the time average in eq. (2) is performed over all straight paths ${\cal P}$ crossing the center of the field overlap. For place fields with identical width $R$, firing rates $a$, and at small lags $s$, we obtain

[TABLE]

with place field distance $r=|x_{n}-x_{n^{\prime}}|$, and $J_{0,1}$ denoting Bessel functions of first kind (see Supplementary Material at [URL] for derivation).

In contrast to 1-d fields where the distance of place fields is reflected by the lag of the correlation peak (theta compression) Dragoi and Buzsaki (2006), correlation functions in 2-d are symmetric because of the symmetry of the path, however, the distance of the place field centers is encoded in the amplitude of the correlation peak at lag 0 (Fig. 1D,E).

Weight Dynamics.

Assuming that the putative grid cells receive inputs from a large number $N\gg 1$ of place cells that sufficiently cover the encoded area, we replace the presynaptic index $n$ by the position of the place field center $x$, i.e., $G_{nn^{\prime}}\to\Gamma(|x_{n}-x_{n^{\prime}}|)$, and thereby translate the learning eq. (3) to continuous coordinates,

[TABLE]

Examples of the convolution kernel $\Gamma(|x|)$ for different learning window functions $W$ are depicted in Fig. 1F, G. The development of the weights follows the pattern formation principles of a lateral inhibition system Swindale (1980). Indeed, the integro-differential equation (IDE) (6) involves non-local interactions effectively implemented through the convolution kernel, inducing a strong close range potentiation and a weaker long range depression of neighbouring synapses, as observed in the typical shape presented in Fig. 1G. A general window-dependent kernel can be obtained for the correlation function from eq. (Model.) as

[TABLE]

[TABLE]

which can take a Mexican-hat type shape for qualitatively different learning window functions $W$, due to the symmetry of the correlation function (Fig. 1F, G). To see this, we can regard Hebbian-like windows to be modelled as the product of a Gaussian and a polynomial of some order $m$, $W(s)=\exp[-s^{2}/(2\rho^{2}\mu^{2})]\,P_{m}(s/\rho)$. The functionals defined in eq. (8) then inherit the symmetries from the cross-correlation, since all of the odd terms in the polynomial cancel out during integration. Thus, in subsequent numerical investigations we focus on windows up to second polynomial order

[TABLE]

whose free parameters $\rho$ and $\mu$ determine its zeroes $s_{0}=\pm\rho$ and negativity $\int W=W_{0}(1-\mu^{2})$. In Figure 1F,G we used $\rho=23$ ms, $\mu=1.025$ and added a linear term $s/\rho$ to the polynomial to get the asymmetric window (grey lines). The $W$-dependent functions $c,\alpha$, and $\beta$ defining the kernel $\Gamma$ are given in the Supplementary Material [URL].

A numerical evaluation of the learning IDE (6) with periodic boundary conditions reveals that the spatially isotropic kernel $\Gamma$ can result in hexagonal packing structures (Fig. 2A, B). Simulations of spiking Poisson neurons confirm these predictions of the meanfield theory (Fig. 2A). As indicated by the kernel function in eq. (7), the grid spacing only depends on the spatial scale $\sigma$ of the place fields in the input (Fig. 2C).

An experimentally accessible quantity to compare our model results to is the ratio of grid field radius $R_{g}$ to grid spacing $l$ as indicated in Fig. 2A. In experiments, the ratio $R_{g}/l$ for grid cells has been determined to be about $0.3$ (Hafting et al., 2005), and for a perfectly hexagonal grid $R_{g}/l$ relates to the fraction $\nu$ of field size per area as $\nu=\pi\,R_{g}^{2}/(\frac{\sqrt{3}}{2}\,l^{2})$. The field fraction $\nu$ can be readily accessed from our numerics as the fraction of non-zero synaptic weights, and for the used learning window fits the experimentally obtained $R_{g}/l$ (Fig. 2D) for all choices of $\sigma$. Finally, learning converges faster for large spacing (large $\sigma$) consistent with a larger amplitude of $\Gamma$ [eq. (7) and Fig. 2B,E].

Linear theory at early times.

Some analytical understanding of the weight dynamics from eq. (6) can be gained from a neural field theory approach Swindale (1980); Linsker (1986); Ermentrout et al. (1986); Miller et al. (1989); Shouval et al. (1996). In this framework, we can neglect the effect of the non-linearities at early times, and focus only on the convolution term $\Gamma\ast J$. The emerging dynamics can be readily understood by looking at the evolution of the weights in Fourier space $\partial_{t}\hat{J}(k)=\hat{\Gamma}(k)\,\hat{J}(k)$, which makes evident that the wave number $k_{m}$ maximizing the kernel Fourier transform $\hat{\Gamma}$ (see Supplementary Material [URL]) will exponentially overgrow all other modes (if $\hat{\Gamma}(k_{m}){>}0$), thus setting the initial periodicity of the pattern.

In a finite region of the parameter space $(\mu,\rho)$ of the learning window, the bimodal (Mexican-hat) shape of $\Gamma$ ensures the existence of a Turing instability, i.e., a transition to single maximum of $\hat{\Gamma}$ (with $\hat{\Gamma}(k_{m}){>}0$) at a non-zero $k_{m}{>}0$ (Fig. 3A-C). Similar to previous work on pattern formation in lateral inhibition systems (e.g. Ermentrout and Cowan (1979)), the permitted parameter region ($k_{m}$ exists and is positive) gives rise to stripe-like and hexagonal patterns (Fig. 3D). In the Supplementary Material [URL], we also provide a complementary description of the pattern formation process based on an approximation of eq. (6) by a partial differential equation (Swift-Hohenberg equation Swift and Hohenberg (1977)), which corroborates the results from the linear theory.

To connect the resulting patterns to other feed-forward models of grid field formation, we parameterize $\Gamma$ by the shape factor $r_{0}/r_{m}$ (Fig. 3C), which is the fraction between the zero and the minimum of $\Gamma$. The shape factor $r_{0}/r_{m}$ reduces the two-parameter learning window to a single qualitatively descriptive parameter, which can be used to describe the bimodal kernel $\Gamma$ independently of the hypothesized biological mechanism. If $r_{0}/r_{m}$ is large ($\sim 0.8$), $\Gamma$ shows only little negativity and the emerging pattern is stripe-like (Fig. 3C, D), if $r_{0}/r_{m}$ is small, $\Gamma$ exhibits strong negativity, the firing fields become dispersed and the pattern looses hexagonality. Hexagonal patterns arise for $r_{0}/r_{m}$ roughly between 0.65 and 0.75 (Fig. 3D). In this region, the shape factor virtually completely determines the geometrical properties of the steady state (Fig. 3E, F). Values of $r_{0}/r_{m}$ that give rise to hexagonal grids can also be identified via the ratio of field width per grid spacing $R_{g}/l$. According to our theory, the experimentally observed value 0.3 Hafting et al. (2005) is achieved with a shape factor of about $r_{0}/r_{m}=0.7$ (Fig. 3F). For higher values of $r_{0}/r_{m}$, $R_{g}/l$ increases to a point where a periodic pattern cannot further dissociate into disjoint fields and the stable pattern becomes stripe-like. For lower $r_{0}/r_{m}$, $R_{g}/l$ decreases, and at some point, the small fields no longer repel each other strongly enough to produce a symmetrical arrangement.

Positivity constraint.

The grid spacing $l$ predicted by the linear theory, however, consistently underestimates the spacing derived from the numerical solution of the mean field dynamics (Fig. 2C). The reason for this error is that, after the initial growth phase, the synaptic weights are influenced by the non-linearities, most importantly the constraint that they cannot become negative.

The impact of this positivity constraint can be intuitively understood if we interpret the convolution $\Gamma\ast J$ as an operation that detects the best overlap of an oscillatory pattern $J\propto\cos(k\,x)$ with a given kernel $\Gamma$. However, after the lowest weights reach zero they stop contributing to the convolution and a slightly lower wave number $\tilde{k}_{m}$ maximizing

[TABLE]

will be favored as the fastest growing mode ($\Theta$ denoting the Heaviside function). Similarly, a particular field size $R_{g}$ maximizing $\tilde{\Gamma}$ will be selected. In the experimentally relevant case $|k|\,R_{g}=2\pi R_{g}/l=2\pi\times 0.3$, numerical maximization of eq. (10) yielded the predicted wave number $\tilde{k}_{m}$ (solid line in Fig. 2C), which excellently agrees with the numerical solutions of the meanfield dynamics.

Conclusion.

For a large variety of STDP windows, the spike-timing correlations of 2-d place cells can account for a feed-forward learning of hexagonal grid patterns. Synaptic plasticity thereby averages over running trajectories of tens of minutes, hence, translating the temporal correlations into a dense code for space. Our model thus predicts that grid cells are generated in the output structures of the hippocampus, e.g., the deep layers of the medial entorhinal cortex Witter (1993) or the parasubiculum van Groen and Wyss (1990). While our linear theory provides a good prediction of grid spacing as well as for conditions that permit structure formation, determining the boundary between hexagonal and stripe-like patterns is less straight-forward and has to take into account the non-linearities. The standard approach, non-linear bifurcation analysis Ermentrout and Cowan (1979, 1980a, 1980b), is difficult because of the strong non-linearity introduced via the positivity constraint, which strongly influences the selection of the final pattern. Despite this drawback, our model provides a universal framework in that it encompasses current models of grid field formation that can be mapped to convolutions with Mexican hat-type kernels that give rise to a Turing instability.

The authors are grateful to Andreas Herz and Anton Sirota for discussions and Martin Stemmler for comments on the manuscript. This work was funded by the German Research Association (DFG), Grant No. LE2250/5-1.