Model of neural activity

We consider a single cell that receives synaptic input from N spatially-tuned excitatory neurons. Input spike trains for i = 1, 2, …, N are generated by independent inhomogeneous Poisson processes with instantaneous rates where δ(t) is the Dirac delta function, and is the timing of the k^th input spike at synapse i. Similarly, the output spike train is generated by an inhomogeneous Poisson process with instantaneous rate r^out(t) where 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, (1) where r₀ 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: (2) where τ_S and τ_L are the short and long filter time constants (0 < τ_S < τ_L), and the parameter μ > 0 sets the filter integral (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–71]. From a signal-processing perspective, the adaptation kernel K performs a temporal band-pass 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. The temporal adaptation kernel K.

(A) Impulse response of the filter (Eq 2). A positive peak with amplitude K(0) = 1/τ_S − μ/τ_L ≈ 3.4 spikes/s is followed by a slow negative response. Note that the kernel is small for t > τ_max, i.e., |K(t)| < 0.01|K(0)| for t > τ_max, with τ_max = 5τ_L = 0.8 s. (B) Frequency response of the filter. The dashed vertical line indicates the filter’s resonance frequency k_res = 1.23 s⁻¹. Parameter values: τ_S = 0.1 s, τ_L = 0.16 s, μ = 1.06. The integral of the filter is 1 − μ = −0.06.

https://doi.org/10.1371/journal.pcbi.1005782.g001

Model of synaptic plasticity

We assume spike-timing dependent plasticity (STDP) at the input synapses [e.g. 72–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 (3)
2. For each pre-synaptic spike at synapse i, we set (4)

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] (5) where the time constant τ_W > 0 controls the maximal time lag at which plasticity occurs, and 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 and prevents the individual weights to grow unbounded. This normalization term mimics local homoeostatic processes observed experimentally [78–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 (6)

Model of input spatial tuning

We consider excitatory inputs with firing rates that are tuned to the spatial position of a virtual rat exploring a square arena of side-length L, i.e., (7) where x_t is the position of the virtual rat at time t, and is a spatial tuning curve. We characterize the spatial tuning curves 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.

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 (11) where the angle θ_t sets the direction of motion and 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 .

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Example virtual-rat trajectories.

Colored lines denote example virtual-rat trajectories obtained by integrating Eq 11 starting at the center of the gray disk. Filled dots indicate the position of the virtual rat at time τ_max = 5τ_L = 0.8 s. Note that the trajectories are smooth within time stretches shorter than τ_max. Parameter values: v = 0.25 m/s, θ_σ = 0.7. The disk radius is vτ_max = 20 cm.

https://doi.org/10.1371/journal.pcbi.1005782.g002

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] (12) with . The functions 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., . Following Kempter et al. [74] we derive (13) where the ensemble averages read (14) (15) Finally, from Eqs 12–15 we obtain (16) where we defined (17) (18) (19) Note that in deriving Eq 16 we approximated the temporal average of the input rates with the spatial average r_av of the input tuning curves . 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 , 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 and , 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 (Eq 17), which holds for τ_W ≪ σ/v. In this case, we obtain (20) 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 and the adaptation kernel K (Sec Input correlation for general inputs, Eq 47): (21) where L² 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 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 .

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 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: (22) where 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², the sum in Eq 22 can be replaced by an integral over all the receptive-field centers r′: (23) where the correlation function C(r, r′) 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′) = C(r − r′, 0) = C(r − r′), where we omit the second argument 0 ≔ (0, 0) for readability. In this case, the integral in Eq 23 can be expressed as a two-dimensional convolution in space: (24)

Fig 3 shows the correlation C as a function of the input receptive-field distance |r − r′| 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-field distances (attraction domain), negative for intermediate distances (repulsion domain), 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].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Input correlation function C for spatially-regular inputs.

The function is circularly symmetric, i.e., it depends only on the distance |r − r′| between the receptive-field centers r and r′ (Eq 54). In the attraction domain (red shaded area) the correlation is positive and the synaptic weights grow in the same direction. In the repulsion domain (blue shaded area) the correlation is negative and the synaptic weights grow in opposite directions. Parameter values: σ = 6.25 cm, r_av = 0.4 s⁻¹, τ_S = 0.1 s, τ_L = 0.16 s, μ = 1.06, W_tot = 1 s, L = 1 m, v = 0.25 m/s.

https://doi.org/10.1371/journal.pcbi.1005782.g003

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: (25) where we defined the Fourier transform pair (26) k is a two-dimensional wave vector, and is the imaginary unit. Solving Eq 25 for k ≠ (0, 0), we obtain (27) where denotes the weight spectrum at time t = 0, and we defined (28) 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 grow or decay exponentially with rates proportional to the eigenvalues λ(k). Therefore, a structure in the synaptic weights emerges on a time scale (29) 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 input-correlation function (Eq 28), which is circularly-symmetric for circularly-symmetric inputs. In this case, in Sec Input correlation for spatially-regular inputs (Eq 62) we derive (30) where and (Eqs 63 and 64) are the zeroth-order Hankel transforms (Eq 59) of the input tuning curve (Eq 9) and of the equivalent adaptation kernel in space (31) Finally, by plugging Eqs 30 into 28, we obtain (32)

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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Impact of the adaptation kernel on grid-pattern formation.

(A1-A2) Critical spatial frequency k_max (A1) and largest eigenvalue λ_max (A2) as a function of the kernel integral 1 − μ and the long kernel time constant τ_L. The short time constant is τ_S = 0.1 s. The black lines are iso-levels (see annotated values). Regions enclosed by two adjacent iso-lines are colored uniformly (darker colors denote larger values). Within the black region in A1 we obtain λ_max ≤ 0 s⁻¹ (see white region in A2). Within the black region in A2 we obtain k_max = 0 m⁻¹ (see white region in A1). The dashed horizontal line indicates zero-integral kernels. The star denotes the parameter values τ_S = 0.1 s, τ_L = 0.16 s, μ = 1.06 of the kernel in Fig 1. (B1-B2) Same as in A but varying the short kernel time constant τ_S. The long time constant is τ_L = 0.16 s. The eigenvalue spectrum is estimated from Eq 32. Further parameter values: σ = 6.25 cm, r_av = 0.4 s⁻¹, W_tot = 1 s, ρ = 900 m⁻², L = 1 m, v = 0.25 m/s, a = 1.1 s⁻¹.

https://doi.org/10.1371/journal.pcbi.1005782.g004

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 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 5A shows the eigenvalue spectrum λ(k) for a choice of the parameter values such that this condition is satisfied. With adaptation time constants τ_S = 0.1 s and τ_L = 0.16 s (Eq 2, Fig 1, star in 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⁻¹. In the following, we simulate the emergence of grid-like patterns in this scenario.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Grid-pattern formation with spatially-regular inputs.

(A) Eigenvalue spectrum λ(k) of the averaged weight dynamics (Eq 32). The black solid line shows the continuous spectrum in the limit of infinite-size environments; the red dots show the discrete eigenvalues for a square arena of side length L = 1 m with periodic boundaries. The horizontal dashed line separates positive and negative eigenvalues. The vertical gray line indicates the critical spatial frequency k_max = 3 m⁻¹. The eigenvalue at frequency k = 0 is not shown. Parameter values: τ_S = 0.1 s, τ_L = 0.16 s, σ = 6.25 cm. (B) Time-resolved distribution of N = 900 synaptic weights updated according to the STDP rule in Eqs 3–6. Red triangles indicate the time points shown in C. Inset: fraction of weights close to the lower saturation bound (w_i < 5 ⋅ 10⁻³). (C) Top row: evolution of the synaptic weights over time. Weights are sorted according to the two-dimensional position of the corresponding input receptive-field centers. Note that each panel has a different color scale (maximum weight at the bottom-left corner, see B for distributions). Bottom row: Fourier amplitude of the synaptic weights at the top row. The red circle indicates the frequency k_max = 3 m⁻¹ of the largest eigenvalue (see panel A). (D) Time evolution of weights' Fourier amplitudes for wave vectors k at the critical frequency |k| = k_max. Wave vector angles (color coded) are relative to the largest mode at the end of the simulation (t = 10⁶ s). The black triangles indicate time points in C. (E) Gridness score of the weight pattern over time. The gridness score quantifies the degree of triangular periodicity. See Sec Numerical simulations for further details and parameter values.

https://doi.org/10.1371/journal.pcbi.1005782.g005

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⁴ s (Eq 29) where η = 2 ⋅ 10⁻⁵ is the learning rate and λ_max ≈ 1 s⁻¹ is the largest eigenvalue in the system. Additionally, the weight spectrum is quickly dominated by the critical frequency k_max = 3 m⁻¹ (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⁵ 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., η = 10 ⋅ 10⁻⁵ 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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Time scales of learning.

(A) Median gridness scores of the input synaptic weights for 40 random weight initializations and different learning-rate values, i.e., η = (2, 3, 5, 10) ⋅ 10⁻⁵. The weight development is simulated with the detailed spiking model with spatially-regular inputs and constant virtual-rat speed (see also Fig 5). (B) Median gridness scores of the input synaptic weights simulated with constant (black line) and variable (green line) virtual-rat speeds for 40 random weight initializations. Variable running speeds are obtained by sampling from an Ornstein-Uhlenbeck process with long-term mean m/s, volatility σ_v = 0.1 m ⋅ s^−1.5 and mean-reversion speed θ_v = 10 s⁻¹. The inset shows the distribution of running speeds (mean: 0.25 m/s std: 0.02 m/s). Note that the long-term mean of the process equals the speed v in constant-speed simulations. See Sec Numerical simulations for further details and additional parameter values.

https://doi.org/10.1371/journal.pcbi.1005782.g006

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 , 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–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–97].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Spatial scale of the grid patterns.

Example grid patterns obtained with different adaptation kernels K (Eq 2, top row) and different input tuning curves (Eq 9, left-most column). For each choice of the functions K and , the synaptic weights (left) and their corresponding Fourier spectra (right) at the end of the simulation are shown (t = 10⁶ s). The synaptic-weight maps have different color scales (maximal values at the bottom-left corner). The red circles indicate the spatial frequency k_max of the weight patterns. Synaptic weights were obtained by simulating the average weight dynamics in Eq 16. Note that we used a larger enclosure (L = 2 m) as compared to the one in Figs 5 and 6 (L = 1 m). See Sec Numerical simulations for further details and parameter values.

https://doi.org/10.1371/journal.pcbi.1005782.g007

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Geometric properties of the grid patterns.

(A) Distribution of grid spatial phases (A1) and grid orientations (A2) for patterns at frequency k_max = 3 m⁻¹ in an arena of side-length L = 2 m (σ = 6.25 cm, τ_L = 0.16 s; see also Fig 7, bottom-left panel). Distributions were obtained from the average weight dynamics in Eq 16 for 200 random initializations of the synaptic weights (t = 10⁶ s). Only patterns with gridness scores larger than 0.5 were considered (197/200). Panel A3 shows example weight patterns for the two most common orientations in A2 (maximal values at the bottom-left corner). (B) Same as in A but for patterns at spatial frequency k_max = 2 m⁻¹ in an arena of side-length L = 2 m (σ = 6.25 cm, τ_L = 0.35 s; see also Fig 7, bottom-right panel). A fraction of 182/200 grids had a gridness score larger than 0.5. See Sec Numerical simulations for further details and parameter values.

https://doi.org/10.1371/journal.pcbi.1005782.g008

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 Fig 8). In the examples of Fig 8A3, one of the grid axes tends to align to a border of the arena whereas in the examples of Fig 8B3 one of the grid axes tends to align to a diagonal of the arena. Similar results are obtained by keeping the grid scale fixed and varying the size of the environment, e.g., compare Fig 8A2 (k_max = 3 m⁻¹ and L = 2 m) and Fig 9F (k_max = 3 m⁻¹ and L = 1 m). In general, we expect grid orientations to be uniformly distributed only in infinite-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).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Grid-pattern formation with spatially-irregular inputs.

(A) Four examples of irregular input firing-rate maps (top row) and the corresponding Fourier spectra (bottom row). The maximal firing rate (spikes/s) is reported at the bottom-left corner. The red circles indicate the spatial frequency k_max = 3 m⁻¹. (B) Four examples of output firing-rate maps (top row) and the corresponding Fourier spectra (bottom row). The gridness score is reported at the bottom-right corner. Output firing-rate maps were estimated from the average weight dynamics in Eq 16 (t = 10⁶ s) for four different realizations of the spatial inputs. (C-F) Distribution of gridness scores (C), grid spatial frequencies (D), grid spatial phases (E), and grid orientations (F) for 100 random realizations of the spatial inputs. The red vertical line in C indicates the mean score (0.77). See Sec Numerical simulations for further details and parameter values.

https://doi.org/10.1371/journal.pcbi.1005782.g009

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–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–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–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–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–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² of input receptive fields, on the integral W_tot of the learning window, on the shapes of the input-tuning curves , 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–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 grid-tuning 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 feed-forward 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–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 Mexican-hat 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 spatially-narrow 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 for a neuron with N synapses and temporally-averaged weights . We recall the weight dynamics in Eq 16 (34) By taking the average over the index i at both sides of Eq 34 we obtain (35) where we defined the mean correlation C_av ≔ N⁻²∑_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 (36) to the normalization level (37)

Input correlation for general inputs

In this section we estimate the input correlation matrix (38) for general spatial tuning curves and smooth movement trajectories of the virtual rat (Sec Model of spatial exploration). We start by computing the temporal average of the product between the input activities and the delayed input activities . 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., (39) 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 (40) Note that Eq 40 is only valid in an approximate sense because Eq 39 assumes T → ∞, but the averaging time window has finite length T ≪ τ_str where τ_str is structure-formation time constant (Eq 29). From Eq 40 follows (41) (42) (43) (44) 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′, t − τ) is the joint probability density of the virtual rat being at position x at time t and at position x′ 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².

Eq 44 shows that the temporal average can be estimated from the input tuning curves and , and the conditional probability density p(x′, 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 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′, 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′| = τv from the initial position x′, that is (45) where v is the speed of the virtual rat (Eq 11), and the denominator ensures that ∫ dx′p(x′, 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′ − x: (46) From Eq 46, the temporal average is approximated by the integral of the spatial cross-correlation over a circle of radius τv. Finally, by using Eq 46 in Eq 20, we obtain (47)

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′) by labeling neurons according to their receptive-field centers r and r′: (48) where 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′: (49) (50) where 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: (51) It is easy to see that the auto-correlation of a Gaussian is still a Gaussian: (52) from which we derive (53) 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 (54) where is the zeroth-order modified Bessel function of the first kind.

Eigenvalue spectrum for spatially-irregular inputs

In this section we estimate the expected eigenvalue spectrum 〈λ_irr(k)〉 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): (65) where and are the zeroth-order Hankel transforms of the input tuning curve (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 linearly-related to the input power spectrum where 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 〈λ_irr(k)〉 for spatially-irregular inputs is linearly-related to the expected input power , that is, (66) where is the two-dimensional Fourier transform of the spatially-irregular tuning curve , 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 . We recall that the input maps are obtained by the superimposing M Gaussian receptive fields (Eq 10) (67) with (68) 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 (69) where is the zeroth-order Hankel transform of the Gaussian function . 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 (70) Therefore, for spatially-irregular inputs, the expected power spectrum is proportional to the power spectrum of a single Gaussian with scale factor Φ ≥ 0. Note that for |k| = 0 we obtain Φ = 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 (71) Finally, from Eqs 65 and 71 we find (Eq 33) (72) In the next section we estimate the scale factor Φ for |k| > 0.

Approximation of the scale factor Φ.

The scale factor (73) is the second moment of the ratio of the random variables (74) 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 side-length 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 μ ≔ (〈x〉, 〈y〉) is (75) where z ≔(x, y), Δ_x ≔ x − 〈x〉, Δ_y ≔ y − 〈y〉, and f_x and f_y are the derivatives of f with respect to x and y. Therefore (76) By neglecting the higher-order joint moments and substituting f_x(μ) = 1/〈y〉 and f_y(μ) = −〈x〉/〈y〉² we obtain (77) and (78) In the following, we use Eq 78 to approximate the scale factor Φ (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 Φ 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] (79) where for A_pm uniformly distributed in interval (0, 1). Therefore, the first two moments of α_p read (80)

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 (81)

Finally, by using Eqs 80 and 81 in Eq 73 we obtain (82) (83) Fig 10A shows the scale factor Φ as a function of the number M of superimposed Gaussian fields (Eq 82). Note that the approximation is more accurate for large values of M, which correspond to lower values of |Cov(α_p, β_p)|. Fig 10B shows the largest eigenvalue in the system as a function of M. The good match between the theoretical curve and the numerical estimations supports the validity of Eq 66. Additionally, Eq 66 predicts that, irrespectively of the value of M, the largest eigenvalue λ_max = λ(k_max) is always at the critical frequency of k_max = 3 m⁻¹ for σ = 6.25 cm and τ_L = 0.16 s, which matches the numerical results in Fig 9.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Scale factor Φ and largest eigenvalue λ_max for spatially-irregular inputs.

(A) The scale factor Φ for M > 1 superimposed fields (Eq 70). The black dots are obtained by estimating the power spectrum at frequency |k| = 1 m⁻¹ for 3600 input realizations. The red line is the theoretical curve in Eq 82. (B) The largest eigenvalue λ_max as a function of the number of superimposed fields M. The black dots are obtained by computing the eigenvalues of the correlation matrix C_ij − aδ_ij for N = 3600 inputs, where δ_ij is the Kronecker delta (Eq 21). The red line is obtained from Eqs 71 and 82. Note that, according to Eq 71, the largest eigenvalue is always at the critical frequency k_max = 3 m⁻¹ for any value of M. Parameter values as in Fig 9 (see Sec Numerical simulations).

https://doi.org/10.1371/journal.pcbi.1005782.g010

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ⁱⁿ applied to the input spike trains , and we model output after-spike hyperpolarizing potentials (AHPs) with a kernel K^out applied to the output spike train S^out: (84) where r₀ ≥ 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: (85) And by taking the Fourier transform (86) at both sides of Eq 85 we obtain (87) 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., . Finally, by defining the equivalent filter (88) the inverse Fourier transform of Eq 87 reads (89) 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 (90) and (91) where τ_in, τ_out > 0 are decay time constants, and the parameter μ_out > 0 scales the integral of the output kernel . We assume that the input kernel Kⁱⁿ (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 (92) where we used (93) Finally, the inverse Fourier transform of Eq 92 reads (94) 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⁻¹) 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 band-pass filter that is already present at the input.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Grid scale with after-spike hyperpolarizing potentials.

The critical spatial frequency k_max is plotted as a function of the output-kernel integral −μ_out and the output-kernel time constant τ_out (Eqs 31 and 32 with K = K^eq). The black lines are iso-levels (see annotated values). Regions enclosed by two adjacent iso-lines are colored uniformly (darker colors denote larger values). The input-kernel time constant is τ_in = 5 ms. Similar results are obtained with different values of τ_in < τ_out. Parameter values: σ = 6.25 cm, v = 0.25 m/s, L = 1 m. r_av = 0.4 s⁻¹.

https://doi.org/10.1371/journal.pcbi.1005782.g011

Numerical simulations

Model parameters and derived quantities are summarized in Tables 1 and 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Model parameters.

https://doi.org/10.1371/journal.pcbi.1005782.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Default parameter values for the numerical simulations.

See also Table 1 for short descriptions of the parameters. TL: top-left, TR: top-right, BL: bottom-left, BR: bottom-right. Note that in Fig 6A the learning rate η is varied from 2 ⋅ 10⁻⁵ to 10 ⋅ 10⁻⁵ and that in Fig 6B the virtual-rat running speed is sampled from an Ornstein-Uhlenbeck process with long-term mean .

https://doi.org/10.1371/journal.pcbi.1005782.t002

Data analysis