Inferring circuit mechanisms from sparse neural recording and global perturbation in grid cells
Abstract
A goal of systems neuroscience is to discover the circuit mechanisms underlying brain function. Despite experimental advances that enable circuitwide neural recording, the problem remains open in part because solving the ‘inverse problem’ of inferring circuity and mechanism by merely observing activity is hard. In the grid cell system, we show through modeling that a technique based on global circuit perturbation and examination of a novel theoretical object called the distribution of relative phase shifts (DRPS) could reveal the mechanisms of a cortical circuit at unprecedented detail using extremely sparse neural recordings. We establish feasibility, showing that the method can discriminate between recurrent versus feedforward mechanisms and amongst various recurrent mechanisms using recordings from a handful of cells. The proposed strategy demonstrates that sparse recording coupled with simple perturbation can reveal more about circuit mechanism than can full knowledge of network activity or the synaptic connectivity matrix.
https://doi.org/10.7554/eLife.33503.001Introduction
In systems neuroscience we seek to discover how neural responses and complex functionality can emerge from the dynamical interactions of neurons in circuits. For instance, the circuit mechanisms that give rise to orientation tuning in primary visual cortex have been closely studied for the better part of a century (Hubel and Wiesel, 1959). Despite these efforts, arbitrating between between different candidate mechanisms has been difficult. Our experimental tools are typically observational: Neurons are recorded, often during a behavior, in increasing numbers today (Dombeck et al., 2010; Ahrens et al., 2012; Ziv et al., 2013; Dunn et al., 2016). Our theoretical models usually run in the ‘forward’ direction: We build hypothesized circuits to reproduce the observed activity data. Because there often is a manytoone mapping from plausible models to neural activity, it is difficult to know which model more accurately describes the underlying system. For this reason, it remains unsettled whether – to return to a familiar example – orientation tuning arises mostly from selective feedforward summation of inputs or lateral interactions (RivlinEtzion et al., 2012; Kim et al., 2014; Takemura et al., 2013; Ferster and Miller, 2000; Sompolinsky and Shapley, 1997).
Here, we show that grid cells (Hafting et al., 2005) provide a unique opportunity to understand cortical circuit mechanism, when coupled with a novel approach for doing so. The promise of our approach lies in the fact that (1) it is not merely observational but rather relies on perturbation, and (2) it provides a novel theoretical measure (the ‘distribution of relative phase shifts’ or DRPS) along which several competing feedforward and recurrent grid cell models can be distinguished with the perturbative experiments.
The structure of grid cell responses – with their periodic tuning to 2D space – makes the system particularly amenable to dissection, as we will see below. Grid cells have already yielded insight into their underpinnings: All cells with a common spatial tuning period remain confined to a single 2D manifold in activity space, and this manifold is invariant over time even when grid cell tuning curves deform as the animals are moved between novel and familiar environments (Yoon et al., 2013; Fyhn et al., 2007), as well as during REM and nonREM sleep (Gardner et al., 2017; Trettel et al., 2017). These findings imply the existence of a 2D continuous attractor dynamics within or feeding into the grid cell circuit.
Many models reproduce the spatially periodic responses of individual grid cells or groups of cells (Fuhs and Touretzky, 2006; Burak and Fiete, 2006; McNaughton et al., 2006; Hasselmo et al., 2007; Burgess et al., 2007; Kropff and Treves, 2008; Guanella et al., 2007; Burak and Fiete, 2009; Welday et al., 2011; Dordek et al., 2016). These include models in which the mechanism of grid tuning is a selective feedforward summation of spatially tuned responses (Kropff and Treves, 2008; Dordek et al., 2016; Stachenfeld et al., 2017), recurrent network architectures that lead to the stabilization of certain population patterns (Fuhs and Touretzky, 2006; Burak and Fiete, 2006; Guanella et al., 2007; Burak and Fiete, 2009; Pastoll et al., 2013; Brecht et al., 2014; Widloski and Fiete, 2014), the interference of temporally periodic signals in single cells (Hasselmo et al., 2007; Burgess et al., 2007), or a combination of some of these mechanisms (Welday et al., 2011; Bush and Burgess, 2014). They employ varying levels of mechanistic detail and make different assumptions about the inputs to the circuit. Because exclusively singlecell models lack the lowdimensional networklevel dynamical constraints observed in grid cell modules (Yoon et al., 2013), and are further challenged by constraints from biophysical considerations (Welinder et al., 2008; Remme et al., 2010) and intracellular responses (Domnisoru et al., 2013; SchmidtHieber and Häusser, 2013), we do not further consider them here. The various recurrent network models (Fuhs and Touretzky, 2006; Burak and Fiete, 2006; McNaughton et al., 2006; Guanella et al., 2007; Burak and Fiete, 2009; Brecht et al., 2014) produce single neuron responses consistent with data and further predict the longterm, acrossenvironment, and acrossbehavioral state cell–cell relationships found in the data (Yoon et al., 2013; Fyhn et al., 2007; Gardner et al., 2017; Trettel et al., 2017), but are indistinguishable on the basis of existing data and analyses. Here we examine ways to distinguish between a subset of grid cell models, specifically between the recurrent and feedforward models, and also between various recurrent network models. We call this subset of models our candidate models. Our goal is not to provide new models of grid cell activity, but rather to show, through theory and modeling, how the candidate models could be feasibly distinguished through experiment.
The candidate models form a diverse set, with differences that carry important implications for mechanism and for how the network could have developed from plasticity mechanisms. The candidates first broadly partition into recurrent and feedforward models, depending on whether the dynamics that originate spatial tuning and velocity integration are within (recurrent) or upstream (feedforward) of the grid cell layer. Recurrent models further partition on the basis of two key features: topology (periodic or not) and locality of connectivity (from local to global).
Among recurrent models, the first candidate models are aperiodic networks (Figure 1a) (Burak and Fiete, 2009; Widloski and Fiete, 2014): Network connectivity has no periodicity (flat, holefree topology) and it is purely local (with respect to an appropriate or ‘topographic’ rearrangement of neurons only nearby neurons connect to each other). Despite the aperiodic and local structure of the network, activity in the cortical sheet is periodically patterned (under the same topographic arrangement). In this model, coactive cells in different activity bumps in the cortical sheet are not connected, implying that periodic activity is not mirrored by any periodicity in connectivity. Interestingly, this aperiodic network can generate spatially periodic tuning in single cells because, as the animal runs, the population pattern can flow in a corresponding direction and as existing bumps flow off the sheet, new bumps form at the network edges, their locations dictated by inhibitory influences from active neurons in other bumps (Figure 1e). From a developmental perspective, associative learning rules can create an aperiodic network (Widloski and Fiete, 2014), but only with the addition of a second constraint: Either that associative learning is halted as soon as the periodic pattern emerges, so that strongly correlated neurons in different activity neurons do not end up coupled to each other, or that the lateral coupling in the network is physically local, so that grid cells in the same network cannot become strongly coupled through associative learning even if they are highly correlated, because they are physically separated. In the latter case, the network would have to be topographically organized, a strong prediction.
Second are fully periodic networks (Figure 1c) (Guanella et al., 2007; Fuhs and Touretzky, 2006; Pastoll et al., 2013; Brecht et al., 2014). The network is topologically a torus with periodic boundary conditions between the pairs of opposite edges, and connectivity is global: There is no neural rearrangement under which network connectivity will be local. It is mathematically equivalent to view this network as having a single activity bump (Burak and Fiete, 2006; Guanella et al., 2007) or having multiple periodically arranged bumps with interbump connections (Burak and Fiete, 2009). In this network, periodic connectivity underlies periodic activity. Developmentally, a fully periodic network would naturally arise if associative plasticity continued postpattern formation, so that the topology of activity and connectivity would come to mirror each other (Widloski and Fiete, 2014).
Third are partially periodic networks (Burak and Fiete, 2009) (Figure 1b) with periodic boundary conditions (torus topology) but only local connectivity on the torus after appropriate rearrangement of neurons. In this model, neural activity on the cortical sheet is multipeaked and periodic (under appropriate rearrangement). Conceptually and developmentally, these networks are the strangest: None of the connectivity in the bulk of the network reflects the periodic nature of activity within it, except for the connectivity necessary to connect together neurons across the edges of an initially aperiodic sheet of cells. The wiring of this ‘edge’ subset of neurons must, unlike the rest of the cells, depend on details of the periodic activity pattern to make sure that opposite edge bumps are ‘aligned’ before joining (Figure 1—figure supplement 2). It is unclear how activitydependent plasticity rules, which could wire together faraway edge neurons based on activity, would refrain from doing the same for the rest, to maintain otherwise local connectivity.
The fourth potential combination of topology and locality is not permitted: it is not possible to obtain gridlike activity from neurons with global connectivity (and singlebump activity) but aperiodic boundaries (topologically flat holefree networks), Figure 1d.
Feedforward models of grid cell activity form a robust and growing set. In these models, grid cells merely sum and transform with a pointwise nonlinearity inputs that are already spatially tuned with roughly uniform coverage and resolution across the environment (Figure 1f) (Kropff and Treves, 2008; Welday et al., 2011; Mhatre et al., 2012; Bush and Burgess, 2014; Hasselmo and Brandon, 2012; Dordek et al., 2016; Stachenfeld et al., 2017); thus, it is implicitly assumed that the upstream inputs to grid cells have already performed path integration. These feedforward models, which we propose could be distinguished from recurrent models with the proposed perturbative approach, themselves segment into two major varieties. One type (Welday et al., 2011; Mhatre et al., 2012; Bush and Burgess, 2014; Hasselmo and Brandon, 2012) generates lowdimensional grid cell population activity across environments (e.g., in Welday et al., 2011, three upstream circuits, each a 1D continuous attractor network, integrate one component of animal velocity aligned to each of the three primary directions of a triangular lattice; the combined response is 2dimensional, and preserved across environments; other models of this type differ in details but are similar in this regard), as predicted also by the recurrent models and found in the data (Yoon et al., 2013). In the second type (Kropff and Treves, 2008; Dordek et al., 2016; Stachenfeld et al., 2017), the grid cell pattern for an environment depends on the place cell pattern for that environment. Thus, when the place cell representations remap across environments, the model grid cells will not preserve their spatial relationships.
Our candidate models are the set of recurrent and feedforward models described and cited above. They are architecturally and mechanistically distinct in ways both large and subtle: they differ in whether grid cells or their upstream inputs are performing velocitytolocation integration, in whether spatial patterning originates wholly or only partly within grid cells, and in the structure of their recurrent circuitry. As already noted, some of the subtleseeming structural differences have important implications for circuit development: different connectivity profiles and topologies require distinct models of plasticity and experience during circuit formation (Widloski and Fiete, 2014). Nevertheless, candidate recurrent and feedforward models that exhibit approximate 2D continuous attractor dynamics are difficult to distinguish on the basis of existing data.
As we discuss at the end, neither complete single neuronresolution activity records nor complete single synapseresolution weight matrices (connectomes) will fully suffice to distinguish between the candidate models because they are observational or correlative techniques: they do not probe the causal origin of the observed responses.
We show how it is nevertheless possible to gain surprisingly detailed information about the grid cell circuit from a feasible perturbationbased experimental strategy, enough to discriminate between the candidate models.
Results
A perturbationbased probe of circuit architecture
The question of mechanism is focused on a prespecified set of neurons or local circuit: Is the observed lowdimensional grid cell activity primarily based on recurrent interactions within the set, and how, or is it inherited from feedforward drive originating outside this set? We refer to a perturbation as simple, lowdimensional and global in this context if it affects most cells within this set without regard to their individual functional identities, and does not affect those outside. In what follows, we consider the set to consist of all grid cells and conjunctive cells in one (or more) grid modules (Stensola et al., 2012), as well as the interneurons that surround them; toward the end we discuss the effects of perturbing subpopulations or bigger sets.
The central idea is as follows: Globally perturbing either the timeconstant of neurons or the gain of recurrent inhibition is predicted to affect cell–cell spatial tuning relationships in candidate models in a specific way that can be robustly observed and characterized from ultrasparse sampling of neurons in the network, and the predicted effects differ across candidate mechanisms.
To present the idea, we consider a thought experiment on the aperiodic recurrent network models. We will retake the larger perspective, of discriminating between the various model categories, immediately afterward. In aperiodic models, perturbing the gain of recurrent inhibition or the timeconstant of neurons will induce a shift in the period of the internal population pattern (Figure 2—figure supplement 1). Let us quantify the change in period by the population period stretch factor, $\alpha =\frac{{\lambda}_{pop,post}}{{\lambda}_{pop,pre}}1$ (where ${\lambda}_{pop,pre}$ is the preperturbation population period). Without loss of generality, suppose that the focus of pattern expansion is at the left edge, Figure 2a (blue: original pattern, red: expanded pattern). Each neuron can be assigned a population phase with respect to the period of the population pattern: If the phase at the left edge is called 0 (again without loss of generality), neurons lying at integer multiples of the original period also originally had a phase of 0 (Figure 2b). However postexpansion, the population phase of a neuron originally one period away from the left edge is no longer zero (Figure 2a,b). Let us call the shift in the population phase of this neuron one ‘quantum’ (Figure 2b), and denote it by $\mathrm{\Delta}$. The quantum of shift must be $\mathrm{\Delta}=1\frac{{\lambda}_{pop,pre}}{{\lambda}_{pop,post}}=\frac{\alpha}{1+\alpha}$ ($\approx \alpha $ for small perturbations). A neuron $K$ periods away will shift in phase by $K$ quanta, Figure 2b. If there are $M$ bumps in the population pattern, the largest shift will be $M$ quanta, or $M\mathrm{\Delta}$ (modulo 1). If we construct a distribution of shifts in population phase pre to postexpansion for cells across the network, the distribution will be quantal, with $2M$ peaks (assuming the biggest phase shift, $M\mathrm{\Delta}$, is less than $1/2$, because phase is a periodic variable that we parameterize as running between $1/2$ to $1/2$; this condition can be met by keeping the perturbation small, such that $\mathrm{\Delta}<1/(2M)$), Figure 2—figure supplement 2 and Figure 2c. In other words, for small perturbations, the number of peaks in this distribution is predicted to be twice the number of bumps in the original population pattern. We will call this distribution of relative phase shifts the DRPS.
Practically, however, the grid cell network might not be topographically wellordered on a sufficiently fine scale (Heys et al., 2014), and one cannot simply image the population response and expect to read off pattern phases for each cell as in Figure 2a,b.
Fortunately, the distribution of shifts in the difficulttoobserve population phases of cells, based on instantaneous and topographically ordered population activity snapshots, is mirrored in the distribution of shifts in the relative phase of the straightforwardtoobserve and timeaveraged spatial tuning curves of cells (Figure 2d). Consider a pair of cells one population period apart preperturbation, so they have the same population phase (circle, square or square, triangle in Figure 2a). These cells are coactive and have the same spatial tuning curves, and thus a relative spatial tuning phase of 0 (circle, square or square, triangle in Figure 2d, top). Postperturbation, their spatial tuning curves will be shifted relative to each other by the same amount as the shift in their individual population phases (circle, square or square, triangle in Figure 2d, bottom). In other words, cells one bump apart in the original population pattern will exhibit one quantum of shift in their relative spatial tuning. The relative phase of spatial tuning for a pair originally separated by $K$ periods will shift by $K$ quanta postperturbation (e.g., circle, triangle in Figure 2d, bottom: spatial tuning curves shift by two quanta in phase because these cells were two periods apart in the original population pattern).
This series of theoretical observations leads us to construct a predicted distribution of relative phase shifts (DRPS) from all pairs of neurons, Figure 2e. The DRPS is quantal and has the same number of peaks as the distribution of shifts in population phase (Figure 2c). Indeed, multiplying the number of peaks in the multimodal DRPS by $1/2$ gives the number of bumps in the original population pattern, if the quantal shift size is sufficiently small. The DRPS is a property of patterning in an abstract space, independent of how neurons are actually arranged in the cortical sheet. It can be obtained from the spatial tuning curves of cells recorded simultaneously through either conventional electrophysiology or imaging. As we show later, a robust estimate of the full DRPS can be obtained from only a handful of cells.
In 2D, relative phase is a vector. The two components are each computed simply as in 1D, along each of the two principal axes of the spatial tuning grid. For an aperiodic network, for small enough perturbations, the total number of bumps in the population pattern can be inferred to be a quarter of the product of the number of peaks in the two relative phase shift distributions from the two components of the relative phase (Figure 2f–h).
Experimental knobs to modulate the population pattern
To generate the DRPS in experiment and use it to distinguish between grid cell models requires an experimental knob that can be turned to change the period of the population activity pattern. Temperature is one potential knob: Cooling a biological system reduces reaction rates and increases timeconstants through the Arrhenius effect (Katz and Miledi, 1965; Thompson et al., 1985; Moser and Andersen, 1994; Long and Fee, 2008). However, existing models of grid cells are based on simplified ratebased or linearnonlinear Poisson (LNP) spiking units, and it is unclear which parameters to modify to correctly predict the effects of cooling the neural circuit: Varying a ‘neural’ timeconstant parameter in a recurrent network of simple units may or may not change the population pattern period, depending on whether PSP height is scaled together with the timeconstant change (Widloski and Fiete, 2014; Beed et al., 2013) or not. To better predict the effects of cooling on grid cell period, we constructed more detailed grid cell models using cortical HodgkinHuxley model neurons (Pospischil et al., 2008) whose parameters accommodate thermal effects (Hodgkin et al., 1952; Katz and Miledi, 1965).
The population period in aperiodic grid cell models built from HodgkinHuxley neurons is pulled in opposing directions by temperature modulations in ionchannel biophysics and synaptic signalling (Figure 2—figure supplement 4). However, the dominant influence on network response comes from the growth in the PSP timeconstant with cooling and results in an overall expansion of the population period (Figure 2—figure supplement 4).
This result allows us to conclude that the net effect of cooling the biological circuit should be an expansion in the period, if the circuit is recurrently connected and aperiodic. It also allows us to continue using simple ratebased and LNP spiking models because we can now interpret how to scale parameters as a function of temperature: It is most appropriate to scale the timeconstant inversely with temperature, while essentially keeping the PSP height fixed (Figure 2—figure supplement 4).
The strength of recurrent inhibition is another experimental knob. Unlike temperature, manipulating the gain of inhibitory synaptic conductances has a relatively unambiguous interpretation in grid cell models. Experimentally, the strength of inhibition might be modulated by locally infusing allosteric modulators that increase inhibitory channel conductances (e.g. benzodiazipines; Rudolph and Möhler, 2004 and personal communication with C. Barry). In both cortical HodgkinHuxley based models grid cell models (Figure 2—figure supplement 4) and ratebased models, a gain change in inhibitory conductances predicts a change in the period of the population pattern (Figure 2—figure supplement 4 and Moser et al., 2014, Widloski and Fiete, 2014).
To summarize, thermal perturbation (cortical cooling) and biochemical perturbation (drug infusions to alter the gain of recurrent inhibition) are two experimental manipulations that could, according to the models, alter the period of a recurrently formed population pattern and thus may act as appropriate knobs to enable the construction of the DRPS.
Discriminating among recurrent architectures
In dynamical simulations of the various plausible candidates (Materials and methods), the same global perturbations have different effects, resulting in distinct predicted DRPS’s. To generate maximally robust and easytointerpret predictions, we focus on how the candidate models differ with respect to one simple property of the DRPS: the overall width of its envelope.
In aperiodic networks (Figure 1a) with smooth boundaries for accurate integration (Burak and Fiete, 2009), an incremental increase in the strength of global perturbation results in incremental expansion of the population activity pattern (Figure 3a, red, and Figure 2—figure supplement 1) (see Widloski, 2015 for an analysis of boundary conditions and permitted number of peaks). Thus, the peaks in the DRPS will incrementally spread out, producing a DRPS envelope that gradually and smoothly widens with perturbation strength (Figure 3b–c, red). In addition, because the change in period is incremental when the perturbation strength is gradually increased, it may be possible to estimate the number of bumps in the population pattern by counting peaks in the DRPS.
Partially periodic networks (Figure 1c), unlike aperiodic networks, must because of their symmetry accommodate an integer number of complete activity bumps in a way that is perfectly periodic (Widloski, 2015). The bumps and spacings within a partially periodic network are identical and the population pattern (if the geometry of the 2D pattern is fixed) is characterized simply by the number of bumps, which is constrained to be an integer. Incrementally increasing the perturbation strength is thus predicted to first result in no change, followed by a sudden change when the network can accommodate an entire additional bump, Figure 3a (purple) (or an additional row of bumps in 2D, assuming the pattern does not rotate as a result of the perturbation; see Discussion). As soon as a new bump has been inserted into the population pattern, the phase shifts will be large even for cells in adjacent bumps, and the DRPS will be wide. To summarize, for partially periodic networks, incremental changes in perturbation strength are therefore predicted to result in a stepwise (stepping to maximal width) change in the DRPS (Figure 3b–d, purple).
Counting peaks to estimate the number of bumps in the underlying population pattern after a stepwise change in the DRPS will likely result in substantial underestimation: because the phases shift by a large step when a change occurs, if a shift of $M$ quanta already exceeds one cycle, the DRPS will not distinguish between an $M$bump and a $K$bump network ($K>M$; Figure 2—figure supplement 2 and e.g., Figure 3b: compare peaks in the solid and dashed lines for small and large perturbations, respectively).
Finally, in the fully periodic network (Figure 1b) the globally periodic connectivity completely determines the population period of the pattern, and changes in the neural timeconstants or network inhibition strength do not alter it (Figure 3a, blue). Thus, the same global perturbations that effected changes in the population period in the other recurrent models (Figure 3a, red and purple) have no effect in the fully periodic network. The DRPS is consequently predicted to remain narrow, unimodal, and peaked at zero (Figure 3b–c, blue).
Discriminating feedforward from recurrent architectures
If lowdimensional dynamics and spatially tuned responses first originate upstream of the perturbed set, then the perturbations will leave unchanged the spatial tuning phases of grid cells, preserving grid cell–grid cell relationships. This prediction holds even if grid cells play a role in constructing their particular patterns of spatial tuning, for instance by combining elements that are already spatially tuned as when stripetuned inputs are combined to generate 2D lattice responses (Mhatre et al., 2012; Welday et al., 2011; Bush and Burgess, 2014) (Figure 1f). Thus, for feedforward models, as for fully periodic (recurrent) networks, the DRPS is predicted to remain narrowly peaked at zero across a range of perturbation strengths (Figure 3c, dashed green line).
Further, perturbing grid cells but not their spatially patterned feedforward inputs will not affect their spatial tuning. By contrast, in all recurrent models (Figure 1a–c), perturbing the grid cell network induces a change in the efficacy with which feedforward velocity inputs drive the population phase over time, thus the spatial tuning period of cells is predicted to change even if the population period does not (as in fully periodic networks – see Figure 3—figure supplement 1), Figure 3d. This expansion in spatial tuning period with global perturbation strength is predicted to hold for all three recurrent network classes, and distinguishes fully periodic recurrent networks from feedforward ones.
Finally, in both feedforward and recurrent neural network models, the amplitude of the grid cell response will change in response to perturbation (Figure 3e). This universal prediction of amplitude change with perturbation can be used as an assay of whether the attempted global perturbation is in effect.
Data limitations and robustness
We consider two key data limitations. First, it is not yet experimentally feasible to record from all or even a large fraction of cells in a grid module. Interestingly, the proposed method is tolerant to extreme subsampling of the population: a tiny random fraction grid cells from the population (10 out of e.g. 1600 cells, or 0.6%) can capture the essential structure of the full DRPS, Figure 4a, including its overall width and the detailed locations of its multiple peaks. This robustness to subsampling is dramatically better than in statistical inference methods, where even ‘sparse’ methods can require $\sim 2$ orders of magnitude denser data (Soudry et al., 2015).
The second limitation arises from the limited accuracy with which spatial tuning and relative phase can be estimated from finite data. In tests that depend only on the width of the DRPS (e.g. Figure 3), this phase uncertainty is not a serious limitation.
Resolving the relative phase accurately becomes important when counting DRPS peaks to estimate how many bumps are in the underlying population pattern of a recurrent network. The spacing between DRPS peaks determines the required tolerance in relative phase (Figure 4—figure supplement 1). DRPS peak spacing (in the aperiodic network) increases with the stretch factor at small stretch factors (Figure 4b and Figure 4—figure supplement 1), but the stretch factor must still obey $\alpha \approx \mathrm{\Delta}<1/(2M)$ (where $M$ equals the larger of the number of bumps along the two dimensions of the population pattern; Figure 4b) to avoid underestimating the number of bumps in the population pattern.
Fortunately, it is possible to gain progressively better estimates of relative phase over time even if there is substantial drift in the spatial responses of cells, because relative phases remain stable in a fixed network (Yoon et al., 2013) (here ‘fixed’ means that a given perturbation strength is stably maintained). Many estimates of relative phase may be made from short pieces of the trajectory, and these estimates averaged together (similar to the methods used in Yoon et al., 2013 and Bonnevie et al., 2013).
To distinguish $M=5$ bumps per dimension based on structure within the DRPS requires a stretch factor $\alpha \approx \mathrm{\Delta}<1/(2M)=0.1$, and a phase noise of 0.02 or smaller (Figure 4—figure supplement 1), which would require an approximately 8 min recording (estimated from grid cell and trajectory data, http://www.ntnu.edu/kavli/research.gridcelldata), Figure 4c. Distinguishing seven bumps would require $\alpha \le 0.07$, phase noise less than $0.01$, and a 35 min recording.
In summary, the proposed method has high tolerance to subsampling and more limited tolerance to phase uncertainty, which can be reduced by averaging estimates over time.
A decision tree for experimental design
We lay out a decision tree with an experimental workflow for discriminating between disparate feedforward and recurrent grid cell mechanisms, all of which exhibit approximate 2D continuous attractor dynamics at the population level (Figure 5).
We start with the ‘specific’ approach, which, according to our model, has more discriminatory power than the ‘nonspecific approach’ described later. The experimental demands of this approach are to be able to stably induce a global perturbation in at least one grid module, and to do so at 2–3 different strengths. Critically, the perturbation must be one of the two specific types discussed above: a perturbation of the strength (gain) of inhibition in the network, or of the network time constants. The data to be collected are simultaneous recordings from several grid cells as the animal explores novel enclosures with no proximal spatial cues, over a $\ge 20$ minute trajectory.
First, before applying a perturbation, characterize spatial tuning (periods) and cell–cell relationships (relative spatial phases). Next, apply a series of 2–3 global perturbations of increasing strength. At each perturbation strength, characterize the spatial tuning of cells and cellcell relationships.
A change in the amplitude of the grid cells’ response across the different perturbations should signal that the perturbation is having an effect, regardless of underlying mechanism (Figure 5, first triangle on left).
If the different perturbation strengths do not cause a change in the spatial tuning periods of single cells (but the response amplitudes do change), it follows that velocity integration and spatial patterning are originating elsewhere, consistent with some feedforward mechanism (Figure 5, green). To confirm, verify that cell–cell relationships remain unchanged across perturbations, as also predicted for feedforward networks.
If there is a change in the spatial tuning period, characterize the cell–cell relationships in each perturbation condition. Plot the DRPS from each perturbed condition relative to the preperturbation condition, and quantify its width and if possible the separation between its peaks. If the DRPS width or peak separation increases steadily and smoothly with perturbation strength, that implies an aperiodic recurrent architecture (Figure 5, red). If the DRPS peak separation or width exhibits a steplike change, it is consistent with a partially periodic recurrent network (Figure 5, purple). Together with a change in the spatial tuning period, a DRPS that remains narrowly peaked at zero, with no change in width with perturbation strength, is consistent with a fully periodic network (Figure 5, blue).
Finally, if the network is either an aperiodic or partially periodic recurrent network, the number of peaks in the DRPS for each relative phase dimension is a lower bound on the quantity $2M$, where $M$ is the number of bumps in the population pattern along that dimension. If the stretch factor $\alpha $ times the number of bumps is smaller than 1/2 and the DRPS is multiply peaked the number of DRPS peaks should equal twice the number of population activity bumps along the corresponding dimension (Figure 5, final triangle and gray oval).
The ‘specific’ approach above should provide insight into the underlying dynamics of the system with respect to the candidate models, regardless of outcome. By contrast, a ‘nonspecific’ approach (Figure 5, dashed box) could do the same, but only for certain outcomes. Suppose that after a number of any type of perturbations to the system, with known or unknown underlying mechanisms and at a local or systemic scale, one measures the DRPS. If the DRPS does not exhibit multiple peaks then, because this outcome is consistent with many possibilities and the nature of the perturbation is not precisely known or controlled to change the inhibitory gain or neural timeconstant (the specific perturbations that provide higher discriminatory power), one cannot conclude anything about circuit architecture. On the other hand, if the DRPS after nonspecific perturbation does exhibit multiple, equispaced peaks, one can conclude with high confidence that the brain generates an underlying multibump population pattern through recurrent mechanisms with partially periodic or aperiodic structure. This is because a multipeaked DRPS is a highly specific outcome of recurrent patternformation dynamics.
Questions about experimental contingencies
Will it be possible to distinguish a ‘no effect’ result from an ‘experiment is not working’ result? Specific inhibitory gain or timeconstant perturbations are predicted to change the amplitude of the neural tuning curves relative to the unperturbed case in all models of Figure 1. An amplitude change from a specific perturbation is the signal that the experiment is working (Figure 5, first triangle on left).
If some circuit perturbation results in an amplitude change, can we learn something about the circuit from it? If a perturbation affects response amplitudes without primarily affecting inhibitory gains or neural time constants, it qualifies as a ‘nonspecific’ perturbation. The ability to learn about circuit architecture then depends on the outcome: a multipeaked DRPS is informative, but a nonmultipeaked DPRS is not (see Figure 3 and Figure 5 and above). For instance, direct perturbation of the amplitude of neural responses primarily through a change in the activation threshold of neurons (putatively the mechanism in Kanter et al., 2017 through the action of DREADDs (Sternson and Roth, 2014; SánchezRodríguez et al., 2017)) is predicted to result in amplitude changes but not a pattern period change in any of the candidate models (cf. Figure 5—figure supplement 3), and therefore cannot discriminate between feedforward and recurrent models unless the perturbation also affects gains or time constants. Computation of the DRPS would help resolve the question.
What if the perturbation targets only a (random) subset of grid cells in the circuit? The qualitative predictions are unchanged if only a (random) subset of the neurons in the grid cell circuit are targeted by the perturbation. Quantitatively, the size of the effect (change in period for a given perturbation strength) will be scaled in recurrent network models by an amount proportional to the fraction, Figure 5—figure supplement 1.
What if the perturbation affects the timeconstants of only excitatory (E) or only inhibitory (I) cells? In an implementation of recurrent grid cell models with separate E and I populations (Widloski and Fiete, 2014), if the perturbation is specific to only E or only I cells (so long as the I cells are the ones that mediate inhibitory interactions between grid cells), the qualitative predictions regarding the direction of change in the spatial tuning period and in the DRPS remain unchanged (Figure 5—figure supplement 2).
What if the perturbation affects the gain of only E or only I cells? The effect is the same as a timeconstant change in the particular populations (see response above). On the other hand, changing the bias rather than the slope of the E cells, I cells, or both is a nonspecific perturbation (see question above); the population period in recurrent models is predicted to not vary with a bias change (Figure 5—figure supplement 3).
What if the perturbed population includes speed cells (Kropff et al., 2015) and pure head direction cells? Changing the amplitude or timeconstant of speed cells and head direction cells could further change the spatial tuning periods of grid cells by affecting the gain of velocity inputs to the system. However, perturbation of these cells should not by itself induce a shift in relative phases in recurrent and feedforward network models, so DRPS predictions remain as if the perturbed population did not include these cells.
What if the perturbed population includes conjunctive gridhead direction cells (Sargolini et al., 2006)? The experiment would answer the question of whether grid activity patterns are generated through velocity integration by the targeted set of cells, without distinguishing the relative contributions of these two populations. Narrowing the set of targeted cells to only one of these populations would provide finergrained answers about mechanism. Similarly, if the experiment targets only grid cells and the mechanism is found to be feedfoward, the same experiments can be repeated one level upstream.
Would recurrent network models in which I cells are spatially untuned generate different predictions? The model for spatially tuning in E cells and untuned I responses is in Widloski and Fiete, 2014. In simulation, the population period in the model can move in the opposite direction as a function of perturbation compared to models with spatially tuned I cells. However, the DRPS is sensitive to the magnitude but not the sign of phase shifts hence predictions about the DRPS remain qualitatively the same in both types of models (Figure 5—figure supplement 4).
What if errors in spatial phase estimation exceed the uncertainty assumed in this paper? Most predictions (Figure 3) for differentiating between candidate mechanisms depend only on the DRPS envelope and its width and not on its detailed multipeaked structure. Thus, they are fairly robust to phase estimation uncertainty (left two columns in (Figure 4b). Better phase estimation is only required if the goal is to determine the number of bumps in an underlying recurrently patterned network. As described earlier, this uncertainty can be reduced by averaging over longer trajectories.
Will it be possible to discover if grid cell responses are based on selective feedforward summation of the intermingled nongrid cells in MEC? If the perturbation can be confined to exclude the intermingled nongrid cells in MEC, then it should be possible to use the perturbation and DRPS approach to tell apart a feedforward mechanism of spatial tuning inherited from nongrid MEC cells from a recurrent mechanism.
What if corrective inputs from place cells or external cues override the perturbationinduced change in the grid cell response? This is a real possibility in familiar environments, and can mask a multipeaked DRPS even if it would exist in the absence of corrective cues. For this reason, postperturbation experiments should be performed in novel, featureless environments (Figure 5—figure supplement 5; more discussion below).
Discussion
It is interesting to compare the potential of the present approach for discovering mechanism with other approaches. A highquality, fullcircuit connectome (Seung, 2009; Briggman et al., 2011) can specify the topology and locality of the network architecture. In other words, with appropriate analysis of the obtained data it should be possible to learn whether the connectivity matrix is ‘local’ (Widloski and Fiete, 2014) (Figure 1a), partially periodic (Figure 1b), or fully periodic (Figure 1c).
Network topology is, however, but one ingredient in circuit mechanism: Determining whether the observed connectivity actually accounts for the activity still requires inference (for instance, given a set of connections and weights, it is unknown whether they are strong enough to drive pattern formation in neural activity; determining this involves writing down a model of neural dynamics with the observed coupling). Even with further inference steps, whether the network originates certain functions like velocitytoposition integration and spatial tuning de novo (as in Figure 1a–c) or only amplifies or alters spatial tuning inherited from elsewhere (as in Figure 1f) cannot be answered by connectomics data. Despite their functional differences, feedforward and recurrent network models may exhibit similar lateral connectivity between grid cells. By contrast, the perturbative approach outlined here has the potential to reveal whether the function of path integration and spatial tuning originates in the perturbed set. The same approach can be sequentially applied to candidate areas progressively upstream of the grid cells.
Next, fullcircuit activity data at singleneuron resolution can reveal much about the dynamics and dimensionality of the population response in the circuit. But without perturbation, inferring mechanism from activity alone is problematic: Materials and methods to estimate connectivity from activity (Pillow et al., 2008; Roudi et al., 2009; Honey et al., 2009) yield only effective couplings that reflect collective and externally driven correlations in addition to the true couplings. In other words, activity data alone without perturbation does not indicate where the observed activity arises or its mechanisms.
In summary, while connectomics and largescale recording will provide vast amounts of valuable information, they are by themselves fundamentally correlative and thus not sufficient for discriminating between the candidate models discussed here. As we have shown, they may also not be immediately necessary: a lowdimensional or ‘global’ perturbative approach which does not require targeting specific individual neurons according to their responses can yield rich information about mechanism, and can do so with a far sparser dataset.
Interestingly, cooling and similar perturbation experiments have been performed in V1 (Michalski et al., 1993; Ferster and Miller, 2000) but were not as revealing about underlying mechanism as they promise to be in grid cells. Why is this? Unfortunately, the candidate models of orientation tuning in V1 are ring networks (fully periodic, singlebump) or a feedforward mechanism, and as we have seen, these two models do not differ in their predictions for the DRPS (Figure 3). The multibump spatial tuning of grid cells derived from velocity integration at some stage offers a way to distinguish feedforward from recurrent models because perturbation at the integration stage is predicted produce a change in the spatial tuning curve period, an opportunity that does not exist in in V1. Thus, grid cells offer a unique opportunity to uncover the circuit mechanisms that support tuning curves and computation in the cortex, and our modeling work shows how to do so.
Assumptions
We have assumed that the population pattern is stable against rotations (but the spatial tuning curves of cells are permitted to rotate) because a rotation would induce large changes in the DRPS and obscure the predicted effects of pattern expansion. Our assumption is supported by the observation that cell–cell phase relationships between grid cells are conserved across time and environments (Yoon et al., 2013), which can only hold if the underlying population pattern does not rotate.
The simplification that relative phases in the population pattern can be obtained from relative spatial tuning phases is valid if the intrinsically determined relative spatial phases of cells are not overridden by external spatial inputs. For instance, if an external cue (landmark or boundary) is associated with a specific configuration of grid cell phases, with the association acquired preperturbation, then the cue could activate the same configuration of grid cells postperturbation, which can interfere with the perturbationinduced shifts in the intrinsic relative phases between these cells. To avoid this possibility it is important, postperturbation, to assess spatial tuning relationships between cells only in novel environments, where there are no previously learned associations between external cues and the grid cell circuit. Ideally, these novel environments will be relatively free of spatial cues that resemble previously encountered cues and boundaries. Thus, the best environments for postperturbation testing would be circular 2D arenas, differently colored, patterned, and scented, and with minimal distal cues beyond a global orienting cue; or virtual environments with visually textured but landmarkfree walls (Yoon et al., 2016).
Even in novel environments, intrinsic error correcting mechanisms hypothesized in Sreenivasan and Fiete, 2011 might trigger preperturbation grid cell configurations: a configuration of grid cells, after it is associated with a specific place field, can be triggered simply by activation of that place cell by another but similar grid cell configuration in the novel environment. We explore this possibility in a model and show that even after constructing associations of grid configurations with place fields at every location in two familiar environments, grid cell activations in a novel environment do not trigger activation of the learned place fields and their associated grid configurations from the familiar environments Figure 5—figure supplement 5. Based on this result, we believe that postperturbation relative phases in grid cells may be relatively unaffected by intrinsic errorcorrection mechanisms in relatively featureless novel environments.
An interesting corollary to the possibility that previously learned reset or corrective inputs may coactivate cells that are outofphase cells postperturbation (as is possible for partially periodic and aperiodic recurrent mechanisms) in familiar environments is that such resets should degrade rather than improve the quality of grid cell spatial tuning postperturbation in previously learned environments.
Finally, it is important to note that if, in feedforward models, one were to include strong, continuous (rather than punctate, landmarkbased) feedback from the grid cell layer to the spatially tuned inputs (as in Bush and Burgess, 2014), the network would effectively become a recurrent circuit that we have not included as a candidate. Similarly, we have excluded from our analysis recurrent network models of the spatial circuit with heterogeneous tuning and connectivity (Cueva and Wei, 2018; Banino et al., 2018; Kanitscheider and Fiete, 2016); these models do not yet capture the modular dynamics of the grid cell system, in which cells cluster in spatial period and those with similar period have the same orientation without the help of external aligning cues. When these models are refined, and if the result is a distinct mechanism for modular grid cell dynamics than the candidate models considered here, it will be interesting to perform our proposed perturbations in them to obtain their predictions for experiment.
Materials and methods
Figure 1 and Figure 1—figure supplement 1 are schematic. In Figure 2, Figure 4a–b, Figure 2—figure supplement 2 and Figure 4—figure supplement 1, relative phase is computed from the population phases using idealized (handdrawn) periodic population patterns that expand (${\delta}^{ij}={\varphi}_{pop}^{i}{\varphi}_{pop}^{j}$), without the use of neural network simulations. Figure 3, Figure 1—figure supplement 1, Figure 2—figure supplement 1, Figure 3—figure supplement 1, Figure 2—figure supplement 4, Figure 5—figure supplement 1, Figure 5—figure supplement 2, Figure 5—figure supplement 3, and Figure 5—figure supplement 4, which distinguish between different recurrent architectures, are obtained by simulating the grid cell system in a neural network. Briefly, the network consists of excitatory and inhibitory neurons (except in Figure 1—figure supplement 1 – see figure caption for details) with linearnonlinear Poisson (LNP) spiking dynamics (Burak and Fiete, 2009; Widloski and Fiete, 2014) (except for Figure 2—figure supplement 4, where we use HodgkinHuxley dynamics). Structured lateral interactions between neurons pattern the neural population responses. Relative spatial tuning phases are computed from the tuning curves of different neurons, obtained by simulating the network response over 1 min long simulated quasirandom trajectories. The analysis of relative phase shifts, tuning amplitude and period in a network includes all cells with sufficiently good spatial tuning profiles: this set includes all cells in the fully and partially periodic networks and $3/4$ of the cells in aperiodic networks (from the central part of the network). Since the inhibitory and excitatory populations share similar population patterning and spatial tuning in these simulations (except in Figure 5—figure supplement 4), we arbitrarily display results from the inhibitory population.
Neural network simulations
We use two different neuron models in our network simulations: LNP and HodgkinHuxley neurons, described below. Roman subscripts (e.g. $i,j$) refer to individual cells within population $P$. The population index $P$ can take the values $\{$I, E${}^{R}$, E${}^{L}\}$, designating inhibitory cells or excitatory cells that receive rightward or leftward velocity input, respectively. Integration is by the Euler method with timestep $dt$.
Linearnonlinearpoisson (LNP) neurons
Request a detailed protocolThe timevarying firing rate ${r}_{i}^{P}(t)=f({G}_{i}^{P}(t))$ of the $(P,i)$th cell is an instantaneous function of its timevarying summed input ${G}_{i}^{P}(t)$ with thresholdlinear transfer function $f$:
Neurons emit spikes according to an inhomogeneous point process with rate ${r}_{i}^{P}(t)$ and coefficient of variance of CV = 0.5 (see Burak and Fiete, 2009 and Widloski and Fiete, 2014 for details on generating a subPoisson point process). LNP dynamics were used in all simulations except for Figure 2—figure supplement 4.
Cortical HodgkinHuxley (CHH) neurons
Request a detailed protocolThe membrane potential of the $(P,i)$th neuron is given by:
where ${C}_{m}$ is the capacitance of the membrane, ${I}^{ion}(V)$ is the sum of the cell’s intrinsic ionic currents, and ${I}^{syn}(V)$ is the current from recurrent and feedforward synaptic inputs to the cell. The ionic current is modeled as (Pospischil et al., 2008):
where the $\overline{g}$’s represent maximal conductance values and the $\overline{V}$’s are the reversal potentials of the leak conductance (L), the fast (K) and slow (M) potassium conductances, and the sodium conductance (Na). The dynamics and parameter settings of $n,m,q,h$ are as in Pospischil et al., 2008 (we have replaced the ‘p’ gating variable in Pospischil et al., 2008 with the notation ‘q’). For CHH neurons, the time of a spike is defined as the timestep when the voltage crosses 0 mV from below. CHH dynamics were used in Figure 2—figure supplement 4.
Synaptic activation
Request a detailed protocolFor both LNP and CHH neurons, spikes by the $(P,i)$th neuron activate all its outgoing synapses according to:
where ${t}_{i,b}^{P}$ is the time of the $b$th spike and $\delta (t)$ is the Dirac delta function. The sum is over all spikes of the cell.
Network inputs and interactions
Request a detailed protocolWe based our grid cell network models on the connectivity and weights that emerge from plasticity rules over a plausible developmental process, given in Widloski and Fiete, 2014, and thus might better represent the grid cell system than a model fully wired by hand. Moreover, the network contains both inhibitory and excitatory units (with the number of inhibitory units equalling 1/5 the number of excitatory units, like in cortex).
Synaptic input to LNP cells
The total synaptic input ${G}_{i}^{syn,P}(t)$ into the $(P,i)$th LNP cell is given by
where ${\alpha}^{vel,P}$ is the velocity input (described below), in multiplicative form; ${G}_{i}^{rec,P}={\sum}_{{P}^{\prime}}{\sum}_{j=1}^{{N}^{{P}^{\prime}}}{W}_{ij}^{P{P}^{\prime}}{s}_{j}^{{P}^{\prime}}$ is the recurrent network input; ${G}^{0},{G}^{{0}^{\prime},P}$ are (small, positive) constant bias terms (${G}^{0}=50,{G}^{{0}^{\prime},I}=0,{G}^{{0}^{\prime},{E}^{L}}={G}^{{0}^{\prime},{E}^{R}}=15$); and ${A}_{i}^{P}$ is a smooth envelope that modulates neural activity magnitudes across the network (described below).
To model additive velocity input, as in Figure 2—figure supplement 1c, we replace Equation 5 with the following:
where ${G}^{vel,P}={W}^{vel}{\alpha}^{vel,P}$ and ${W}^{vel}=200$ (${\alpha}^{vel,P}$ described below).
Synaptic input to CHH cells
The total synaptic current ${I}_{i}^{syn,P}$ into the $(P,i)$th CHH neuron is given by
where ${\overline{V}}^{{P}^{\prime}}$ is the reversal potential for synaptic inputs from population ${P}^{\prime}$ (${\overline{V}}^{E}=0$ mV and ${\overline{V}}^{I}=80$ mV), ${g}_{i}^{rec,P}={\sum}_{j}^{{N}^{{P}^{\prime}}}{W}_{ij}^{P{P}^{\prime}}{s}_{j}^{{P}^{\prime}}$ is the recurrent network input, ${I}^{0}$ is a constant bias and ${\alpha}^{vel,P},{A}_{i}^{P}$ are the same velocity and envelope terms mentioned above.
Velocity input
The cells in the $P$th population receive a common motionrelated input proportional to animal velocity along preferred direction ${\widehat{e}}^{P}$:
where $\overrightarrow{v}$ is the instantaneous velocity of the animal and ${\beta}^{vel}$ is a scalar gain parameter. ${\widehat{e}}^{P}=$ (0,0), (0,1), (0,–1) for the I, E${}^{R}$ and E${}^{L}$ populations, respectively. Unless otherwise noted, the velocity input is derived from a 1 min quasirandom trajectory (Widloski and Fiete, 2014).
Recurrent weights
We based the recurrent weights ${W}_{ij}^{P{P}^{\prime}}$ from cell $j$ in population ${P}^{\prime}$ to $i$ in $P$ on those from the mature network of Widloski and Fiete, 2014.
We first describe weights in their periodic form. Let ${x}_{ij}=i\gamma j$, where $\gamma =\frac{{N}_{P}}{{N}_{{P}^{\prime}}}$ (${N}_{P}$ is the number of neurons in population $P$). We also define the norm, ${x}_{{N}_{P}}\equiv \mathrm{min}({N}_{P}x,x)$. The E$\to $I weights (i.e., $P=I$ and ${P}^{\prime}={E}^{L},{E}^{R}$) are written as
where $\eta $ controls the overall weight strength, $\mathrm{\Delta}$ and $\sigma $ control the shift and width, respectively, of the Gaussian profile, and $\rho $ is a scale factor that is used to shift from partially periodic ($\rho =1$) to fully periodic ($\rho =11$). The parameter $\rho $ takes the same values for the IE and II weights (which are described below). The parameters are set as follows:
Weight  $\eta $  $\mathrm{\Delta}$  $\sigma $ 

E${}^{L}\mathrm{\to}$ I  11.5  2  4 
E${}^{R}\mathrm{\to}$ I  11.5  2  4 
The I$\to $E weights are written as
where $\mathrm{\Theta}$ is the Heaviside function ($\mathrm{\Theta}(x)=1$ for $x\ge 0$ and 0 otherwise). The first Heaviside function cuts out weights along the diagonal (the width of which is controlled by the parameter $\delta $), while the second, third, and fourth Heaviside functions together act as a windowing function to set to zero portions of the matrix to make the weights qualitatively resemble the developmental weights from Widloski and Fiete, 2014. The parameters are set as follows:
Weight  $\eta $  $\mathrm{\Delta}$  $\sigma $  $\mu $  $\delta $ 

I $\mathrm{\to}$ E${}^{L}$  4  8  10  1  3 
I $\mathrm{\to}$ E${}^{R}$  4  8  10  1  3 
Finally, the I$\to $I weights are written as
which is essentially a sumofGaussians with the central portion removed (the width of which is controlled by $\delta $). The parameters are set as follows:
Weight  $\eta $  $\mathrm{\Delta}$  $\sigma $  $\delta $ 

I $\mathrm{\to}$ I  12  4  6  3 
For the aperiodic network, the weights have the same form and parameter values as above (with $\rho =1$), except with the following replacements:
where $.$ is the absolute value and ${A}_{ij}^{P{P}^{\prime}}={A}_{i}^{P}{A}_{j}^{{P}^{\prime}}$ is an envelope function used to enforce a tapered profile on the weights, similar to Burak and Fiete, 2009:
where ${r}_{i}^{P}=i\frac{{N}_{P}}{2}$, and$\kappa =0.3$ determines the range of the taper while ${a}_{0}=30$ controls its steepness.
Changing pattern period by varying the ‘neural’ timeconstant and the gain of recurrent inhibition in a network of LNP neurons
The period of the population pattern can be varied by rescaling the synaptic activation time constant, ${\tau}_{syn}$. It can also be varied by changing a gain parameter ${\gamma}_{inh}$ that controls the strength of synaptic weights from the inhibitory neurons: we set ${W}^{PI}\to {\gamma}_{inh}{W}^{PI}$, and allow ${\gamma}_{inh}$ to be varied away from unity.
The effect of timeconstant on period in the different networks is quite nontrivial: It cannot be derived from a linear stability analysis on the network equations since it depends strongly on nonlinear interactions within the network bulk and with the network boundaries (Widloski, 2015). Instead, we study the effect though simulation of the nonlinear dynamics of the networks.
As noted in the main manuscript, neuromodulators can drive the requisite gain changes in recurrent weights. We show, through the more detailed HodgkinHuxley neuron simulations described below, that temperature may be used in experiments to cause similar changes in period as can be affected by changing recurrent weight strength, and that the effects of temperature change resemble the effects of changing the timeconstant in the LNP model.
We study HodgkinHuxley (HH) dynamics to predict, with the help of more biophysically detailed neuron models and the documented variation of their parameters with temperature, the effects of cooling on population activity in grid cells. Specifically, we use a ‘regular spiking’ HH model of cortical neurons (Pospischil et al., 2008), which we supplement with models that describe temperatureinduced changes in the parameters (Hodgkin et al., 1952; Katz and Miledi, 1965).
Effects of temperature and neuromodulation on HH dynamics
Some HH models include modifications that capture the effects of temperature variation (Hodgkin et al., 1952; Katz and Miledi, 1965). These temperature effects are modeled by ${Q}_{10}$ factors that multiply the timeconstants (${Q}_{10}^{\tau}=$ 3) and amplitudes (${Q}_{10}^{a}=$ 1.3) of the ionic conductances. At temperature $T$ (in ${}^{\circ}$C), the conductance amplitudes $\overline{g}(T)$ and time constants $\tau (T)$ have the following form:
where ${T}_{0}$ is 36${}^{\circ}$. We applied the ${Q}_{10}$ factor for $\overline{g}$ to the ionic conductance amplitudes ${\overline{g}}_{L},{\overline{g}}_{K},{\overline{g}}_{M},{\overline{g}}_{Na}$ as well as to the synaptic conductance amplitudes ${W}_{ij}^{P{P}^{\prime}}$. We also simultaneously applied the ${Q}_{10}$ factor for $\tau $ to the conductance and synaptic timeconstants ${\tau}_{n},{\tau}_{q},{\tau}_{m},{\tau}_{h}$ and ${\tau}_{syn}$. (For gating variable $x$, the time constant ${\tau}_{x}$ is defined as ${\tau}_{x}=1/({\alpha}_{x}+{\beta}_{x})$, where ${\alpha}_{x}$ and ${\beta}_{x}$ are the rate constants governing the gating variable’s dynamics (Pospischil et al., 2008).)
Finally, to isolate which parameters drove the strongest thermal effects on population patterning and the direction of these effects (so that we could extract lessons for how to vary parameters in grid cell models with simpler neuron dynamics) we applied thermal changes to the ionic conductances only (changing ${\overline{g}}_{L},{\overline{g}}_{K},{\overline{g}}_{M},{\overline{g}}_{Na},{\tau}_{n},{\tau}_{q},{\tau}_{m},{\tau}_{h}$ according to the ${Q}_{10}$ factors while holding ${W}_{ij}^{P{P}^{\prime}}$ and ${\tau}_{syn}$ constant), or to the synaptic conductances only (changing ${W}_{ij}^{P{P}^{\prime}}$ and ${\tau}_{syn}$ according to the ${Q}_{10}$ factors while holding the ionic conductance parameters fixed).
To simulate the effects of a neuromodulatory gain change in inhibitory synapses, we set ${W}^{PI}$ to ${\gamma}_{inh}{W}^{PI}$, where ${\gamma}_{inh}$ is the prefactor modulating the strength of inhibition.
Simulation parameters
Request a detailed protocolLNP dynamics
${N}_{{E}_{L}}={N}_{{E}_{R}}=$ 400 neurons; ${N}_{I}=$ 160 neurons; CV = 0.5; $dt=$0.5 ms; ${\tau}_{syn}=$30 ms*; ${\beta}^{vel}=$ 1; ${\gamma}_{inh}$ = 1*. (*: Indicates that parameters can change through perturbation.)
Aperiodic network with CHH dynamics
All ionic conductance parameters are identical to those described in Pospischil et al., 2008 for the RS model; as noted there, the parameters are set to values corresponding to a temperature of ${T}_{0}={36}^{\circ}C$. Synaptic weight definitions and parameter values same as LNP dynamics for aperiodic network (above), except that all $\eta $ values are scaled by the factor $0.0015$ 400 neurons; ${N}_{I}=$ 160 neurons; $dt=$0.025 ms; ${\tau}_{syn}=$15 ms*; ${\beta}^{vel}=$ 0.8; ${C}_{m}$ = 1 $\mu $ F/cm$}^{2$; ${\overline{g}}_{L}=$0.1 ms/cm$}^{2$*; ${\overline{g}}_{K}=$5 ms/cm${}^{2}$*; ${\overline{g}}_{M}=$0.07 ms/cm${}^{2}$*; ${\overline{g}}_{Na}=$50 ms/cm${}^{2}$*; ${\overline{V}}_{L}=$−70 mV; ${\overline{V}}_{K}=$−90 mV; ${\overline{V}}_{Na}=$50 mV; ${I}^{0}=$ 3 $\mu $A/cm$}^{2$; $\rho =$1. ${\gamma}_{inh}$ = 1*. (*: Indicates that parameters can change through perturbation).
LNP dynamics with EE connections
All network parameters and synaptic weight definitions same as for LNP network (see above) with the addition of EE connections (see below), except that the weight parameters have the following changes:
Weight  $\eta $  $\mathrm{\Delta}$  $\sigma $  $\mu $  $\delta $ 

E${}^{L}\mathrm{\to}$ I  3  2  8  
E${}^{R}\mathrm{\to}$ I  3  2  8  
I $\mathrm{\to}$ E${}^{L}$  3.25  8  8  1  3 
I $\mathrm{\to}$ E${}^{R}$  3.25  8  8  1  3 
I $\mathrm{\to}$ I  4  4  6  3 
The EE weights for the periodic networks are written similar to the EI weights as
and have the following parameters:
Weight  $\eta $  $\mathrm{\Delta}$  $\sigma $ 

E${}^{L}\mathrm{\to}$ E${}^{L}$  5.5  4  4 
E${}^{R}\mathrm{\to}$ E${}^{R}$  5.5  4  4 
E${}^{L}\mathrm{\to}$ E${}^{R}$  5.5  0  4 
E${}^{R}\mathrm{\to}$ E${}^{L}$  5.5  0  4 
As in the LNP case, to get the aperiodic version of the EE weights, replace
where the envelope function ${A}_{ij}$ is described above.
Alternative formulation of the DRPS
Request a detailed protocolAs before, the $i$th cell’s firing phase within the periodic population activity pattern, defined as the cell’s population phase, is ${\varphi}_{pop}^{i}=((i1)mod{\lambda}_{pop})/{\lambda}_{pop}$ (with the arbitrary choice, made without loss of generality, that neuron 1 has phase 0) (Figure 2—figure supplement 3b, blue curve). For each cell in the population, plotting the preperturbation phase against the postperturbation phase (red vs. blue curves in Figure 2—figure supplement 3b) shows that the data is quantized and lies on a series of parallel manifolds, Figure 2—figure supplement 3c. This quantization is captured via the following transformation to the phase shifts:
(we have assumed that the true stretch factor, $\alpha $, is known – later, we will show how $\alpha $ can be inferred from the data) followed by a modulo operation
and then reflecting about the midpoint of the interval
The distribution of these phase shift values, Figure 2—figure supplement 3d, has three special properties: (1) The distribution is quantized, due to the fact that population activity pattern itself is quantized. (2) The number of peaks in the distribution is exactly equal to the number of bumps in the population activity pattern (this holds only for sufficiently small perturbations, such that $\mathrm{\Delta}M<0.5$, where $M$ is the number of bumps in the preperturbation population activity pattern – see Figure 2—figure supplement 2 for explanation). (3) The peak separation in the distribution is exactly equal to the stretch factor, $\alpha $. The transformations described in Equations 2022 require knowledge of the stretch factor, $\alpha $, a quantity that is not directly observable. However, it can be inferred from the data, because the desired $\alpha $ value is the one that makes the distribution the most peaky. This is equivalent to projecting the data onto its orthogonal axis, Figure 2—figure supplement 3c. Peakyness is quantified as the Pearson’s correlation coefficient between the DRPS and a comblike function defined over the same interval. The comb function is a series of deltafunctions laid out with a spacing equal to $\alpha $. The desired $\alpha $ stretch factor is the one that maximizes this correlation (not shown).
Correction of grid cell phases by grid celldriven place cells
Request a detailed protocolThe model described below is based on work in Sreenivasan and Fiete, 2011. We assume $M$ modules, each with $N$ grid cells. The $i$th grid cell from the ${m}^{th}$ module in the ${k}^{th}$ environment has the following simplified tuning curve response:
where ${\lambda}_{m}$ is the spatial period of the ${m}^{th}$ module, ${\varphi}_{i}=2\pi i/N$ is the cell’s phase relative to others within the module (fixed across environments), and ${\stackrel{~}{\varphi}}_{m,k}$ is a random modulewide phase shift that is specific to each module and each environment, Figure 5—figure supplement 5a. The synaptic projections from grid cells to place cells are set as follows: Assume a population of $P$ place cells. For the $i$th place cell in the ${k}^{th}$ familiar environment, assign a random place preference, ${x}_{i,k}^{PC}$. The synaptic weight from the ${(j,m)}^{th}$ grid cell onto the ${(i)}^{th}$ place cell is incremented based on experience in environment $k$. The increment is Hebbian, given by the amplitude of the grid cell tuning curve at that place cell’s preferred location:
The total weight from grid cell $(j,m)$ to place cell $i$ is given by the sum of increments over all $L$ familiar environments:
Given these weights, the ${i}^{th}$ place cell’s full subthreshold activity in environment $k$ is simply a weighted sum over the activities of the grid cells across modules, based on its weights:
This description of place cell subthreshold activations holds for both familiar and novel environments; the only difference between familiar and novel environments is that in the latter there has been no increment of the grid cellplace cell weights based on coincident grid cellplace cell activity, Figure 5—figure supplement 5a–b. In the current implementation, we allowed every cell to have a field in every familiar environment. We see that even in this case, the subthreshold activations of PCs in the novel environment are far lower than at place fields in familiar environment; in other words, they will not be activated and drive correction or resetting of the grid cell phases in the novel environment. Including the measured degrees of sparseness in PCs should lead to even less interference than seen in simulated novel environment conditions.
Measures used in main text
DRPS in 1D
Request a detailed protocolThe relative phase of cell $i$ and $j$ is defined as ${\delta}^{ij}={d}^{ij}mod\lambda /\lambda $, where ${d}^{ij}$ is the offset in the central peak of the crosscorrelation in their spatial tuning curves, and $\lambda $ is their common spatial period (in the main text, for Figures 2, 4ab, Figure 2—figure supplement 2, and Figure 4—figure supplement 1, the relative phase is computed directly from the population phases, that is, ${\delta}^{ij}={\varphi}_{pop}^{i}{\varphi}_{pop}^{j}$. The relative phase magnitude is given by $\delta =\mathrm{min}(\delta ,1\delta )$, where $\cdot $ is the absolute value norm. The DRPS is computed by making a distribution of phase magnitude shifts, ${\delta}_{pre}{\delta}_{post}$, where ${\delta}_{pre}$ and ${\delta}_{post}$ are the relative phases measured pre and postperturbation.
2D relative phase
Request a detailed protocolFor two cells $i$ and $j$, let $\overrightarrow{d}$ be the displacement vector which measures the 2D offset in the central peak of the crosscorrelation in their spatial tuning curves. The displacement vector is converted into a 2D phase $\overrightarrow{\delta}$ according to $\overrightarrow{\delta}=({\delta}_{1},{\delta}_{2})=f({d}_{1}^{proj}/{\lambda}_{1}mod1,{d}_{2}^{proj}/{\lambda}_{2}mod1)$, where ${\overrightarrow{d}}^{proj}=({d}_{1}^{proj},{d}_{2}^{proj})$ is the oblique projection of $\overrightarrow{d}$ onto the principal vectors ${\lambda}_{1}{\widehat{e}}_{1}$ and ${\lambda}_{2}{\widehat{e}}_{2}$, and
DRPS in 2D
Request a detailed protocolThe DRPS in 2D is computed separately for the two components of the 2D relative phase. That is, given the relative phase vector $\overrightarrow{\delta}=({\delta}_{1},{\delta}_{2})$, the DRPS is computed by making a distribution of phase magnitude shifts for each component: ${\delta}_{1,pre}{\delta}_{1,post}$ and ${\delta}_{2,pre}{\delta}_{2,post}$, where the magnitude is defined as the absolute value norm: $\cdot =\cdot .$
Bootstrap resampling and phase uncertainty
Request a detailed protocolGiven an original spike map of $M$ total spikes (with locations) from one cell, we created a new spike map of $N$ ($N<M$) total spikes, by picking spikes (with their corresponding location coordinates) from the original map one at a time, at random, and with replacement. The same was done for a second, simultaneously recorded cell. From these sampled spike trains for a pair of cells, we estimated relative phase (by computing the location of the peak closest to the origin in the crosscorrelation of the spatial maps of the two cells, as in Yoon et al., 2013). The procedure was performed 100 times, generating 100 bootstrapped relative phase estimates per cell pair. Phase uncertainty was measured as the peak location of the Rayleigh distribution that best fit the distribution of magnitudes of the bootstrapped relative phase estimates.
Spatial tuning curves
Request a detailed protocolFor a given cell and trajectory, we build a histogram of spike counts at each location (bin size = 1 cm), then normalize the count in each bin by the amount of time spent in it. The normalized histogram is smoothed by convolution with a boxcar filter (width = 5 bins) to yield a spatial tuning curve.
Spatial tuning period and amplitude
Request a detailed protocolThe spatial tuning period is measured as the inverse of the spatial frequency with the highest peak in the power spectrum of the spatial tuning curve (excluding the peak at 0 frequency). Likewise, the spatial tuning amplitude is measured as the mean spike rate density across the bins of the spatial tuning curve. The quantities reported in Figure 3 and Figure 3—figure supplement 1 are averaged over all cells in the population.
Population activity period and gridness
Request a detailed protocolThe population activity gridness is taken to be the power of the largest frequency component of the power spectrum measured from a normalized snapshot (frame) of the population activity (normalized = mean subtracted, followed by division by standard deviation). The power spectrum is rescaled by the factor 2/L${}^{2}$, where L is the number of bins in the population activity vector from which the power spectrum was computed. The population activity vector is shortened to include only the middle onehalf of the population, so that for the E${}^{L}$ population, L is 100. From the power spectrum, the population activity period is taken to be the wavelength at which the power spectrum has the largest peak. Throughout the paper, both the reported population period and gridness are averaged over the last 10000 snapshots of the population activity pattern from a given trial.
Velocity response
Request a detailed protocolVelocity response is measured as the translation speed (neurons/sec) of the network pattern to fixed input velocity, computed by tracking the displacement of the pattern for 10 s, smoothing the resulting trajectory with an 4 s moving average filter, and then measuring the average speed of the middlehalf of the trajectory.
Periodicity score for the DRPS
Request a detailed protocolWe smooth the histogram of relative phase shifts (by convolution with a 2bin Gaussian kernel) and normalize it (by mean subtraction and division by the standard deviation). Next, we compute the power spectrum, rescaling the result by $2/{L}^{2}$, where $L$ is the number of bins in the histogram ($L=200$). The periodicity score is set to be the power of the largest amplitude nonzero frequency component in the scaled power spectrum. This score returns 1 if the DRPS is a pure sinusoid. It returns 0 if the DRPS is flat and returns an average value of $<0.2$ if the DRPS were constructed bin by bin by taking independent, identically distributed (iid) samples from a uniform distribution on the unit interval.
References

Grid cells require excitatory drive from the hippocampusNature Neuroscience 16:309–317.https://doi.org/10.1038/nn.3311

An isomorphic mapping hypothesis of the grid representationPhilosophical Transactions of the Royal Society B: Biological Sciences 369:20120521.https://doi.org/10.1098/rstb.2012.0521

Do we understand the emergent dynamics of grid cell activity?Journal of Neuroscience 26:9352–9354.https://doi.org/10.1523/JNEUROSCI.285706.2006

Accurate path integration in continuous attractor network models of grid cellsPLoS Computational Biology 5:e1000291.https://doi.org/10.1371/journal.pcbi.1000291

An oscillatory interference model of grid cell firingHippocampus 17:801–812.https://doi.org/10.1002/hipo.20327

A hybrid oscillatory interference/continuous attractor network model of grid cell firingJournal of Neuroscience 34:5065–5079.https://doi.org/10.1523/JNEUROSCI.401713.2014

Functional imaging of hippocampal place cells at cellular resolution during virtual navigationNature Neuroscience 13:1433–1440.https://doi.org/10.1038/nn.2648

Neural mechanisms of orientation selectivity in the visual cortexAnnual Review of Neuroscience 23:441–471.https://doi.org/10.1146/annurev.neuro.23.1.441

A spin glass model of path integration in rat medial entorhinal cortexJournal of Neuroscience 26:4266–4276.https://doi.org/10.1523/JNEUROSCI.435305.2006

A model of grid cells based on a twisted torus topologyInternational Journal of Neural Systems 17:231–240.https://doi.org/10.1142/S0129065707001093

A model combining oscillations and attractor dynamics for generation of grid cell firingFrontiers in Neural Circuits 6:30.https://doi.org/10.3389/fncir.2012.00030

Measurement of currentvoltage relations in the membrane of the giant axon of LoligoThe Journal of Physiology 116:424–448.https://doi.org/10.1113/jphysiol.1952.sp004716

Receptive fields of single neurones in the cat's striate cortexThe Journal of Physiology 148:574–591.https://doi.org/10.1113/jphysiol.1959.sp006308

The effect of temperature on the synaptic delay at the neuromuscular junctionThe Journal of Physiology 181:656–670.https://doi.org/10.1113/jphysiol.1965.sp007790

Path integration and the neural basis of the 'cognitive map'Nature Reviews Neuroscience 7:663–678.https://doi.org/10.1038/nrn1932

The effect of reversible cooling of cat's primary visual cortex on the responses of area 21a neuronsThe Journal of Physiology 466:133–156.

Conserved spatial learning in cooled rats in spite of slowing of dentate field potentialsThe Journal of Neuroscience 14:4458–4466.https://doi.org/10.1523/JNEUROSCI.140704458.1994

Grid cells and cortical representationNature Reviews Neuroscience 15:466–481.https://doi.org/10.1038/nrn3766

Minimal HodgkinHuxley type models for different classes of cortical and thalamic neuronsBiological Cybernetics 99:427–441.https://doi.org/10.1007/s0042200802638

Analysis of GABAA receptor function and dissection of the pharmacology of benzodiazepines and general anesthetics through mouse geneticsAnnual review of pharmacology and toxicology 44:475–498.https://doi.org/10.1146/annurev.pharmtox.44.101802.121429

Cellular mechanisms of spatial navigation in the medial entorhinal cortexNature Neuroscience 16:325–331.https://doi.org/10.1038/nn.3340

New perspectives on the mechanisms for orientation selectivityCurrent Opinion in Neurobiology 7:514–522.https://doi.org/10.1016/S09594388(97)800311

Efficient "Shotgun" inference of neural connectivity from highly subsampled activity dataPLoS Computational Biology 11:e1004464.https://doi.org/10.1371/journal.pcbi.1004464

The hippocampus as a predictive mapNature Neuroscience 20:1643–1653.https://doi.org/10.1038/nn.4650

Chemogenetic tools to interrogate brain functionsAnnual Review of Neuroscience 37:387–407.https://doi.org/10.1146/annurevneuro071013014048

ThesisGrid cell attractor networks: development and implications. PhD thesisUniversity of Texas, Austin.

Specific evidence of lowdimensional continuous attractor dynamics in grid cellsNature Neuroscience 16:1077–1084.https://doi.org/10.1038/nn.3450

Longterm dynamics of CA1 hippocampal place codesNature Neuroscience 16:264–266.https://doi.org/10.1038/nn.3329
Article and author information
Author details
Funding
Human Frontier Science Program (HFSPRGP0062/2014)
 Ila R Fiete
National Science Foundation (NSFCRCNS IIS1311213)
 Ila R Fiete
Howard Hughes Medical Institute
 Ila R Fiete
Simons Foundation
 Ila R Fiete
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We are grateful to Dori Derdikman, Caswell Barry, Laura Colgin, and Lisa Giocomo for useful discussions and Rishidev Chaudhuri and Ingmar Kanitscheider for comments on the manuscript. This work was supported in part by grants from the Human Frontiers in Science Program (HFSPRGP0062/2014), the National Science Foundation (NSFCRCNS IIS1311213), the Howard Hughes Medical Institute through the Faculty Scholars Program, and the Simons Foundation through the Simons Collaboration on the Global Brain.
Version history
 Received: November 11, 2017
 Accepted: July 7, 2018
 Accepted Manuscript published: July 9, 2018 (version 1)
 Version of Record published: August 6, 2018 (version 2)
Copyright
© 2018, Widloski et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,052
 views

 335
 downloads

 14
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Immunology and Inflammation
 Neuroscience
CD4^{+}CD25^{+}Foxp3^{+} regulatory T cells (Treg) have been implicated in pain modulation in various inflammatory conditions. However, whether Treg cells hamper pain at steady state and by which mechanism is still unclear. From a metaanalysis of the transcriptomes of murine Treg and conventional T cells (Tconv), we observe that the proenkephalin gene (Penk), encoding the precursor of analgesic opioid peptides, ranks among the top 25 genes most enriched in Treg cells. We then present various evidence suggesting that Penk is regulated in part by members of the Tumor Necrosis Factor Receptor (TNFR) family and the transcription factor Basic leucine zipper transcription faatflike (BATF). Using mice in which the promoter activity of Penk can be tracked with a fluorescent reporter, we also show that Penk expression is mostly detected in Treg and activated Tconv in noninflammatory conditions in the colon and skin. Functionally, Treg cells proficient or deficient for Penk suppress equally well the proliferation of effector T cells in vitro and autoimmune colitis in vivo. In contrast, inducible ablation of Penk in Treg leads to heat hyperalgesia in both male and female mice. Overall, our results indicate that Treg might play a key role at modulating basal somatic sensitivity in mice through the production of analgesic opioid peptides.

 Neuroscience
Biological synaptic transmission is unreliable, and this unreliability likely degrades neural circuit performance. While there are biophysical mechanisms that can increase reliability, for instance by increasing vesicle release probability, these mechanisms cost energy. We examined four such mechanisms along with the associated scaling of the energetic costs. We then embedded these energetic costs for reliability in artificial neural networks (ANNs) with trainable stochastic synapses, and trained these networks on standard image classification tasks. The resulting networks revealed a tradeoff between circuit performance and the energetic cost of synaptic reliability. Additionally, the optimised networks exhibited two testable predictions consistent with preexisting experimental data. Specifically, synapses with lower variability tended to have (1) higher input firing rates and (2) lower learning rates. Surprisingly, these predictions also arise when synapse statistics are inferred through Bayesian inference. Indeed, we were able to find a formal, theoretical link between the performancereliability cost tradeoff and Bayesian inference. This connection suggests two incompatible possibilities: evolution may have chanced upon a scheme for implementing Bayesian inference by optimising energy efficiency, or alternatively, energyefficient synapses may display signatures of Bayesian inference without actually using Bayes to reason about uncertainty.