RESEARCH ARTICLE

Replay as wavefronts and theta
sequences as bump oscillations in a gridcell attractor network
Louis Kang1,2*, Michael R DeWeese1,2
1

Redwood Center for Theoretical Neuroscience, Helen Wills Neuroscience Institute,
University of California, Berkeley, Berkeley, United States; 2Department of Physics,
University of California, Berkeley, Berkeley, United States

Abstract Grid cells fire in sequences that represent rapid trajectories in space. Duringlocomotion, theta sequences encode sweeps in position starting slightly behind the animal andending ahead of it. During quiescence and slow wave sleep, bouts of synchronized activity
represent long trajectories called replays, which are well-established in place cells and have been
recently reported in grid cells. Theta sequences and replay are hypothesized to facilitate many
cognitive functions, but their underlying mechanisms are unknown. One mechanism proposed forgrid cell formation is the continuous attractor network. We demonstrate that this establishedarchitecture naturally produces theta sequences and replay as distinct consequences of modulatingexternal input. Driving inhibitory interneurons at the theta frequency causes attractor bumps tooscillate in speed and size, which gives rise to theta sequences and phase precession, respectively.
Decreasing input drive to all neurons produces traveling wavefronts of activity that are decoded as
replays.

*For correspondence:louis.kang@berkeley.edu
Competing interests: The
authors declare that nocompeting interests exist.
Funding: See page 20
Received: 01 March 2019
Accepted: 15 November 2019
Published: 18 November 2019
Reviewing editor: Sachin
Deshmukh, Indian Institute of
Science Bangalore, India
Copyright Kang and
DeWeese. This article isdistributed under the terms of
the Creative Commons
Attribution License, which
permits unrestricted use andredistribution provided that theoriginal author and source arecredited.

Introduction
The hippocampal region contains spatially tuned cells that generally encode an animal’s position inits spatial environment. Place cells in the hippocampal formation fire at one or a few locations(O’Keefe and Dostrovsky, 1971), and grid cells in the medial entorhinal cortex (MEC) fire at many
locations that form a triangular lattice in space (Hafting et al., 2005). However, at short timescales,
these neurons fire in coordinated sequences that represent trajectories away from the current animal
position.
The first type of sequence occurs when the animal is moving and the local field potential (LFP) of
the hippocampal region is dominated by an oscillation of 5–11 Hz, the theta range (Vanderwolf, 1969; O’Keefe, 1976; Vinogradova, 1995). During every cycle of this oscillation, neurons corresponding to locations slightly behind the animal fire first, followed by those corresponding to thecurrent location and finally locations ahead of the animal (Skaggs et al., 1996; O’Neill et al., 2017).
These so-called theta sequences involving many neurons are related to a single-neuron phenomenon
called phase precession (O’Keefe and Recce, 1993; Hafting et al., 2008). When an animal firstenters the firing field of a place or grid cell, the neuron spikes late in the theta cycle. As the animal
moves through the field, subsequent spikes tend to arrive at smaller theta phase, or earlier within acycle. Thus, activity within each theta cycle starts with neurons whose peak firing occurs behind theanimal and ends with neurons whose peak firing occurs ahead of the animal, which is a theta
sequence (Skaggs et al., 1996).
The second type of sequence occurs when the animal is idle and not moving, a state called quiet
wakefulness or quiescence; it also occurs during slow wave sleep, which will not concern us here.
During idle periods, the LFP of the hippocampal region loses theta oscillations but is instead

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

1 of 34

Research article

Neuroscience Physics of Living Systems

intermittently punctuated by sharp-wave ripples with power across 140–200 Hz (Buzsáki, 1986;
Chrobak and Buzsáki, 1996). During these events, spatially tuned neurons fire in coordinated bursts
that represent rapid, long-distance trajectories. These replays are well-established in place cells(Lee and Wilson, 2002; Foster and Wilson, 2006) and have recently been observed in grid cells
within superficial layers of the MEC (O’Neill et al., 2017). Note that replays can also involve gridcells within deep layers of the MEC (Ólafsdóttir et al., 2016); these findings, which have been disputed (Trimper et al., 2017), will not be addressed by this work.
Hippocampal experiments on theta sequences and replay have established a rich phenomenology, which includes their properties near branch points of a maze (Wu and Foster, 2014;
Belchior et al., 2014), their modulation by reward history (Johnson and Redish, 2007;
Wikenheiser and Redish, 2015a; Ambrose et al., 2016; Wu et al., 2017), their ability to predict
future actions (Pfeiffer and Foster, 2013; Singer et al., 2013; Wikenheiser and Redish, 2015b;
Xu et al., 2019), and their potential role in disease pathogenesis (Middleton et al., 2018). Thesefindings suggest that hippocampal sequences facilitate many cognitive functions, such as memory
consolidation, credit assignment, and action planning (Foster, 2017; Ólafsdóttir et al., 2018;
Zielinski et al., 2018). However, the mechanisms that produce theta sequences and replay amongplace and grid cells are still unclear.
Meanwhile, a mechanism for producing grid cells themselves, the continuous attractor model, hasbeen accumulating experimental support (Couey et al., 2013; Yoon et al., 2013; Heys et al., 2014;
Dunn et al., 2015; Fuchs et al., 2016; Zutshi et al., 2018; Gu et al., 2018). This class of networks
posits particular configurations of synapses within the MEC that are local and symmetric (Burak andFiete, 2009; Widloski and Fiete, 2014). By incorporating key biological features such as fully spiking neural dynamics, we find that new phenomena emerge in this recurrent architecture: grid cells
will exhibit either theta sequences or replays depending on the external input provided by otherbrain regions. These inputs correspond to experimentally measured changes between active andinactive behavior (Steriade et al., 2001; Simon et al., 2006) and do not involve changes in connectivity. From this perspective, the same architecture that produces grid cells causes them to participate in sequences that may be leveraged by other brain regions, such as the hippocampus, topursue desired cognitive goals.

Results
Our fully spiking continuous attractor model produces grid cells
We first describe our implementation of a continuous attractor network and demonstrate that it
reproduces classic grid responses when a simulated animal explores a 2D open field. Model details
are provided in Materials and methods and Table 1. We simulate a single grid module—that is, gridcells of one scale and orientation (Stensola et al., 2012)—with a continuous attractor network basedon Burak and Fiete (2009) and Widloski and Fiete (2014). Neurons are assembled in a 2D sheetwith five overlapping populations, of which four represent excitatory principal cells and one represents inhibitory interneurons (Figure 1A). Each excitatory neuron excites its neighbors across all
populations, and each inhibitory neuron inhibits only excitatory neurons located a certain distanceaway. We will address the purpose of having four slightly different excitatory populations shortly.
To accurately simulate rapid sequences of spikes, we implement fully spiking grid cells that obey
leaky integrate-and-fire dynamics. Each neuron has a membrane potential that tends toward a
steady-state value called the drive, which represents a combination of resting potential and broad
input from other brain regions (Figure 1B). If the potential exceeds a threshold of 1 in arbitrary units,a spike is emitted, the potential is reset to 0, and the neuron’s targets experience a postsynaptic
jump in potential after a brief synaptic delay. During locomotion, excitatory principal cells are driven
to fire by various cortical and subcortical inputs (Kerr et al., 2007; Bonnevie et al., 2013). Thus,
their drive exceeds threshold at the center of the neural sheet; it decays towards the edges to avoid
edge effects that disrupt continuous attractor dynamics (Burak and Fiete, 2009). Meanwhile, inhibitory interneurons have subthreshold drive that oscillates at a theta frequency of 8 Hz due to input
from the medial septum (Brandon et al., 2011; Koenig et al., 2011; Gonzalez-Sulser et al., 2014;
Unal et al., 2015). With this architecture and random initial membrane potentials, the neural sheet
naturally develops local regions of activity, called bumps, that are arranged in a triangular lattice andKang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

2 of 34

Research article

Neuroscience Physics of Living Systems

Table 1. Main model parameters and their values unless otherwise noted.
Values that change between runs, idle periods, and allocentric corrections are indicated accordingly.
Parameter

Variable

Value

Neurons per populationnn

232  232

Simulation timestep

Dt

1 ms

Exc. membrane time constant

tþ
m

40 ms

Inh. membrane time constant

tm

20 ms

Exc.-to-exc. synaptic delay

tþþ
s

5 ms

Exc.-to-inh. synaptic delay

ts þ

2 ms

Inh. synaptic delayts

Exc. drive maximum

aþ
max

2 ms
(
2:0
1:6
2:0

Exc. drive minimumaþ
minExc. drive scaled spreadaþInh. drive magnitude

amag

Inh. drive theta amplitude

athInh. drive theta frequency

f

8 Hz

Exc. synaptic strength

wþ
mag

0.2

Exc. synaptic spreadrwþ

6Inh. synaptic strength

wmag

2.8

Inh. synaptic distancerw

12Exc. synaptic shift3

Exc. velocity gain

a

runs
idle
allo.

0.8
n
1:2 runs
0:9 idle(
0:72 runs
0:0
idle
0:72 allo.
(
0:2 runs
0
idle
0
allo.

0.25 s/m

Exc. noise magnitudeP

Var½z ðr; tÞ

0.0022Inh. noise magnitude

Var½z ðr; tÞ

0.0022are coherent across populations (Figure 1C and Figure 1—video 1). This self-organized grid is an
attractor state of the network.
How do we produce neurons with grid-like spatial tuning from the grid-like activity pattern on theneural sheet? This transformation is performed by the four excitatory populations, each of which corresponds to a perpendicular direction along the neural sheet and a perpendicular direction in thespatial environment. Excitatory neurons have output synapses biased in their preferred sheet direction (Figure 1A), and they receive more input when the animal moves along their preferred spatial
direction. When the animal moves along, say, the ’North’ direction, neurons in the N population
have increased drive; since their outputs are biased in the ’up’ direction on the neural sheet, theactivity pattern moves up. In this way, the activity pattern on the sheet moves synchronously with theanimal’s 2D trajectory (Figure 1D), and the grid pattern is projected onto single neuron spatial
responses (Figure 1E and Figure 1—video 2). Each active neuron in the sheet is a grid cell with thesame scale and orientation (Figure 1—figure supplement 1).

Grid cells are spatially tuned and theta-modulated along a linear track
For the rest of the paper, we simulate 1D trajectories consisting of runs along a linear track separated by idle periods at either end (Figure 2A). As in the 2D case, attractor bumps on the neural
sheet move synchronously with animal motion. But over the course of a simulation, the grid-like pattern on the neural sheet drifts (Burak and Fiete, 2009) such that for a given track position, theKang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

3 of 34

Research article

S
E

inh.

S

inh.

D

inhibitory

threshold
drive

reset

membrane potential

excitatory

N
W
S
E

W

C

run

B

N

membrane potential

A

Neuroscience Physics of Living Systems

reset

E

spatial coords.
North
West

threshold
drive

East

South

spikes
40 ms ago

0

now

1
corr.

Figure 1. Model architecture and generation of 2D grid cells. (A) Our model consists of a neural sheet with five overlapping populations, four of them
excitatory—N, W, S, and E—and one inhibitory. Each density plot depicts the synaptic outputs in the sheet of a neuron at the origin. (B) Each neuron isdriven to a particular membrane potential, which exceeds the spiking threshold for excitatory neurons at the center of the sheet and oscillates at 8 Hz
for inhibitory neurons while the animal is running. (C) Snapshot of neural activity showing S and inhibitory populations separately; other excitatory
populations have activity patterns similar to that of the S population. Each pixel is a neuron and dark colors indicate recent spiking. (D) Left, segment of
a 2D open field trajectory. Right, neural activity over the course of the segment with each neuron colored according to the position at which it attained
its maximum firing rate. Each attractor bump moves in synchrony with animal motion. (E) Left, two sample grid cells with spikes shown as colored dots
superimposed on the animal’s trajectory. Each neuron’s location in the sheet is indicated by a circle of corresponding color in D. Right, autocorrelationof rate maps calculated from spikes at left. Black scale bars, 50 neurons. White scale bars, 50 cm.
The online version of this article includes the following video and figure supplement(s) for figure 1:
Figure supplement 1. 2D firing maps for many grid cells.
Figure 1—video 1. Emergence of a grid-like pattern on the neural sheet from randomly initialized membrane potentials.
https://elifesciences.org/articles/46351#fig1video1
Figure 1—video 2. Path integration over an open field trajectory produces 2D grid cells.
https://elifesciences.org/articles/46351#fig1video2

corresponding bump locations on the neural
sheet slowly change. This drift introduces errors
in path-integration (Hardcastle et al., 2015),
which are believed to be corrected by allocentric
input from border cells (Ocko et al., 2018;
Keinath et al., 2018), boundary vector cells(Evans et al., 2016), or landmark cells(Raudies and Hasselmo, 2015; Savelli et al.,2017). Thus, we implement brief allocentric corrections in our model. Attractor bump locations
on the neural sheet that correspond to either
end of the track are learned during simulation
setup. They are periodically reintroduced during
the main simulation as excitatory input between
idle periods and runs (Figure 2A,B; seeMaterials and methods). In Appendix 2, we present a version of our model without allocentric

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

Video 1. Neural activity during runs, idle periods, andallocentric corrections. Left, position of the animal(black square) along a 1D track. Right, neural activity of
the E population. Each pixel is a neuron, with blackcorresponding to current spikes and lightest graycorresponding to spikes 40 ms ago. Red circles
indicate regions of recording.
https://elifesciences.org/articles/46351#video1

4 of 34

Research article

Neuroscience Physics of Living Systems

corrections that still demonstrates many of our results, though to a weaker degree.
Grid cells at different locations in the sheet generally correspond to different positions along the
track (Figure 2C). From four regions in the sheet chosen for recording, we select up to 150 excitatory neurons whose spatial firing fields are stable from lap to lap and collectively tile the track(Figure 2D; Materials and methods). Over replicate simulations with different random initial membrane potentials, many simulations contain only grid cells with single fields as in Figure 2D. In otherreplicates, the attractor grid self-organizes on the neural sheet with an orientation such that two
attractor bumps pass over recorded neurons and produce two fields (Figure 2—figure supplement
1; Appendix 1). All of our single-simulation examples presented below in the main text correspondto the simulation depicted in Figure 2D that exhibits single fields. Similarly, the grid cells in
O’Neill et al. (2017), with which we will primarily compare our results, often possess a single dominant field along an arm of a T-maze. However, both types of simulations contribute to our results,
and we will provide examples for simulations exhibiting two fields in our figure supplements.
At the theta timescale, grid cell activity is highest when the oscillating inhibitory drive is lowest
and excitatory cells are most disinhibited (Figure 2E). This observation allows us identify theta
phases in our simulations with those defined through LFP measurements. Experimentally, grid cells(Hafting et al., 2008) and place cells (Buzsáki, 2002; Feng et al., 2015) are most active at troughsof the theta-filtered LFP, which can be explained in terms of extracellular currents generated by neural activity (Buzsáki et al., 2012). O’Neill et al. (2017) define theta cycles to span from trough totrough of the LFP. Accordingly, we align our theta cycles to start with phase 0˚ at troughs of theinhibitory drive (Figure 2E).

Theta phase precession arises from oscillations in attractor bump size
Phase precession results
During runs, activity passes through recorded neurons as attractor bumps move along the neural
sheet (Figure 3A, Figure 3—figure supplement 1, and Video 1). However, spike trains from single
neurons exhibit finer temporal organization with respect to the theta oscillation in a variety of ways
(Figure 3B, Figure 3—figure supplement 2, and Figure 3—figure supplement 3). Many neurons
are phase-independent and fire throughout the theta cycle. Some are phase-locking and strongly
prefer to fire at 360˚ (equivalently, 0˚ or 720˚). Finally, some exhibit phase precession with a preferred
phase that starts around 360˚ when the animal enters a grid field but decreases as the animal progresses through it.
We first quantify these relationships using circular-linear correlation (Equation 13), which indicates whether theta phase changes with position (Kempter et al., 2012). Our correlation magnitudes 0.17 ± 0.11 (mean ± s.d.; Figure 3C) are low compared to experiments, which report grid cellswith mean » 0.3 and place cells with mean » 0.4 (O’Neill et al., 2017). This difference arises fromtwo sources. Figure S10 of O’Neill et al. (2017) shows some highly correlated neurons with magnitudes up to 0.8; these are absent from our simulations. It also shows fewer neurons with correlation
value close to 0, which corresponds to either phase-independent or phase-locking neurons. Sinceboth subgroups lack a preferred theta phase that consistently changes with position, circular-linear
correlation cannot distinguish between them (Figure 3B, top two rows).
To further characterize phase behavior, we use circular-linear regression (Kempter et al., 2012),
which can differentiate between all three subgroups of phase relationships presented in Figure 3B.
Phase-independent neurons are defined to have a regression fit score less than a cutoff of 0.4(Figure 3D; Equation 11). Neurons with high fit score are then distinguished as either phase-locking
or precessing by the absolute value of their regression slope in units of field size, which is called thephase precession range (Figure 3E). We use a range of 60˚ as the cutoff. Experimentally, phase-locking and precessing neurons are found in Layers III and II (LIII and LII), respectively (Hafting et al.,2008). A wide variety of phase precession ranges are reported across subtypes of LII neurons withmedians spanning from 50˚ for pyramidal cells (Ebbesen et al., 2016) to 170˚ for stellate cells(Reifenstein et al., 2016). Our cutoff is close to the lower limit of this span. Yet, regardless of thecutoff value, regression slopes are biased towards negative values for rightward laps and positive
values for leftward laps, which means that the preferred firing phase tends to decrease as an animal
moves through a grid field. This directional tendency is a key biological feature of phase precession

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

5 of 34

Research article

Neuroscience Physics of Living Systems

AB

brief allocentric correction

C

threshold
drive

reset

D

60

membrane potential

inhibitory

membrane potential

excitatory

rate [spk/s]
0
20

threshold
drive
reset

E

inhibitory drive

cm

m

0c

Figure 2. Grid cells along a 1D track. (A) Trajectory consisting of runs along a track separated by idle periods ateither end. Between the end of an idle period and the start of a run, the network receives brief allocentric input.
(B) Allocentric input corrects the location of attractor bumps on the neural sheet (Materials and methods). (C) Left,
track diagram. Right, neural activity over runs with each neuron colored according to the track position at which it
attained maximum firing rate. Red circles indicate regions of recording. Scale bar, 50 neurons. (D) Firing fields ofrecorded grid cells sorted by position of maximum rate. (E) Multiunit activity of neurons in D averaged over theta
cycles, which span from one trough of the oscillating inhibitory drive to the next. Data is repeated over two cycles
for clarity.
The online version of this article includes the following figure supplement(s) for figure 2:
Figure supplement 1. Firing fields for multiple simulations.

seen experimentally (Hafting et al., 2008), and it is also maintained for other cutoff values applied
to the fit score (Figure 3–Figure Supplement 4A,B).
We also construct spike density plots using all neurons in each subgroup (Figure 3F). These plotsshow that phase-locking neurons have a preferred phase around 360˚ throughout the field and thatphase-precessing neurons have a preferred phase also around 360˚ at field entry. This value corresponds to theta troughs when excitatory neurons are most disinhibited (Figure 2E) and matches
measurements for LIII and LII neurons (Hafting et al., 2008; note that their phase convention differsfrom ours by 180˚). More detailed features also seem to agree between phase-precessing neurons in
our simulations and LII neurons in experiments. First, the peak density decreases by » 75˚ over the
first 50% of field progress, and second, there is a smaller density at around 60% field progress and140˚ theta phase (Figure 3F; Hafting et al., 2008). This second density has been separately related
to bursting (Climer et al., 2013) and a second source of neural input (Yamaguchi et al., 2002;
Chance, 2012); further investigation is required to reconcile these assertions with its origin in ourmodel.
The dependence of phase behavior on simulation features and parameters is explored in detail in
the figure supplements of Figure 3. We will briefly list a few key results interspersed with experimental findings that support them. Phase precession statistics do not depend on the direction of animal
motion (Figure 3–Figure Supplement 5E,F); in experiments, Climer et al. (2013) report omnidirectional precession in a 2D environment. Precession range and correlation magnitude increase withanimal speed (Figure 3–Figure Supplement 6A, p = 0.002 and 0.0006 for precession range and correlation magnitude, respectively, by pooling data across track orientations and applying the KrushalWallis H test); in experiments, Jeewajee et al. (2014) report that faster motion is associated withsteeper precession slopes and higher correlation magnitudes (Figure 3a, results for pdcd, which corresponds to field progress defined in this paper).

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

6 of 34

Research article

A

Neuroscience Physics of Living Systems

Binhibitory drivephaseindependent

phaselocking

phaseprecessing

C

D

F

all
3264 neurons

phaseindependentEphaselocking or
precessing

phase-independent
2123 neurons

phaseprecessing

phase-locking
653 neurons

phaselocking

phaseprecessing

phase-precessing
488 neurons

densitymax
0

Figure 3. During runs, certain neurons exhibit theta phase locking and precession. (A) Grid cell spike rasters for four laps along the track. Vertical bluelines indicate theta cycle boundaries. (B) Relationship between animal position along the track and theta phase for representative neurons. Dots
represent spikes during runs in the directions indicated by arrows, with each spike repeated at two equivalent phases for clarity. Lines indicate fit bycircular-linear regression. Numbers in each panel from top left to top right indicate magnitudes of correlation coefficients, regression fit scores, andregression slopes. (C–E) Data across all replicate simulations. (C) Magnitudes of circular-linear correlation coefficients. Mean ± s.d.: 0.17 ± 0.11. (D) Fit
scores for circular-linear regression. (E) Regression slopes in units of field size for neurons with fit score >0.4. The predominance of negative values forrightward runs and positive for leftward runs indicates decreasing phase as the animal traverses grid fields in either direction. (F) Spike densities fordifferent subgroups. An animal enters a grid field at progress 0% and exits it at 100%.
Figure 3 continued on next page

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

7 of 34

Research article

Neuroscience Physics of Living Systems

Figure 3 continued
The online version of this article includes the following figure supplement(s) for figure 3:
Figure supplement 1. Spike rasters over multiple runs.
Figure supplement 2. Position-phase relationships during rightward runs for many neurons.
Figure supplement 3. Position-phase relationships during leftward runs for many neurons.
Figure supplement 4. Further analysis of subgroups with different phase relationships.
Figure supplement 5. Dependence of phase relationships on excitatory population and number of fields.
Figure supplement 6. Phase relationship properties over different simulation parameters.
Figure supplement 7. Phase relationship properties over different simulation parameters, continued.

To summarize our phase precession results, our simulation produces a subgroup of phase-locking
neurons akin to MEC LIII neurons and a subgroup of phase-precessing neurons akin to MEC LII neurons. Further consideration of correlation coefficients and precession ranges, with connections tomeasurements, will be provided in the Discussion.

Attractor bump oscillations explain phase precession
We can understand the emergence of phase precession through the effect of oscillating inhibitory
drive on attractor bumps. Higher inhibitory drive suppresses the activity of the excitatory grid cells
and decreases the size of bumps (Figure 4A,B). Imagine we record from a single neuron as a bumpof activity moves through it on the neural sheet. As an attractor bump moves, its size oscillates with
inhibitory drive (Figure 4C, top, and Figure 4—video 1). The activity bump will most likely first reachthe grid cell when it is large and growing in size, which corresponds to the end of the theta cycle
when the inhibitory drive is decreasing towards its lowest point. Thus, as an animal enters a gridfield, the first spike will likely have large theta phase slightly less than 360˚. During the next theta
cycle, the center of the bump has moved closer to the neuron, so the edge of the bump will reachthe neuron when it is growing in size but not as large. Thus, the next spike will have slightly smaller
theta phase. This is the origin of phase precession in our model. We continue generating spikes withsimplified dynamics in which the neuron fires when it is contained within the bump, but there is a
40 ms refractory period corresponding to membrane potential building up towards threshold after itis reset. This procedure represents one lap along the track at constant velocity.
We represent additional laps by performing the same simplified procedure as above, but with different shifts in theta phase, since the animal does not always enter the grid field at the same initial
phase. Spikes accumulated across laps show phase precession (Figure 4C, bottom), though theeffect is strongest at the beginning of the field. In the middle of the field, spikes still precess, butthey cluster around values of both 360˚ and 180˚. At the end of the field, a few spikes even precessin the opposite direction, increasing in phase with time. Thus, bump oscillations alone can explain
the direction of phase precession and the preferred phase at field entry, but spike phases do notdecrease perfectly linearly with progress through the field.
To analyze attractor bump oscillations more precisely, we calculate the average bump shape as afunction of theta phase (Figure 4D and Figure 4—video 2; Appendix 1). We now imagine that thisaverage bump passes multiple times through a recorded neuron whose spiking probability is proportional to bump activity (Figure 4E). Each instance of this stochastic process produces a phase-locking, phase-precessing, or phase-independent neuron. By explicitly rescaling the activity of theaverage bump without otherwise changing its time-varying shape, we can generate neurons
enriched in one of these phase relationships (Figure 4F and Figure 4—figure supplement 1). Underhigh activity, neurons are driven to fire across all theta phases and exhibit phase independence.
Under low activity, neurons only fire during the most permissive theta phase and exhibit phase locking at 360˚. Under intermediate levels of activity, neurons can respond to oscillations in bump size as
demonstrated in Figure 4C and exhibit phase precession.
This finding on how activity level influences phase behavior is supported by the experimental
report that higher spike counts are associated with steeper precession slopes (Jeewajee et al.,2014, Figure 3b, results for pdcd, which corresponds to field progress defined in this paper). Our
main model also reflects this finding: neurons with low, medium, and high firing rates have increased
likelihoods for phase locking, precession, and independence (Figure 3–Figure Supplement 4D).

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

8 of 34

Research article

Neuroscience Physics of Living Systems

A inhibitory drive = 0.35

Binhibitory drive = 0.85
N

72°

0°

inhibitory driveE

inhibitory drive

N

D
theta
phase

C

144°

216°

288°

average
bump
shape
0

F

low bump activity

max
activity

direction of
bump motion

medium bump activity

high bump activity

Figure 4. Bump size oscillations explain phase precession in models with simplified dynamics. (A,B) Data from simulations with fixed inhibitory driveand constant animal velocity. (A) Snapshots of neural activity. Scale bar, 50 neurons. (B) The diameter of bumps on the neural sheet decreases linearlywith inhibitory drive (linear regression R2 ¼ 0:90, ANOVA p ~ 10 9 ). (C) Phase precession in a conceptual model with bump size oscillations. We imagine
an attractor bump, with size oscillations described by B, passing through a recorded grid cell. Top, a single lap. The recorded neuron is at location 0
and fires a spike (black dot) whenever contained within the bump (gray area), subject to a 40 ms refractory period. Bottom, relationship between theta
phase and time across multiple laps with different initial phases. Spikes occur around 360˚ (equivalently, 0˚) at the start of the field, and their phase
generally decreases with time within in the grid field. (D) Attractor bump shape at different theta phases averaged over theta cycles and individual
bumps (Appendix 1). Grays scaled separately for each theta phase. Scale bar, three neurons. (E,F) Phase behavior in a simplified model using averagebump dynamics. (E) We imagine the average attractor bump passing through a recorded grid cell. Top, a single lap. The recorded neuron is at location
0 and stochastically fires a spike (black dot) with rate proportional to bump activity, subject to a 40 ms refractory period. Bottom, relationship between
time and theta phase across multiple laps with different initial phases. (F) Relationship between time and theta phase using average bumps whose
activity has been rescaled to different maximum values: 40, 50, and 100 spikes/s. Dots represent spikes generated according to E. Lines indicate fit bycircular-linear regression. Numbers in each panel from top left to top right indicate magnitude of correlation coefficient, regression fit score, andregression slope.
The online version of this article includes the following video and figure supplement(s) for figure 4:
Figure supplement 1. Simplified model using average bump dynamics produces different phase behaviors for different levels of rescaled activity.
Figure 4—video 1. Phase precession in a simplified conceptual model.
https://elifesciences.org/articles/46351#fig4video1
Figure 4—video 2. Average attractor bump shape as a function of theta phase.
https://elifesciences.org/articles/46351#fig4video2

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

9 of 34

Research article

Neuroscience Physics of Living Systems

Thus, the heterogeneous phase behaviors observed in our main simulations arise, at least in part,from heterogeneous levels of bump activity. Three factors that determine the level of activity experienced by a neuron are its distance to the center of the neural sheet, its sheet location relative toattractor bumps, and its preferred firing direction relative to animal motion. Accordingly, the prevalence of different phase behaviors varies with these factors (Figure 3–Figure Supplement 4E,G and
Figure 3–Figure Supplement 5B,D).

Theta sequences arise from oscillations in attractor bump speed
Next, we use the firing fields illustrated in Figure 2D to decode the animal’s position from the population activity presented in Figure 3A (Materials and methods). The decoded position generally
matches the animal’s actual position (Figure 5A and Figure 5—figure supplement 1), but at thetheta timescale, there are deviations whose regularities are revealed by averaging over theta cycles.
To do so, we take quadruplets of consecutive theta cycles, ignoring those whose decoded positions are close to the ends of the track and reversing those corresponding to right-to-left motion(Materials and methods). We align the decoded positions with respect to the actual position of theanimal midway through each theta quadruplet, and we average these shifted decoded positions
over theta quadruplets. This average decoded trajectory generally increases across the quadruplet,corresponding to the animal’s actual forward motion (Figure 5B and Figure 5—figure supplement2). However, within each theta cycle, the decoded position increases rapidly, before retreating at
cycle boundaries. We identify these forward sweeps as theta sequences, which represent motion at
0.88 ± 0.15 m/s, approximately twice the actual speed of 0.50 ± 0.05 m/s (mean ± s.d.; Figure 5C).

A

inhibitory driveBinhibitory drive

prob.
0.5
0

prob.
max
0

C

D

inhibitory drive

Figure 5. During runs, the decoded position exhibits theta sequences. (A) Decoded position for four laps along
the track corresponding to Figure 3A. Vertical blue lines indicate theta cycle boundaries. (B) Decoded position
shifted by actual position and averaged over theta cycles. The thick black fit line indicates the time window withthe steepest increase in decoded position. Vertical thin black lines indicate theta cycle boundaries. (C) Mean
actual run speed (error bars indicate s.d. over time) and mean theta sequence speed, as indicated by the slope ofthe thick black line in B (error bars indicate s.d. over replicate simulations). (D) Difference between the maximum
likelihood decoded position from B and the actual animal position as a function of time. Within each theta cycle,
the segment from the smallest to the largest value is emphasized as a theta sequence.
The online version of this article includes the following figure supplement(s) for figure 5:
Figure supplement 1. Decoded position over many runs.
Figure supplement 2. Decoded position shifted by actual position and averaged over theta cycles for many
simulations.
Figure supplement 3. Dependence of theta sequence speed on number of fields and simulation parameters.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

10 of 34

Research article

Neuroscience Physics of Living Systems

This matches experimental observations of average theta sequence speed » 1 m/s compared to theanimal’s speed of » 0.4–0.5 m/s (O’Neill et al., 2017). We also illustrate these sequences by comparing the maximum likelihood decoded position and the actual position as a function of time
(Figure 5D). Averaged over theta cycles, the decoded position lags behind the actual position at
cycle boundaries and then advances ahead of it midcycle.
Theta sequences, like phase precession, arise naturally in our model from an effect of oscillating
inhibitory drive on attractor bumps. In addition to determining bump size, inhibitory drive also
affects the speed at which bumps move along the neural sheet—higher inhibitory drive produces
faster bump motion for a given animal speed (Figure 6A,B). At the middle of theta cycles, inhibitory
drive is highest, so attractor bumps move quickly and represent fast animal trajectories. These theta
sequences are faster than the animal’s actual speed, which is approximately constant within a theta
cycle and thus encoded as the average bump speed over the entire cycle. We can also see the originof theta sequences directly from the average dynamics of attractor bumps (Figure 4D, Figure 6C,
and Figure 4—video 2). The location of peak activity moves with a sawtooth-like pattern that, overa theta cycle, lags behind the overall motion of the bump and advances ahead of it (Figure 6D).
These forward surges produce the sequences in decoded position (Figure 5D).
Our model’s mechanism for theta sequences can explain why their speed is roughly twice theactual speed, as observed experimentally (O’Neill et al., 2017). Suppose the theta oscillation
spends roughly the same amount of time at its peaks, when attractor bumps move at maximal
speed, and at its troughs, when bumps do not move at all. Then, the average bump speed, whichA

inhibitory drive = 0.35

1 s agoN

now
spikesD
C

B

inhibitory drive = 0.85
Ninhibitory drive

inhibitory drive

activitymax
0

Figure 6. Bump speed oscillations explain theta sequences. (A,B) Data from simulations with fixed inhibitory driveand constant animal velocity. (A) Snapshots of neural activity. Scale bar, 50 neurons. (B) The speed of bumps on
the neural sheet increases linearly with inhibitory drive (linear regression R2 ¼ 0:91, ANOVA p ~ 10 9 ). Thus, theaverage decoded position in Figure 5B advances most rapidly around the middle of the theta cycle. (C) Activity of
the average attractor bump depicted in Figure 4D along its central axis. (D) Top, peak activity location from C asa function of time. Bottom, peak activity position relative to constant bump motion as a function of time. Withineach theta cycle, the segment from the smallest to the largest value is emphasized and corresponds to a theta
sequence in Figure 5D.
The online version of this article includes the following figure supplement(s) for figure 6:
Figure supplement 1. Relationship between attractor bump size and speed.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

11 of 34

Research article

Neuroscience Physics of Living Systems

encodes the actual animal speed, is half the maximal bump speed, which encodes the theta
sequence speed. Indeed, simulations whose grid cells contain two fields, whose track is oriented differently, whose trajectory is faster or slower, or whose parameters are chosen differently all havetheta sequence speeds between 1.2 and 2.5 times the actual speed (Figure 5—figure supplement
3).

Replay arises from traveling wavefronts of activity
We now simulate idle periods at either end of the track. We would like to know how external inputs
to the MEC change at transitions between active and quiescent states, but these measurements
have not been performed to our knowledge. Thus, we instead adapt experimental findings taken
during slow wave sleep, which resembles the quiescent state in many behavioral and electrophysiological ways (Buzsáki, 1996). The neocortex, which activates MEC principal cells, is broadly lessactive during slow wave sleep (Steriade et al., 2001), and the medial septum, which inhibits MEC
interneurons, increases its firing rate while losing most of its oscillatory nature (Simon et al., 2006).
We model these effects respectively as decreased excitatory drive and decreased inhibitory drive
without oscillations (Figure 7A). The network rapidly switches between these inputs and those of
Figure 1B, which represent the active state, as the animal transitions between quiescence and runs.
During quiescence, attractor bumps still form a lattice at the center of the neural sheet, where
excitatory drive still exceeds threshold (Figure 7B and Video 1). However, these bumps vanish
toward the edge of the sheet, where they are replaced by traveling wavefronts of activity. These
wavefronts cannot be solely sustained by excitatory drive, which is subthreshold in this outer region;
instead, they must be nucleated and then self-sustained through excitatory-to-excitatory connections. In our model, nucleation occurs at activity bumps toward the edge of the sheet (Video 1).
Once nucleated, the wavefront propagates freely through the outer region because inhibitory neurons, with low drive, are not activated by the wavefront. In contrast, wavefronts cannot propagate in
the center region, because excitatory neurons there receive enough drive to spike vigorously andactivate nearby inhibitory neurons (Figure 7B, inhibitory population). These interneurons, with theirsurround distribution of synaptic outputs (Figure 1A), constrain neural activity to localized bumps
and prevent wavefront propagation.
When a wavefront passes through a region of recorded neurons, it can appear as a sequence of
spikes (Figure 7C and Figure 7—figure supplement 1) that is decoded as a rapid trajectory(Figure 7D and Figure 7—figure supplement 2). We identify these events as replays, which traverse
much of the track over ~100 ms (Figure 7E and Figure 7—figure supplement 3). Due to the stochastic nature of wavefront generation, the number of detected replays can vary considerably acrosslaps (Figure 7—figure supplement 2) and across replicate simulations (Figure 7F). However, replays
have a characteristic speed determined by that of wavefront propagation; it is 7.0 ± 2.6 m/s, which isabout 14 times faster than the actual speed of 0.50 ± 0.05 m/s (mean ± s.d.; Figure 7G). This agrees
with experiments, which measure replays at » 4–6 m/s and actual motion at » 0.4–0.5 m/s (O’Neill et al., 2017). Replay speed depends on the number of grid fields and the values of simulation parameters (Figure 7—figure supplement 4). In fact, grid cells with two fields encode two different track positions, so replays among them can have their distribution of decoded position split
between parallel lines (Figure 7—figure supplement 3). We have implemented two different detection methods that yield slightly different replay speeds: one that is identical to the single-field case,and another in which multiple lines participate in the fit (Figure 7—figure supplement 4A-D ;Appendix 1). Nevertheless, across all simulation parameters and detection methods, replays arealways much faster than both the actual animal speed and the theta sequence speed, as observed
experimentally (O’Neill et al., 2017). This happens because traveling wavefronts lack inhibition andshould be the network phenomenon that propagates fastest.

Replay properties under different network configurations
Lower velocity gain produces replays with directional bias
Under the main parameter values in Table 1, attractor bumps travel long distances on the neural
sheet when the animal moves from one end of the track to another (Figure 2C). These distances are
longer than the grid scale, so different positions on the track can correspond to the similar locations
on the neural sheet, preventing a consistent correspondence between neural sheet location andKang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

12 of 34

Research article

Neuroscience Physics of Living Systems

A

B

idle

S

inh.

inhibitory

threshold
drive

reset

membrane potential

membrane potential

excitatory

threshold

drive
reset

spikes
40 ms ago

C

now

D

prob.
0.5
0

E

F

G

Figure 7. Lower drive during quiescence produces traveling wavefronts of activity that are decoded as replays. (A) During idle periods at either end of
the track, excitatory and inhibitory drives decrease, and the latter no longer oscillates. (B) Snapshot of neural activity showing S and inhibitory
populations separately. Red circles indicate regions of recording. Scale bar, 50 neurons. (C) Grid cell spike rasters for four idle periods. (D) Decoded
position corresponding to C. Rapid trajectory sweeps arise from traveling wavefronts as depicted by the S population in B. (E) Replays from one
simulation with cyan fit lines. (F) Number of replays across simulations. (G) Mean actual run speed (error bars indicate s.d. over time) and mean replay
speed, as indicated by the slopes of cyan lines in E (error bars indicate s.d. over replays).
The online version of this article includes the following figure supplement(s) for figure 7:
Figure supplement 1. Spike rasters over multiple idle periods.
Figure supplement 2. Decoded position over many idle periods.
Figure supplement 3. Many replays.
Figure supplement 4. Dependence of replay count and speed on number of fields and simulation parameters.

track position (yellow and purple areas in Figure 2C). With lower velocity gain, attractor bumps
travel shorter distances, and such a correspondence is apparent (yellow and purple areas in
Figure 8A).
When an animal is idle at either end of the track, attractor bumps are still located at neurons that
represent the current position (Figure 8A and Figure 8—video 1). Recall that replays are nucleated
at attractor bumps (Video 1 and Figure 8A). As a wavefront radiates from a bump, it will likely firstactivate neurons whose firing fields represent nearby positions, followed by those with increasingly
distant fields, due to the correspondence between neural sheet locations and track positions

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

13 of 34

Research article

Neuroscience Physics of Living Systems

lower velocity gain

A

S

60

B

C

cm

m

0c

spikes within last 10 ms

F

E
drive
threshold

reset

S

inh.

G

membrane potential

idle
excitatory

membrane potential

periodic boundary conditions

Drun
excitatory

threshold
drive

reset

S

inh.

Figure 8. Replays acquire additional properties in networks with different configurations of attractor bumps. (A–C) With lower velocity gain, neural
sheet locations have a clear one-to-one correspondence with track positions, which allows replays to exhibit a directional bias. (A) Left, track diagramwith animal idle at 0 cm. Right, corresponding snapshot of neural activity superimposed on a background colored according to the position at which
each neuron attains its maximum firing rate (cf. Figure 2C, right panel). Red circles indicate regions of recording. Attractor bumps still represent the
animal position, and a wavefront can be seen emanating from one bump at the top right. (B) Direction of replay propagation across simulations. For
each actual position and replay direction, width of gray bars indicates the number of simulations with a particular replay count. Means indicated withlines. p-values from the Mann-Whitney U test show that a significantly greater number of replays propagate away from either actual position. (C)Distances between actual position and decoded position at the start or end of replays. Paired medians compared by the Wilcoxon rank-sum test withindicated p-value. (D–G) With periodic boundary conditions and uniform excitatory drive, attractor bumps and activity wavefronts appear throughout
the neural sheet. (D,E) Excitatory drive during runs and idle periods. (F,G) Snapshots of neural activity during runs and idle periods showing S andinhibitory populations. Each pixel is a neuron and dark colors indicate recent spiking. Red circles indicate regions of recording. Scale bars, 50 neurons.
The online version of this article includes the following video and figure supplement(s) for figure 8:
Figure supplement 1. Replay properties in simulations with lower velocity gain.
Figure supplement 2. Replay properties in simulations with periodic boundary conditions and uniform excitatory drive.
Figure supplement 3. Phase precession and theta sequence properties in simulations with periodic boundary conditions and uniform excitatory drive.
Figure 8—video 1. Replays preferentially radiate from attractor bumps encoding the current animal position.
https://elifesciences.org/articles/46351#fig8video1
Figure 8—video 2. Neural activity during runs, idle periods, and allocentric corrections for a simulation with periodic boundary conditions and uniform
excitatory drive.
https://elifesciences.org/articles/46351#fig8video2(Figure 8—video 1). Thus, the decoded position tends to start near the actual position of the animal
and travel away from it (Figure 8B). Accordingly, the decoded starting position of replays are closer
to the actual position than the decoded ending position (Figure 8C).

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

14 of 34

Research article

Neuroscience Physics of Living Systems

However, replays are not all detected immediately after wavefront nucleation; wavefronts traveling for some time may pass through a region of recorded neurons from any direction. Thus, despite
the predominance of replays directed away from the actual animal position, some start at the opposite side of the track and travel toward the actual position (Figure 8B). This directional bias has notbeen directly assessed for grid cells, but it appears in experiments on hippocampal replay (Diba andBuzsáki, 2007; Davidson et al., 2009).

Periodic boundary conditions and uniform drive produces replays
throughout the neural sheet
Our main model has nonperiodic boundary conditions, which requires excitatory drive to decrease
towards the edge of the neural sheet to prevent boundary effects that disrupt path integration(Burak and Fiete, 2009). During idle periods, the persistence of attractor bumps at the center ofthe sheet maintain a representation of the animal’s position (Figure 7B). With periodic boundaries
that eliminate edges, we can implement uniform excitatory drive during both runs and idle periods(Figure 8D,E). This allows the entire neural sheet to exhibit attractor bumps during runs and activity
wavefronts during idle periods (Figure 8F,G and Figure 8—video 2). Thus, replays can be recorded
from anywhere on the neural sheet. Note that allocentric input at the end of each idle period isessential in this case because positional information in the form of attractor bumps completely disappears. Simulations with periodic boundary conditions and uniform drive produce all the essential
replay and theta-related phenomena observed in our main model (Figure 8—figure supplement 2
and Figure 8—figure supplement 3).

Discussion
Summary of our results
Decades of hippocampal experiments have established theta phase precession, theta sequences,and replay as central features of place cell dynamics (O’Keefe, 1976; Buzsáki, 1986; Skaggs et al.,
1996; Lee and Wilson, 2002; Foster and Wilson, 2006). Evidence has been accumulating that gridcells, discovered more recently, also exhibit these phenomena (Hafting et al., 2008;
Ólafsdóttir et al., 2016; O’Neill et al., 2017). These three phenomena appear to be similar to one
another in many ways, but experiments have also demonstrated important discordances. Theta
sequences have been understood as the direct consequence of individually phase-precessing neurons (O’Keefe and Recce, 1993; Skaggs et al., 1996). However, place cell experiments show thatphase precession can appear without theta sequences (Feng et al., 2015). Theta sequences andreplays both represent temporally compressed trajectories, and both seem important for memory
consolidation and planning (Foster and Knierim, 2012; Zielinski et al., 2018). Yet, replay trajectories are significantly faster than theta sequences, and theta sequences only propagate in the direction of motion whereas replays can propagate in any direction relative to the resting animal(Davidson et al., 2009; Feng et al., 2015; Zheng et al., 2016; O’Neill et al., 2017).
We unify phase precession, theta sequences, and replay under a single model while maintaining
distinctions among them. Phase precession and theta sequences arise naturally from oscillations in
attractor bump size and speed, respectively, which are both consequences of medial septum input.
However, bump size and speed are not generally correlated across changes in other model parameters (Figure 6—figure supplement 1). Thus, it is conceivable that one theta-related phenomenon
can exist without the other under certain conditions. Replays in our model arise from the disappearance of attractor bumps and their replacement by traveling wavefronts of activity. Wavefront propagation involves dynamical processes that are fundamentally different from attractor bump motion,and thus has led to different decoded speeds and directions for replays and theta sequences.
Crucially, we produce these phenomena in a continuous attractor network that can still generate
2D grid cells through path integration. Our implementation shares similarities with an early versionof the model (Burak and Fiete, 2009) and a version whose connectivity can be learned(Widloski and Fiete, 2014). However, there are three key differences. First, we use fully spiking neurons, which are rarely used in grid cell attractor models—we have only found them utilized in
English (2017) and Bush and Burgess (2014), the latter of which combines an attractor networkwith the oscillatory interference model for grid cells (Burgess et al., 2007; Hasselmo et al., 2007).

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

15 of 34

Research article

Neuroscience Physics of Living Systems

Typically, neural dynamics are founded on firing rates, which can be used to generate spike trains
through a memoryless Poisson process (Burak and Fiete, 2009; Widloski and Fiete, 2014); thesemodels do not capture the strong short-time anticorrelation in spiking due to the reset of membranepotential. Second, we include excitatory-to-excitatory connections in accordance with their recent
experimental discovery (Fuchs et al., 2016; Winterer et al., 2017; Schmidt et al., 2017;
Zutshi et al., 2018). These connections allow wavefronts to propagate among excitatory neurons
while nearby inhibitory neurons remain silent. Third, we vary the external input to the network
between active and quiescent behavioral states in accordance with observations (Steriade et al.,
2001; Simon et al., 2006). Consequently, the network can instantly switch between operating
regimes that exhibit either theta-related phenomena or replays.
Our model presents a distinct conception of entorhinal replay that can be tested experimentally.
Direct observation of propagating wavefronts can be pursued by calcium imaging of the MEC(Heys et al., 2014; Diehl et al., 2019). Otherwise, tetrode recordings of grid cell replay in open fieldenvironments may be analyzed for wavefront dynamics. In two dimensions, wavefronts would not
represent a single trajectory, but instead a manifold of positions that sweep out a plane throughspace. New decoding techniques that permit the grid cell population to simultaneously represent
inconsistent locations would be required to distinguish between these two possible representationsof motion.

Relationships to experiments and other models
We chose to build a continuous attractor model that can produce theta-related behavior and replay
with as simple an architecture as possible, but we imagine that more complex network geometries
are possible. For example, during quiescence, our model maintains activity bumps at the center of
network and only exhibits replays at the edges; with alternative implementations of excitatory drive,
multiple areas may contain persistent attractor bumps that nucleate wavefronts. We have also
ignored other known components of the grid system. We do not model multiple grid modules(Stensola et al., 2012; Urdapilleta et al., 2017; Kang and Balasubramanian, 2019). During runs,
spike frequency adaptation (Yoshida et al., 2013) would reduce the number of spikes late in thegrid field, which could enhance phase precession according to the conceptual model presented in
Figure 4 and increase correlation magnitude. During quiescence, the MEC exhibits Up and Down
states marked by high and low levels of activity; perhaps transitions between these states, which is
associated with CA1 ripples (Hahn et al., 2012), can nucleate multiple wavefronts that propagate insynchrony. Moreover, simulation features used to study gamma oscillations in grid cells(Nolan, 2018; Chrobak and Buzsáki, 1998; Colgin et al., 2009) can be introduced to elucidate therole of gamma during theta-related phenomena and replays (Pfeiffer and Foster, 2015; RamirezVillegas et al., 2018).
Our model does not explicitly implement biophysical distinctions between LII and LIII neurons or
between stellate and pyramidal cells within LII. These neural subpopulations have different characteristic ranges of phase precession that even vary between experimental reports. Ebbesen et al.
(2016) found the median precession range of LII stellate and pyramidal cells to be » 130˚ and » 50˚,
respectively, which is smaller than the values of » 110˚ and » 170˚ found by Reifenstein et al.
(2016). Meanwhile, LIII grid cells, which often possess conjunctive head direction tuning(Sargolini et al., 2006) like the neurons in our model, predominantly exhibit phase-locking withphase precession range close to 0˚ according to Hafting et al. (2008). However, Reifenstein et al.
(2014) reports small, but significant precession in LIII, and single-lap analyses by Ebbesen et al.
(2016) finds that LIII neurons precess even more than LII pyramidal cells. O’Neill et al. (2017) does
not distinguish between these cell types and observes all of these phase behaviors; this heterogeneity is dynamically captured by the excitatory principal cells in our model. Simulations implementing
multiple types of principal cells would be required to precisely assign the distinctions in theta phase
properties to the appropriate cell type.
Other models have been proposed for theta phase precession (Tsodyks et al., 1996;
Thurley et al., 2008; Navratilova et al., 2012; Thurley et al., 2013) and replay (Ponulak and Hopfield, 2013; Azizi et al., 2013; Bayati et al., 2015; Chenkov et al., 2017; Gönner et al., 2017;
Haga and Fukai, 2018). For example, one grid cell model uses after-spike depolarization within a
1D continuous attractor network to generate phase precession and theta sequences(Navratilova et al., 2012). This experimentally observed feature pulls attractor bumps back toKang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

16 of 34

Research article

Neuroscience Physics of Living Systems

recently active locations every theta cycle. A place cell model combines an oscillating membranepotential, as found in our model, with synaptic facilitation (Thurley et al., 2008). These features
could be added to our model to strengthen theta-related phenomena and potentially increase correlation magnitude, but they not required. Previously proposed replay models were all intended toaddress place cells, not grid cells. Several of them encode replay trajectories into synaptic weights,
either through hard-wiring or a learning mechanism (Chenkov et al., 2017; Gönner et al., 2017;
Haga and Fukai, 2018). Two models have suggested that replays originate from wavefronts of activity propagating through networks of place cells (Ponulak and Hopfield, 2013; Bayati et al.,2015; Palmer et al., 2017). These wavefronts then enhance connections through the network
though spike-timing-dependent plasticity rules, which could account for the ability of reward tomodulate hippocampal replay (Pfeiffer and Foster, 2013; Ambrose et al., 2016; Palmer andGong, 2013).

Possible implications
Although phase precession, theta sequences, and replay arise without learning in our model, they
are not necessarily epiphenomenal in the sense that they serve no cognitive purpose. Grid cells arethought to provide a spatial representation that is immediately available upon entering an environment and that remains stable over days of experience (Hafting et al., 2005; Leutgeb et al., 2007;
Diehl et al., 2019). They provide spatial information to the hippocampus (Cheng and Frank, 2011;
Hales et al., 2014), whose place cells associate location with environmental context (Wilson and
McNaughton, 1993; Anderson and Jeffery, 2003) and reward (Hollup et al., 2001; Mamad et al.,2017) to aid the animal in pursuing its goals. In this vein, we envision that superficial MEC layers act
as a pattern generator for phase precession, theta sequences, and replay, especially in new environments. These phenomena are then presented to the hippocampus to be refined by experience andmodulated by reward.
The concept that sequences in hippocampus arise from superficial MEC layers has been experimentally substantiated for theta phenomena: lesioning the MEC disrupts phase precession and theta
sequences in place cells (Schlesiger et al., 2015), but grid cell phase precession is not affected byhippocampal inactivation (Hafting et al., 2008). Driving place cell replays with grid cell replays
would require temporal coordination between the two, for which experimental support is mixed.
O’Neill et al. (2017) does not find coordination between grid cell replays in superficial MEC layers
and sharp-wave-ripples (SWRs) in CA1, during which place cell replays appear. However,Ólafsdóttir et al. (2016) and Ólafsdóttir et al. (2017) report that hippocampal replays represent
trajectories coherently with grid cell activity in deep MEC layers, and Yamamoto and Tonegawa(2017) reports that SWRs in CA1 are preceded by SWRs in superficial MEC layers. Further investigation would help to elucidate the temporal relationships between hippocampus and MEC. Also, notethat we do not exclude the possibility that place cell replays can also be initiated within the hippocampus or by subcortical nuclei (Stark et al., 2014; Schlingloff et al., 2014; Oliva et al., 2016;
Angulo-Garcia et al., 2018; Sasaki et al., 2018).
Place cell activity appears to contain an abundance of certain short sequences that represent contiguous locations in novel environments and participate frequently in replays (Liu et al., 2018). Thisfinding is related to the disputed (Silva et al., 2015) observation of hippocampal preplay(Dragoi and Tonegawa, 2011; Grosmark and Buzsáki, 2016; Liu et al., 2019). From the perspective of our model, such preexisting sequences may arise from wavefronts generated intrinsically by
grid cells, which then drive corresponding place cells. After experience, hippocampal plasticity supplements the sequences with additional place cells to enhance their representation of meaningfultrajectories (Grosmark and Buzsáki, 2016; Drieu et al., 2018). In this way, wavefronts in the MEC
can accelerate learning by providing a scaffold for the hippocampus to modify and improve. Meanwhile, the MEC maintains a symmetric representation of space, so the next experience, possibly withdifferent spatial properties, can be learned equally well.
Yet, we emphasize that our model focuses on the mechanistic origins of grid cell sequences andits results hold regardless of whether the suggested relationships between MEC and hippocampus
are borne out by future experiments.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

17 of 34

Research article

Neuroscience Physics of Living Systems

Materials and methods
This section provides an overview of our simulation methods. Full details are provided in Appendix
1.

Model architecture and dynamics
Our architecture is inspired by Burak and Fiete (2009) and Widloski and Fiete (2014). A 2D neural
sheet contains five overlapping populations, each with neurons at r ¼ ðx; yÞ, where x ¼ 1; . . . ; n andy ¼ 1; . . . ; n. There are four excitatory populations: N, S, W, and E; we index them by P and use thesymbol þ for parameters common to all of them. There is one inhibitory population .
Each excitatory neuron is described by a membrane potential fP ðr; tÞ and a spike indicatorsP ðr; tÞ, which equals either 1 if neuron r in population P spiked at time t or 0 otherwise. When thepotential exceeds a threshold of 1 in arbitrary units, the neuron fires a spike and its potential is reset
to 0:
fP ðr; tÞ

0 and sP ðr; tÞ

(1)

1:

We prevent fP ðr; tÞ from decreasing past –1 as a limit on hyperpolarization. The same definitions
and dynamics apply to the inhibitory population .
Membrane potentials follow leaky integrator dynamics. The excitatory neurons obey
tþ
m
¼

X
P0 ;r0wþ ðr

r0

fP ðr; t þ DtÞ fP ðr; tÞþ fP ðr; tÞDtX00w ðr r0 Þs ðr0 ; t^eP ÞsP ðr0 ; t tþþ
s Þþ
r0

ts Þ

(2)



^ P  VðtÞ þ zP ðr; tÞ:þ aþ ðrÞ 1 þ aE

And the inhibitory neurons obeyf ðr; t þ DtÞ f ðr; tÞþ f ðr; tÞDtX00wþ ðr r0 ^
eP ÞsP ðr0 ; t ts þ Þ þ a ðtÞ þ z ðr; tÞ:
¼
tm

(3)

P0 ;r0

Neural drive is given by a, which is modulated by the animal’s velocity VðtÞ for excitatory neurons.
^ P , which we call North,
Each excitatory population P prefers a perpendicular direction in space E
South, West, and East for the populations N, S, W, and E. Neural connectivity is given by w. Each
excitatory population P has its synaptic outputs shifted by a small amount  in a preferred direction^
eP on the neural sheet, which we call up, down, left, and right for the populations N, S, W, and E.
Spikes produce an instantaneous change in postsynaptic membrane potential after a synaptic delay
ts , which is longer for excitatory-to-excitatory connections (tþþs ) than for excitatory-to-inhibitory conþnections (ts ) based on axon morphology (Schmidt et al., 2017). Independent, zero-mean, normally
distributed noise z is added to each neuron at each timestep. See Table 1 for complete variable definitions and values.

Drive and connectivity
The drive to excitatory populations is(þaþþmin þ ðamax
a ðrÞ ¼aþminwhere  ¼

aþ
min Þ

1þcosðp=aþ Þ2

<a+a+,(4)

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2
2 n
x nþ1
þ y nþ1= 2 is a scaled radial distance for the neuron at r ¼ ðx; yÞ. It equals 02
2

at the center of the neural sheet and approximately one at the midpoint of an edge. The drive toinhibitory subpopulations isa ðtÞ ¼ amag

where

ath cosð2pft þ

0 Þ;

(5)

0 is a phase offset chosen randomly at the start of each lap.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

18 of 34

Research article

Neuroscience Physics of Living Systems

The synaptic connectivity from excitatory neurons to all neurons is based on the symmetricfunction(þ 1þcosðpr=rwþ Þþr <rw+
2w ðrÞ ¼ wmag
(6)0rrw+,
where r ¼ jrj. The synaptic connectivity from inhibitory neurons to excitatory neurons is(w Þ
wmag 1 cosðpr=r
r <2rw2
w ðrÞ ¼
0r2rw-.

(7)

There is no connectivity from inhibitory neurons to other inhibitory neurons.
See Table 1 for complete variable definitions and values.

Brief allocentric correction
The grid-like pattern of excitatory drive during brief allocentric corrections (Figure 2B) represents
inputs from neurons with allocentric responses, such as border, boundary vector, or landmark cells(Evans et al., 2016; Raudies and Hasselmo, 2015; Savelli et al., 2017). Interactions between allocentric neurons and grid cells have also been implemented in other continuous attractor models(Ocko et al., 2018; Keinath et al., 2018). In our model, allocentric input is learned in an accelerated
manner during simulation setup. Immediately before the main simulation phase, we simulate one runin each track direction followed by a learning period at either end of the track. The locations of
attractor bumps during learning periods are stored by tallying the number of spikes produced ateach neural sheet location. These tallies are linearly rescaled to serve as the allocentric drive to excitatory populations.

Single neuron recording and decoding
We simulate our model 10 times with slightly different run trajectories randomly generated by fractional Brownian motion processes, followed by 10 more times with the same trajectories but
reversed velocities. These simulations produce spike trains for all neurons in the network. To match
multiunit tetrode recordings, we select two sets of up to 150 excitatory neurons in each simulation.
Each of these sets is an independent recording to be analyzed separately, which we call a ‘simulation’ in our results.
We select neurons for recording as follows. For the first set, we choose four points on the neural
sheet whose distance from the center is 95 neurons and whose angular positions, or clock positions,are random but equally spaced. From circular areas of radius 12 centered at these points, we randomly select up to 150 excitatory neurons whose firing fields are stable across laps. For the second
set, we choose the four circular areas with the same distance from the center of the sheet and withangular positions offset 45˚ from the first four areas. We again select up to 150 neurons with stablefiring fields. All single-simulation results reported in the main text come from the same simulation.
We use grid cell firing fields accrued from all laps in both directions to decode the animal’s position from the population activity of its recorded neurons. We use the standard Bayesian decoding
procedure that assumes independent, Poisson-distributed spike counts and a uniform prior(Zhang et al., 1998; O’Neill et al., 2017).

Acknowledgements
The authors thank John Widloski, David Foster, Fritz Sommer, Zengyi Li, Guy Isely, Kate Jeffery,
Hugo Spiers, and Srdjan Ostojic for stimulating discussions and helpful ideas. Part of this work was
performed at the Simons Institute for the Theory of Computing.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

19 of 34

Research article

Neuroscience Physics of Living Systems

Additional information
Funding
FunderGrant reference number

Author

Adolph C. and Mary Sprague Postdoctoral fellowship
Miller Institute for Basic Research in Science, University ofCalifornia Berkeley

Louis Kang

Army Research Office

Michael R DeWeese

W911NF-13-1-0390

The funders had no role in study design, data collection and interpretation, or thedecision to submit the work for publication.
Author contributions
Louis Kang, Conceptualization, Software, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing; Michael R DeWeese, Resources, Supervision, Funding acquisition, Visualization, Writing—original draft, Writing—review and editing
Author ORCIDs
Louis Kang
https://orcid.org/0000-0002-5702-2740
Decision letter and Author response
Decision letter https://doi.org/10.7554/eLife.46351.sa1
Author response https://doi.org/10.7554/eLife.46351.sa2

Additional files
Supplementary files. Source code 1. Source code for the main simulation written in C.
. Transparent reporting form

Data availability
Source code for the simulations have been included as supporting files.

References
Ambrose RE, Pfeiffer BE, Foster DJ. 2016. Reverse replay of hippocampal place cells is uniquely modulated bychanging reward. Neuron 91:1124–1136. DOI: https://doi.org/10.1016/j.neuron.2016.07.047, PMID: 27568518
Anderson MI, Jeffery KJ. 2003. Heterogeneous modulation of place cell firing by changes in context. The Journalof Neuroscience 23:8827–8835. DOI: https://doi.org/10.1523/JNEUROSCI.23-26-08827.2003, PMID: 14523083
Angulo-Garcia D, Ferraris M, Ghestem A. 2018. Spatio-temporal organization of cell assemblies in nucleus
reuniens during slow oscillations. bioRxiv. DOI: https://doi.org/10.1101/474973
Azizi AH, Wiskott L, Cheng S. 2013. A computational model for preplay in the Hippocampus. Frontiers in
Computational Neuroscience 7:161. DOI: https://doi.org/10.3389/fncom.2013.00161, PMID: 24282402
Bayati M, Valizadeh A, Abbassian A, Cheng S. 2015. Self-organization of synchronous activity propagation inneuronal networks driven by local excitation. Frontiers in Computational Neuroscience 9:e00069. DOI: https://doi.org/10.3389/fncom.2015.00069
Belchior H, Lopes-Dos-Santos V, Tort AB, Ribeiro S. 2014. Increase in hippocampal theta oscillations during
spatial decision making. Hippocampus 24:693–702. DOI: https://doi.org/10.1002/hipo.22260, PMID: 24520011
Bonnevie T, Dunn B, Fyhn M, Hafting T, Derdikman D, Kubie JL, Roudi Y, Moser EI, Moser MB. 2013. Grid cells
require excitatory drive from the Hippocampus. Nature Neuroscience 16:309–317. DOI: https://doi.org/10.
1038/nn.3311, PMID: 23334581
Brandon MP, Bogaard AR, Libby CP, Connerney MA, Gupta K, Hasselmo ME. 2011. Reduction of theta rhythm
dissociates grid cell spatial periodicity from directional tuning. Science 332:595–599. DOI: https://doi.org/10.
1126/science.1201652, PMID: 21527714
Burak Y, Fiete IR. 2009. Accurate path integration in continuous attractor network models of grid cells. PLOS
Computational Biology 5:e1000291. DOI: https://doi.org/10.1371/journal.pcbi.1000291, PMID: 19229307

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

20 of 34

Research article

Neuroscience Physics of Living Systems
Burgess N, Barry C, O’Keefe J. 2007. An oscillatory interference model of grid cell firing. Hippocampus 17:801–
812. DOI: https://doi.org/10.1002/hipo.20327, PMID: 17598147
Bush D, Burgess N. 2014. A hybrid oscillatory interference/continuous attractor network model of grid cell firing.
The Journal of Neuroscience 34:5065–5079. DOI: https://doi.org/10.1523/JNEUROSCI.4017-13.2014,
PMID: 24695724
Buzsáki G. 1986. Hippocampal sharp waves: their origin and significance. Brain Research 398:242–252.
DOI: https://doi.org/10.1016/0006-8993(86)91483-6, PMID: 3026567
Buzsáki G. 1996. The hippocampo-neocortical dialogue. Cerebral Cortex 6:81–92. DOI: https://doi.org/10.1093/
cercor/6.2.81, PMID: 8670641
Buzsáki G. 2002. Theta oscillations in the Hippocampus. Neuron 33:325–340. DOI: https://doi.org/10.1016/
S0896-6273(02)00586-X, PMID: 11832222
Buzsáki G, Anastassiou CA, Koch C. 2012. The origin of extracellular fields and currents–EEG, ECoG, LFP and
spikes. Nature Reviews Neuroscience 13:407–420. DOI: https://doi.org/10.1038/nrn3241, PMID: 22595786
Chance FS. 2012. Hippocampal phase precession from dual input components. Journal of Neuroscience 32:
16693–16703. DOI: https://doi.org/10.1523/JNEUROSCI.2786-12.2012, PMID: 23175823
Cheng S, Frank LM. 2011. The structure of networks that produce the transformation from grid cells to place
cells. Neuroscience 197:293–306. DOI: https://doi.org/10.1016/j.neuroscience.2011.09.002, PMID: 21963867
Chenkov N, Sprekeler H, Kempter R. 2017. Memory replay in balanced recurrent networks. PLOS Computational
Biology 13:e1005359. DOI: https://doi.org/10.1371/journal.pcbi.1005359, PMID: 28135266
Chrobak JJ, Buzsáki G. 1996. High-frequency oscillations in the output networks of the hippocampal-entorhinal
Axis of the freely behaving rat. The Journal of Neuroscience 16:3056–3066. DOI: https://doi.org/10.1523/
JNEUROSCI.16-09-03056.1996, PMID: 8622135
Chrobak JJ, Buzsáki G. 1998. Gamma oscillations in the entorhinal cortex of the freely behaving rat. The Journalof Neuroscience 18:388–398. DOI: https://doi.org/10.1523/JNEUROSCI.18-01-00388.1998, PMID: 9412515
Climer JR, Newman EL, Hasselmo ME. 2013. Phase coding by grid cells in unconstrained environments: twodimensional phase precession. European Journal of Neuroscience 38:2526–2541. DOI: https://doi.org/10.1111/
ejn.12256, PMID: 23718553
Colgin LL, Denninger T, Fyhn M, Hafting T, Bonnevie T, Jensen O, Moser MB, Moser EI. 2009. Frequency of
gamma oscillations routes flow of information in the Hippocampus. Nature 462:353–357. DOI: https://doi.org/
10.1038/nature08573, PMID: 19924214
Couey JJ, Witoelar A, Zhang SJ, Zheng K, Ye J, Dunn B, Czajkowski R, Moser MB, Moser EI, Roudi Y, Witter MP.
2013. Recurrent inhibitory circuitry as a mechanism for grid formation. Nature Neuroscience 16:318–324.
DOI: https://doi.org/10.1038/nn.3310, PMID: 23334580
Davidson TJ, Kloosterman F, Wilson MA. 2009. Hippocampal replay of extended experience. Neuron 63:497–
507. DOI: https://doi.org/10.1016/j.neuron.2009.07.027, PMID: 19709631
Diba K, Buzsáki G. 2007. Forward and reverse hippocampal place-cell sequences during ripples. Nature
Neuroscience 10:1241–1242. DOI: https://doi.org/10.1038/nn1961, PMID: 17828259
Diehl GW, Hon OJ, Leutgeb S, Leutgeb JK. 2019. Stability of medial entorhinal cortex representations over time.
Hippocampus 29:284–302. DOI: https://doi.org/10.1002/hipo.23017, PMID: 30175425
Dragoi G, Tonegawa S. 2011. Preplay of future place cell sequences by hippocampal cellular assemblies. Nature
469:397–401. DOI: https://doi.org/10.1038/nature09633, PMID: 21179088
Drieu C, Todorova R, Zugaro M. 2018. Nested sequences of hippocampal assemblies during behavior support
subsequent sleep replay. Science 362:675–679. DOI: https://doi.org/10.1126/science.aat2952, PMID: 30409880
Dunn B, Mørreaunet M, Roudi Y. 2015. Correlations and functional connections in a population of grid cells.
PLOS Computational Biology 11:e1004052. DOI: https://doi.org/10.1371/journal.pcbi.1004052, PMID: 2571490
8Ebbesen CL, Reifenstein ET, Tang Q, Burgalossi A, Ray S, Schreiber S, Kempter R, Brecht M. 2016. Cell TypeSpecific differences in spike timing and spike shape in the rat parasubiculum and superficial medial entorhinal
cortex. Cell Reports 16:1005–1015. DOI: https://doi.org/10.1016/j.celrep.2016.06.057, PMID: 27425616
English G. 2017. Spiking Neural Network Models for the Emergence of Patterned Activity in Grid Cell
Populations. Switzerland: ETH Zurich and University of Zurich.
Evans T, Bicanski A, Bush D, Burgess N. 2016. How environment and self-motion combine in neuralrepresentations of space. The Journal of Physiology 594:6535–6546. DOI: https://doi.org/10.1113/JP270666,
PMID: 26607203
Feng T, Silva D, Foster DJ. 2015. Dissociation between the experience-dependent development of hippocampal
theta sequences and single-trial phase precession. Journal of Neuroscience 35:4890–4902. DOI: https://doi.
org/10.1523/JNEUROSCI.2614-14.2015, PMID: 25810520
Foster DJ. 2017. Replay comes of age. Annual Review of Neuroscience 40:581–602. DOI: https://doi.org/10.
1146/annurev-neuro-072116-031538
Foster DJ, Knierim JJ. 2012. Sequence learning and the role of the Hippocampus in rodent navigation. Current
Opinion in Neurobiology 22:294–300. DOI: https://doi.org/10.1016/j.conb.2011.12.005, PMID: 22226994
Foster DJ, Wilson MA. 2006. Reverse replay of behavioural sequences in hippocampal place cells during theawake state. Nature 440:680–683. DOI: https://doi.org/10.1038/nature04587, PMID: 16474382
Fuchs EC, Neitz A, Pinna R, Melzer S, Caputi A, Monyer H. 2016. Local and distant input controlling excitation inlayer II of the medial entorhinal cortex. Neuron 89:194–208. DOI: https://doi.org/10.1016/j.neuron.2015.11.
029, PMID: 26711115

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

21 of 34

Research article

Neuroscience Physics of Living Systems
Gönner L, Vitay J, Hamker FH. 2017. Predictive Place-Cell sequences for Goal-Finding emerge from goal
memory and the cognitive map: a computational model. Frontiers in Computational Neuroscience 11:84.
DOI: https://doi.org/10.3389/fncom.2017.00084, PMID: 29075187
Gonzalez-Sulser A, Parthier D, Candela A, McClure C, Pastoll H, Garden D, Sürmeli G, Nolan MF. 2014.
GABAergic projections from the medial septum selectively inhibit interneurons in the medial entorhinal cortex.
Journal of Neuroscience 34:16739–16743. DOI: https://doi.org/10.1523/JNEUROSCI.1612-14.2014,
PMID: 25505326
Grosmark AD, Buzsáki G. 2016. Diversity in neural firing dynamics supports both rigid and learned hippocampal
sequences. Science 351:1440–1443. DOI: https://doi.org/10.1126/science.aad1935, PMID: 27013730
Gu Y, Lewallen S, Kinkhabwala AA, Domnisoru C, Yoon K, Gauthier JL, Fiete IR, Tank DW. 2018. A Map-like
Micro-Organization of grid cells in the medial entorhinal cortex. Cell 175:736–750. DOI: https://doi.org/10.
1016/j.cell.2018.08.066, PMID: 30270041
Hafting T, Fyhn M, Molden S, Moser MB, Moser EI. 2005. Microstructure of a spatial map in the entorhinal
cortex. Nature 436:801–806. DOI: https://doi.org/10.1038/nature03721, PMID: 15965463
Hafting T, Fyhn M, Bonnevie T, Moser MB, Moser EI. 2008. Hippocampus-independent phase precession inentorhinal grid cells. Nature 453:1248–1252. DOI: https://doi.org/10.1038/nature06957, PMID: 18480753
Haga T, Fukai T. 2018. Recurrent network model for learning goal-directed sequences through reverse replay.
eLife 7:e34171. DOI: https://doi.org/10.7554/eLife.34171, PMID: 29969098
Hahn TT, McFarland JM, Berberich S, Sakmann B, Mehta MR. 2012. Spontaneous persistent activity in entorhinal
cortex modulates cortico-hippocampal interaction in vivo. Nature Neuroscience 15:1531–1538. DOI: https://doi.org/10.1038/nn.3236, PMID: 23042081
Hales JB, Schlesiger MI, Leutgeb JK, Squire LR, Leutgeb S, Clark RE. 2014. Medial entorhinal cortex lesions only
partially disrupt hippocampal place cells and hippocampus-dependent place memory. Cell Reports 9:893–901.
DOI: https://doi.org/10.1016/j.celrep.2014.10.009, PMID: 25437546
Hardcastle K, Ganguli S, Giocomo LM. 2015. Environmental boundaries as an error correction mechanism forgrid cells. Neuron 86:827–839. DOI: https://doi.org/10.1016/j.neuron.2015.03.039, PMID: 25892299
Hasselmo ME, Giocomo LM, Zilli EA. 2007. Grid cell firing may arise from interference of theta frequency
membrane potential oscillations in single neurons. Hippocampus 17:1252–1271. DOI: https://doi.org/10.1002/
hipo.20374, PMID: 17924530
Heys JG, Rangarajan KV, Dombeck DA. 2014. The functional micro-organization of grid cells revealed by cellularresolution imaging. Neuron 84:1079–1090. DOI: https://doi.org/10.1016/j.neuron.2014.10.048, PMID: 25467
986
Hollup SA, Molden S, Donnett JG, Moser MB, Moser EI. 2001. Accumulation of hippocampal place fields at thegoal location in an annular watermaze task. The Journal of Neuroscience 21:1635–1644. DOI: https://doi.org/
10.1523/JNEUROSCI.21-05-01635.2001, PMID: 11222654
Jeewajee A, Barry C, Douchamps V, Manson D, Lever C, Burgess N. 2014. Theta phase precession of grid andplace cell firing in open environments. Philosophical Transactions of the Royal Society B: Biological Sciences
369:20120532. DOI: https://doi.org/10.1098/rstb.2012.0532
Johnson A, Redish AD. 2007. Neural ensembles in CA3 transiently encode paths forward of the animal at adecision point. Journal of Neuroscience 27:12176–12189. DOI: https://doi.org/10.1523/JNEUROSCI.3761-07.
2007, PMID: 17989284
Kang L, Balasubramanian V. 2019. A geometric attractor mechanism for self-organization of entorhinal grid
modules. eLife 8:e46687. DOI: https://doi.org/10.7554/eLife.46687, PMID: 31373556
Keinath AT, Epstein RA, Balasubramanian V. 2018. Environmental deformations dynamically shift the grid cell
spatial metric. eLife 7:e38169. DOI: https://doi.org/10.7554/eLife.38169, PMID: 30346272
Kempter R, Leibold C, Buzsáki G, Diba K, Schmidt R. 2012. Quantifying circular-linear associations: hippocampal
phase precession. Journal of Neuroscience Methods 207:113–124. DOI: https://doi.org/10.1016/j.jneumeth.
2012.03.007, PMID: 22487609
Kerr KM, Agster KL, Furtak SC, Burwell RD. 2007. Functional neuroanatomy of the parahippocampal region: thelateral and medial entorhinal Areas. Hippocampus 17:697–708. DOI: https://doi.org/10.1002/hipo.20315,
PMID: 17607757
Koenig J, Linder AN, Leutgeb JK, Leutgeb S. 2011. The spatial periodicity of grid cells is not sustained during
reduced theta oscillations. Science 332:592–595. DOI: https://doi.org/10.1126/science.1201685,
PMID: 21527713
Lee AK, Wilson MA. 2002. Memory of sequential experience in the Hippocampus during slow wave sleep.
Neuron 36:1183–1194. DOI: https://doi.org/10.1016/S0896-6273(02)01096-6, PMID: 12495631
Leutgeb JK, Leutgeb S, Moser MB, Moser EI. 2007. Pattern separation in the dentate gyrus and CA3 of the
Hippocampus. Science 315:961–966. DOI: https://doi.org/10.1126/science.1135801, PMID: 17303747
Liu K, Sibille J, Dragoi G. 2018. Generative predictive codes by multiplexed hippocampal neuronal tuplets.
Neuron 99:1329–1341. DOI: https://doi.org/10.1016/j.neuron.2018.07.047, PMID: 30146305
Liu K, Sibille J, Dragoi G. 2019. Preconfigured patterns are the primary driver of offline multi-neuronal sequence
replay. Hippocampus 29:275–283. DOI: https://doi.org/10.1002/hipo.23034, PMID: 30260526
Mamad O, Stumpp L, McNamara HM, Ramakrishnan C, Deisseroth K, Reilly RB, Tsanov M. 2017. Place field
assembly distribution encodes preferred locations. PLOS Biology 15:e2002365. DOI: https://doi.org/10.1371/
journal.pbio.2002365, PMID: 28898248

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

22 of 34

Research article

Neuroscience Physics of Living Systems
Middleton SJ, Kneller EM, Chen S, Ogiwara I, Montal M, Yamakawa K, McHugh TJ. 2018. Altered hippocampal
replay is associated with memory impairment in mice heterozygous for the Scn2a gene. Nature Neuroscience
21:996–1003. DOI: https://doi.org/10.1038/s41593-018-0163-8, PMID: 29867081
Navratilova Z, Giocomo LM, Fellous JM, Hasselmo ME, McNaughton BL. 2012. Phase precession and variablespatial scaling in a periodic attractor map model of medial entorhinal grid cells with realistic after-spike
dynamics. Hippocampus 22:772–789. DOI: https://doi.org/10.1002/hipo.20939, PMID: 21484936
Nolan M. 2018. A Model for Grid Firing and Theta-Nested Gamma Oscillations in Layer 2 of the Medial
EntorhinalCortex. In: Hippocampal Microcircuits. Springer. p. 567–584. DOI: https://doi.org/10.1007/978-3319-99103-0_15
O’Keefe J. 1976. Place units in the Hippocampus of the freely moving rat. Experimental Neurology 51:78–109.
DOI: https://doi.org/10.1016/0014-4886(76)90055-8, PMID: 1261644
O’Keefe J, Dostrovsky J. 1971. The Hippocampus as a spatial map. Preliminary evidence from unit activity in thefreely-moving rat. Brain Research 34:171–175. DOI: https://doi.org/10.1016/0006-8993(71)90358-1,
PMID: 5124915
O’Keefe J, Recce ML. 1993. Phase relationship between hippocampal place units and the EEG theta rhythm.
Hippocampus 3:317–330. DOI: https://doi.org/10.1002/hipo.450030307, PMID: 8353611
O’Neill J, Boccara CN, Stella F, Schoenenberger P, Csicsvari J. 2017. Superficial layers of the medial entorhinal
cortex replay independently of the Hippocampus. Science 355:184–188. DOI: https://doi.org/10.1126/science.
aag2787, PMID: 28082591
Ocko SA, Hardcastle K, Giocomo LM, Ganguli S. 2018. Emergent elasticity in the neural code for space. PNAS
1154:e11798. DOI: https://doi.org/10.1073/pnas.1805959115
Ólafsdóttir HF, Carpenter F, Barry C. 2016. Coordinated grid and place cell replay during rest. Nature
Neuroscience 19:792–794. DOI: https://doi.org/10.1038/nn.4291, PMID: 27089021
Ólafsdóttir HF, Carpenter F, Barry C. 2017. Task demands predict a dynamic switch in the content of awake
hippocampal replay. Neuron 96:925–935. DOI: https://doi.org/10.1016/j.neuron.2017.09.035, PMID: 29056296
Ólafsdóttir HF, Bush D, Barry C. 2018. The role of hippocampal replay in memory and Planning. Current Biology
28:R37–R50. DOI: https://doi.org/10.1016/j.cub.2017.10.073
Oliva A, Fernández-Ruiz A, Buzsáki G, Berényi A. 2016. Role of hippocampal CA2 region in triggering SharpWave ripples. Neuron 91:1342–1355. DOI: https://doi.org/10.1016/j.neuron.2016.08.008, PMID: 27593179
Palmer J, Keane A, Gong P. 2017. Learning and executing goal-directed choices by internally generated
sequences in spiking neural circuits. PLOS Computational Biology 13:e1005669. DOI: https://doi.org/10.1371/
journal.pcbi.1005669, PMID: 28759562
Palmer JHC, Gong P. 2013. Formation and regulation of dynamic patterns in Two-Dimensional spiking neuralcircuits with Spike-Timing-Dependent plasticity. Neural Computation 25:2833–2857. DOI: https://doi.org/10.
1162/NECO_a_00511
Pfeiffer BE, Foster DJ. 2013. Hippocampal place-cell sequences depict future paths to remembered goals.
Nature 497:74–79. DOI: https://doi.org/10.1038/nature12112, PMID: 23594744
Pfeiffer BE, Foster DJ. 2015. PLACE CELLS. Autoassociative dynamics in the generation of sequences ofhippocampal place cells. Science 349:180–183. DOI: https://doi.org/10.1126/science.aaa9633, PMID: 26160946
Ponulak F, Hopfield JJ. 2013. Rapid, parallel path planning by propagating wavefronts of spiking neural activity.
Frontiers in Computational Neuroscience 7:98. DOI: https://doi.org/10.3389/fncom.2013.00098, PMID: 23882213Ramirez-Villegas JF, Willeke KF, Logothetis NK, Besserve M. 2018. Dissecting the synapse- and FrequencyDependent network mechanisms of in Vivo Hippocampal Sharp Wave-Ripples. Neuron 100:1224–1240.
DOI: https://doi.org/10.1016/j.neuron.2018.09.041, PMID: 30482688
Raudies F, Hasselmo ME. 2015. Differences in Visual-Spatial input may underlie different compression propertiesof firing fields for grid cell modules in medial entorhinal cortex. PLOS Computational Biology 11:e1004596.
DOI: https://doi.org/10.1371/journal.pcbi.1004596, PMID: 26584432
Reifenstein E, Stemmler M, Herz AV, Kempter R, Schreiber S. 2014. Movement dependence and layer specificity
of entorhinal phase precession in two-dimensional environments. PLOS ONE 9:e100638. DOI: https://doi.org/
10.1371/journal.pone.0100638, PMID: 24959748
Reifenstein ET, Ebbesen CL, Tang Q, Brecht M, Schreiber S, Kempter R. 2016. Cell-Type specific phase
precession in layer II of the medial entorhinal cortex. The Journal of Neuroscience 36:2283–2288. DOI: https://doi.org/10.1523/JNEUROSCI.2986-15.2016, PMID: 26888937
Sargolini F, Fyhn M, Hafting T, McNaughton BL, Witter MP, Moser MB, Moser EI. 2006. Conjunctive
representation of position, direction, and velocity in entorhinal cortex. Science 312:758–762. DOI: https://doi.
org/10.1126/science.1125572, PMID: 16675704
Sasaki T, Piatti VC, Hwaun E, Ahmadi S, Lisman JE, Leutgeb S, Leutgeb JK. 2018. Dentate network activity isnecessary for spatial working memory by supporting CA3 sharp-wave ripple generation and prospective firingof CA3 neurons. Nature Neuroscience 21:258–269. DOI: https://doi.org/10.1038/s41593-017-0061-5, PMID: 29335604
Savelli F, Luck JD, Knierim JJ. 2017. Framing of grid cells within and beyond navigation boundaries. eLife 6:
e21354. DOI: https://doi.org/10.7554/eLife.21354, PMID: 28084992
Schlesiger MI, Cannova CC, Boublil BL, Hales JB, Mankin EA, Brandon MP, Leutgeb JK, Leibold C, Leutgeb S.
2015. The medial entorhinal cortex is necessary for temporal organization of hippocampal neuronal activity.
Nature Neuroscience 18:1123–1132. DOI: https://doi.org/10.1038/nn.4056, PMID: 26120964

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

23 of 34

Research article

Neuroscience Physics of Living Systems
Schlingloff D, Káli S, Freund TF, Hájos N, Gulyás AI. 2014. Mechanisms of sharp wave initiation and ripple
generation. Journal of Neuroscience 34:11385–11398. DOI: https://doi.org/10.1523/JNEUROSCI.0867-14.
2014, PMID: 25143618
Schmidt H, Gour A, Straehle J, Boergens KM, Brecht M, Helmstaedter M. 2017. Axonal synapse sorting in medial
entorhinal cortex. Nature 549:469–475. DOI: https://doi.org/10.1038/nature24005, PMID: 28959971
Silva D, Feng T, Foster DJ. 2015. Trajectory events across hippocampal place cells require previous experience.
Nature Neuroscience 18:1772–1779. DOI: https://doi.org/10.1038/nn.4151, PMID: 26502260
Simon AP, Poindessous-Jazat F, Dutar P, Epelbaum J, Bassant MH. 2006. Firing properties of anatomically
identified neurons in the medial septum of anesthetized and unanesthetized restrained rats. Journal of
Neuroscience 26:9038–9046. DOI: https://doi.org/10.1523/JNEUROSCI.1401-06.2006, PMID: 16943562
Singer AC, Carr MF, Karlsson MP, Frank LM. 2013. Hippocampal SWR activity predicts correct decisions during
the initial learning of an alternation task. Neuron 77:1163–1173. DOI: https://doi.org/10.1016/j.neuron.2013.01.
027, PMID: 23522050
Skaggs WE, McNaughton BL, Wilson MA, Barnes CA. 1996. Theta phase precession in hippocampal neuronal
populations and the compression of temporal sequences. Hippocampus 6:149–172. DOI: https://doi.org/10.
1002/(SICI)1098-1063(1996)6:2<149::AID-HIPO6>3.0.CO;2-K, PMID: 8797016
Stark E, Roux L, Eichler R, Senzai Y, Royer S, Buzsáki G. 2014. Pyramidal cell-interneuron interactions underlie
hippocampal ripple oscillations. Neuron 83:467–480. DOI: https://doi.org/10.1016/j.neuron.2014.06.023,
PMID: 25033186
Stensola H, Stensola T, Solstad T, Frøland K, Moser MB, Moser EI. 2012. The entorhinal grid map is discretized.
Nature 492:72–78. DOI: https://doi.org/10.1038/nature11649, PMID: 23222610
Steriade M, Timofeev I, Grenier F. 2001. Natural waking and sleep states: a view from inside neocortical neurons.
Journal of Neurophysiology 85:1969–1985. DOI: https://doi.org/10.1152/jn.2001.85.5.1969, PMID: 11353014
Thurley K, Leibold C, Gundlfinger A, Schmitz D, Kempter R. 2008. Phase precession through synaptic facilitation.
Neural Computation 20:1285–1324. DOI: https://doi.org/10.1162/neco.2008.07-06-292, PMID: 18085985
Thurley K, Hellmundt F, Leibold C. 2013. Phase precession of grid cells in a network model without external
pacemaker. Hippocampus 23:786–796. DOI: https://doi.org/10.1002/hipo.22133, PMID: 23576429
Trimper JB, Trettel SG, Hwaun E, Colgin LL. 2017. Methodological caveats in the detection of coordinatedreplay between place cells and grid cells. Frontiers in Systems Neuroscience 11:57. DOI: https://doi.org/10.
3389/fnsys.2017.00057, PMID: 28824388
Tsodyks MV, Skaggs WE, Sejnowski TJ, McNaughton BL. 1996. Population dynamics and theta rhythm phase
precession of hippocampal place cell firing: a spiking neuron model. Hippocampus 6:271–280. DOI: https://doi.
org/10.1002/(SICI)1098-1063(1996)6:3<271::AID-HIPO5>3.0.CO;2-Q, PMID: 8841826
Unal G, Joshi A, Viney TJ, Kis V, Somogyi P. 2015. Synaptic targets of medial septal projections in the
Hippocampus and extrahippocampal cortices of the mouse. Journal of Neuroscience 35:15812–15826.
DOI: https://doi.org/10.1523/JNEUROSCI.2639-15.2015, PMID: 26631464
Urdapilleta E, Si B, Treves A. 2017. Selforganization of modular activity of grid cells. Hippocampus 27:1204–
1213. DOI: https://doi.org/10.1002/hipo.22765, PMID: 28768062
Vanderwolf CH. 1969. Hippocampal electrical activity and voluntary movement in the rat.
Electroencephalography and Clinical Neurophysiology 26:407–418. DOI: https://doi.org/10.1016/0013-4694(69)90092-3
Vinogradova OS. 1995. Expression, control, and probable functional significance of the neuronal theta-rhythm.
Progress in Neurobiology 45:523–583. DOI: https://doi.org/10.1016/0301-0082(94)00051-I, PMID: 7624485
Widloski J, Fiete IR. 2014. A model of grid cell development through spatial exploration and spike timedependent plasticity. Neuron 83:481–495. DOI: https://doi.org/10.1016/j.neuron.2014.06.018, PMID: 25033187
Wikenheiser AM, Redish AD. 2015a. Hippocampal theta sequences reflect current goals. Nature Neuroscience
18:289–294. DOI: https://doi.org/10.1038/nn.3909, PMID: 25559082
Wikenheiser AM, Redish AD. 2015b. Decoding the cognitive map: ensemble hippocampal sequences anddecision making. Current Opinion in Neurobiology 32:8–15. DOI: https://doi.org/10.1016/j.conb.2014.10.002,
PMID: 25463559
Wilson MA, McNaughton BL. 1993. Dynamics of the hippocampal ensemble code for space. Science 261:1055–
1058. DOI: https://doi.org/10.1126/science.8351520, PMID: 8351520
Winterer J, Maier N, Wozny C, Beed P, Breustedt J, Evangelista R, Peng Y, D’Albis T, Kempter R, Schmitz D.
2017. Excitatory microcircuits within superficial layers of the medial entorhinal cortex. Cell Reports 19:1110–
1116. DOI: https://doi.org/10.1016/j.celrep.2017.04.041, PMID: 28494861
Wu C-T, Haggerty D, Kemere C, Ji D. 2017. Hippocampal awake replay in fear memory retrieval. Nature
Neuroscience 20:571–580. DOI: https://doi.org/10.1038/nn.4507
Wu X, Foster DJ. 2014. Hippocampal replay captures the unique topological structure of a novel environment.
Journal of Neuroscience 34:6459–6469. DOI: https://doi.org/10.1523/JNEUROSCI.3414-13.2014
Xu H, Baracskay P, O’Neill J, Csicsvari J. 2019. Assembly responses of hippocampal CA1 place cells predictlearned behavior in Goal-Directed spatial tasks on the radial Eight-Arm maze. Neuron 101:119–132.
DOI: https://doi.org/10.1016/j.neuron.2018.11.015, PMID: 30503645
Yamaguchi Y, Aota Y, McNaughton BL, Lipa P. 2002. Bimodality of theta phase precession in hippocampal place
cells in freely running rats. Journal of Neurophysiology 87:2629–2642. DOI: https://doi.org/10.1152/jn.2002.87.
6.2629, PMID: 12037166

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

24 of 34

Research article

Neuroscience Physics of Living Systems
Yamamoto J, Tonegawa S. 2017. Direct medial entorhinal cortex input to hippocampal CA1 is crucial for
extended quiet awake replay. Neuron 96:217–227. DOI: https://doi.org/10.1016/j.neuron.2017.09.017, PMID: 28957670
Yoon K, Buice MA, Barry C, Hayman R, Burgess N, Fiete IR. 2013. Specific evidence of low-dimensional
continuous attractor dynamics in grid cells. Nature Neuroscience 16:1077–1084. DOI: https://doi.org/10.1038/
nn.3450, PMID: 23852111
Yoshida M, Jochems A, Hasselmo ME. 2013. Comparison of properties of medial entorhinal cortex layer II
neurons in two anatomical dimensions with and without cholinergic activation. PLOS ONE 8:e73904.
DOI: https://doi.org/10.1371/journal.pone.0073904, PMID: 24069244
Zhang K, Ginzburg I, McNaughton BL, Sejnowski TJ. 1998. Interpreting neuronal population activity byreconstruction: unified framework with application to hippocampal place cells. Journal of Neurophysiology 79:
1017–1044. DOI: https://doi.org/10.1152/jn.1998.79.2.1017, PMID: 9463459
Zheng C, Bieri KW, Hsiao YT, Colgin LL. 2016. Spatial sequence coding differs during slow and fast gamma
rhythms in the Hippocampus. Neuron 89:398–408. DOI: https://doi.org/10.1016/j.neuron.2015.12.005,
PMID: 26774162
Zielinski MC, Tang W, Jadhav SP. 2018. The role of replay and theta sequences in mediating hippocampalprefrontal interactions for memory and cognition. Hippocampus 91:1–13. DOI: https://doi.org/10.1002/hipo.
22821
Zutshi I, Fu ML, Lilascharoen V, Leutgeb JK, Lim BK, Leutgeb S. 2018. Recurrent circuits within medial entorhinal
cortex superficial layers support grid cell firing. Nature Communications 9:3701. DOI: https://doi.org/10.1038/
s41467-018-06104-5, PMID: 30209250

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

25 of 34

Research article

Neuroscience Physics of Living Systems

Appendix 1
Additional methods
Model details
Setup
We initialize each neuron on the neural sheet to a random membrane potential between 0 and
1. We evolve the network according to Equation 2 and Equation 3 for 500 timesteps without
velocity input to generate a rough grid-like pattern of activity. To make the grid more regularand to anneal grid defects, we perform three evolutions of timesteps with constant velocity
input V of magnitude 0.5 m/s and angles, p/2 – p/5, 2p/5, and p/4 successively. We thenevolve the network for 4 laps without idle periods.
Finally, we learn allocentric inputs as follows. We simulate 2 additional laps with 1000
timesteps of idle periods at each end of this track. During this time, both the excitatory andinhibitory drives are maintained at their run values without theta oscillations. Over the final 500
timesteps at each end, the number of spikes produced at each neural sheet location acrossthe 4 excitatory populations are tallied. These values are linearly rescaled to serve as
allocentric inputs to excitatory neurons, with the highest spike count corresponding toþaþ
max and anything below the 5th percentile of nonzero activity corresponding to amin . Then themain simulation phase begins.

Simulations with periodic boundary conditions
Since annealing grid defects is difficult with periodic boundary conditions, we initialize
membrane potentials with a grid-like configuration on the neural sheet as in Burak and Fiete(2009). The grid scale is chosen to be similar to that self-organized without periodic boundary
conditions. The rest of simulation setup proceeds similarly as above, except the three
evolutions with constant velocity input last only 500 timesteps each.
Parameter values that differ from those in Table 1 are n = 192, aþ
mag ¼ 1:2 during runs and0.92 when idle for the uniform excitatory drive, aþ
max ¼ 1:4 during allocentriccorrections, amag ¼ 0:75 during runs and –0.1 when idle, and ath ¼ 0:15 during runs. Results areshown in Figure 8F, G.

Trajectory
For 2D open field simulations, we use velocity inputs extracted by Burak and Fiete (2009)
from a real animal trajectory reported in Hafting et al. (2005).
For 1D track simulations, we generate trajectories as follows. Each trajectory consists of 36
runs alternating in direction, each lasting 1.5 s. Each run is followed by an idle period of
1.5 s and an allocentric correction period of 0.5 s. The velocity angle is p/5, chosen to avoid
directions of symmetry in the neural sheet.
In each run, the speed linearly increases from 0 m/s to 0.5 m/s over the first 0.3 s, remains
at 0.5 m/s over the next 0.9 s, and linearly decreases from 0.5 m/s to 0 m/s over the last 0.3 s.
To the middle 0.9 s of constant speed, we add noise. The noise is generated by a fractional
Brownian motion process with drift 0, volatility 1, and Hurst index 0.8. We then add an overall
constant at each timestep to make the noise sum to 0 (and thus, to leading order, contribute 0
additional displacement) and introduce an overall scale factor to make its maximum
magnitude 0.1 m/s. We add this processed noise to the constant velocity of 0.5 m/s. The initial
theta phase of each run is chosen randomly, and the phase increases during runs at a constantrate to span 360˚ over 125 ms.
Excitatory and inhibitory drives change between runs and idle periods according toTable 1. They do so by linearly interpolating between corresponding values over the first
0.3 s of each idle period.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

26 of 34

Research article

Neuroscience Physics of Living Systems

Single neuron recording and position decoding
Selection of neurons for recording
As described in Materials and methods, we select up to 150 neurons total for recording fromfour equally spaced circular areas on the neural sheet (up to 200 neurons total for simulations
with periodic boundary conditions). For the main simulation, these areas are centered 95
neurons from the center of the sheet; for the simulations with lower velocity gain in Figure 8AC, the distance is 75 neurons; for the simulations with periodic boundary conditions inFigure 8D-G, the distance is 60 neurons. In each area, we randomly choose 30–50 neurons
from each excitatory population and compute their firing fields for each run as follows. We
partition space into 5 cm bins and compute the firing rate of each bin as the number of spikes
divided by the occupancy time; these raw fields are smoothed by a Gaussian filter of width 4bins. Firing fields for leftward and rightward runs are considered separately. In at least one
direction, recorded neurons must have maximum firing rate above 0.5 Hz for all runs, and thenormalized correlation, defined as u  v=jujjvj for two vectors u and v, must be above 0.6 forall pairs of runs.

Firing fields
For 2D open field simulations, we partition space into bins indexed by R = (X, Y) and computethe firing rate of each bin as the number of spikes divided by the occupancy time. These raw
fields are smoothed by a Gaussian filter of width 5 bins to produce the field map F(R). We
define the normalized autocorrelation map asX1NðRÞ
FðR0 ÞFðR0 þ RÞCðRÞ ¼

R0

1Nð0Þ

X
R0

FðR0 ÞFðR0 Þ(8);where N(R) is the number of pairs of positions separated by R. Autocorrelation maps areshown in Figure 1E.
For 1D track simulations, we partition space into 2 cm bins indexed by X and compute thefiring rate of each bin as the number of spikes divided by the occupancy time over all laps.
These raw fields are smoothed by a Gaussian filter of width 1 bin to produce the field map F(X). Field maps are shown in Figure 2D.
To categorize a simulation as exhibiting ’one field’ or ’two fields,’ we first order all neurons
by the track position at which they are most active. We exclude the middle 40% of all neurons,since we do not expect neurons with central field positions to exhibit two fields. We then shifteach field such that its position of peak activity is at 0 cm and calculate the fraction of field
activity located more than 18 cm away. If this fraction is greater than 0.1, we consider thesimulation as exhibiting two fields. The center of mass of this distant activity, using theabsolute value of position, is the distance between fields for the simulation.

Position decoding
We use the standard Bayesian techniques to decode position from the population activity ofgrid cells (Zhang et al., 1998; O’Neill et al., 2017). Let sk(t) count the number of spikes firedby the recorded neuron k within a time window of DT around t. Each neuron has a field map
Fk(X) as described above. For every time t, we wish to calculate the probability that the animal
is at position X given the population activity s(t). We do so through Bayes’s rule, assuming aflat prior and a spiking likelihood that obeys independent Poisson statistics:
!
!
Y
Y
X
sk ðtÞ
pðXjsðtÞÞ / pðsðtÞjXÞ ¼
pðsk ðtÞjXÞ /Fk ðXÞ
exp
DT
Fk ðXÞ :
(9)
k

k

k

We increment time in steps of 5 ms. For theta-related analysis, we use a sliding decoding
window of width DT = 20 ms. For analysis of replays, whose timescales are shorter and spikes
are more concentrated, we use DT = 10 ms and require each window to contain at least two
spikes.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

27 of 34

Research article

Neuroscience Physics of Living Systems

Theta-related analysis
Multiunit activity
To calculate the multiunit activity averaged over theta cycles, we sum spike trains overrecorded neurons during each run, extract the multiunit activity of each complete theta cycle,
and average these extracts. This theta-modulated multiunit activity is shown in Figure 2E.

Phase precession
We use circular-linear regression and correlation as described in Kempter et al. (2012) toquantify the relationship between theta phase and position within a grid field. For eachneuron, we combine spikes across all rightward runs and across all leftward runs, considering
each direction separately. For each direction, we identify a separation between grid fields as a
10 cm span without spikes. For each neuron and each direction, we only consider the mostcentral field. We exclude fields that spike within 3 cm from the track ends to avoid edge
effects, and also exclude fields containing fewer than 30 spikes or spanning less than
12 cm wide.
For each field, we record the position Xj and theta phase j at which each of its spikes
j occurred. We hypothesize that the phase of each spike is linearly influenced by its position:
j ¼ 2pqXj þ

0;

(10)

where q is the phase precession slope and 0 is the phase offset. We wish to fit theseparameters by minimizing the total circular distance between true phases and corresponding
predicted phases. The circular distance between two phases and  is defined to be 2[1 – cos( – )]. As explained in Kempter et al. (2012), we can transform the minimization over into
the following maximization:v"
ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
#2 "
#2ﬃ
u
u 1X
X
1
^q ¼ arg maxq tcosð j 2pqXj Þ þsinð j 2pqXj Þ ;
(11)
S j
S j
where is the number of spikes. We perform the maximization numerically over q. The
maximum value is called the mean resultant length (or vector strength), and it is used as our fit
score. The minimizing value of 0 can be then calculated directly:Psinð j 2p^
qXj Þ^0 ¼ arctan P j
:
(12)
cosð
2p^
qXj Þ
j
j
Regression results are shown in Figure 3D, E.
We can also estimate the circular-linear correlation coefficient, which is analogous to thePearson correlation coefficient for two linear variables. Note that since the marginal
distribution of Q may be uniform, we use the following expression (Kempter et al., 2012):Pjeið j Qj Þið j þQj ÞjeP

qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ;^ ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
P2Þ P sin2 ðQj QÞ2j sin ð j
j

(13)

P
P
where Qj ¼ 2pj^qjXj is a circular variable related to Xj, and  ¼ arctanð j sin j = j cos j Þ andPP ¼ arctanð sin Qj =Q
cos Qj Þ are respective circular means. This correlation coefficient isjjshown in Figure 3C.

Normalized bump coordinate
Over the course of a simulation, attractor bumps trace out trails of activity on the neural sheet(Figure 2C). For a neuron at a given sheet location, we quantify how centrally an attractor
bump passes through it with a metric called the normalized bump coordinate. We first tally
the total number of spikes occurring at each location on the neural sheet across all excitatory

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

28 of 34

Research article

Neuroscience Physics of Living Systems

populations. This yields an overall firing rate for each sheet location. For a neuron at a given
location, we look for the central axis of the nearest attractor bump trail as follows. We find thehighest overall firing rate in a window of length 10 neurons and width 4 neurons centered at
the given location and oriented with long side perpendicular to the track orientation. The
normalized bump coordinate is the ratio between the given location’s overall firing rate andthis highest overall firing rate, which should correspond to an attractor bump’s central axis.
Bump coordinate results are shown in Figure 3—figure supplement 4G.

Simplified phase precession models
For our conceptual model with bump size oscillations, we imagine an attractor bump moving
through a recorded neuron. Although the bump naturally oscillates in speed, we advance thebump with constant speed in the neural sheet; similar results (not shown) are obtained if the
speed oscillations depicted in Figure 6B are included. We stipulate that the middle of thebump passes through the recorded neuron at position 0 and time 0. Thus, the 1D profile ofthe bump is given by the region betweenb ðtÞ ¼ vt  ½b þ Db cosð2pft

0 Þ;

(14)

where b Db and b þ Db are the smallest and largest radii of the bump, v is the constant bumpvelocity, f = 8 Hz, and the phase offset 0 is incremented by p/5 for every new lap. We
choose v ¼ 17 neurons=s, b ¼ 3:3 neurons, and Db = 1.4 neurons, which are close to the values
presented in Figure 4B and Figure 6B. Whenever b– (t)  0  b+(t), the neuron fires a spike,subject to a refractory period of 40 ms. Results are shown in Figure 4C.
For our simplified model with average bump dynamics, we first determine the averagebump shape as a function of theta phase as described later. As in the conceptual model
described immediately above, we imagine the central axis of this average bump moving
through a recorded neuron. We simulate 18 laps with evenly spaced initial theta phases. The
recorded neuron’s spiking rate is proportional to bump activity. We can linearly rescale bumpactivity to a desired maximum level across all theta phases; that is, the bump would still exhibit
different levels of maximum activity at each theta phase. Spikes are generated with subPoisson statistics using the decimation procedure described in Burak and Fiete (2009); we
pﬃﬃﬃ
choose the coefficient of variation to be 1= 8. We also use a refractory period of 40 ms.
Results are shown in Figure 4E, F.

Theta sequences
We visualize theta sequences by averaging the decoded position over theta cycles. From allruns along the track, we take overlapping quadruplets of consecutive theta cycles, discarding
those during which the animal comes within of the track ends to avoid edge effects. We assign
time 0 to the middle of each theta quadruplet. We then shift the decoded positions of eachquadruplet by the actual position of the animal at time 0, such that the origin corresponds tothe actual animal position at the middle of the theta quadruplet. Before averaging, we
proceed as in O’Neill et al. (2017) and divide the probability at each timestep by its
maximum value. We then separately average each timestep of the quadruplet, ignoring cycles
without spikes within that timestep.
To find the slope of the fastest decoded motion within the theta cycle, we use a similar
approach to O’Neill et al. (2017). For each time window whose midpoint lies within the third
cycle of the averaged theta quadruplet, we find the fit line that captures the most decoded
probability. In O’Neill et al. (2017), this involves summing all the probability that lies within acertain number of position bins away from each possible line. Instead of setting such a hard
cutoff, we pass the averaged probability distribution at each timestep through a Gaussian
filter of width. For each candidate line, we sum over the filtered probability of the closest
position bin at each timestep; this process is analogous to weighting each bin according to a
Gaussian function of its distance to the line and summing the unfiltered probability over all
bins. The fit line is the one with the largest total probability. Finally, across all time windows,
we choose the fit line with highest slope to represent the averaged theta sequence. These
sequences and the Gaussian-filtered, averaged probability distribution is shown in Figure 5B.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

29 of 34

Research article

Neuroscience Physics of Living Systems

Replay analysis
Replay detection
As in O’Neill et al. (2017), we look for replays within highly synchronous events (HSEs). To
detect HSEs, we first filter the multiunit activity during idle periods by a Gaussian of width andfilter size, and we remove all zero-valued timesteps. HSEs are time intervals lasting at least
during which the processed multiunit activity exceeds its 20th percentile and contains at least
one timestep that exceeds its 80th percentile (or 60th percentile for simulations with periodic
boundary conditions).
For each HSE, we look for linear replays by finding the fit line that captures the mostdecoded probability as we did for theta sequences. We first pass the probability distribution
at each timestep during the HSE through a Gaussian filter whose width is and whose maximum
value is 1. For each candidate line, we sum over the filtered probability of the closest position
bin at each timestep and divide by the number of timesteps that the line spends within track
limits. The result is an average probability captured by the line, which we call a replay score.
For each HSE, the fit line is the one with the highest replay score. Finally, we exclude HSEs
whose fit lines have a replay score less than 0.6 or spend less than 30 ms or 30 cm within track
limits. The remaining fit lines are replays. Extracted replays are shown in Figure 7E.
For simulations whose neurons exhibit two fields, the probability distribution of thedecoded position may be spread over both fields. Thus, we supplement each candidate linewith two more lines of the same slope. These three lines are separated by the distance
between grid fields. We sum over the filtered probability of the closest position bins to eachof these lines.

Replay properties
The speed of a replay is the magnitude of the slope of its fit line, as reported in Figure 7G.
The start position of a replay is the fit line position at the beginning of its HSE, which isclipped to 0 cm or 60 cm if it lies outside track limits. Similarly, the end position of a replay is
the fit line position at the end of its HSE, which is similarly clipped. These positions arereported in Figure 8C.

Simulations with fixed inhibitory drive and constant velocity
Setup and trajectory
These simulations are structured similarly to the main simulations on a 1D track. However, wechoose a constant inhibitory drive amag and eliminate theta oscillations with ath ¼ 0. During
setup, instead of evolving the network with an animal trajectory, we evolve 1000 timesteps
with constant velocity of magnitude 0.5 m/s and angle p/5. For the main simulation, we evolve
another 1000 timesteps with the same velocity.

Attractor bump speed
To extract the speed of attractor bumps, we first divide time into 40 ms bins. Within each bin,
we sum the spikes at each position on the neural sheet across all excitatory populations, andwe pass the result through a Gaussian filter of width. We thus obtain smoothed activity maps
at 40 ms intervals.
We find the seven activity peaks that remain closest to the center of the neural sheet
throughout the simulation. For each time bin, we take the average position of all seven peaks,
and fit this average motion with a linear model. The norm of the resulting velocity is the bumpspeed, as reported in Figure 6B.

Attractor bump size
To extract the size of attractor bumps, we return to the raw spike trains, which we register
across the course of the simulation by compensating for the overall bump motion. That is, toeach spike emitted by an excitatory neuron, we assign a shifted neural sheet position equal toKang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

30 of 34

Research article

Neuroscience Physics of Living Systems

the true position of the neuron minus the bump velocity obtained above multiplied by its
spiking time.
Around each of the seven activity peaks identified above, which are stationary after
registration, we draw a circle of radius 12 neurons and calculate the center of mass of all
enclosed spikes over the entire simulation. This point is the center of the registered activity
bump. We then calculate the registered distance of each spike to the center; the 90th
percentile is taken to be the radius of the activity bump. Bump diameters averaged over theseven activity peaks are reported in Figure 4B.

Attractor bump shape
To determine the average attractor bump shape during runs, we use periodic simulations withuniform excitatory drive to generate a set of attractor bumps that experience the same drive.
We reintroduce oscillations into the inhibitory drive. These simulations use aþ
mag ¼ 1:2,
amag ¼ 0:7, and ath ¼ 0:15.
We first perform registration, as described above, for 12 attractor bumps. We count the
number of spikes across all bumps at each registered location on the neural sheet (using
0.2 neuron bins) as a function of time within the theta cycle (using 1 ms bins). We then apply
Gaussian filters of width 3 time bins and 5 position bins to obtain average bump shape as afunction of theta phase, as shown in Figure 4D and Figure 4—video 2.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

31 of 34

Research article

Neuroscience Physics of Living Systems

Appendix 2
Simulations without brief allocentric correction
Model details
Simulation setup proceeds similarly to simulations with allocentric correction, except weevolve the network for 6 laps including idle periods after the three evolutions with constantvelocity input. We also do not learn allocentric connections.
Without allocentric correction, the main simulation can only be run for 12 laps instead of
36. Each lap consists of a 1.5 s run and a 2.5 s idle period. Excitatory and inhibitory drives
linearly interpolate between run and idle values over the first and last 0.3 s of each idle period.
We perform 20 replicate simulations in each direction instead of 10. Parameter values thatdiffer from those in Table 1 are amag ¼ 0:65 during runs and 0.1 when idle,  ¼ 2, a ¼ 0:12 s=m,
and Var½zP ðr; tÞ ¼ Var½z ðr; tÞ ¼ 0:0122 .

Analysis details
With fewer laps compared to simulations with allocentric correction, we use more lenient
cutoffs for both theta-related and replay analysis. During phase relationship analysis, we
require a minimum of 15 spikes, instead of 30. We also use a regression fit score cutoff of 0.3,
instead of 0.4, below which neurons are considered phase-independent. During replay
detection, replays must span a minimum of 20 cm, instead of 30 cm.

Results
In simulations without brief allocentric corrections, many theta-related and replay properties
are present a weaker level. We can only run the simulations for 12 laps, and we must use lower
velocity gain so that the attractor bumps do not travel far on the neural sheet (Appendix 2—
figure 1A and Appendix 2—video 1). Otherwise, attractor drift will cause the firing fields tochange significantly between the start and the end of the simulation.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

32 of 34

Research article

Neuroscience Physics of Living Systems

A

B

60

cm

m

0c

C

D
F

inhibitory drive

E
prob.
max
0

G

H

I

J

Appendix 2—figure 1. Theta-related and replay properties in simulations without brief allocentric corrections. (A) Left, track diagram. Right, neural activity over runs with each neuron coloredaccording to the track position at which it attained maximum firing rate. Red circles indicate
regions of recording. (B) Relationship between animal position and theta phase forrepresentative neurons. Clockwise from top left, relationships depicted are phaseindependence, phase locking, phase precession, and phase precession. Dots represent spikes
during runs in the directions indicated by arrows. Lines indicate fit by circular-linear regression.
Numbers in each panel from top left to top right indicate magnitudes of correlation
coefficients, regression fit scores, and precession ranges. (C–E) Data across all replicate
simulations. (C) Magnitudes of circular-linear correlation coefficients. Mean ± s.d.: 0.15 ± 0.11.
(D) Fit scores for circular-linear regression. (E) Regression slopes for neurons with fit score>0.3. (F) Decoded position averaged over theta cycles. Thick black line fit to theta sequence.
(G) Replays with cyan fit lines. (H) Mean actual run speed (error bars indicate s.d. over time),

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

33 of 34

Research article

Neuroscience Physics of Living Systems

mean theta sequence speed (error bars indicate s.d. over replicate simulations), and mean
replay speed (error bars indicate s.d. over replays). (I) Direction of replay propagation acrosssimulations. Means indicated with lines. p-values from the Mann-Whitney U test show that asignificantly greater number of replays propagate away from either actual position. (J)Distances between actual position and decoded position at the start or end of replays. Paired
medians compared by the Wilcoxon rank-sum test with indicated p-value. Scale bars,
50 neurons.

Appendix 2—video 1. Neural activity during runs and idle periods during a simulation without
allocentric corrections.
Left, position of the animal (black square) along a 1D track. Right,
neural activity of the E population. Each pixel is a neuron, with black corresponding to current
spikes and lightest gray corresponding to spikes 40 ms ago. Red circles indicate regions ofrecording.
During runs, our neurons can still be divided into subgroups whose activity exhibits
different phase relationships (Appendix 2—figure 1B). The magnitude of circular-linear
correlation (Appendix 2—figure 1C; 0.15 ± 0.11, mean ± s.d.) is slightly lower than it is withallocentric correction (Figure 3C; 0.17 ± 0.11, mean ± s.d.). Excluding phase-independentneurons with fit score <0.3 (Appendix 2—figure 1D), we see a large population of phaselocking neurons and a small population of phase-precessing neurons (Appendix 2—figure1E). If we use a high phase range cutoff of 120˚ for phase precession, we find that preferred
phase predominantly decreases with progress through the grid field, which is the biologicallyobserved direction. By decoding position from spike trains and averaging over theta cycles,
we see a sawtooth pattern containing theta sequences (Appendix 2—figure 1F), though thepattern is less organized than it is with allocentric correction (Figure 5B).
During idle periods, traveling wavefronts produce replays (Appendix 2—figure 1G andAppendix 2—video 1). The speeds of theta sequences and replays are 1.1 ± 0.3 m/s and
7.0 ± 2.6 m/s respectively (mean ± s.d.; Appendix 2—figure 1H), which agree well withcorresponding speeds in simulations with allocentric correction (0.9 ± 0.2 m/s and 7.0 ± 2.6 m/s), and which also agree well with experiments (O’Neill et al., 2017). As in Figure 8A–C, thelower velocity gain produces a bias in replay direction away from the current animal position(Appendix 2—figure 1I,J).
In summary, after removing allocentric corrections, our model maintains its features
involving theta sequences and replays. It still produces some degree of phase precession, but
for fewer neurons and larger precession ranges.

Kang and DeWeese. eLife 2019;8:e46351. DOI: https://doi.org/10.7554/eLife.46351

34 of 34