A Continuous Attractor Model with Realistic Neural and Synaptic Properties Quantitatively Reproduces Grid Cell Physiology

Computational simulations with data-driven physiological detail can foster a deeper understanding of the neural mechanisms involved in cognition. Here, we utilize the wealth of cellular properties from Hippocampome.org to study neural mechanisms of spatial coding with a spiking continuous attractor network model of medial entorhinal cortex circuit activity. The primary goal is to investigate if adding such realistic constraints could produce firing patterns similar to those measured in real neurons. Biological characteristics included in the work are excitability, connectivity, and synaptic signaling of neuron types defined primarily by their axonal and dendritic morphologies. We investigate the spiking dynamics in specific neuron types and the synaptic activities between groups of neurons. Modeling the rodent hippocampal formation keeps the simulations to a computationally reasonable scale while also anchoring the parameters and results to experimental measurements. Our model generates grid cell activity that well matches the spacing, size, and firing rates of grid fields recorded in live behaving animals from both published datasets and new experiments performed for this study. Our simulations also recreate different scales of those properties, e.g., small and large, as found along the dorsoventral axis of the medial entorhinal cortex. Computational exploration of neuronal and synaptic model parameters reveals that a broad range of neural properties produce grid fields in the simulation. These results demonstrate that the continuous attractor network model of grid cells is compatible with a spiking neural network implementation sourcing data-driven biophysical and anatomical parameters from Hippocampome.org. The software (version 1.0) is released as open source to enable broad community reuse and encourage novel applications.


Introduction
The hippocampus and surrounding circuits, including the entorhinal cortex, play a crucial role in representing and associating space, context, and episodic memory.Despite remarkable advances in characterizing the molecular and cellular properties of hippocampal neurons, a coherent framework to link this wealth of data with network-level computational functions is still missing.Simulating cells involved in spatial cognition with detailed levels of physiological realism may foster mechanistic investigations of their roles in spatial navigation.Grid cells are place-modulated neurons with firing fields at multiple discrete and periodically spaced locations [1].These fields form a hexagonal pattern that tiles the entire space available to the animal.Many computational simulations have helped advance theories on how grid cells and other cellular activities are involved in spatial representation [2][3][4][5].Continuous attractor networks (CANs) are a major class of computational models used to recreate grid cell activities [6].
This study investigates whether the CAN model of grid cells is compatible with measured biophysical details of the underlying neural network, such as neuron type-specific excitability, connectivity, and synaptic signaling.Finding an answer is important because it could add helpful realism to spatial navigation modeling.Prior work successfully modeled selected physiological characteristics (e.g., [7]), but obstacles remain in systematically finding relevant experimental measurements, designing a computationally efficient system that is compatible with those details yet can run at a useful neural scale, and acquiring animal data to compare modeled results.This work tackles those obstacles by taking advantage of existing open-access resources to test a unique combination of biological details in a data-driven spiking neural network (SNN) model of pertinent portions of the entorhinal circuit.The results contribute supporting evidence in favor of the plausibility of the CAN model.
One CAN grid cell model created accurate path integration based on inputs that encode only a rat's velocity using a firing rate coding [2].Attractor states in the model occur as bumps of firing that are greater than areas outside of those regions.The study stipulates that grid cells have preferred directions, i.e., they are each specialized to respond more to certain direction-modulated input.An offset in the outgoing synaptic connections of grid cells influences the directional movement of the bumps when activating those grid cells.A key property that helps generate bump activities is synaptic connectivity (and related weights) in center-surround (CS) structures.The CS structure supports bump firing in center regions and inhibits firing in surrounding regions.Another CAN grid cell model, using exponential integrate and fire neurons, included additional components such as place cells, signaling between grid cells and interneurons (INs), and three different synaptic receptor types [4].
The present work uses neuron types and their properties from Hippocampome.org[8] to test whether previously proposed models can reproduce key physiological details observed in animal recordings when incorporating data-driven biophysical details.Hippocampome.org is a knowledge base of rodent hippocampal formation circuitry and offers a wide variety of neural properties that are ready to be integrated into simulations [9].Parameters available from Hippocampome.orgalong with the corresponding experimental measurements include the Izhikevich Model (IM) of neuron firing and excitability [10][11][12], Tsodyks-Markram (TM) model of synaptic signaling and short-term plasticity [13][14][15][16], neuron counts [17], and connection probabilities [18].Hippocampome.orgalso provides a wealth of neuron type-specific details about hippocampal activity and molecular expression that can add to the biological realism of computational simulations [19][20][21][22].We recently created a large-scale model of resting-state dynamics in CA3 using Hippocampome.orgproperties [23] and showed its suitability to simulate pattern completion via cell assemblies [24].The present work complements those studies by also demonstrating the application of Hippocampome.orgresources to modeling complex network interactions efficiently with reusable code.This study's software is released as an open source resource to explore how varying property parameters within reported biological ranges affect the reproduction of electrophysiological data from mice.Comparing experimental recordings to simulated recordings assesses the ability of a physiologically realistic CAN model to recreate real neural activities.Testing of this model's neural parameter predictions in rodents can verify their accuracy and be used to further refine the model.

Results
Our research design consisted of utilizing experimentally recorded mouse locomotion to activate conjunctive (EC LI-II multipolar pyramidal) cells based on speed and direction, and place (CA1 pyramidal) cells based on position, as described further in Supplementary Materials V. We then compared the corresponding simulated and real grid cell plots to gauge whether a spiking neural network implementation of the continuous attractor model with parameters based on experimental measurements could produce firing dynamics consistent with those observed in animal studies.
Simulations recreated well several properties of real mouse recordings, including grid field spacings, sizes, and orientation, as well as firing rates (Figure 1).Multiple studies have reported grid field size and spacing increase from the dorsal to the ventral end of the MEC [25,26].This consistency of scale change along the dorsoventral axis suggests a purpose for this pattern.Our model could create grid fields with different scales qualitatively similar to those observed experimentally, such as large, intermediate, and small (Figure 1A-C, respectively).
Int. J. Mol.Sci.2023, 24, x FOR PEER REVIEW 3 of 20 cell plots to gauge whether a spiking neural network implementation of the continuous attractor model with parameters based on experimental measurements could produce firing dynamics consistent with those observed in animal studies.Simulations recreated well several properties of real mouse recordings, including grid field spacings, sizes, and orientation, as well as firing rates (Figure 1).Multiple studies have reported grid field size and spacing increase from the dorsal to the ventral end of the MEC [25,26].This consistency of scale change along the dorsoventral axis suggests a purpose for this pattern.Our model could create grid fields with different scales qualitatively similar to those observed experimentally, such as large, intermediate, and small (Figure 1A-C, respectively).The simulations could generate specific rotations of grid patterns (see Supplementary Materials V for methods) qualitatively similar to those observed in real recordings (Figure 1A).When considering the whole analyzed time period, identical in real experiments and simulations (142.5, 141.4, and 60.2 min in Figure 1A-C, respectively), the mean firing rate of the grid cells in Figure 1 matched to a one decimal place accuracy.We thus compared the mean firing rates across an equal number of simulated and real grid cells (N = 29).For this purpose, we selected six simulated experiments to approximate the range of firing observed in real cells, including large-grid scales and low-and high-firing rates in smallgrid scales.Scale sizes and firing rates were grouped based on the mean levels of the real cells (<=50% for small or low and >50% for large or high), and scale measurement methods are in Supplementary Materials VI.In particular, the simulated small-scale cells included eleven low-and eleven high-mean-firing rate cells, and large-scale cells included seven cells, which were similar to that observed in real cells.Because we could not assume the normality of the mean firing rates, we opted for a non-parametric Wilcoxon rank sum test.The firing rates of the real cells had a median of 1.76 spikes/s and interquartile range IQR of 1.55 spikes/s, while the corresponding values for the simulated cells were 1.75 spikes/s (median) and 2.02 spikes/s IQR (Figure 2).The p-value of 0.93 indicates a lack of support The simulations could generate specific rotations of grid patterns (see Supplementary Materials V for methods) qualitatively similar to those observed in real recordings (Figure 1A).When considering the whole analyzed time period, identical in real experiments and simulations (142.5, 141.4, and 60.2 min in Figure 1A-C, respectively), the mean firing rate of the grid cells in Figure 1 matched to a one decimal place accuracy.We thus compared the mean firing rates across an equal number of simulated and real grid cells (N = 29).For this purpose, we selected six simulated experiments to approximate the range of firing observed in real cells, including large-grid scales and low-and high-firing rates in small-grid scales.Scale sizes and firing rates were grouped based on the mean levels of the real cells (<=50% for small or low and >50% for large or high), and scale measurement methods are in Supplementary Materials VI.In particular, the simulated small-scale cells included eleven low-and eleven high-mean-firing rate cells, and large-scale cells included seven cells, which were similar to that observed in real cells.Because we could not assume the normality of the mean firing rates, we opted for a non-parametric Wilcoxon rank sum test.The firing rates of the real cells had a median of 1.76 spikes/s and interquartile range IQR of 1.55 spikes/s, while the corresponding values for the simulated cells were 1.75 spikes/s (median) and 2.02 spikes/s IQR (Figure 2).The p-value of 0.93 indicates a lack of support for the hypothesis that the firing rates of simulated and real cells have different distributions.for the hypothesis that the firing rates of simulated and real cells have different distributions.We also measured and compared grid field sizes and spacing between real and simulated recordings (see Supplementary Materials VI for methods).These comparisons utilize the same cells used in the firing rates analyses and the Wilcoxon test for the same reasons described for firing rates.The grid field sizes of real cells had a median of 11.16 cm and IQR of 2.79 cm, and simulated cells had a median of 11.23 cm and IQR of 4.74 cm (Figure 2).Grid field spacing in real cells had a median of 44.87 cm and IQR of 16.03 cm, and simulated cells had a median of 43.56 cm and IQR of 15.53 cm.The p-value for sizes was 0.70 and for spacing 0.71, both indicating a lack of support for the hypothesis that field sizes and spacings have different distributions between simulated and real cells.
Next, we investigated the robustness of the SNN implementation of the CAN model, as quantified by the grid score, relative to variations of neuronal and synaptic model parameters within their biologically realistic ranges.For IM models, focusing on the stellate cells, we considered parameter ranges to be biologically realistic if the resultant firing pattern had a fitness error up to 200% of the minimum, corresponding to the Hippocampome.orgbest ranking parameter value (see Section 4 for details).For TM models, we defined biological realism as two standard deviations from the mean of each parameter listed on Hippocampome.org.When testing the TM parameters, we focused on the synaptic connectivity from stellate cells to INs, varying the values of all three sets of glutamatergic connections (onto axo-axonic, basket, and basket multipolar cells) at the same time.
This analysis demonstrated that the emergence of grid score dynamics in the CAN model is more sensitive to certain model parameters than others (Figures 3 and S3).In particular, several IM and TM parameters have wide ranges of variability that generate above-threshold grid scores.These include stellate cell IM parameters k, a (Figure 3A), Vt, Vpeak (Figure S3B), and Vr (Figure S3C), as well as stellate-to-IN TM parameters g, taud (Figure 3B), and tauu (Figure 3D).This indicates that grid field functionality could be preserved amid a broad diversity of physiological modulatory factors that could affect the intrinsic excitability and synaptic signaling of MEC LII stellate cells.In contrast, the range of viable values for TM parameter taux producing sufficient grid scores was substantially more limited (Figure 3D).This indicates that taux may be specialized toward accommodating grid field activities.The situation is more complex for IM parameters b and d.While each of them has fairly broad ranges of values yielding strong grid patterns, they also depend on each other via negative correlation: the highest values of parameter b only allow high grid scores when paired with the lowest values of parameter d, and vice versa (Figure 3C).A similar relation is present between IM parameters C and Vmin (Figure S3A).We also measured and compared grid field sizes and spacing between real and simulated recordings (see Supplementary Materials VI for methods).These comparisons utilize the same cells used in the firing rates analyses and the Wilcoxon test for the same reasons described for firing rates.The grid field sizes of real cells had a median of 11.16 cm and IQR of 2.79 cm, and simulated cells had a median of 11.23 cm and IQR of 4.74 cm (Figure 2).Grid field spacing in real cells had a median of 44.87 cm and IQR of 16.03 cm, and simulated cells had a median of 43.56 cm and IQR of 15.53 cm.The p-value for sizes was 0.70 and for spacing 0.71, both indicating a lack of support for the hypothesis that field sizes and spacings have different distributions between simulated and real cells.
Next, we investigated the robustness of the SNN implementation of the CAN model, as quantified by the grid score, relative to variations of neuronal and synaptic model parameters within their biologically realistic ranges.For IM models, focusing on the stellate cells, we considered parameter ranges to be biologically realistic if the resultant firing pattern had a fitness error up to 200% of the minimum, corresponding to the Hippocampome.orgbest ranking parameter value (see Section 4 for details).For TM models, we defined biological realism as two standard deviations from the mean of each parameter listed on Hippocampome.org.When testing the TM parameters, we focused on the synaptic connectivity from stellate cells to INs, varying the values of all three sets of glutamatergic connections (onto axo-axonic, basket, and basket multipolar cells) at the same time.
This analysis demonstrated that the emergence of grid score dynamics in the CAN model is more sensitive to certain model parameters than others (Figures 3 and S3).In particular, several IM and TM parameters have wide ranges of variability that generate above-threshold grid scores.These include stellate cell IM parameters k, a (Figure 3A), V t , V peak (Figure S3B), and V r (Figure S3C), as well as stellate-to-IN TM parameters g, tau d (Figure 3B), and tau u (Figure 3D).This indicates that grid field functionality could be preserved amid a broad diversity of physiological modulatory factors that could affect the intrinsic excitability and synaptic signaling of MEC LII stellate cells.In contrast, the range of viable values for TM parameter tau x producing sufficient grid scores was substantially more limited (Figure 3D).This indicates that tau x may be specialized toward accommodating grid field activities.The situation is more complex for IM parameters b and d.While each of them has fairly broad ranges of values yielding strong grid patterns, they also depend on each other via negative correlation: the highest values of parameter b only allow high grid scores when paired with the lowest values of parameter d, and vice versa (Figure 3C).
A similar relation is present between IM parameters C and V min (Figure S3A).The areas within thresholds in these plots help reveal the sensitivity of each score to the parameters.The areas within thresholds in these plots help reveal the sensitivity of each score to the parameters.To illustrate the effect of these biophysical parameters on grid cell activity, it is useful to examine examples of the firing rate maps.In particular, values from areas of the parameter space giving rise to a high grid score yield strong grid patterns (Figure 3E), while values outside of those areas cause the firing to lose its grid pattern (Figure 3F).In this example, the average firing rate is greater in the low grid score case compared to the high grid score one (~2.6 spikes/s in panel E and ~15.0 spikes/s in panel F).The firing activity is largely focused on field areas in the high-grid score case (Figure 3E), while it is widespread in the environment in the low-grid score case (Figure 3F).Of note, the simulation settings used to produce the grid scores and other results in Figure 3, aside from the parameters explored, remained constant and were the same as those used in Figure 1B, corresponding to an intermediate grid scale cell.The results remained reasonably consistent when extending the parameter search analysis to large-and (to a narrower extent) small-field cells (Figure S4).
One of the advantages of an SNN implementation of the CAN grid cell model is the ability to inspect firing dynamics from single cells to the entire circuit.Throughout the simulations, network activity yielded realistic, sustained, non-repetitive spiking patterns (Figure 4 and Supplementary Video at Hippocampome.org/RTPlots2).Specific activity patterns emerged in each collection of neuron types.At any given time, only a subset of grid cells is active, and the constellation of this subset changes smoothly over time generating the appearance of curved lines in the raster plot, consistent with an activity bump moving along a continuous attractor space.Stellate cell firing stimulated INs, increasing their spiking immediately after the activation of grid cells.In turn, INs inhibited the grid cell layer in different locations due to the CS connectivity.Conjunctive cell spiking reflects a combination of a baseline firing level and activity targeting grid cells with specific preferred directions.This results in nearly all conjunctive cells firing together at times with only specific groups of cells spiking in between bouts of synchronous activity.
CA1 pyramidal cells showed relatively sparser firing compared to other neuron types.This is consistent with the known physiology of these neurons and with the function of place cells, which only respond to restricted areas of an environment.The single bump of CA1 PC place fields limits their firing to a relatively narrow group of neurons during the course of a few seconds.In contrast, grid cells, and the INs connected to them, fired in multiple places due to the multi-bump structure of the CAN.The intended movement velocity of the animal also affects the level of firing at each time.A speed control function differentially activates conjunctive cells depending on the speed of the animal, and that signal propagates throughout connected cells (Supplementary Materials V).Conjunctive cells also more strongly stimulate grid cells with preferred directions that match the direction of animal movement at each timestep, which contributes to signaling changes.
Next, we investigated the rhythms of the network activity displayed in Figure 4 by power spectrum density analysis using the combined spiking activity of each neuron type group.When considered all together, MEC LII neurons had the highest power in the beta band (12-25 Hz), followed by the gamma band (25-100 Hz), with 41.7% and 25.5% of total power, respectively.Taken individually, stellate, axo-axonic, basket, and basket multipolar cells had their highest powers in the beta band (55.9-44.6%),and second highest in the gamma band (33.4-24.8%; Figure 4C).In the hippocampal formation, phase coupling is especially informative relative to the slow network oscillations [21,23,27].We selected the peak firing of stellate cells within the beta rhythm as a common frequency to evaluate mean resultant vector length (MRVL) values of all neuron types during the first minute of the simulation.We based this choice on the stellate cells being directly connected to all other neuron types and the beta rhythm having the highest power across neuron types (Figure 4C).Stellate, axo-axonic, basket, and basket multipolar cells had an MRVL value of 0.04, and multipolar pyramidal cells had an MRVL value of 0.05.These neuron types had Rayleigh test p-values smaller than 0.001.These results indicate these neuron types had significant but weak spike-phase coupling to the beta rhythmicity in grid cell firing [28].CA1 pyramidal cells showed relatively sparser firing compared to other neuron types.This is consistent with the known physiology of these neurons and with the function of place cells, which only respond to restricted areas of an environment.The single bump of CA1 PC place fields limits their firing to a relatively narrow group of neurons during the course of a few seconds.In contrast, grid cells, and the INs connected to them, fired in multiple places due to the multi-bump structure of the CAN.The intended movement velocity of the animal also affects the level of firing at each time.A speed control function differentially activates conjunctive cells depending on the speed of the animal, and that signal propagates throughout connected cells (Supplementary Materials V).Conjunctive cells also more strongly stimulate grid cells with preferred directions that match the direction of animal movement at each timestep, which contributes to signaling changes.
Next, we investigated the rhythms of the network activity displayed in Figure 4 by The activity of the grid cells in their corresponding neural space reflected the multibump structure of the CAN model (Figure 5A and Supplementary Video).At the same time, the spiking neural network implementation of the CAN model also allows the virtual recording of voltage and current in a selection of its simulated neurons (Figure 5B-D and Supplementary Video).The current plots recorded in the stellate cells combine the contributions from both fast (AMPA and GABA A ) and slow (NMDA and GABA B ) receptors from all sources, namely CA1 pyramidal, EC LI-II multipolar pyramidal, EC LII axo-axonic, time, the spiking neural network implementation of the CAN model also allows the virtual recording of voltage and current in a selection of its simulated neurons (Figure 5B-D and Supplementary Video).The current plots recorded in the stellate cells combine the contributions from both fast (AMPA and GABAA) and slow (NMDA and GABAB) receptors from all sources, namely CA1 pyramidal, EC LI-II multipolar pyramidal, EC LII axoaxonic, MEC LII basket, and ECLII basket multipolar cells.The currents recorded in the axo-axonic cells correspond to the excitatory input (AMPA and NMDA) from stellate cells.In order to facilitate future related research, we released the simulation software pertaining to this work as open source at Github.com/Hippocampome-Org/spatial_nav.We also release as open access all mouse recording data (Osf.io/ud27p)from [29] and five new cells recorded for this work at Osf.io/jd6rv.Among several functions, the released software automates parameter exploration and generates recording scores and plots.Researchers can use this free resource to reproduce the simulation results, further investigate the spiking network model, expand on this work, and repurpose parts or the entirety of the code for studies of different subjects.The software includes usage instructions and documentation on its components, covering a wide range of topics to guide users from beginner to advanced levels of work with these simulations.Additional details on hardware and software versions used in this study are also available in the documentation.

Discussion
Using a biophysically detailed simulation, we recreated well several properties of grid cells observed in animals, including field spacing and size as well as rotation of grid patterns.Furthermore, the simulations recreated the increase in the spacing and size of individual grid fields similar to those observed in rodents along the dorso-ventral axis of the MEC [30,31].Grid cells can be grouped into discretized modules that share similar field sizes and spacings, grid pattern orientation, and other characteristics [1,32,33].These modules have been theorized to represent different spatial scales of an environment in a computationally efficient manner [34][35][36] and could work together in spatial cognition to inform an animal of its location [37,38].Controlling field spacing and size in the simulation helped recreate properties of different grid cell modules throughout a realistically broad range of scales.This opens the prospect of utilizing this spiking neural network implementation of the CAN in future explorations of the computational function of multi-scale grid representation.
Our simulation also reproduced with good accuracy several characteristics of the grid cell firing rates observed in real animals, including the average, the variability, and the differences between in-field and between-field levels.Prior CAN model implementations typically created homogeneous firing fields, contrasting the natural heterogeneity of field firing and shape in experimental recordings.The firing variability of our simulation was an emergent property of the data-driven neural and synaptic properties modeled based on empirical evidence sourced from Hippocampome.org.Given the potential importance of variable neuronal firing in network function [39], it may be interesting to investigate computationally the role of heterogeneous grid cell activity in spatial navigation and pathfinding.
We deemed it valuable to recreate grid pattern rotations relative to the environment borders because studies indicate that those rotations could be non-random and functionally relevant.In particular, the offset angle of the grid pattern axes from the walls minimized symmetry with the area borders [40].Moreover, moving a cue card along the wall of a circular environment coincided with grid cell pattern rotations [30,41].Interestingly, rotation of geometrically polarized enclosures (e.g., square walled areas) affected grid pattern rotation more than distal cues (e.g., rectangular cue card), but in geometrically symmetrical enclosures (e.g., circular areas) the rotation of distal cues had greater effects on grid pattern orientations [42].The data-driven SNN implementation of the CAN released with this work could be useful in the future to investigate the biophysical and circuit mechanisms underlying the geometric relationship between the spatial representation of environmental stimuli and grid patterns.
Parameter exploration identified ranges of biophysical values supporting grid cell functionality, and those ranges constitute testable hypotheses bridging quantitative physiology and circuit function.For instance, future studies could attempt to fit the firing patterns of putative MEC LII stellate cells recorded in vivo with IM models.This would allow a statistical comparison of the distributions of IM parameters corresponding to grid and nongrid cells against the predictions of our simulations.In principle, similar experiments could be designed to fit synaptic activity between grid cells and INs with TM model parameters, although in vivo recording of unitary synaptic signals remains technically challenging in live behaving animals.
The synaptic connectivity in the simulation was constrained by values reported for each directional pair of neuron types on Hippocampome.organd the peer-reviewed literature.For each presynaptic and postsynaptic neuron type, Hippocampome.orgestimates the overall probability of connection for the whole MEC based on axonal and dendritic morphology.In contrast, most experimental studies quantifying connectivity, typically based on paired neural recordings, rely on a limited sample of neurons in relative local proximity to each other.We designed our simulation with intermediate connectivity between these two extremes, consistent with the modeled neuron counts in each type falling between full and local scale.For instance, because experimental studies reported 25.9% connectivity between stellate cells (potential grid cells) and INs in mice, we constrained the connectivity in our model so as not to exceed that amount.The effective simulation of cells with biologically realistic grid patterns and high grid scores using this connectivity indicates potential matches to real cell properties, but more animal recordings are needed to fully characterize circuit connectivity.
In rats, more connections are present between LII basket cells and principal cells in the dorsal region as compared to the ventral region of the MEC [43], including stellate cells that could be grid cells.The CAN model in this work followed a similar pattern of connectivity with grid cell to IN connections relative to grid scales.When simulating large-grid scales, found in ventral MEC, only a limited number of inhibitory CS rings fit within the neural layer of grid cells before notably overlapping with each other, which caused unwanted bump movements.In our SNN simulation, this CS inhibition is generated by fast-spiking INs, which are activated by excitatory synaptic connections from grid cells.Smaller grid scales, found in dorsal MEC, allowed many more rings without overlap, and therefore a greater number of grid cell to IN connections with effective bump movements.The approximate average number of five, twelve, and sixteen IN connections from each grid cell, in large-, intermediate-, and small-grid scale simulations, respectively, follows the same connectivity trend reported along the dorsoventral axis.A more in-depth analysis of circuit connectivity in our simulation is available in Supplementary Materials IV.
Future rodent recording can test if similar features of the firing patterns observed in the simulated raster plots also emerge in vivo.Since the simulation only includes a subset of MEC neuron types, the comparison would require techniques to distinguish those kinds of neurons, such as stellate cells, while recording thousands at a time with a 1-millisecond temporal resolution.When available, we predict that such experiments will reveal grid cell activity occurring in curved line patterns as observed in our simulated raster plots, given that the cell firing reflects animal movement.Such firing is predicted to occur along a toroidal manifold [44].To display that pattern, it will be sufficient to arrange the grid cells in the raster plot in sequential order of their response to the animal positions in the environment.Moreover, the periodicity of grid cells emerging from a sequence code of trajectories [36,45] in an environment could produce various patterns of how the activity bump moves in time, where the pattern is informative of the trajectory of the animal, as observed in our data.The firing of grid cells in the simulation within activity bumps often follows closely the movement of the mouse (Figure 4A).
One of the brain structures targeted earliest by Alzheimer's disease is the hippocampal formation, including the entorhinal cortex [46,47], as characterized by altered neural excitability [48] and synaptic dysfunction [49].Spiking neural network simulations exploring how such perturbations could affect grid fields, in parallel with real recordings of animals engaged in a spatial task, could help bridge the gap between data from animal experiments and human relevance, and lead to early-detection or treatment advancements.In addition to investigating disease, biologically realistic computational models such as those developed in this work could more generally foster a deeper understanding of the role of attractor dynamics in multiple forms of learning and memory [50].
In our simulations, similar to prior CAN models, a single layer of excitatory neurons provides direction and speed input signals to grid cells.However, speed signals could also come from inhibitory INs [51], and a near real-time accurate speed code by firing rate may require fast-spiking INs [52,53].Another alternative or complementary source of speed signal could come from the correlation of speed with theta rhythms [29,53] and related cholinergic [23] or glutamatergic [54,55] input from the medial septum.Expanding this SNN implementation of the CAN to include those additional mechanisms might help shed light on their computational plausibility.
In the simulation, grid cells rely on place cell input to assist with drift correction.Boundary cells, visual landmarks, and other mechanisms may also contribute to drift correction [56].The place cell input used in this simulation causes extra firing to appear in one grid field relative to others.While this is often observed in real animal data [29], the current model might have challenges in exactly reproducing the activity of grid cells with uniform firing fields.Well-designed automated speed control functions in this model can help reduce the need for drift correction, but incorporating additional drift correction methods might aid in capturing more types of grid field firing.Tightly coordinated computational and experimental research could reveal major insights into the mechanisms of drift correction.
In rare instances (<2% of recorded time), the animals traveled especially fast, e.g., 90 cm/s.We considered those speeds as outliers and did not tune the model to match those cases.Attempting to optimize the simulations for outlier speeds with the current design decreased the accuracy of the results at slower speeds.This was the case when an animal rapidly transitioned from extremely fast movements to very slow speeds, causing a spurious neural signal to linger during slower speeds.It would be interesting for future research to explore the circuit dynamics underlying these transitions in parallel animal experiments and computational simulations.

Methods
A recent review of publications on simulations of the hippocampal formation identified CAN grid cell models as a notable case relevant to spatial navigation [57].Here, we used Hippocampome.orgparameters to build a new SNN of grid cells adopting multiple methods from previous CAN models [2,4].In particular, we preserved the original model realism to only include inhibitory signaling between grid cells (an optional choice in [2]) as opposed to other models that rely on excitatory signaling between grid cells [58].Also, as in the original model [2], preferred directions in grid cells managed bump movement, and each cell had a north, south, east, or west preferred direction.We adapted the model to convert its firing rate representation into the Izhikevich dynamical system formalism, which captures individual neuron spikes in a biologically realistic manner while also accommodating the integration of Hippocampome.orgparameters.We also used place cell input to help correct for drift in bump attractors as in a published CAN model based on exponential integrate-and-fire neurons [4].We selected the TM synapse model using neuron type-specific parameters from Hippocampome.orgfor use in our work.In addition to factors similar to those included in previous models (e.g., [4]), such as receptor-specific signal decay constants and reversal potentials, the TM model also includes variables for facilitation and depression, accounting for short-term plasticity.

Mapping Functional Neuron Types to Morphological Neuron Types
Grid cells are most abundant in layer 2 (LII) of the medial entorhinal cortex (MEC), but also exist in other brain subregions such as the pre-and parasubiculum [59,60].Grid cells in MEC LII can be either glutamatergic stellate cells or pyramidal cells [61].In rats, the primary population of excitatory cells in MEC LII are stellate cells [62,63].Those stellate cells express the large glycoprotein reelin [64,65] and largely influence each other through indirect connections via inhibitory interneurons.In contrast, pyramidal cells in LII may have direct connections with each other [18].Given the different interaction mechanisms and the empirical evidence on relative abundance, we chose to map grid cells to MEC LII stellate cells.
Anatomical, electrophysiological, and molecular information regarding GABAergic INs in the entorhinal cortex is notoriously scarce [9,66].Nevertheless, the identification of grid cells with MEC LII stellate cells provides definitive constraints as to the nature of the INs that could provide CS inhibition.In particular, because MEC LII stellate cells have soma and dendrites limited to layers 1 and 2, the INs that grid cells can use to communicate with each other must have axons extending into those layers.Moreover, parvalbumin-expressing (typically perisomatic-targeting) INs affect grid cell behaviors more than other IN neuron types [67].We thus chose to include in our model the three Hippocampome.orgneuron types with those characteristics, namely LII axo-axonic, LII basket, and LII basket-multipolar cells, to capture a diversity of neural properties that may be present in animals (Figure 6).cate with each other must have axons extending into those layers.Moreover, parvalbumin-expressing (typically perisomatic-targeting) INs affect grid cell behaviors more than other IN neuron types [67].We thus chose to include in our model the three Hippocampome.orgneuron types with those characteristics, namely LII axo-axonic, LII basket, and LII basket-multipolar cells, to capture a diversity of neural properties that may be present in animals (Figure 6).The neuron type chosen to represent place cells is CA1 pyramidal cells because those are well-characterized neurons with place cell characteristics [68].For simplicity, we simulated direct connections from place cells to grid cells.Indirect connectivity through deeplayer MEC cells may be more realistic, though recent evidence also supports potential direct connections from CA1 to MEC LII [69].Our simplified abstract connection represents both the indirect pathway and the possible direct pathway, and is similar to the design adopted by prior models [4].That similarity in network architecture helps test the effect of adopting data-driven parameters from Hippocampome.orgwhile otherwise minimizing changes from previous methods.Our model does not include connections that may occur from grid cells to place cells as modeled in [70].The neuron type chosen to represent place cells is CA1 pyramidal cells because those are well-characterized neurons with place cell characteristics [68].For simplicity, we simulated direct connections from place cells to grid cells.Indirect connectivity through deep-layer MEC cells may be more realistic, though recent evidence also supports potential direct connections from CA1 to MEC LII [69].Our simplified abstract connection represents both the indirect pathway and the possible direct pathway, and is similar to the design adopted by prior models [4].That similarity in network architecture helps test the effect of adopting data-driven parameters from Hippocampome.orgwhile otherwise minimizing changes from previous methods.Our model does not include connections that may occur from grid cells to place cells as modeled in [70].
Animal studies helped inform how to model the sources of signals that provide movement velocity information to grid cells.Coding for speed and direction contributing to velocity signals may come from separate neuron types or from conjunctive cells.One potential source of speed signal is firing rates in a subpopulation of MEC neurons [71].Head-direction cells (HDCs) could provide running direction signals to grid cells, although HDCs can code more accurately for head direction than running direction [72].HDCs have been identified in 11 brain regions including MEC LII [26,73,74].Based on prior models and rodent studies, we included a cell group that provides velocity-modulated excitation to grid cells [2,4].We mapped these direction-and speed-modulated conjunctive cells to LI-II Multipolar pyramidal cells, the most abundant glutamatergic neurons in EC based on Hippocampome.org[17].These neurons have potential synaptic connections to stellate cells [18,75] and can thus supply direct excitatory input to grid cells (see Section 3 for an alternative based on GABAergic signaling).

Computational Modeling: Neurons and Synapses
We utilized CARLsim [76] as the simulation platform because of its efficiency at high-performance computing of spiking neural networks using graphical processing units (GPU).The adoption of IM to describe neuronal excitability and firing allows simulations to recreate biologically realistic diversity and complexity of neuron type-specific input-output dynamics.This modeling formalism utilizes 4 parameters (K, a, b, and d) to control the intrinsic dynamics of the neuron, e.g., affecting bursting, rheobase, and spike frequency adaptation.Five additional IM parameters scale the neuronal response by setting the cell capacitance (C), resting voltage (Vr), threshold (Vt), spike cutoff value (Vpeak), and postspike reset value (c).The 9 IM parameters (Table 1) were taken from Hippocampome.orgfor all glutamatergic neuron types.The available empirical evidence for MEC inhibitory cells is insufficient to constrain data-driven IM.Therefore, we used standard fast-spiking IM parameters suitable for these types of perisomatic, parvalbumin-expressing neurons [77].The numbers of neurons in each type included in the simulation (Table 1) are also consistent with the Hippocampome.orgcell census ranges for the mouse.
Table 1.IM parameters.a These values were utilized in the intermediate and small grid-scale simulations (e.g., Figure 1B,C).The number of INs was increased to 1200 in each of the three types for large grid-scale simulations (Figure 1A) to help accommodate control of the larger bump sizes.We also strived to maximize the simulation biological realism in setting the properties of synaptic connections for each pair of neuron types, including the numbers of connections and the TM model parameters (Table 2).We considered both Hippocampome.orgprobabilities [18] and additional animal data-driven values to select the connectivity ratios between neuron types.In particular, a study in mouse MEC LII reported 35.7% connectivity from INs to stellate cells and 25.9% from stellate cells to INs [78], while another mouse study reported a 26.9% connection probability in both directions [79].Additional studies have found higher connectivity between stellate cells and INs in rats [43,80].The synaptic connections between INs and grid cells in the simulation follow a CS distribution of connection weights, with CS circumferences configured with the aim to match grid scales found in animal data, as detailed in Supplementary Materials IV (the provided software helps automate construction of these distributions).All other connections had unitary weights and the model did not include long-term plasticity.

Neuron
The TM model calculates not only the intensity and duration of synaptic signals via conductance (g) and kinetic decay (tau d ) parameters, but also short-term plasticity with parameters (U, tau u , and tau x ) controlling facilitation and depression [15].Hippocampome.orgprovides data-driven ranges of TM parameters for fast (AMPA and GABA A ) currents (Hippocampome.org/synapse).For each pair of neuron types, all TM parameters in our simulations were within two standard deviations of the Hippocampome.orgmeans for 56-day-old male mice at 32 degrees Celsius in voltage clamp.This range represents 95% of the observed values under the assumption of a Gaussian distribution [81].We set the conductance level of slow currents (NMDA and GABA B ) based on available evidence as described in Supplementary Materials VII, and used the CarlSim default parameters for the NMDA and GABA B time decays (tau d ).We created software programs to automatically test how simulation parameters within (and beyond) the ranges of Hippocampome.orgIM neuron and TM synapse model values affected grid cell activity.Specifically, we adopted the fitness error from prior work [11] to quantify the deviation of firing patterns obtained with a given selection of IM parameters from the Hippocampome.orgexperimental data for the same neuron type [10].Parameter exploration started from Hippocampome.orgIM parameters and used an evolutionary algorithm (EA) software (as in [11]) to test higher and lower values.We selected a firing pattern fitness error threshold of 200% of the minimum, which corresponded to the Hippocampome.orgparameter value.This choice balanced accuracy and variability in reproducing the experimental firing patterns.
To evaluate how well the firing activity of a stellate cell reflected a grid pattern, we utilized a numerical score (grid score) previously introduced to quantify experimental data [29].We chose a grid score threshold of 0.2 as an indicator of an acceptable grid pattern after a review of many plots with different scores.In these parameter exploration analyses, for computational practicality, we used a shorter simulation time than the full animal experiment time (see Supplementary Materials I for additional details).
We analyzed spiking from the simulation to detect if rhythm-related properties were present.We created multi-taper power spectral density plots (to suppress noise) with Gaussian smoothing using a 3-ms moving window size (to minimize harmonic artifacts).We used a Rayleigh test to determine statistical significance in spike-phase coupling [82].We also computed MRVL values with the CircStat MATLAB toolbox [83] to quantify the strength of spike-to-phase coupling in the 0-1 range [84].

Subjects and Neural Recording
We compared our SNN continuous attractor model of grid cells with benchmark experimental data from live-behaving mice.Of 33 experimentally recorded grid cells utilized here, 28 came from a re-analysis of our earlier published work [29], and 5 came from new recordings collected for this study.The earlier experiments investigated differences in grid cell firing between sessions recorded in light and darkness [29].We only utilized the data from normal lighting conditions, i.e., the time periods that the lights were turned on.
The additional data were collected from 3 male mice (Mus musculus), 4 months old (1 wild type, C57BL/6J, and 2 heterozygous transgenic ChAT-Cre tg/wt , B6;129S6-Chat tm1(cre)Lowl /MwarJ; The Jackson Laboratory).All experimental procedures were performed in accordance with the regulations stated in the Guide for Care and Use of Laboratory Animals by the National Research Council and approved by the George Mason University Institutional Animal Care and Use Committee (Protocol No.: 0501, approved on 14 January 2022).Mice were implanted with a microdrive containing 4 movable tetrodes targeting the superficial layers of the left MEC at the following stereotaxic coordinates: 3.4 mm lateral to the midline and 0.5 mm anterior to the transverse sinus angled at a polar angle of six degrees.For an additional experimental purpose not reported in this study, a virus injection of calcium indicator GCaMP8s and a chronic implant of an optical fiber into the medial septum was performed.After surgery, mice were individually housed on a reversed 12 h light/dark cycle.The surgery for microdrive implantation and the electrophysiological recordings were performed as described previously [29].Mice were given one week to recover from the surgery.Recordings of spikes and local field potentials were carried out while the animals randomly foraged for scattered pieces of cereal bits (Froot Loops, Kellogg Company, Battle Creek, MI, USA) in a wood open field environment of 45 × 45 cm 2 with 30 cm high walls until animals sampled a sufficiently large area of the open field.One visual cue card (a white triangle) was placed on one of the walls of the open field.Tetrodes were turned in increments of 50 µm until spiking activity was detected and theta oscillations were observed indicating the tetrodes had entered the superficial layers of the MEC.Tracking of the position and head direction of the mouse was achieved by using head-mounted infrared LEDs.After data collection, animals were anesthetized deeply with isoflurane and transcranial perfusion with DPBS 1X + 10% buffered formalin was performed.To visualize tetrode track and MEC recording sites, sagittal brain sections (40-50 µm) were sliced using a vibratome and stained with a Cresyl violet solution.
The data analysis, including the classification of putative principal neurons, computation of spatial firing rate maps, and grid score, was performed as described previously [29].

Hardware and Software
Results from this work were produced by running CARLsim6 on Linux with Ubuntu 21.10 and CentOS 8.5.We utilized a stand-alone machine with a Nvidia GeForce RTX 3090 GPU and an Intel i9-13900KS CPU.We also used a supercomputer at George Mason University Office of Research Computing with an Nvidia A100 GPU and AMD Milan Core 4.0 GHz CPU.Typically, the simulation used 5 GB of RAM and 10 GB of GPU RAM.The experiment in Figure 1B was used to assess simulation speed performance.One second of simulated time took 1.4 s of real time on the stand-alone machine and 2.4 s of real time on the supercomputer.Additional software and hardware details are listed in the GitHub documentation.

Conclusions
This work integrated detailed neural and synaptic properties from Hippocampome.org and the experimental literature to enhance the realism of a grid cell CAN model.The results showed many qualitative and quantitative similarities to real recordings.Exploration of the effects of multiple physiological parameters on grid patterns narrowed down the biologically realistic ranges that may enable proper spatial navigation, generating testable hypotheses for future animal studies.Reproduction of key animal recording features with good accuracy, while integrating neuron type-specific empirical data on cell excitability and synapse signaling, lends support for the CAN model of grid cell activity.

Figure 1 .
Figure 1.Comparing simulated and real results.Panels (A-C) show results with different grid scales.Values reported at the top of the panels are grid score (G), mean firing rate in spikes/s (m), and peak firing rate, also in spikes/s (p).Top row: trajectory of animal in the environment (black lines) and spike location (red markers).Middle row: firing rate map of the neuron (color scale bars in spikes/s).Bottom row: autocorrelogram of the rate map (color bar ranges: (A), −74% [percent is relative to maximum correlation] to 100%; (B), −51% to 100%; (C), −60% to 100%).

Figure 1 .
Figure 1.Comparing simulated and real results.Panels (A-C) show results with different grid scales.Values reported at the top of the panels are grid score (G), mean firing rate in spikes/s (m), and peak firing rate, also in spikes/s (p).Top row: trajectory of animal in the environment (black lines) and spike location (red markers).Middle row: firing rate map of the neuron (color scale bars in spikes/s).Bottom row: autocorrelogram of the rate map (color bar ranges: (A), −74% [percent is relative to maximum correlation] to 100%; (B), −51% to 100%; (C), −60% to 100%).

Figure 2 .
Figure 2. Distributions of grid cell metrics in real and simulated recordings.The data points analyzed are each cell's mean of every value type (e.g., 7 field size values from a cell).Red lines are at median data points, blue rectangles represent data between the 25th and 75th percentiles, and black whiskers extend to the first and last values in the sets of data points.

Figure 2 .
Figure 2. Distributions of grid cell metrics in real and simulated recordings.The data points analyzed are each cell's mean of every value type (e.g., 7 field size values from a cell).Red lines are at median data points, blue rectangles represent data between the 25th and 75th percentiles, and black whiskers extend to the first and last values in the sets of data points.

Figure 3 .
Figure 3. Effects of neuronal and synaptic model parameters on stellate cell grid scores.The plots in panels (A-D) were obtained by testing nine values per parameter and interpolating scores between those values.Black outlines encompass regions with grid scores above the 0.2 threshold (cf.color scale bars on the right of the plots).Teal outlines indicate the biologically realistic parameter ranges.Green stripes represent intersection areas with both acceptable grid scores and biologically realistic

Figure 3 .
Figure 3. Effects of neuronal and synaptic model parameters on stellate cell grid scores.The plots in panels (A-D) were obtained by testing nine values per parameter and interpolating scores between Int. J. Mol.Sci.2023, 24, x FOR PEER REVIEW 7 of 20

Figure 4 .
Figure 4. Neuron type-specific spiking dynamics.(A) Raster plot of the spiking of all neurons during two seconds of simulation of animal data.This activity was from the simulation in Figure 1B.The red line represents the movement path of a mouse that was used as velocity input into the simulation: the corresponding right y-axis label represents two-dimensional environment position converted to one-dimension with the formula ( − 1) *   +  , where y is the y-axis location position, x is the x-axis location position, and   is the length of the x-axis in cm.(B) Population activity from the combined spiking of all simulated neurons in MEC LII, namely grid cells, INs, and CjCs.(C) Power spectral densities of neural spiking during the first minute of the simulation.The spiking in each plot is from all neurons within each color-coded neuron type group.

Figure 4 .
Figure 4. Neuron type-specific spiking dynamics.(A) Raster plot of the spiking of all neurons during two seconds of simulation of animal data.This activity was from the simulation in Figure 1B.The red line represents the movement path of a mouse that was used as velocity input into the simulation: the corresponding right y-axis label represents two-dimensional environment position converted to one-dimension with the formula (y − 1) * x size + x, where y is the y-axis location position, x is the x-axis location position, and x size is the length of the x-axis in cm.(B) Population activity from the combined spiking of all simulated neurons in MEC LII, namely grid cells, INs, and CjCs.(C) Power spectral densities of neural spiking during the first minute of the simulation.The spiking in each plot is from all neurons within each color-coded neuron type group.

Figure 5 .
Figure 5. Simulated voltage and current.The simulation settings were those used to create results in Figure 1B.(A) Rate map of all grid cell spiking within a 400 ms time span.The color bar on the right is the count of spikes in that time period.(B-D) Recordings during that same time span of a stellate cell and IN with x-and y-axis position (2, 1) in their respective neural layers.

Figure 5 .
Figure 5. Simulated voltage and current.The simulation settings were those used to create results in Figure 1B.(A) Rate map of all grid cell spiking within a 400 ms time span.The color bar on the right is the count of spikes in that time period.(B-D) Recordings during that same time span of a stellate cell and IN with xand y-axis position (2, 1) in their respective neural layers.

Figure 6 .
Figure 6.Overview of simulation components and methods.Green arrows represent components utilized in the simulation.Blue arrows represent directions of connections between the neural layers.The striped blue arrow indicates a connection between layers that may only exist indirectly in nature.Abbreviations: MEC, medial entorhinal cortex; LII, layer 2; CA1, cornu ammonis 1.

Figure 6 .
Figure 6.Overview of simulation components and methods.Green arrows represent components utilized in the simulation.Blue arrows represent directions of connections between the neural layers.The striped blue arrow indicates a connection between layers that may only exist indirectly in nature.Abbreviations: MEC, medial entorhinal cortex; LII, layer 2; CA1, cornu ammonis 1.

Table 2 .
Connection parameters.The values here are those used to produce the results in Figure1B.Values are rounded to the 3rd decimal place.The number of connections column has numbers rounded to the nearest interval.a The connectivity was specific to each grid cell and varied within the standard deviation listed.b The individual number of connections reflects the CS distribution (details in Supplementary Materials IV).c These values were adjusted for the scale of the simulation (details in Supplementary Materials IV).