A single-cell spiking model for the origin of grid-cell patterns

Spatial cognition in mammals is thought to rely on the activity of grid cells in the entorhinal cortex, yet the fundamental principles underlying the origin of grid-cell firing are still debated. Grid-like patterns could emerge via Hebbian learning and neuronal adaptation, but current computational models remained too abstract to allow direct confrontation with experimental data. Here, we propose a single-cell spiking model that generates grid firing fields via spike-rate adaptation and spike-timing dependent plasticity. Through rigorous mathematical analysis applicable in the linear limit, we quantitatively predict the requirements for grid-pattern formation, and we establish a direct link to classical pattern-forming systems of the Turing type. Our study lays the groundwork for biophysically-realistic models of grid-cell activity.


Introduction
Grid cells are neurons of the medial entorhinal cortex (mEC) tuned to the position of the animal in the environment [1,2]. Unlike place cells, which typically fire in a single spatial location [3,4], grid cells have multiple receptive fields that form a strikingly-regular triangular pattern in space. Since their discovery, grid cells have been the object of a great number of experimental and theoretical studies, and they are thought to support high-level cognitive functions such as self-location [e.g. 5,6], spatial navigation [e.g. [7][8][9], and spatial memory [10,11]. Nevertheless, to date, the mechanisms underlying the formation of grid spatial patterns are yet to be understood [12,13]. The attractor-network theory proposes that grid fields could arise from a path-integrating process, where bumps of neural activity are displaced across a low-dimensional continuous attractor by self-motion cues [14][15][16][17][18][19][20][21]. The idea that self-motion inputs could drive spatial firing is motivated by the fact that mammals can use path integration for navigation [22], that speed and head-direction signals have been recorded within the mEC [23,24], and that, in the rat [1,25] but not in the mouse [26,27], grid firing fields tend to persist in darkness. However, grid-cell activity may rely also on non-visual sensory inputs-such as olfactory or tactile cues -even in complete darkness [28]. Additionally, the attractor theory alone cannot explain how grid fields are anchored to the physical space, and how the properties of the grid patterns relate to the geometry of the enclosure [29][30][31].
Here we focus on the idea that grid-cell activity does not originate from self-motion cues, but rather from a learning process driven by external sensory inputs. In particular, it was proposed that grid patterns could arise from a competition between persistent excitation by spatially-selective inputs and the reluctance of a neuron to fire for long stretches of time [49][50][51][52][53]. In this case, Hebbian plasticity at the input synapses could imprint a periodic pattern in the output activity of a single neuron. Spatially-selective inputs, i.e., inputs with significant spatial information, are indeed abundant within the mEC [54-56] and its afferent structures [57-61] And spike-rate adaptation, which is ubiquitous in the brain [62], could hinder neuronal firing in response to persistent excitation.
Kropff and Treves [49] explored this hypothesis by means of a computational model; see also [63][64][65][66][67] and Sec Related models for similar works. The emergence of grid-like patterns was demonstrated with theoretical arguments and with numerical simulations of a rate-based network. However, because of a relatively abstract level of description, the outcomes of the model could not be easily confronted with experimental data. Specifically, the simulations included a network-level normalization mechanism that constrained the mean and the sparseness of the output activity, and it remained unsettled whether grid patterns could emerge in a single-cell scenario. Additionally, the synaptic weights did not obey Dale's law. And the robustness of the model was not tested against shot noise due to stochastic spiking. Finally, the link between the numerical simulations and the underlying mathematical theory remained rather loose.
To overcome these issues, we propose here a single-cell spiking model based on similar principles as the model by Kropff and Treves [49], but that is, on the one hand, more biologically realistic, and on the other hand, better suited for mathematical treatment. Importantly, we show that grid patterns can emerge from a single-cell feed-forward mechanism needless of any network-level interaction (although recurrent dynamics may be still required to explain the coherent alignment of grid patterns [1]). To increase biological plausibility, we consider a stochastic spiking neuron model, and we constrain the synaptic weights to non-negative values (Dale's law). Finally, by studying the model analytically, we quantitatively predict the requirements for grid-pattern formation, and we establish a direct link to classical pattern-forming systems via the Turing instability [68].

Model of neural activity
We consider a single cell that receives synaptic input from N spatially-tuned excitatory neurons. Input spike trains S in i ðtÞ ≔ P k dðt À t in i;k Þ for i = 1, 2, . . ., N are generated by independent inhomogeneous Poisson processes with instantaneous rates r in i ðtÞ where δ(t) is the Dirac delta function, and t in i;k is the timing of the k th input spike at synapse i. Similarly, the output spike train S out ðtÞ ≔ P k dðt À t out k Þ is generated by an inhomogeneous Poisson process with instantaneous rate r out (t) where t out k denotes the timing of the k th output spike. We assume that inputs are integrated linearly at the output, and that the output neuron is equipped with an intrinsic spike-rate adaptation mechanism, that is, where r 0 is a baseline rate, w i is the synaptic weight of input neuron i, and the function K is a temporal filter modeling the spike-rate adaptation dynamics. Note that the instantaneous output rate r out depends only on the temporal history of the input spikes and that there is no reset mechanism after the emission of an output spike. The impulse response of the adaptation kernel K is the sum of two exponential functions: where τ S and τ L are the short and long filter time constants (0 < τ S < τ L ), and the parameter μ > 0 sets the filter integral R 1 0 dt KðtÞ ¼ 1 À m ( Fig 1A). Intuitively, at the arrival of an input spike, the firing probability of the output neuron is first increased for a short time that is controlled by the time constant τ S , and then decreased for a longer time that is controlled by the time constant τ L . This second hyper-polarization dynamics effectively hinders the neuron to fire at high rates for long stretches of time, mimicking a spike-rate adaptation mechanism [69][70][71]. From a signal-processing perspective, the adaptation kernel K performs a temporal bandpass filtering of the input activity (Fig 1B), and the two time constants τ S and τ L control the resonance frequency k res at which the filter response is maximal. Note that in Sec Pattern formation with after-spike potentials we study a variant of the present model where neuronal adaptation is obtained though after-spike hyperpolarizing potentials associated to the output activity of the neuron.

Model of synaptic plasticity
We assume spike-timing dependent plasticity (STDP) at the input synapses [e.g. [72][73][74][75][76]. Input and output spikes trigger weight changes Δw i according to the following rule: 1. For each pair of a post-synaptic spike and a pre-synaptic spike at synapse i, we set 2. For each pre-synaptic spike at synapse i, we set where η ( 1 is a small learning rate, and the STDP learning window W(Δt) sets the weight change as a function of the time difference Δt ≔ t pre − t post between pre-and post-synaptic spikes. We consider a symmetric STDP learning window [77] WðDtÞ where the time constant τ W > 0 controls the maximal time lag at which plasticity occurs, and W tot ¼ R 1 À 1 dt WðtÞ is the integral of the learning window. The first part of the learning rule (Eq 3) is the classical Hebbian term whereas the second part (Eq 4) is a local normalization term that stabilizes the average synaptic strength w av ¼ N À 1 P N i¼1 w i and prevents the individual weights to grow unbounded. This normalization term mimics local homoeostatic processes observed experimentally [78][79][80]; see also [81] for a review. The parameters α > 0 and β set, respectively, the rate of weight decay and the target average weight w av (Sec Weight normalization). Importantly, the synaptic weights are constrained to non-negative values by imposing the hard bounds

Model of input spatial tuning
We consider excitatory inputs with firing rates r in i that are tuned to the spatial position of a virtual rat exploring a square arena of side-length L, i.e., where x t is the position of the virtual rat at time t, and C in i is a spatial tuning curve. We characterize the spatial tuning curves C in i in two alternative scenarios: 1. spatially-regular inputs, i.e., each input has a single spatial receptive field; 2. spatially-irregular inputs, i.e., each input has multiple spatial receptive fields at random locations.
The first scenario, which is reminiscent of hippocampal place-cell activity [3,82,83], is easier to study analytically and cheaper to simulate numerically. The second scenario, which is reminiscent of parasubicular activity [57-61], is motivated by the anatomy of the entorhinal circuit (Sec Input spatial tuning and the origin of grid-cell patterns). In both cases, we consider circularly-symmetric receptive fields that cover the arena evenly. Indeed, place fields in open environments do not show systematic shape biases, and, in the absence of specific reward or goal locations, their centres are roughly homogeneously distributed [3, 57-61, 82, 83]. Note, however, that border-like inputs [84,85]-which are not radially-symmetric-are present in the real system, but not explicitly modeled here. Finally, for simplicity, we assume periodic boundaries at the edges of the arena. Spatially-regular inputs. In the case of spatially-regular inputs, we assume tuning curves of the form where r i is the receptive-field center of neuron i and G is a Gaussian function: The parameter σ > 0 sets the width of the receptive field, and r av is the average firing rate in the environment. We assume that the input receptive-field centers r i cover the entire arena evenly. This input scenario is considered for the mathematical derivations in Sec Weight dynamics for spatially-regular inputs and for the numerical simulations in Secs Emergence of grid spatial patterns and Geometrical properties of the grid patterns.
Spatially-irregular inputs. In the case of spatially-irregular inputs, each tuning curve C in i is the sum of M > 1 Gaussian receptive fields with random amplitudes A ij and random receptive-field centers r ij with i = 1, 2, . . ., N and j = 1, 2, . . ., M, that is, The scaling factors b i ¼ P M j¼1 A ij normalize the inputs C in i to the same average rate r av , and all the superimposed fields share the same field size σ (Eq 9). The field amplitudes A ij are uniformly distributed in the range (0, 1), and the receptive-field centers r ij are uniformly distributed in the environment. This input scenario is considered for the mathematical derivations in Sec Eigenvalue spectrum for spatially-irregular inputs and for the numerical simulations in Sec Pattern formation with spatially-irregular inputs.

Model of spatial exploration
The movement of the virtual rat follows a smooth random walk that satisfies the following three assumptions: (i) the movement speed v is constant in time; (ii) the random walk is isotropic and ergodic with respect to the auto-covariance; (iii) the virtual-rat trajectories are smooth within time stretches shorter than the time length τ max = 5τ L of the adaptation kernel K ( Fig  1A). Note that assumption (i) is obviously not valid in general. However, because synaptic plasticity acts on a time scale that is much slower than behaviour, the relevant variable for pattern formation is the rat running speed averaged over long stretches of time (e.g. minutes), which can be considered approximately constant. We assume an average running speed of 25 cm/s, which is experimentally plausible [86]. Assumptions (ii) and (iii) hold by ignoring directional anisotropies deriving from the geometry of the environment, and by observing that experimental rat trajectories are approximately straight over short running distances (e.g, over distances shorter than 25 cm) [86].
Mathematically, the two-dimensional virtual-rat trajectories x t are sampled from the stochastic process where the angle θ t sets the direction of motion and W t is a standard Wiener process (Fig 2). The parameters v and σ θ control the speed of motion and the tortuosity of the trajectory. Note that we also perform simulations with variable running speeds. In this case, the speed is sampled from an Ornstein-Uhlenbeck process with long-term mean " v ¼ v.

Mathematical results on grid-pattern formation
The grid-cell model presented above is studied both analytically and numerically. In this section, we obtain an equation for the average dynamics of the synaptic weights, and we derive the requirements for spatial pattern formation. In Sec Numerical results on grid-pattern formation we demonstrate the emergence of grid-like activity by simulating both the detailed spiking model and the averaged system. The analytical results presented here may be skipped by the less mathematically-inclined reader.
We study structure formation in the activity of an output cell by averaging the weight dynamics resulting from the stochastic activation of input and output neurons (Sec Model of neural activity) and the STDP learning rule (Sec Model of synaptic plasticity), while a virtual rat explores a two-dimensional enclosure and the inputs are spatially tuned (Secs Model of input spatial tuning and Model of spatial exploration). We take both ensemble averages across spike-train realizations and temporal averages within a time window of length T. The averaging time length T separates the time scale of neural activation (of the order of the width τ W of the learning window W) from the time scale τ str of structure formation, i.e., τ W ( T ( τ str . Because τ str is inversely proportional to the learning rate η (Eq 29), such averaging is always possible provided that the learning rate η is small enough. In other words, we assume that within a time T, the virtual rat has roughly explored the entire environment, but the synaptic weights did not change considerably. In this case, the dynamics of the synaptic weights w i is approximated by a drift-diffusion process, where the deterministic drift term reads [74] with " w i ! 0. The functions S in i and S out denote input and output spike trains (Sec Model of neural activity), the angular brackets denote ensemble averages over input and output spike trains, and the overbars denote temporal averages, i.e., " f ðtÞ ≔ T À 1 R t tÀ T ds f ðsÞ. Following Kempter et al. [74] we derive where the ensemble averages read Finally, from Eqs 12-15 we obtain where we defined a ≔ r av a À Note that in deriving Eq 16 we approximated the temporal average of the input rates r in i with the spatial average r av of the input tuning curves C in i . This approximation holds with the assumption that in a time T the virtual rat roughly covers the entire space evenly.
By ignoring the non-linear weight constraints " w i ! 0, the average weight dynamics is described by a linear system with coupling terms C ij (Eq 16). The coefficients C ij are given by the temporal correlations of the input rates r in i and r in j , filtered by the adaptation kernel K and the STDP learning window W (Eq 17).
To further simplify the calculations, we assume that the low-pass filtering introduced by the STDP learning window can be neglected for the purpose of studying pattern formation. In particular, we assume that the learning window W decays much faster than the changes in the input correlations r in i ðt þ sÞr in j ðt À tÞ(Eq 17), which holds for τ W ( σ/v. In this case, we obtain where W tot is the integral of the learning window (Eq 5). Finally, by assuming smooth virtual-rat trajectories at constant speed v, the correlation matrix C ij can be estimated solely from the input tuning curves C in i and the adaptation kernel K (Sec Input correlation for general inputs, Eq 47): where L 2 is the area explored by the virtual rat. In Eq 21, the matrix element C ij is obtained by integrating the spatial cross-correlation of the input tuning curves C in i ? C in j over circles of radius τv, and by weighting each integral with the amplitude of the adaptation kernel K at time τ. Note that Eq 21 holds for generic spatial tuning curves C in i . Weight dynamics for spatially-regular inputs. To study the emergence of spatial patterns, we now consider the simplified scenario of spatially-regular inputs (Sec Spatially-regular inputs). That is, the input tuning curves C in i are circularly-symmetric Gaussian functions that cover the entire space evenly (Eqs 8 and 9). This input representation is particularly useful because it establishes a direct mapping between the neuron identity (the index i) and a position in physical space (the receptive-field center r i ). Therefore, studying pattern formation in the activity of the output neuron is reduced to studying pattern formation in the space of the synaptic weights. Note, however, that such a simple input scenario is not necessary for the formation of grid-cell patterns in general, as shown in Sec Pattern formation with spatially-irregular inputs.
With spatially-regular inputs, the average weight dynamics in Eq 16 can be rewritten by labeling the synaptic weights according to the corresponding receptive-field centers r i : where " wðr i Þ ¼ " w i and C(r i , r j ) = C ij . Additionally, in the limit of a large number N ) 1 of input neurons and receptive fields that cover the environment with constant density ρ = N/L 2 , the sum in Eq 22 can be replaced by an integral over all the receptive-field centers r 0 : where the correlation function C(r, r 0 ) is the continuous extension of the correlation matrix C ij = C(r i , r j ). Because the inputs are translation invariant (Eq 8), the correlation function C is also translation invariant, i.e., C(r, r 0 ) = C(r − r 0 , 0) = C(r − r 0 ), where we omit the second argument 0 ≔ (0, 0) for readability. In this case, the integral in Eq 23 can be expressed as a twodimensional convolution in space: Fig 3 shows the correlation C as a function of the input receptive-field distance |r − r 0 | for the adaptation kernel K in Fig 1 and Gaussian input fields with size σ = 6.25 cm (Eq 9). The function C has the shape of a typical Mexican-hat kernel, i.e., it is positive for short receptive-  , and zero otherwise. In this case, the synaptic weights of close-by input fields grow together whereas the synaptic weights of input fields that are further apart are repelled from each other (Eq 24). Such a competitive Mexican-hat interaction is at the basis of many pattern-forming systems found in nature, and it is directly related to diffusion-driven instabilities of the Turing type [see e.g. 87].
Eigenvalue spectrum for spatially-regular inputs. To study spatially-periodic solutions, we take the two-dimensional Fourier transform with respect to r at both sides of Eq 24: where we defined the Fourier transform pair k is a two-dimensional wave vector, and j ¼ ffiffiffiffiffiffi ffi À 1 p is the imaginary unit. Solving Eq 25 for k 6 ¼ (0, 0), we obtainŵ ðkÞ ¼ŵ 0 ðkÞexpðZlðkÞtÞ ð27Þ whereŵ 0 ðkÞ denotes the weight spectrum at time t = 0, and we defined The function λ(k) defines the eigenvalue spectrum of the dynamical system in Eq 24, and the corresponding eigenfunctions are the elements of the Fourier basis exp(2πj k Á r). Eq 28 is also called the dispersion relation of the system. Note that solving Eq 25 for k = (0, 0) one obtains the dynamics of the total synaptic weight, which is kept normalized by the learning rule (Sec Weight normalization). From Eq 27, the Fourier modes of the synaptic weightsŵðkÞ grow or decay exponentially with rates proportional to the eigenvalues λ(k). Therefore, a structure in the synaptic weights emerges on a time scale where λ max ≔ max k [λ(k)] is the largest eigenvalue in the system. Importantly, the eigenvalues λ(k) are linearly related to the Fourier transform of the inputcorrelation functionĈðkÞ (Eq 28), which is circularly-symmetric for circularly-symmetric inputs. In this case, in Sec Input correlation for spatially-regular inputs (Eq 62) we derivê whereG andK sp (Eqs 63 and 64) are the zeroth-order Hankel transforms (Eq 59) of the input tuning curve G (Eq 9) and of the equivalent adaptation kernel in space Finally, by plugging Eqs 30 into 28, we obtain From Eqs 27 and 32 we recognize a necessary condition for spatial patterns to emerge: the eigenvalue spectrum λ(k) = λ(k) shall have a global maximum λ max > 0 at a frequency k max > 0. In this case, all the Fourier modes k at the critical frequency |k| = k max are unstable (Eq 25), and spatially-periodic patterns could emerge. Fig 4 shows the critical frequency k max (panels A1 and B1) and the largest eigenvalue λ max (panels A2 and B2) as a function of the parameters of the adaptation kernel K, i.e., the short time constant τ S , the long time constant τ L , and the kernel integral 1 − μ (Eq 2). The input receptive-field width σ is kept constant. In panels A1 and B1, the green-shaded regions correspond to parameter values where periodic grid-like patterns could emerge (k max > 0). Conversely, the white regions denote parameter values where place-cell-like receptive fields could emerge (k max = 0) [88]. We note that the spatial scale of the periodic patterns depends on the long adaptation time constant τ L (panel A1), but is largely unaffected by the short time constant τ S (panel B1). Additionally, the largest spatial frequencies are obtained for small values of τ L and negative kernel integrals (panel A1). This leads us to the following predictions: the grid scale shall depend on the long temporal dynamics of the adaptation kernel, and the smallest grid scales require adaptation kernels with an overall inhibitory effect on the activity of the output neuron. We also note that larger values of τ L correspond to larger values of λ max (panel A2). Thus, we predict that grids at larger scales shall develop faster than grids at smaller scales (Eq 29).
Importantly, the formation of grid-like patterns also requires a nonlinearity in the system. Indeed, for triangular lattices to emerge, only three wave vectors k of the same length |k| shall survive. But this cannot be achieved in a linear system where all Fourier modes develop independently from each other (Eq 27). Yet the non-linear weight constraints imposed in our model (Eq 6) are sufficient to generate triangular patterns (Sec Emergence of grid spatial patterns).
In summary, the theory presented here gives necessary conditions for spatial pattern formation, and it predicts how the shape of the adaptation kernel K influences the scale of the grids and the relative time required for their formation. The theory remains however agnostic about the specific two-dimensional periodicity of the resulting patterns, i.e., it cannot predict whether the final solutions are, e.g., planar waves, square, rhomboidal, or triangular lattices. Further mathematical insights on this topic could be obtained by using perturbation methods [see e.g. 89], but this is beyond the scope of the current manuscript.

Numerical results on grid-pattern formation
In Sec Mathematical results on grid-pattern formation we derived an equation for the average dynamics of the synaptic weights w i , under the STDP learning rule and the stochastic activation of input and output neurons (Eq 16). In the case of spatially-regular inputs, we then computed the systems eigenvalue spectrum λ(k) in terms of the Gaussian input tuning curve G and the temporal adaptation kernel K (Eq 32). We showed that periodic spatial patterns could emerge if the eigenvalue spectrum λ(k) had a global maximum λ max > 0 at a frequency k max > 0 (Fig 4).  Fig 4), and Gaussian input receptive fields of size σ = 6.25 cm (Eq 9), the eigenvalue spectrum peaks at the critical frequency k max = 3 m −1 . In the following, we simulate the emergence of grid-like patterns in this scenario.
Emergence of grid spatial patterns. First, we simulate the detailed spiking model with spatially-regular inputs (Sec Spatially-regular inputs). The results are shown in Fig 5B-5E. In line with the theory, a structure emerges in the synaptic weights (Fig 5B and 5C) on a time scale of τ str = 1/(ηλ max ) % 5 Á 10 4 s (Eq 29) where η = 2 Á 10 −5 is the learning rate and λ max % 1 s −1 is the largest eigenvalue in the system. Additionally, the weight spectrum is quickly dominated by the critical frequency k max = 3 m −1 (Fig 5C, bottom row) at which the eigenvalue spectrum has a global maximum (Fig 5A).
Importantly, the synaptic weights also develop a periodic triangular symmetry, which is reminiscent of grid-cell patterns. Such triangular symmetry emerges after a substantial fraction of weights has hit the low saturation bound (Eq 6, Fig 5B, inset). Periodic pattern formation is indeed a strictly non-linear phenomenon, and excluding the spike generation process, weight saturation is the only non-linearity present in the system. In the linear regime, all Fourier modes k with frequency |k| = k max exponentially grow with equal rate ηλ max and independently from each other (Eq 25). In this case, the random weight pattern at time t = 0 s is amplified at the frequency k max , but no periodic structure emerges. In the non-linear regime, instead, the exponentially growing modes are mutually coupled, and a spontaneous symmetry breaking occurs: only three Fourier modes with wave vectors that are 60 degrees apart survive in our simulations (see Fig 5C and 5D).
In the example of Fig 5, a triangular symmetry starts to emerge after 2 Á 10 5 s, i.e., about 50 hours of exploration of the virtual rat. Yet the time scale of learning depends on the learning rate η and on the largest eigenvalue λ max (Eq 29), which are under-constrained by experimental data (see Eq 32 for the dependence of λ max on other model parameters). From a theoretical standpoint, the speed of learning is limited by the noise in the system, which is due to the virtual-rat random walk and the stochastic spiking of the neurons. To theoretically explore this limit and test the robustness of the model against noisy initial conditions, we simulated the development of the synaptic weights for different values of the learning rate and multiple random initializations of the synaptic weights. The results are reported in Fig 6A. With larger learning rates, grid-like patterns emerge faster. However, if the learning rate is too large, e.g.,  https://doi.org/10.1371/journal.pcbi.1005782.g006 η = 10 Á 10 −5 in our simulations, the gridness score fluctuates at low levels and no stable grid pattern emerges (yellow line in Fig 6A). Therefore, our results suggest that tens of hours of spatial exploration are required for stable grid patterns to emerge. Finally, the model is robust to random initializations of the synaptic weights, and to variations of the running speed of the virtual rat (Fig 6B).
Geometrical properties of the grid patterns. We now discuss the geometrical properties of the simulated grid patterns. A periodic triangular grid is characterized by three fundamental properties: i) the grid scale, i.e., the distance between two neighboring peaks; ii) the grid spatial phase, i.e., the spatial offset of the grid peaks with respect to a reference point; and iii) the grid orientation, i.e, the angle between one of the three grid axes and a reference direction.
Grid scale. In our model, the grid scale is set by the critical frequency k max at which the eigenvalue spectrum has a global maximum (Eq 32 and Fig 5). This critical frequency depends only on the movement speed v of the virtual rat, the width σ of the input tuning curve G, and the temporal dynamics of the adaptation kernel K (Fig 4). Therefore, grid patterns at different scales are obtained, for example, by varying the width σ of the input receptive fields or the long time scale τ L of the adaptation kernel (Fig 7, see also Fig 4). This theoretical result is consistent with the facts that spatial tuning in the hippocampal formation is typically broader ventrally than dorsally [90][91][92], and that grid scales vary in the same direction [1,93]. Additionally, we predict that the adaptation time scale may also have a dorso-ventral gradient, similarly to other intrinsic cellular properties in the mEC [e.g. 71, [94][95][96][97].
Grid spatial phase. With evenly-distributed input fields and periodic boundaries, the spatial phases of the grid patterns depend only on the initial condition of the synaptic weights, i.e., random weight initializations result in uniformly-distributed grid phases (Fig 8A1 and 8B1). This result is in line with the phases of nearby grid cells being roughly evenly distributed in experimental data [1], but see also [98]. Yet it remains unclear whether the same results would be obtained in the case of non-periodic boundaries.
Grid orientation. With periodic boundary conditions, our model produces grid orientations that are distributed non-uniformly. Precisely, the distribution of grid orientations depends on the scale of the pattern relative to the size of the environment, e.g., in the same environment patterns at different scales tend to align differently (compare panels A2 and B2 in sized environments, or in environments that are much larger than the pattern size. Nevertheless, because grid orientation depends on the boundary conditions, it remains difficult to compare the distributions obtained here with the ones observed experimentally [1,93,99,100]. Finally, in order to explain grid alignment across cells and/or environments [1,93], collateral interactions between developing grid cells may be required [50, 51, 101] (see also Sec Recurrent dynamics).
Pattern formation with spatially-irregular inputs. We demonstrated the emergence of grid-like patterns in the case of spatially-regular inputs, i.e., for each input cell having a single Gaussian receptive field in space (Sec Spatially-regular inputs). We now show that similar results are obtained in the case of spatially-irregular inputs (Sec Spatially-irregular inputs). We generate spatially-irregular inputs C in i by superimposing M > 1 Gaussian receptive fields with equal width σ, but random centers and random amplitudes (Eq 10, see Fig 9A for examples). The functions C in i are normalized such that their average firing rate r av is constant for all input neurons and independent from the number M of superimposed receptive fields.
We test grid-pattern formation in this scenario by simulating the average dynamics of the synaptic weights (Eqs 16 and 21) for random realizations of the input tuning curves C in i , with N = 3600 input neurons and M = 10 receptive fields per neuron. We then estimate output firing-rate maps from the synaptic weights at the end of the simulations (t = 10 6 s). The results are shown in Fig 9B and 9C. In the majority of the cases (73/100) a regular grid-like pattern emerges at the output.
Like in the case of spatially-regular inputs, the spatial scale of the output patterns depends on the long adaptation time constant τ L and on the width σ of the input receptive fields. Indeed, for σ = 6.25 cm and τ L = 0.16 s, we obtain output grid patterns with spatial frequency k max = 3 m −1 (Fig 9D), which is equal to the one obtained for spatially-regular inputs with the same parameter values (Figs 5 and 7, bottom-left panel). This can be understood by the fact that the expected eigenvalue spectrum for spatially-irregular inputs hλ irr (k)i is qualitatively similar to the eigenvalue spectrum λ(k) for spatially-regular inputs (Sec Eigenvalue spectrum for spatially-irregular inputs, Eq 72): where 0 F 1 is a scale factor. We also find that the scale factor F depends on the number M of superimposed fields, i.e., F % 4/(3M) for M > 3 (Eq 82), meaning that structure formation is slower for larger numbers of superimposed fields (Eq 29). Finally, like in the case of spatially-regular inputs, with periodic boundary conditions the spatial phases of the simulated grids distribute evenly in the arena (Fig 9E), and the grid orientations tend to cluster according to the grid scale and the size of the environment (Fig 9F,

Discussion
We studied the origin of grid-cell patterns in a single-cell spiking model relying solely on 1) spatially-tuned feed-forward inputs, 2) spike-rate adaptation, 3) and synaptic plasticity at the input synapses. We considered two input scenarios: spatially-regular inputs (reminiscent of place-cell activity), and spatially-irregular inputs (reminiscent of parasubicular activity). First, we studied the average dynamics of the system analytically, and we derived necessary conditions for the emergence of spatially-periodic solutions (Sec Mathematical results on grid-pattern formation). We then simulated the model numerically, and showed that grid-like patterns emerge both with spatially-regular and spatially irregular inputs (Sec Numerical results on grid-pattern formation). In the following, we discuss the main assumptions and predictions of our model.

Input spatial tuning and the origin of grid-cell patterns
We assumed that the feed-forward input activity is spatially tuned. Such spatial tuning could be provided by hippocampal place cells, or by other cortical or sub-cortical structures with less regular spatial firing. From a theoretical point of view, we find that grid patterns emerge faster with place-cell-like inputs, i.e., with inputs having a single receptive field in space. From an anatomical point of view, both scenarios seem plausible. On the one hand, grid-cell activity requires excitatory drive from the hippocampus [102], which projects to the deep layers of the mEC [103,104] where grid cells are found [23, 60]. On the other hand, parasubicular inputs target layer II of the mEC [61, [105][106][107][108] where grid cells are most abundant [23,60]. Although a small fraction of parasubicular cells already shows grid-like tuning [60, 61], the activity in parasubiculum is often characterized by multiple spatially-irregular fields [57-61] similar to those assumed in our model (Fig 9A).
That grid-cell activity could originate from parasubicular inputs is further supported by the detailed layout of the entorhinal circuit. Layer II principal neurons segregate into stellate and pyramidal cells, which are distinguished by their morphology, intrinsic properties [69], and immunoreactivity [109][110][111]. Interestingly, pyramidal-cell somata cluster into anatomical patches [110,111], which are preferentially targeted by parasubicular axons [61]; and the spiking activity in parasubiculum precedes the activity of layer II pyramidal cells by a few degrees in the theta cycle [61]. Such a network configuration suggests that grid patterns may originate in the layer II pyramidal cells via parasubicular inputs, and be inherited by the stellate cells via feed-forward projections. Consistent with this view is that both stellate and pyramidal cells show grid spatial tuning [55], and that direct intra-laminar connections are found from pyramidal onto stellate cells and not vice-versa [112,113]; but see [114].
In summary, our model is consistent with entorhinal grid-cell activity originating either in the superficial layers via parasubicular input or in the deep layers via hippocampal input. It is also possible that multiple sites of origin exist, and that grid-like tuning is inherited-and even sharpened-via feed-forward projections from the deep to the superfical layers [115][116][117][118][119] or from the superficial to the deep layers [104].

Spike-rate adaptation
Our grid-cell model relies on the presence of a spike-rate adaptation mechanism. Spike-rate adaptation has been observed throughout the cortex [62], and is prominent in layer II of the mEC, in both stellate and pyramidal neurons [69,70]. Yoshida et al. [71] also reported a dorso-ventral gradient in the adaptation strength of layer II entorhinal cells. However, because adaptation was found to be stronger ventrally than dorsally, Yoshida et al. [71] interpreted their results as evidence against grid-cell models based on adaptation. Yet the critical variable controlling the grid scale is not the strength of adaptation, but rather its temporal dynamics (Fig 7), which was not systematically analyzed [71]; see also [101] for a similar discussion on this point.
We modeled spike-rate adaptation by applying a temporal kernel K to the input spike trains (Eq 1). The kernel K, was composed of a brief depolarization peak and a slower hyper-polarizing potential (on a time scale of hundreds of milliseconds). Such a slow hyper-polarizing potential reduced the output firing rate in response to persistent excitation, and it filtered the input activity in a low-frequency band (i.e. with a resonance frequency of about 1 Hz, see Fig  1). The shape of the kernel was motivated by long-lasting hyper-polarizing potentials following excitatory post-synaptic potentials found in hippocampal CA1 pyramidal neurons [120], although similar responses have not been observed in the mEC yet.
However, the formation of grid-cell patterns could rely on any other cellular or synaptic mechanism that effectively acts as a band-pass filter on the input activity. A candidate mechanism is the after-spike hyperpolarizing potential (AHP). AHPs are indeed observed in the superficial layers of the mEC where single action potentials are followed by both a fast (2-5 ms) and a medium AHP (20-100 ms) [69,97,121]. To assess whether such hyperpolarizing potentials could underlie grid-pattern formation, we extended our model to account for AHPs (Sec Pattern formation with after-spike potentials). However, we found that grids at typical spatial scales cannot be obtained by AHPs alone. Yet after-spike potentials could amplify the effects of a band-pass filtering mechanism that is already present at the input.
Spike-rate adaptation could also rely on hyperpolarization-activated cation currents (I h ), which depend on HCN channels [122,123]. Fast I h currents (mediated by HCN1 channels) have been shown to control the theta-frequency resonance of entorhinal stellate cells in vitro [95,97,[124][125][126]. Instead, slower I h currents (mediated by HCN2-4 channels) could generate in entorhinal cells the low-frequency resonance assumed by our model (Fig 1B).

Synaptic plasticity
We propose that grid-cell patterns emerge from a synaptic reorganization of the mEC network, which is assumed to be plastic. This is in line with both LTP and LTD being reported in the entorhinal cortex [121,[127][128][129][130], but see also [131]. Additionally, asymmetric STDP was observed in the mEC [76]. Although we used a symmetric learning window in our model, the exact window shape has little effect on grid-pattern formation, provided that its temporal width (on the order of tens of milliseconds) is much shorter than the correlation length of the input activities (on the order of hundreds of milliseconds).
Structure formation via Hebbian learning is typically a slow process. In our model, grid-like patterns emerge on a time scale that is inversely proportional to the learning rate η and to the maximal eigenvalue λ max (Eq 29). The latter depends on the spatial density ρ = N/L 2 of input receptive fields, on the integral W tot of the learning window, on the shapes of the input-tuning curves G, and on the dynamics of the adaptation kernel K (Eq 32). Because most of these quantities are under-constrained by empirical data, a direct comparison with experimental time scales remains difficult. Yet learning shall be slow enough such that the input correlations that drive structure formation dominate over random fluctuations of the synaptic weights, which are due to the random walk of the virtual rat and the shot noise of the stochastic spiking. In our simulations, we find that this requires tens of hours of spatial exploration (Fig 6A).
Such slow process may seem in contrast with grid-cell activity appearing immediately in a novel environment [1,132]. However, grid-like tuning may not need to be learned in each environment anew, but rather recalled-and possibly refined-from the experience of similar environments explored in the past. Although hippocampal place cells [133,134] and entorhinal non-grid spatial cells [56] seem to remap completely in novel spaces, pattern formation could still leverage on residual correlations across environments that are hardly observable from the simultaneous recordings of only a few tens of neurons. Additionally, grid-cell learning could generalize across spatial contexts through border and boundary-vector inputs [84,85], which are invariant across environments.
We suggest that a structure in the synaptic weights may be formed during the animal's ontogenetic development, i.e., within a two-week period after the animal leaves the nest [135][136][137]. Consistent with this hypothesis is that stable spatial firing is observed before grid-cell maturation, e.g., hippocampal place cells develop prior to grid cells [135,136] and irregular spatial cells are present before grid cells [137].

Recurrent dynamics
We studied the emergence of grid patterns in a purely single-cell model, ignoring any network-level interaction between the neurons. However, because excitatory and inhibitory recurrent circuits have been described in the mEC [19,20,112,113,138], grid cells are likely to be mutually coupled [139,140]. Such recurrent connections could explain the modular organization of grid-cell properties [93,101] and their coherent responses to environmental changes [139]. Feedback interactions within a module may also amplify an initially broad gridtuning given by the feed-forward inputs, similarly to the sharpening of receptive fields in visual cortex [141,142]. Finally, recurrent dynamics may sustain grid-like activity when the feed-forward inputs are temporally untuned, like in attractor models [14]. Still, spatially-tuned feedforward inputs could be required for the initial formation of grid-like patterns [see e.g. 21].

Related models
Our work-and the one by Kropff and Treves [49]-belong to a broad category of grid-cell models based on spatially-tuned feed-forward inputs and Hebbian synaptic plasticity [63][64][65][66][67]. In all these models, periodic spatial patterns arise via a common underlying principle: the input correlations that drive the dynamics of the synaptic weights have the form of a Mexicanhat kernel (Fig 3). What distinguishes the models among each other-and generates distinct predictions-is the specific mechanism by which such Mexican-hat interactions are obtained.
In our model, a Mexican-hat kernel results from the intrinsic adaptation dynamics of the output neuron, which controls the grid scale directly (Fig 7).
By contrast, in the models by Castro and Aguiar [63] and Stepanyuk [64], Mexican-hat correlations arise from the learning rule itself, i.e., by assuming that synaptic plasticity switches between LTP and LTD based on pre-and post-synaptic activities [143]. In this case, the grid spatial scale shall be affected by interfering with the learning rule.
In a different model, Dordek et al. [65] obtain Mexican-hat correlations by constraining the input activity to be effectively zero-mean. The authors discuss that such a zero-mean constraint could originate either from lateral inhibition or from a zero-mean temporal filter controlling the output activity of the neuron. In the latter case, the model by Dordek et al. [65] is analogous to the present one. We note, however, that effectively zero-mean inputs are neither necessary nor sufficient for grid patterns to emerge. Instead, pattern formation depends on the dynamics of the temporal filter and on the shape of the input tuning curves, but not on their means. This can be easily understood by considering the system's eigenvalue spectrum in Fourier space (Eq 30), where the zero-frequency mode (k = 0) is not relevant for the emergence of spatially-periodic patterns. Also note that the smallest grid scales in our model are obtained with negative-mean temporal filters (Fig 4). Yet our results agree with the ones of Dordek et al. [65] in that the non-linearity introduced by imposing non-negative synaptic weights is sufficient for a triangular symmetry to emerge.
Alternatively, Mexican-hat correlations could emerge from phase-precessing feed-forward inputs [66]. In this case, grid-cell activity shall be impaired when phase precession is disrupted.
Finally, Weber and Sprekeler [67] proposed a model where the interplay between spatiallynarrow feed-forward excitation and spatially-broad feed-forward inhibition generates a Mexican-hat kernel. This model predicts that the grid scale shall be affected by manipulating inhibitory inputs to the mEC.

Model predictions and conclusion
We presented a single-cell model for the origin of grid-cell activity based on Hebbian synaptic plasticity and spike-rate adaptation. Our work builds upon the model by Kropff and Treves [49] and improves its original formulation in several aspects: 1) grid-like patterns emerge form a purely single-cell mechanism independently of any network-level interaction; 2) neuronal activities are spike-based and stochastic; 3) the input synaptic weights are purely excitatory; 4) the dynamics of the synaptic weights is studied analytically and linked to classical Turing-like patterns.
The present model makes the following experimental predictions. First, grid-cell patterns shall be affected by disrupting synaptic plasticity during ontogenetic development, which is consistent with preliminary data from Dagslott et al. [144]. Second, adult grid-cell activity shall be influenced by systematic behavioral or environmental biases in the first weeks of spatial exploration, e.g., by rising animals in environments without boundaries or with non-zero surface curvature [52,145]. Third, the grid scale shall be affected by three factors: 1) the spatial tuning-width of the feed-forward inputs; 2) the average speed of the rat during ontogenetic development; 3) the time constant of the recovery from spike-rate adaptation. Fourth, grids at larger scales shall develop faster as compared to grids at smaller scales (Fig 4).
We believe that manipulations of the intrinsic adaptation properties of single cells are key to distinguish our model from other feed-forward models based on Hebbian learning (Sec Related models). To this end, further experimental work shall be devoted to pinpoint the biophysical mechanisms underlying adaptation in the mEC. Extensions of the present model could also explain how the geometry of the enclosure affects grid-cell symmetry [99], and how grid-like tuning emerges in non-spatial contexts [146,147].
To conclude, our study contributes to a better understanding of the fundamental principles governing grid-cell activity, and lays the groundwork for more biophysically-realistic grid-cell models.

Weight normalization
Here we derive the dynamics of the mean synaptic weight w av ¼ N À 1 P N i¼1 " w i for a neuron with N synapses and temporally-averaged weights " w i . We recall the weight dynamics in Eq 16 By taking the average over the index i at both sides of Eq 34 we obtain where we defined the mean correlation C av ≔ N −2 ∑ ij C ij . Note that we used the property ∑ j C ij = NC av for all i, which holds true for translation-invariant inputs. Therefore, for NC av < a, the mean weight w av decays exponentially with time constant to the normalization level

Input correlation for general inputs
In this section we estimate the input correlation matrix for general spatial tuning curves C in i and smooth movement trajectories of the virtual rat (Sec Model of spatial exploration). We start by computing the temporal average r in i ðtÞr in j ðt À tÞof the product between the input activities r in i ðtÞ and the delayed input activities r in j ðt À tÞ. We assume that the stochastic process X t controlling the virtual-rat trajectory (Eq 11) is ergodic with respect to the auto-covariance, i.e., where the angular brackets denote statistical expectation. By using this ergodicity property (Eq 39) and the spatial tuning of the inputs (Eq 7), we derive Note that Eq 40 is only valid in an approximate sense because Eq 39 assumes T ! 1, but the averaging time window has finite length T ( τ str where τ str is structure-formation time constant (Eq 29). From Eq 40 follows where the integrals in Eqs 42-44 run over all positions in the environment (a square arena of side-length L), and p(x, t, x 0 , t − τ) is the joint probability density of the virtual rat being at position x at time t and at position x 0 at time t − τ. From Eqs 43 to 44, we used the fact that, for large times t, the virtual rat has equal probability of being in any position x, i.e., p(x, t) = 1/L 2 . Eq 44 shows that the temporal average r in i ðtÞr in j ðt À tÞcan be estimated from the input tuning curves C in i and C in j , and the conditional probability density p(x 0 , t − τ|x, t). This conditional probability density has not yet been solved for correlated random walks in two dimensions [148]. Nevertheless, an additional approximation is possible. Because the temporal average r in i ðtÞr in j ðt À tÞis weighted by the adaptation adaptation kernel K(τ) (Eq 38), and K(τ) is negligible for τ > τ max % 5τ L (Eq 2), we are interested in the conditional probability p(x 0 , t − τ|x, t) only at lags τ < τ max . In this case, for movement trajectories that are sufficiently smooth, we can assume that in a time τ the virtual rat has moved to a position x at distance |x − x 0 | = τv from the initial position x 0 , that is where v is the speed of the virtual rat (Eq 11), and the denominator ensures that R dx 0 p(x 0 , t − τ|x, t) = 1; see also Fig 2 for exemplary virtual-rat trajectories in this scenario. We now use Eq 45 in Eq 44, and let z ≔ x 0 − x: From Eq 46, the temporal average r in i ðtÞr in j ðt À tÞis approximated by the integral of the spatial cross-correlation C in i ? C in j over a circle of radius τv. Finally, by using Eq 46 in Eq 20, we obtain Input correlation for spatially-regular inputs In this section we compute the input correlation function C and its Fourier spectrumĈ in the case of spatially-regular inputs (see Sec Weight dynamics for spatially-regular inputs). First, we rewrite the input correlation matrix C ij in Eq 21 as a continuous function C (r, r 0 ) by labeling neurons according to their receptive-field centers r and r 0 : where C in r ðxÞ ≔ Gðjx À rjÞ is a Gaussian input tuning curve centered at position r (Eq 9). Because the inputs are translation invariant, the correlation function C depends only on the translation vector u ≔ r − r 0 : where C in 0 ðxÞ ≔ GðjxjÞ is the tuning curve centered at the origin 0 = (0, 0). Next, we substitute in Eq 50 the definition of the integral operator in Eq 46: It is easy to see that the auto-correlation of a Gaussian is still a Gaussian: from which we derive where φ is the angle between the vectors u and z. Finally, by expressing in polar coordinates the vector z ≔ |z|[cos(φ), sin(φ)], from Eqs 51 and 53 we obtain where I 0 ðxÞ ≔ 1=ð2pÞ R 2p 0 dφ expðx cos ðφÞÞ is the zeroth-order modified Bessel function of the first kind.

Fourier spectrum of the input correlation function.
Here we compute the Fourier spectrum of the correlation function C in Eq 51. First, we observe that the second integral in Eq 51 is a two-dimensional cross-correlation in the variable z between the functions δ(|z| − τv) and C in 0 ? C in 0 j z evaluated at point u. Therefore, by taking the two-dimensional Fourier transform with respect to u at both sides of Eq 51 yieldŝ where we defined the Fourier transform pair: with k Á u = |k||u| cos(θ), and we used the definition of the zeroth-order Bessel function Because the tuning function C in 0 ðxÞ ≔ GðjxjÞ is circularly symmetric, its two-dimensional Fourier transformĈ in 0 ðkÞ is proportional to the zeroth-order Hankel transform of G: where we defined the zeroth-order Hankel transform pair: GðkÞ ≔ Z 1 0 dr r GðrÞJ 0 ðkrÞ and GðrÞ ¼ By using Eq 58 in Eq 55 we obtain and by defining the equivalent adaptation kernel in space we findĈ Finally, the zeroth-order Hankel transforms of the Gaussian tuning curve G (Eq 9) and of the adaptation kernel in space K sp (Eqs 61 and 2) read Eigenvalue spectrum for spatially-irregular inputs In this section we estimate the expected eigenvalue spectrum hλ irr (k)i for spatially-irregular inputs (Secs Spatially-irregular inputs and Pattern formation with spatially-irregular inputs). We recall that, for spatially-regular inputs, in Sec Mathematical results on grid-pattern formation we obtained (Eq 32): whereG andK sp are the zeroth-order Hankel transforms of the input tuning curve G (Eq 9) and of the equivalent adaptation kernel in space K sp (Eqs 31 and 61). Note that the parameters ρ, L, W tot , and a do not depend on k. From Eq 65, the eigenvalue spectrum λ(k) is linearlyrelated to the input power spectrum jĈ in 0 ðkÞj 2 where C in 0 ðxÞ ≔ GðjxjÞ is an input tuning curve centered at the origin 0 ≔ (0, 0) (Sec Input correlation for spatially-regular inputs).
Here, in analogy to Eq 65, we assume that the expected eigenvalue spectrum hλ irr (k)i for spatially-irregular inputs is linearly-related to the expected input power hjĈ in p ðkÞj 2 i, that is, whereĈ in p ðkÞ is the two-dimensional Fourier transform of the spatially-irregular tuning curve C in p ðxÞ, and the angular brackets denote statistical expectation across input realizations (see Eq 56 for a definition of the two-dimensional Fourier transform). The validity of this assumption is confirmed numerically at the end of this section.
Let us compute the expected input power spectrum hjĈ in p ðkÞj 2 i. We recall that the input maps C in p ðxÞ are obtained by the superimposing M Gaussian receptive fields (Eq 10) with GðrÞ ¼ ð9Þ L 2 r av 2ps 2 exp À r 2 2s 2 and b p ≔ The field amplitudes A pm ! 0 are uniformly distributed in the range (0, 1), and the receptive field centers r pm are uniformly distributed in the environment (see Fig 9A for examples). From Eq 67 we derive whereGðkÞ is the zeroth-order Hankel transform of the Gaussian function GðrÞ. In deriving Eq 69, we used the shift property of the Fourier transform and the equivalence between the Fourier and the zeroth-order Hankel transforms for circularly-symmetric functions (Eq 58).
Finally, from Eq 69 we obtain Therefore, for spatially-irregular inputs, the expected power spectrum hjĈ in p ðkÞj 2 i is proportional to the power spectrum 4p 2G 2 ðkÞ of a single Gaussian G with scale factor F ! 0. Note that for |k| = 0 we obtain F = 1 (Eqs 69 and 70), which means that the average rate r av is independent of the number M of input receptive fields and their specific spatial arrangement. Using Eq 70 in Eq 66 yields Finally, from Eqs 65 and 71 we find (Eq 33) In the next section we estimate the scale factor F for |k| > 0. Approximation of the scale factor F. The scale factor is the second moment of the ratio of the random variables where the field amplitudes A pm ! 0 are independently and uniformly distributed in the range (0, 1) and the field centers r pm are independently and uniformly distributed in a square of sidelength L.
In general, for two random variables x and y, the first order Taylor expansion of the ratio f(z) = x/y around the expected value μ ≔ (hxi, hyi) is where z ≔(x, y), Δ x ≔ x − hxi, Δ y ≔ y − hyi, and f x and f y are the derivatives of f with respect to x and y. Therefore In the following, we use Eq 78 to approximate the scale factor F (Eq 73). We start by giving an intuitive interpretation of the random variables α p and β p . Consider a M-steps random walk on the complex plane with random directions r pm Á k and random step sizes A pm . The coefficients α p measure the total distance traveled by the random walker, and the coefficients β p measure the total length of the path (Eq 74). Note that the larger the number of steps M, the smaller is the correlation between the distance traveled α p and the total path length β p , i.e., |Cov(α p , β p )| ( 1 for M ) 1. In this case we can neglect the covariance term in Eq 78, and the factor F is approximated by knowing only the first two moments of the distributions of α p and β p .
For |k| > 1/L, the random directions r pm Á k (mod 1) are approximately uniformly distributed in the range (0, 1). In this case, the traveled distance α p follows a Rayleigh distribution with density [149] where hA 2 pm i ¼ 1=3 for A pm uniformly distributed in interval (0, 1). Therefore, the first two moments of α p read The total path length β p is the sum of M random variables uniformly distributed in (0, 1), which follows an Irwin-Hall distribution. Therefore, the first two moments of β p are Finally, by using Eqs 80 and 81 in Eq 73 we obtain Pattern formation with after-spike potentials Here we study whether grid-like patterns could emerge by means of after-spike hyperpolarizing potentials (see discussion in Sec Spike-rate adaptation). To this end, we consider a model of the output neural activity that is alternative to the one presented in the main text (Sec Model of neural activity, Eq 1). We model input post-synaptic potentials (PSPs) with a kernel K in applied to the input spike trains S in j , and we model output after-spike hyperpolarizing potentials (AHPs) with a kernel K out applied to the output spike train S out : where r 0 ! 0 is a baseline firing rate. First, we show that the average dynamics of Eq 84 can be rewritten in terms of an equivalent kernel K eq applied to the input spikes only. We average Eq 84 across input and output spike train realizations: From Eqs 85 to 87 we assumed that the input and the output kernels are causal, i.e., K in,out (t) = 0 for t < 0, and that the output kernel has integral different from 1, i.e., which is equivalent to Eq 15 with K eq = K. Next, we compute the equivalent filter K eq for a simple choice of the input and output kernels and where τ in , τ out > 0 are decay time constants, and the parameter μ out > 0 scales the integral of the output kernel R 1 0 dt K out ðtÞ ¼ À m out . We assume that the input kernel K in (modeling an incoming PSP) decays faster than the output kernel K out (modeling an output AHP), i.e., τ in < τ out . From the definition of the filter K eq in Eq 88 we obtain where we usedK Finally, the inverse Fourier transform of Eq 92 reads for t ! 0 and K eq (t) = 0 for t < 0. Eq 94 shows that the equivalent filter K eq is a difference of two exponentials, similarly to the kernel K in Eq 2. Note however that the two exponentials are scaled differently as compared to the original filter K. Additionally, if the integral of the output kernel is negative, the integral of the equivalent filter is always positive (Eq 88 with ω = 0). To test whether spatially-periodic patterns could still emerge in this scenario, we compute the eigenvalue spectrum λ(k) and the critical spatial frequency k max by using Eqs 31 and 32 with K = K eq . Surprisingly, we find that typical grid scales (e.g., k max > 2 m −1 ) are obtained for output-kernel time constants of the order of seconds, which seem biologically unrealistic ( Fig  11). Therefore, we conclude that AHPs alone are not sufficient to generate grid-like patterns.
Nevertheless, AHPs could still support structure formation by amplifying the effects of a bandpass filter that is already present at the input.

Numerical simulations
Model parameters and derived quantities are summarized in Tables 1 and 2.
Simulation of the detailed spiking model. The detailed spiking model (Figs 5 and 6) is simulated using the Brian2 simulation software [150]. Neural and synaptic variables are integrated with a time step of 1 ms. The random walk of the virtual rat that is updated every 10 ms. The physical space explored by the virtual rat is discretized in 200 2 square bins.
Simulation of the averaged weight dynamics. The average weight dynamics (Eq 16) is integrated by using the forward Euler method with integration time step of 50 s (Figs 7-9). The input correlation matrix C is computed using Eq 54 for spatially-regular inputs, and using Eq 21 for spatially-irregular inputs.
Initialization of the synaptic weights. At the initial condition the synaptic weights are normally distributed around the target normalization level w 1 av ¼ 5 Á 10 À 3 . The standard deviation is 10 −4 for the spiking simulations and 10 −3 for the average weight dynamics.

Data analysis
Grid properties. We compute the grid spatial scale from the two-dimensional Fourier amplitude of the grid pattern. We estimate the radial amplitude profile by averaging over the angular dimension. We then define the grid scale as the frequency where the amplitude profile has a global maximum.
The grid orientation is estimated from the spatial auto-correlogram of the grid pattern. We detect the peak closest to the center in the first quadrant of the auto-correlogram. We then define the grid orientation as the angle between the detected peak and the horizontal axis. We define the grid spatial phase as the position of the closest peak to the center in the crosscorrelation between the grid pattern and a reference grid at the same scale.
Gridness score. We estimate the gridness score similarly to Langston et al. [135]. First, we compute the spatial auto-correlogram of the weight (or firing-rate) pattern and we retain only points within a ring of outer radius R i and inner radius R i /2. We then compute the gridness score g i as where ρ i (φ) is the Pearson's correlation coefficient between the original ring (of outer radius R i ) and the same ring rotated by φ degrees. The final gridness score is defined as the maximum g i by varying the outer radius R i between 0.7/k max and 2.5/k max where k max is the spatial frequency of the pattern. Estimation of output firing-rate maps. The output firing-rate maps C out in Fig 9B are computed as follows: where r 0 is the baseline firing rate, w i are the synaptic weights at the end of the simulation, C in i are the input spatial maps, and K sp is the equivalent adaptation kernel in space (Eq 31). The convolution with the filter K sp accounts for the average effect of the temporal kernel K on the output firing rate.  Fig 6A the learning rate η is varied from 2 Á 10 −5 to 10 Á 10 −5 and that in Fig 6B the virtual-rat running speed is sampled from an Ornstein-Uhlenbeck process with long-term mean " v ¼ v.