Selforganization of modular activity of grid cells

Abstract A unique topographical representation of space is found in the concerted activity of grid cells in the rodent medial entorhinal cortex. Many among the principal cells in this region exhibit a hexagonal firing pattern, in which each cell expresses its own set of place fields (spatial phases) at the vertices of a triangular grid, the spacing and orientation of which are typically shared with neighboring cells. Grid spacing, in particular, has been found to increase along the dorso‐ventral axis of the entorhinal cortex but in discrete steps, that is, with a modular structure. In this study, we show that such a modular activity may result from the self‐organization of interacting units, which individually would not show discrete but rather continuously varying grid spacing. Within our “adaptation” network model, the effect of a continuously varying time constant, which determines grid spacing in the isolated cell model, is modulated by recurrent collateral connections, which tend to produce a few subnetworks, akin to magnetic domains, each with its own grid spacing. In agreement with experimental evidence, the modular structure is tightly defined by grid spacing, but also involves grid orientation and distortion, due to interactions across modules. Thus, our study sheds light onto a possible mechanism, other than simply assuming separate networks a priori, underlying the formation of modular grid representations.


| I N TR ODU C TI ON
An accurate representation of their own position in space is thought to be a fundamental requirement for animals to perform most of their activities, which entail moving around in the environment. The spatial cognitive map theory outlines such a flexible representation (Tolman, 1948), in which an animal creates schemas that incorporate spatial knowledge. For years, the cognitive map has been associated with the hippocampus (O'Keefe & Nadel, 1978), where, in rodents, principal cells have been seen to represent specific locations in paradigms that require the animal to move in a given environment (O'Keefe & Dostrovsky, 1971). In 2005, the discovery of grid cells, one synapse upstream to the hippocampus, in the medial entorhinal cortex, appeared to reveal a much more geometrical, rigid implementation of the cognitive map concept (Hafting, Fyhn, Molden, Moser, & Moser, 2005;Moser, Kropff, & Moser, 2008). In small flat environments, these grid cells fire at multiple locations with a regular hexagonal structure; unlike place cells, they fire in all environments and maintain rigid coactivity relations (Fyhn, Hafting, Treves, Moser, & Moser, 2007). Reminiscent of a stack of sheets of millimeter paper stapled together, nearby cells display the same hexagonal pattern, although phase-shifted, in one environment and across environments. At the same time, grid spacing increases gradually from the dorsal to the ventral end of entorhinal cortex, suggesting that the same representation is replicated at multiple scales (Brun et al., 2008). At a given scale, the observed rigidity, that is, that phase relationships between nearby grid cells are preserved during remapping (Fyhn et al., 2007), makes one think of strong local networks, which have been hypothesized to approximate a metric system .
Do these parallel representations occur at any possible scale? By systematic analysis of multiple unit recordings in single animals, grid cell activity was found to be organized in discrete modules, with only four to five spatial scales observed in each animal (Stensola et al., 2012). A simple explanation of such modular organization is that it is pre-wired: grid cells may emerge already affiliated with one of four to five networks, within which units interact strongly, while they interact weakly, if at all, across networks. A model that offers such an account, where units learn their affiliation among different pre-defined scales, has already been presented (Pilly & Grossberg, 2014). We want to assess the alternative hypothesis that the modular structure is not prewired, and instead self-organizes. Our approach is based on a model, the self-organization model previously proposed for the formation of grid cell activity (Kropff & Treves, 2008). Here, we consider the concomitant self-organization of the modular structure for an extended network of grid units.

| T HE M OD EL
In the model, three ingredients are necessary for the development of a coherent grid pattern. First, a cellular adaptation mechanism suppresses firing within a certain time window after spells of intense firing; for a rodent moving in an environment at a certain average speed, and single units responding primarily to location-specific inputs, this tends to suppress firing in a ring around each region where the inputs happen to be stronger. In other words, adaptation acts indirectly as a spatial bandpass filtering. Second, the self-organization of multiple fields into the hexagonal pattern requires the gradual reinforcement of those initially random fields that happen not to suppress each other. This can be achieved by Hebbian plasticity on the afferent connections to a single unit, and slowly produces hexagonal grids at the single unit level, with random orientations. Third, different grid units endowed with the same adaptation time constants produce, in a flat environment, distinct but aligned hexagonal patterns. The addition of collateral interactions, modulated by the direction selectivity of each unit (Kropff & Treves 2008), has been shown to have the capability to align the grid patterns into a common orientation for groups of units ("modules"), an alignment which is preserved during remapping (Si, Kropff & Treves, 2012).
Strictly speaking, this network model accounts for the conjunctive (grid x head direction) activity found in layer III and deeper layers of mEC (Sargolini et al., 2006); however, it can be extended to include also, as a downstream population of units, those expressing pure grid activity, as observed in layer II, without any head direction modulation ). In the model, the information about the animal's position may be delivered by different sources, including a major feedback projection from CA1 to deep layers in mEC (Bonnevie et al., 2013;Moser et al., 2014). As the self-organization character of the model requires Hebbian learning to refine and reinforce a specific lay-out of the multiple fields of each unit, it depends on input statistics, whether from CA1 or elsewhere, and ultimately on the exploration of the environment itself. This is what makes the model predict different types of lay-out as the geometry or the dimensionality change, in comparison to a simple planar Euclidean geometry (Stella, Si, Kropff, & Treves, 2013;Urdapilleta, Troiani, Stella, & Treves, 2015;. As previously discussed , the dependence of the model on learning conditions may be critical particularly during development. In this sense, the timescale to develop a population-based grid structure in the grid or conjunctive layer is in the order of days of active exploration.
Detailed information about the model is given in the Appendix.

| Prewired modular organization
In a simple scenario, the self-organization of different modules may reduce to the parallel self-organization of largely separate networks FIGUR E 1 Homogeneous networks with different grid spacing. (a) The parameter b 1 in the self-organization model largely controls this property. Networks with large values of b 1 are tightly organized in a common representation with small grid spacing, whereas networks with a small value of b 1 express a less coherent representation with large grid spacing and more variability. internally homogeneous in terms of the single-unit properties, including grid spacing. A common orientation is produced by collateral interactions within the units of the module, without, or with weak projection to other modules. In Figure 1, we show the simulated emergence of typical grid representations in homogeneous networks, each assigned specific adaptation properties (here restricted to the manipulation of the b 1 inverse time constant). The parameter b 1 , which in the model sets the timescale of the adaptation process, controls the grid spacing developed by the units (see Figure 1a). A common orientation for all units in the ensemble and diverse phase relationships between them are organized by the collaterals and by individual head direction modulation. When isolated units are simulated with directional modulation, there is a slow temporal drift of the firing map along the opposite direction of the associated head-direction selectivity. This is a consequence of the adaptation properties (modeling fatigue), exerting a larger effect in a given direction for each conjunctive unit. The effect is essentially the same as the one reported for the experience-dependent properties of the firing fields in place cells (Mehta, Barnes & McNaughton, 1997;Mehta, Quirk & Wilson, 2000). When conjunctive units interact with each other, drift is reduced because neighboring units "push" each other along all directions, eventually establishing a kind of strained or frustrated fixation, except when the directional modulation is very strong. For units developing a larger grid spacing (small b 1 ), individual fields are also larger and the drift effect lasts longer. To compensate this (assuming it makes sense for the different modules to reach stability at about the same time), as the grid spacing increases, the head-directional modulation should be less tuned. This weaker modulation at large grid spacing has an impact on the variability observed in the geometrical properties of the pattern, which generally display a looser overall alignment. This variability, both in grid spacing and in alignment, is also observed in experimental studies (Stensola et al., 2012). The (putative) computational requirement of variable head direction modulation strength exactly matches what was recently found about the topographical organization of head direction selectivity along the dorso-ventral axis in the mEC (Giocomo et al., 2014). Representative examples of single-unit grid activity displayed by the parallel development of networks are shown in Figure 1b, as well as the prototypical auto-correlograms used to characterize grid patterns.
As expected, therefore, the parallel development of noninteracting networks successfully explains modular formation, by relying on pre-wired compartments.

| Modular self-organization
As an alternative, we consider the hypothesis that the self-organization of the modules may occur by a spontaneous breaking of a continuous distribution of time constants in a single heterogeneous network, via interactions that produce discrete domains of grid spacing. This can be considered analogous, at the level of firing maps, to the synchronization of interacting oscillators exhibited by the Kuramoto model (Acebron, Bonilla, Perez Vicente, Ritort, & Spigler, 2005;Kuramoto, 1984). In that case, due to the interactions, elements with different natural frequencies acquire a common rotation, if the coupling is strong enough. In our model, instead, we consider a smooth gradient in the parameter b 1 , possibly accompanied by other variable cellular properties along the dorso-ventral axis; and collateral connections that provide for interactions within a given range along the same axis. As shown in Figure 1, each specific value of b 1 in isolation (or in an ensemble of identical elements) would imply the emergence of a specific value of grid spacing. Given a population of N mEC units, labeled by i (i 5 1, . . ., N mEC ), to avoid boundary effects, we find it convenient to take the index i to run along the axis twice from the ventral to the dorsal and back to the ventral end ( Figure 2). Then units lie effectively on a ring, and each interacts via collateral connections with N lat other neighboring units on each side. To simplify, interactions are dependent on distance along this unidimensional circular axis, and thus interactions across the two branches joining ventral to dorsal ends are not considered. We further assume, in our model, that the "natural" grid spacing k (i.e., the grid spacing developed in isolation or in an ensemble of identical units) should be linearly related to the distance D from the dorsal end of mEC (this is a conservative assumption, as will be clear below). Since one can characterize the relationship between the grid spacing, k, and the parameter b 1 as, roughly, valid for the entire range except possibly very large values of k (see by (Giocomo et al., 2014), we impose progressively sharper headdirection selectivity towards the dorsal end, and also decreasing collateral-to-feedforward weights. In detail, Inverse of the width of the tuning function :

| R E SU LTS
The self-organization of modular activity in a heterogeneous network of interacting conjunctive units is evident from the clustering of individual grid patterns in a few common representations. Grid patterns were identified through standard procedures, based on spatial autocorrelograms (Hafting et al., 2005;Sargolini et al., 2006;Stensola et al., 2012). In brief, units with grid scores above a chosen threshold (see Figure 3) were selected for subsequent analysis. From the spatial autocorrelograms, we identified the six fields closest to the origin and selected three as the "axes" of the grid representation. Following (Stensola et al., 2012), only units surpassing fixed criteria on the interaxis angle and the distortion between individual grid distances (distances from the origin to the respective fields) were further analyzed.
After this filter, the grid spacing of each grid unit was measured as the arithmetic mean of grid distances along the three axes. As shown in The ratios of the grid spacing between two consecutive modules are not fixed, but vary between 1.27 and 1.4, and are broadly consistent with experimental data.
The clustering indicated by the distribution of grid spacing around certain values can also be revealed if we rank units according to their grid spacing, and then consider the absolute difference of grid spacing between units. This ordered bi-dimensional map is shown in Figure 3c and clearly exhibits a community structure, where units belonging to the same module have a small difference in grid spacing (white color).
Colored corners indicate the limits of each module. Quantitatively, the URDAPILLETA ET AL.

| 1207
"discreteness" was defined by (Stensola et al., 2012) as the standard deviation of the statistics given by the collection of unit counts in all bins of the (binned) distribution of grid spacing, for a particular bin width. In agreement to that study, we found a consistent large discreteness value, irrespective of the bin width, clearly separated from a control uniform-like distribution (not shown).
In addition to grid spacing, the orientation of the hexagonal pattern is essential to define a common representation and, therefore, where they have found this coherent distortion but also a bi-modal structure of the elliptical organization. Because the distortion does not provide information for the modular architecture in our simulations, we rely on spacing and orientation to define modules with minimal arbitrary manipulation. This is achieved with a k-means clustering selection procedure, with four modules set ad hoc, k 5 4 (although other values were also considered). As the procedure does not produce a unique deterministic outcome, we run the analysis 300 times and evaluated the performance of each repetition by means of the population-based mean of the so-called silhouette coefficient Cs, see (Stensola et al., 2012) for further explanation. Based on the highest mean Cs, the best clustering outcome was selected and the mean properties of each cluster or module, as well as its cardinality, were obtained. A summary is given in Table 1. Properties obtained by the k-means procedure confirm the modular structure already revealed by the analysis based on the different histograms at the population level.
Because all previous results may critically depend on the threshold imposed on the grid score for further consideration of individual units, we analyzed the effects of varying this threshold, which in turn produces smaller/larger datasets. Figure 3g shows the smoothed histograms for different threshold values, whereas Table 2 indicates the number of units and the mean Cs obtained with the (best selected) k-means procedure. Larger values of mean Cs quantifies domains for data points that are better resolved by clusters, with increasingly less degree of overlap between them. A mean Cs of 0.5 may be considered a consistent separation. As noted, the modular structure is reliable, beyond the particular value of the threshold. Naturally, as this value increases, fewer units remain in the dataset but, on the other hand, the quality of the definition of the clusters is enhanced.

| D I SCUSSION
Our simulations point at the possibility that grid modularity may not be hardwired but may result from a simple self-organization process at the local network level. Similarly, the adaptation model (Kropff & Treves, 2008) has earlier demonstrated the possibility that grid structure may emerge from self-organization at the single unit level.
There is a long list of features in the simulated results that could be compared with the experimentally observed modularity. The degree of abstraction of the model from the biophysics and connectivity of real entorhinal networks, and the many simplifications it entails, make a systematic study not so useful. One aspect, however, deserves atten- Note also that in a recent report an inverse progression has been observed in bats (Heys et al., 2016). In this article, the authors performed whole-cell patch recordings along the dorsal-ventral axis of EC in bats and, surprisingly, found that the sag response properties and the resonance properties recorded in layer II neurons of entorhinal cortex in the Egyptian fruit bat demonstrate an inverse relationship along the dorsal-ventral axis compared with the rat.

ACKNOWLEDGMENTS
Extensive discussions with other participants in the GRIDMAP con- where the first term in the brackets determines the feedforward contribution from place units, located one synapse downstream when looking at the feedback branch in the CA1-mEC connectivity layout, and the second term the contribution resulting from collateral interactions. The factor f ui (/) multiplying both contributions is the head direction modulation characterizing each conjunctive unit, specified by its preferred direction u i and the inverse of the tuning width m according to with a basal modulation, c 5 0.2. As mentioned in the main text, the inverse tuning width m varies along the dorso-ventral axis, whereas each individual u i is assigned by randomly sampling the 2p range when establishing unit properties before the simulation starts.

The feedforward contribution in Equation A4
is determined by the activity of "place units" (see Kropff & Treves, 2008, for a discussion of this simplification), each one approximated by a Gaussian function centered in a pre-determined preferred firing location x j0 , with a width r p 5 5 cm, according to  weighted by Ŵ ik . These collateral connections are formulated with the prescription given in (Si et al., 2012), which mimics a long learning process matching individual conjunctive properties, here restricted to a local interaction in the arrangement considered for the heterogeneous population (see Figure 2). In detail, those connections in the range of interaction (|i 2 j| N lat 5 100) are set aŝ whereas those outside this range are zero. In this equation, the head direction selectivity of both pre-and post-synaptic units modulates weight strength, such that it is maximal when their tuning is aligned.
In addition, to mimic learning in coactive conjunctive units, each unit is assigned a spatial selectivity (a place field), in a location randomly selected in the environment (Kropff & Treves, 2008). The "coactivity" , this parameter controls which argument of the threshold function will become negative and therefore, which weights vanish. Thus, only units aligned and close in spatial selectivity between reverberated pre-and post-synaptic activity will interact through collateral connections. In other words, this procedure produces a sparse connectivity where a fraction of about 1/10 connections are effectively instantiated. Importantly, these connections are allowed within the neighborhood of any target unit, |i 2 k| 5 1, . . ., N lat (see main text). The relative contribution of collateral interactions to feedforward processing of "place unit" input is weighted by a factor q, see Equation A4, which here is taken to vary along the dorso-ventral axis (see main text and Figure 2).
The model described above is driven by the random exploration of a virtual agent, constrained in a circular environment of 2 m diameter and moving with a constant speed of 40 cm s 21 . At each time step, the running direction / is randomly updated from a Gaussian distribution with zero mean and standard deviation 0.2 radians, and then used to compute the coordinates of the future position, considering the speed constant. If this position lies beyond the limits of the circular environment, the trajectory is reflected. With this procedure, the arena is covered homogeneously. In a more realistic situation, where the agent covers the environment with a variable running speed, it is the mean velocity that determines the grid spacing, as the slow Hebbian learning process averages out fluctuations (Kropff & Treves, 2008).
The programming code of the model (written in C) and the population results analyzed here are available in a public repository: "ModelDB, accession code 231392", http://senselab.med.yale.edu/ModelDB/ showModel.cshtml?model5231392.