We modeled the dynamics of hippocampal place cells in a strongly recurrent circuit, such as area CA3, as an animal sequentially explored a series of distinct novel ring-like tracks, Figure 1a. The model consisted of recurrently coupled, excitatory stochastic firing rate neurons, each of which received a place-specific external input on any given track, and inhibition was modeled as a global inhibitory feedback, see Materials and methods for model details. To model the global remapping of place fields from one track to another, we randomized the position of the peak input, and as a consequence the place field, of each neuron. In this way, the ordering of cells according to their place field location on one track was random and uncorrelated with their ordering on any other track, see Figure 1b.

Figure 1

Download asset Open asset

We were interested in knowing how the exploration of a novel environment affected the pattern of synaptic connectivity between place cells via a temporally asymmetric plasticity rule, and how this in turn shaped the activity. While the plasticity rule we use is formally referred to as 'Spike-timing dependent’ (Kempter et al., 1999; Song et al., 2000; Pfister and Gerstner, 2006), our model neurons generate spikes stochastically as Poisson processes. Therefore the exact spike timing itself does not play any role but rather only the time variations in the underlying firing rate. Recent theoretical work has shown that the plasticity induced by a number of plasticity rules, including heuristic STDP rules, is in fact dominated by such firing rate variations when irregular, in-vivo like activity is considered (Graupner et al., 2016). We also considered a ‘balanced’ plasticity rule, for which the degree of potentiation and depression were the same, on average, given constant firing rates. Such a rule is a means of avoiding fully potentiating or fully depressing all synapses on long time-scales, that is it allows for the generation and maintenance of structured synaptic connectivity. Alternatively, we could have considered an unbalanced plasticity rule with additional stabilizing mechanisms (Zenke et al., 2015).

In order to study ‘typical’ place-cell dynamics during exploration we first exposed the network to a series of distinct tracks until the matrix of recurrent synaptic weights no longer depended on its initial state. When the trajectory and velocity of the virtual animal was stochastic and distinct on different tracks, the resulting evolution of the synaptic weights also exhibited some variability from track to track, see Figure 2—figure supplement 1. In simulations with constant velocity this variability largely vanished, see Figure 2—figure supplement 2a. After exploring 10 tracks, we exposed the network to another novel track and studied the dynamics in detail; because the exploration process had already become largely stereotyped (despite some variability), the dynamics on the novel track reflected what would be seen qualitatively during the exploration of any other novel track in the future, that is it was dependent only on the learning process itself and not the initial state of the synaptic weight matrix.

In our simulations, we modeled the movement of a virtual animal around the track by varying the external inputs to place cells. Specifically, the external input was maximal for cells with place fields at the animal’s current position and minimal for cells with place fields at the opposite side of the track, with the input decaying smoothly for intermediate locations like a cosine curve, see Materials and methods for details. We modeled the motion of the animal by taking the velocity to be an Ornstein-Uhlenbeck process with a time constant of 10 s; this lead to realistic trajectories with changes of direction, see Figure 2—figure supplement 3. (We will use the metaphor of a virtual animal for conceptual ease, although we really mean that we moved the bump-like external input to our model neurons). Every three minutes of simulation time we stopped the animal at its current location for three seconds and removed the place field input. This led to spontaneous bursting via a synaptic depression-dependent mechanism (Romani and Tsodyks, 2015), reminiscent of sharp-waves seen during awake quiescence and sleep, see Figure 2. Neither the bursting, nor the learning process itself depended qualitatively on the exact amount of time the virtual animal spent moving versus staying still. Synaptic plasticity was allowed to occur during both theta activity and sharp-waves, that is it was never ‘turned-off’.

Figure 2 with 8 supplements see all

Download asset Open asset

As the virtual animal explored the track, the sequential activation of the place cells led to changes in the recurrent connectivity via the plasticity rule, Figure 2a. Whereas the connectivity between cells was initially unstructured on the novel track, it evolved over the span of minutes to tens of minutes such that cells with nearby place fields were more strongly connected than ones with disparate place fields. Furthermore the resulting connectivity was slightly asymmetric, reflecting the bias of the animal’s motion, see Figure 2—figure supplement 3. Although the changes we observed in the recurrent connectivity, see Figure 2a, and concomitant changes in place cell activity were continuous and smooth during the duration of the simulation, there was a dramatic and sharp transition in the structure of the burst activity during awake quiescence. We quantified this transition by measuring the mean pairwise cross-correlation of cells with neighboring place fields on the novel track. Such cells are ‘neighbors’ only on the novel track and not on any other track by virtue of the global remapping or randomization of place-fields. We call this measure the sequential correlation (SC) because it is high when there is a properly ordered sequence of place-cell activations on any given track. The SC during theta-activity (exploration) was already non-zero at the outset of the simulation by virtue of the external place-field input which generates sequential activity, black squares Figure 2b. On the other hand, the SC was initially near zero during spontaneous bursts, red circles Figure 2b. Interestingly the SC abruptly increased during burst activity at a critical time and remained elevated for the duration of the simulation. This transition occurred even in simulations where we included strong heterogeneity in place field tuning, see Figure 2—figure supplement 4. Note that the SC for theta-activity also showed a steady increase leading up to this transition, as did the SC when the total activity was taken into account, without discerning between epochs of movement or awake quiescence, see Figure 2b inset. This steady increase is due to the emergence of so-called theta sequences (Foster and Wilson, 2007). These are bursts of place-cell activity which sweep ahead of the animal’s position on every theta cycle during awake exploration, see Figure 2—figure supplement 5. The emergence of theta sequences here also coincides with a precession of the phase of the underlying firing rate with the ongoing theta rhythm, see Figure 2—figure supplement 5f and discussion for more detail.

Space-time plots of place-cell activity show that the abrupt increase in SC coincides with a transition in the spatio-temporal structure of the bursts: it is initially uncorrelated with the ordering of cells on the novel track, Figure 2c (i), and after the transition there is a clear sequential activation reminiscent of so-called ‘replay’ activity, Figure 2c (ii). In fact, before the transition the bursts are highly structured, but this structure is correlated with a previously explored track and not the novel one. If the learning process was carried on for much longer times all relevant network properties, such as connectivity and mean firing rates, saturated, see Figure 2—figure supplement 8.

We hypothesized that changes in the recurrent connections between place cells in our model during exploration were shaping the spontaneous activity and driving the transition to replay-like activity. The connectivity profile at any point in time could be decomposed into a series of spatial Fourier modes, Figure 3a. We tracked these modes in time and discovered that the transition in burst activity always occurred when the coefficient of the first even mode reached a critical value. Specifically, we conducted additional simulations in which we increased the firing rate of the feedforward inputs which drove place cell activity which in turn led the transition in burst activity to occur at earlier times, see shaded regions in Figure 3b. In fact, a theoretical analysis of our model shows that there is an instability of the spontaneous activity to bursts when the even mode of the connectivity times the gain in the neuronal transfer function is greater than a critical value, see Materials and methods. In our simulations the gain in the transfer function was nearly always the same during bursts because it depended on the mean input to the network during quiescence, which was constant. Therefore it was principally the even mode of the recurrent connectivity which determined the time of the transition. For this reason, the transition to replay occurred when the coefficient of the even mode reached a critical value as seen in Figure 3c. Note that there was no such dependence on the odd mode (dotted lines in Figure 3c).

Figure 3

Download asset Open asset

We then asked how the growth of the even mode depended on the details of the plasticity rule and the frequency of periodic modulation of the place cell activity. We found that the growth of the even mode depended strongly on the frequency, peaking at an intermediate range and decreasing by orders of magnitude at low and high frequencies, Figure 4ai. This occurred independent of the details of the animal’s motion, for example compare Figure 4ai and Figure 2—figure supplement 2e for irregular versus purely clockwise motion respectively. On the other hand the odd mode, which represents a directional bias in the recurrent connectivity, closely tracks the cumulative bias in the animal’s motion, see Figure 2—figure supplement 3, and is independent of the forcing frequency to leading order, see Materials and methods for details. The strong dependence of the growth of the even mode on frequency, meant that, for simulations of up to 1 hr, transitions in the burst activity were observed only when the frequency was in the range of 2–10 Hz, Figure 4aii. Furthermore, the even mode only grew, and transitions only occurred, when the plasticity rule had dominant potentiation at short latencies. For a perfectly anti-symmetric plasticity rule, and for a rule with dominant inhibition at short latencies, the even mode did not change and even decreased, respectively, Figure 4b and c. In the latter case bursts were suppressed entirely (not shown).

Figure 4

Download asset Open asset

We sought to understand the mechanism underlying the evolution of the recurrent connectivity seen in simulations by studying a continuum version of a network of place cells, which can be written

where $r(\theta ,t)$ is the firing rate of a place cell with place field centered at a location $\theta$, $w(\theta ,{\theta }^{\text{ ′}})$ is the synaptic weight from a cell at a position ${\theta }^{\text{ ′}}$ to a cell at a position $\theta$, and $I(\theta ,t)$ is the external input which has the form

The position of the virtual animal is given by $x\left(t\right)$, and the corresponding place-field input is modulated with a frequency $f$. This type of multiplicative theta modulation is seen in intracellular recordings in-vivo, for example see Figures 1 and 5 in (Harvey et al., 2009) and Figure 4 in (Lee et al., 2012).

Figure 5

Download asset Open asset

To model the evolution of the recurrent connectivity we made use of the fact that the step-wise increases in synaptic strength due to the plasticity rule can be approximated as a smooth process as long as plasticity occurs much more slowly than the firing rate dynamics. When this is the case, the rate of change of the synaptic weight from place cell with place field centered at ${\theta }^{\text{ ′}}$ to one with place field at $\theta$ can be written as

where $K\left(T\right)$ is the change in the synaptic weight according to the plasticity rule given a spike pair with latency $T$ (Kempter et al., 1999) and see Materials and methods. This equation reflects the fact that the total change in the synaptic weight is the sum of all the pairwise contributions from the pre- and post-synaptic cells, with each pair of spikes weighted by the plasticity rule with the appropriate latency. (Equations 1–3) represent a self-consistent model for the co-evolution of the firing rates and synaptic weights in the network.

In order to derive an analytical solution we first assume that the neuronal transfer function $\varphi$ is linear. We then make the assumption of slowly evolving synaptic weights explicit by scaling the amplitudes of the potentiations and depressions from the plasticity rule by a small parameter. The upshot is that the connectivity evolves to leading order only on a slow time scale, much slower than the fast neuronal dynamics. Furthermore, we know from numerical simulations that after sufficient exploration the probability of connection between any two cells depends on average only on the difference in place field locations. Therefore, by averaging the connectivity over the fast time $t$ we can write

where the growth rates of the even and odd modes ${\displaystyle {⟨{\dot{w}}_{even}⟩}_t}$ and ${\displaystyle {⟨{\dot{w}}_{odd}⟩}_t}$ are functions of the plasticity rule parameters, the velocity of the animal and the frequency of periodic modulation, see Materials and methods for details. It turns out it is possible to understand these dependencies intuitively and comprehensively without having to study the analytical formulas. Specifically, if we wish to isolate the growth rate of the even mode, which is responsible for driving the emergence of replay in the burst, we can consider place cell pairs where $\theta ={\theta }^{\text{ ′}}$, that is pairs with overlapping place fields. When this is the case we can combine (Equations 3 and 4) to yield

where ${\displaystyle AC\left(T\right)={⟨r\left(\theta ,t\right)r\left(\theta ,t+T\right)⟩}_t}$ is the autocorrelation (AC) of the place-cell activity. Note that despite the similarity in form between (Equation 5) and (Equation 3), the biological interpretation of the two is quite distinct. (Equation 3) describes the changes in the strength of a specific synapse, that from a cell with place-field centered at a position ${\theta }^{\text{ ′}}$ onto a cell with place-field centered at a position $\theta$. The evolution of this single synapse depends only on the sum of the changes from all spike pairs from these two cells. If we averaged this equation over fast fluctuations in time we would find that the slow evolution of the synapse depends on the cross-correlation of the activity of the two neurons times the STDP window. There is an equation like this for every cell pair in the network; clearly for a large network it is not feasible to solve these equations self-consistently. On the other hand, (Equation 5) does not describe the changes in strength of a single synapse. Rather, it describes changes in the strength of a particular pattern of synaptic connectivity in the network. This pattern is one in which cells with highly overlapping place fields have strong and symmetric recurrent connectivity. Furthermore the strength of the synaptic connections decays smoothly with the difference between place field locations. In our theoretical model, in the limit of a large network, the cross-correlation of very nearby cells is simply given by the auto-correlation, which is why it appears in (Equation 5). It is perhaps remarkable that using purely local information, namely the AC of place cell activity, one can infer the growth rate of a network-wide mode of synaptic connectivity. The key assumption that makes this possible is that the synaptic weight between any two cells should depend only on the difference in place-field location and not on the absolute position, (Equation 4). Therefore the growth rate of the even mode is found by multiplying the AC of place-cell activity by the window for plasticity, and integrating. Because the effect of periodic modulation on the AC is straightforward, we can determine graphically how the frequency of modulation interacts with the plasticity rule to drive changes in the burst structure.

We first note that if there is no periodic modulation of place-cell activity then the AC will simply reflect the movement of the animal. This will lead to a very broad AC compared to the time-scale of plasticity. For example, if we assume that the width of the place field is a fraction of the track length (as in our model), then a rat running between 5 and 50 cm/sec on a multi-meter track would have an AC which decays on the order of between seconds and tens of seconds. Therefore, the AC is essentially constant compared to a typical plasticity window of between tens and hundreds of milliseconds, and the integral in (Equation 5) will give nearly zero. Periodically modulating place-cell activity will increase the growth rate as long as potentiation is dominant at short latencies, Figure 5a. Slow modulation will bias the integral in (Equation 5) toward the potentiation lobe of the STDP (Figure 5a left and middle, top), while in an optimal range of frequencies the peaks and troughs of the AC will maximally capture potentiation and flip the sign of the depression lobe (Figure 5a left and middle, middle). Finally, at higher frequencies the plasticity rule undergoes multiple sign flips on a fast time scale, which again gives a near zero value for (Equation 5) (Figure 5a, left and middle, bottom). This means that the maximal growth rate of the even mode occurs for an intermediate range of frequencies: those which modulate place-cell activity on a time scale commensurate with the window for plasticity, Figure 5a right. Note that this has nothing to do with the overall rate of plasticity, which is the same independent of the modulation frequency. That is, even if the AC is flat, recurrent synapses undergo large numbers of potentiations and depressions. Rather, periodic modulation serves to organize the structure of synaptic plasticity at the network level by preferentially strengthening connections between place-cells with overlapping or nearby place field and weakening others. The optimal range of frequencies for growth of the even mode depends only weakly on the ratio of the width of the potentiation to that of the depression lobe, Figure 5c, but significantly on the total width, Figure 5d. Allowing for triplet interactions (as opposed to just pairwise) in the plasticity rule increases the overall growth rate but does not alter the range of optimal frequencies, Figure 5e. On the other hand, the theory predicts that the growth rate of the odd mode is only weakly dependent on the modulation frequency (Figure 5b) as is seen in simulations (Figure 2—figure supplement 2e), and can be understood by considering (Equation 3) with $\theta -{\theta }^{\text{ ′}}=\pi /2$. In this case the growth rate depends on the product of the plasticity rule with the cross-correlation (CC) of cells with disparate place fields. When there is an overall direction bias in motion then the CC will have peak shifted from zero-lag and the product with the STDP rule reliably gives a positive (or negative) growth rate. When there is very little motion bias the CC will be nearly flat, yielding little growth in the odd mode and the resulting connectivity will be highly symmetric.

It is also clear from (Equation 5) that a perfectly anti-symmetric plasticity rule will lead to no growth in the even mode as the product of the rule and the AC will itself be an odd function, see Figure 4b. A rule with dominant depression at short latencies will similarly yield a growth rate which is precisely the inverse of that shown in Figure 5a. This causes a decay in the even mode resulting in a connectivity pattern with locally dominant inhibition, see Figure 4c.

Finally, to go beyond this qualitative treatment and find the exact dependence of the evolution of the synaptic weights on model parameters requires solving (Equations 1–3) self-consistently, see (Equations 43–45) in Materials and methods for the case when $v$ is constant. These are the equations we use to generate the curves in Figure 5. When the AC of the place cell activity is predominantly shaped by the modulation frequency $f$, compared to the motion of the animal (strictly speaking when $v/\left(2\pi f\right)\ll 1$, which is the case already for $f>1$ Hz and for reasonable values of $v$) then the connectivity evolves to leading order as

Equation 6 shows clearly that the growth of the even mode, which drives the emergence of replay, is proportional to strength of oscillatory modulation $I_{theta}$. Additionally it can be seen on inspection that the modulation frequency itself plays a crucial role in setting the growth rate, and peaks at intermediate values. On the other hand, the growth of the odd mode, which represents the degree of asymmetry in the recurrent connectivity, is driven entirely by the motion of the animal to leading order. In the case of constant motion, studied here, it depends only on the velocity of the animal.

The spatio-temporal structure of sharp-wave-like bursts in our model reflected the sequential ordering of place fields on the most recently explored track; correlation with earlier tracks was erased or greatly reduced. In reality, only a fraction of place cells have well-defined place fields in any given environment, that is coding is sparse (Bostock et al., 1991; Fyhn et al., 2007). We incorporated this sparse-coding strategy to our model by providing place-field input to only one half of the total neurons in the network; the other half received constant input. These place cells were chosen randomly from one track to the next and place field locations were assigned randomly as before. Therefore the overlap in the population of place-cells between any two tracks was also fifty percent. We then allowed the network to evolve by having the virtual animal explore thirty distinct tracks, each for 1 hr of simulation time, as before. The resulting matrix of synaptic connections was correlated with the ordering of place-cells in several past explored environments, Figure 6a. The amplitude of the even mode, which is responsible for generating spontaneous replay, decayed roughly exponentially as a function of the recency of the explored track, Figure 6a inset. This suggested that the replay dynamics in the network with sparse coding might be correlated with several previously explored tracks simultaneously (Romani and Tsodyks, 2015). Indeed, when the network was driven with a global, non-selective input, the replay was strongly correlated with the past two explored environments, Figure 6b, whereas the correlation with earlier environments was negligible. On the other hand, when the constant external drive was selective to the subset of place-cells on any given track, replay activity was robustly correlated with activity on that track only, Figure 6c. Such selective input may originate in cortical circuits which store environment-specific sensory information as stable patterns of activity; these patterns correspond to attracting fixed points in network models of long-term memory storage and memory recall (Hopfield, 1982; Tsodyks and Feigel'man, 1988; Recanatesi et al., 2015). Spontaneous switching between cortical representations during slow-wave activity would therefore engage distinct hippocampal replay patterns, which in turn could strengthen and consolidate the corresponding cortical memory trace.

Figure 6

Download asset Open asset

Therefore, sparse coding causes the synaptic connectivity matrix to be simultaneously correlated with place field orderings from multiple tracks and allows for robust replay of activity from those tracks given appropriate inputs. The number of different environments which could be simultaneously encoded in replay depended on model details, particularly the coding level, that is the fraction of active place cells in any given environment (not shown) (Battaglia and Treves, 1998).

Pairwise reactivations of CA1 place-cells with overlapping place fields during awake SWRs improve during exploration; they are stronger during late exploration than early exploration (O'Neill et al., 2006). This holds true not only for pairwise correlations, but also reactivations of entire neuronal ensembles, at least on linear tracks (Jackson et al., 2006). More recent work has shown that significant replay events during awake SWRs in rats running along a three-arm maze emerge abruptly only after a certain number of runs (Wu and Foster, 2014). These results are consistent with our model predictions. We additionally sought to directly test for a time-resolved increase in SC in neuronal data. We looked for this increase in multi-unit recordings of place-cell activity from the hippocampus of rats exploring novel, periodic tracks (Wang et al., 2015).

We first identified cells with well-defined place fields by extracting the coefficients and phase of the first two spatial Fourier modes of their time-averaged activity as a function of the normalized distance along the track (in radians), see Figure 7—figure supplements 1–3. We kept only those cells for which the ratio of the coefficients, the Tuning Index (TI), exceeded the threshold of one, indicating strong spatially selective activity, see methods for details. We then ordered the cells according to their phase (approximately the position of peak firing). The SC of activity over the total duration of the experiment using this ordering was significantly higher than 5000 randomly reshuffled orderings which on average gave SC = 0, see Figure 7a and b.

Figure 7 with 4 supplements see all

Download asset Open asset

When the medial septum (MS) was inactivated via muscimol the SC did not exhibit any dynamics as a function of time, Figure 7c. However, once the animal recovered from the muscimol the SC using the proper phase ordering exhibited an initial ramp over the first ten minutes of exploration and then remained high (significant difference between first point and all others, t-test with multiple-comparison Bonferroni correction, p-values<0.004 and between second and third points), see Figure 7d, solid circles. This is consistent with the model result which showed a similar ramping increase when the total activity (and not just bursts) was considered, see for example Figure 2bii. On the other hand, the SC computed for shuffled phases showed no dynamics and remained close to zero, see Figure 7c,d, open squares. This indicates that there are no global changes in neuronal correlations during exploration, which could occur, for example, due to slow changes in theta modulation or neuronal excitability. Rather, there is a sustained increase in pairwise correlations only between those neurons which encode nearby places, and only when strong theta modulation is present. This finding also held when we considered the more lax criterion $TI\ge 1/2$, see Figure 7—figure supplement 4 (rat A991). In a second animal there was an insufficient number of well-tuned units to repeat the analysis, see Figure 7—figure supplement 4 (rat A992). One possible confound in attributing this increase in correlation to theta modulation alone, is the fact that the firing rates of place-cells during MS inactivation were lower on average than during the post-muscimol experiment. Lower firing rates in the computational model lead to lower rates of plasticity. However, this difference in firing rates was not significant for the well-tuned neurons shown in Figure 7 ($TI\ge 1$), (t-test, p=0.30). An initial, sustained increase in SC was also observed in data from a separate experiment in which an animal was first exposed to a novel hexagonal track, Figure 7e (difference between first point and all others, t-test with correction for multiple comparison, p-values<10–6). Furthermore, no such increase was found on subsequent sessions on the same track, indicating that the change in correlation only occurred when the environment was novel. This result held when we considered the more lax criterion $TI\ge 1/2$ and also when all units were used, regardless of tuning, see Figure 7—figure supplement 3a (rat A992). In a second rat there were too few well-tuned cells to use the criterion $TI\ge 1$, but an initial and sustained increase in SC was also seen using all units, although this correlation was not always significantly different from that calculated from shuffled orderings of the cells (unshaded points in Figure 7—figure supplement 3a, rat A991).

We have presented a computational model of a hippocampal place-cell network in order to investigate how the exploration of novel environments shapes the patterns of recurrent synaptic connectivity. Because place-fields remap randomly from one environment to the next, the recurrent connectivity, shaped by previous learning, is initially uncorrelated with place-cell activity in a novel environment. Our major finding is that the learning rate during the exploration of a novel environment depends almost entirely on the product of the autocorrelation of place-cell activity and the window for plasticity. The integral of this product determines the growth rate of a global, network-wide pattern of synaptic connectivity, which results in strong local recurrence and long-range competition. It is this mode which drives spontaneous replay activity in our model network. The growth rate of this mode is maximum when place-cell activity is periodically modulated on a time-scale commensurate with the plasticity rule, which for realistic time constants yields frequencies in the theta range. Furthermore, lower and higher frequencies than theta lead to learning rates which are orders of magnitude slower. This suggests that the role of theta is to accelerate learning. Note that the overall rate of plasticity is not affected by the presence of oscillations. The number of spike pairs, and hence the number of potentiations and depressions, depends only on the firing rates. Rather, theta oscillations generate repeated pre-post pairings in both directions, which coupled with a plasticity rule with dominant potentiation at short latencies bias plasticity toward potentiating events for neurons with neighboring place fields. One signature of this mechanism is a continuous increase in the pairwise cross-correlation in the activity of neighboring place-cells leading up to a critical time. We have found evidence consistent with this by analyzing the activity of simultaneously recorded hippocampal place cells in a rat during the exploration of a novel track.

In our model we have assumed that plasticity occurs only in the recurrent excitatory synaptic connections, and not in the feed-forward inputs. Therefore we also assume that the place-field input, which peaks at the spatial position of the virtual animal at given moment in time, is itself stable. In fact, consistent with this assumption, it seems most place cells are active from the outset of exploration of a new environment, although see (Hill, 1978; Frank et al., 2004; Monaco et al., 2014). Furthermore cells tend to exhibit only subtle changes in the size and shape of their place fields over time (Mehta et al., 1997; Mehta et al., 2000), also consistent with our model, see Figure 2—figure supplement 2b. On the other hand, it has been shown in area CA1 that some place-cells exhibit place-fields only after several minutes of exploration (Frank et al., 2004; Frank et al., 2006). Recent intracellular recordings indicate that appearance of these ‘hidden’ place-fields requires the coincidence of active dendritic events and synaptic input via the Schaeffer collaterals (Lee et al., 2012; Bittner et al., 2015; Bittner et al., 2017). It may be that this mechanism for place-cell ‘generation’ is particularly salient in cells of area CA1 by virtue of being uniquely positioned to compare and integrate both entorhinal and hippocampal inputs. In any case the strongly recurrent nature of the network we study may make it a more relevant model for circuits in area CA3.

Nonetheless we would expect that changes in spiking activity arising in CA3 due to plasticity in the recurrent connectivity, as in our model, would be reflected in similar changes in the spiking activity of CA1 cells due to the direct inputs via the Schaeffer collaterals. More specifically, in contrast to plasticity mechanisms leading to the formation of place cells themselves, here we have modeled on how plasticity shapes the recurrent connections between already-formed place cells. We find that pairwise correlations between place cells with nearby preferred locations increases during exploration of a novel environment. Assuming such an increase occurs in a strongly recurrent circuit in CA3, we would also expect to observe an increase in correlation in target CA1 pyramidal cells, as long as there exists some mapping from CA3 place cells to CA1 place cells. Such a mapping could be a simple random projection or a more ordered relationship. Recent work suggests that place fields of CA1 pyramidal cells on linear tracks are built up of a weighted sum of CA3 inputs from positions surrounding the relevant place field (Bittner et al., 2017); in such a case the increased correlation in CA3 activity would be shared by CA1 output. Indeed, the data we have analyzed from place cells of area CA1 show increase sequential correlation as predicted by our recurrent model. A recent computational model of STDP-induced formation of place fields in CA1 cells via the feed-forward excitatory connections from CA3 (Schaeffer Collaterals) suggests a role for theta in speeding up place cell formation (D'Albis et al., 2015). This raises the intriguing suggestion that theta may play a key role both in place-cell formation through plasticity of feed-forward inputs, and also in the emergence of replay through plasticity of recurrent synaptic connections as indicated by our work here.

Hippocampal place cells actually exhibit global remapping of their place fields depending on the direction of motion of the animal on linear tracks (McNaughton et al., 1983; Muller et al., 1994). This is perhaps not surprising given that the behaviorally relevant information for the animal is not just the position along the track but also the way in which it is facing; for example this determines how far away a potential reward is, if located at one or both ends of the track. Studies using periodic tracks have shown no such global remapping, but rather some degree of rate remapping (Schwindel et al., 2016), that is the direction of motion affects the shape and amplitude of place fields, but not their position. In the data we have analyzed there is very weak remapping, see Figure 7—figure supplements 1c and 2b, and so the assumption of invariance of place field to direction of motion is a good one. In our model we have made this assumption. The consequence of this is that while exclusively clockwise (CW) or counter-clockwise (CCW) motion will lead to highly asymmetric recurrent connectivity, exploration of both directions will give rise to much more symmetric connectivity, see Figure 2a. In practice any trajectory over a finite amount of time will have a directional bias; in the data we have analyzed the rat spends a slightly larger fraction of the time moving CW compared to CCW, and this will necessarily lead to asymmetries in the connectivity. In linear tracks, due to the global remapping such asymmetries should be even more pronounced.

The inevitable asymmetry in the recurrent connectivity of our model place-cell network due to plasticity during exploration strongly biases spontaneous activity. On a periodic track this replay would be exclusively CW or CCW depending on the corresponding bias in motion during exploration, while learning on a linear track would always produce forward replay. Previous work has shown that perfectly symmetric connectivity can give rise to both forward and backward replay in equal measures, due to spontaneous symmetry breaking of activity (Romani and Tsodyks, 2015). We would argue that such symmetric connectivity is not robust for the reasons given above, although we cannot rule out the existence of homeostatic mechanisms which would conspire to make it so. Rather, we propose here an alternative mechanism for generating backward replay given asymmetric connectivity, based on local sensory input.

Specifically, if global input to our model network is homogeneous then replay occurs only in one direction. This is due to the asymmetry in the connectivity which comes about through a bias in the direction of motion of the animal during exploration. On the other hand, if the input is not homogeneous, but rather spatially localized, then the replay dynamics can be qualitatively distinct. This difference is mediated by the short-term synaptic depression mechanism present in the model. With a strong, spatially localized input, synapses to downstream neurons (in the sense of the asymmetric bias of the connectivity) become rapidly depressed, see Figure 2—figure supplement 6. This forces the activity to travel backward with respect to the bias in the connectivity. In fact, in experiment, when local spatial input is absent, for example when the animal is sleeping in a rest box, forward replay is predominant (Roumis and Frank, 2015). On the other hand, both backward and forward replay are observed when the animal is awake but quiescent on a given track. This is precisely when locally sensory cues are available to the animal, and could potentially shape spontaneous replay events. Recent work shows that some neurons in area CA2 fire more strongly during awake quiescence than during exploration (Kay et al., 2016); they may be providing information regarding local sensory cues. In our scenario the mechanisms leading to forward versus backward replay are distinct and therefore in principle relevant replay statistics such as replay velocity and spiking intensity should also be different.

Here we have considered a simple phenomenological model of plasticity which depends on the timing of spike pairs. Taking into account spike triplets as opposed to only pairs does not alter our findings qualitatively, see Figure 5e, although we are necessarily ensuring a balanced rule through our choice of parameters. Fits of heuristic spike timing-dependent models to data from slice experiments yield parameter values which do not lead to balanced rules such as the ones we have used here (Pfister and Gerstner, 2006). Specifically, when the repetition-rate of spike pairs is high, in-vitro STDP is dominantly potentiating; for such a rule there would be no non-monotonic dependence of the learning rate on modulation frequency at high rates. The mechanism we propose for theta would not be operative in that case. However, in a recurrent network model, unbalanced plasticity rules would quickly lead to synapses either saturating to maximum efficacy or vanishing completely; non-trivial network structure is therefore not possible with such rules. One solution to this problem is to complement the heuristic, unbalanced Hebbian rule, extracted from slice experiments, with a non-Hebbian heterosynaptic rule, as well as additional, slower homeostatic mechanisms (Zenke et al., 2015). It is unclear if the learning rate of this rule would also exhibit a non-monotonic dependence of the modulation frequency, as in our rule. However, heuristic rules describing synaptic plasticity in-vivo are still largely unknown. Therefore one reasonable approach to studying plasticity in recurrent networks in-vivo is to choose the simplest possible rule which allows for the emergence of non-trivial structure; this has been our approach here. It remains to be studied how more realistic voltage- or calcium-based plasticity rules interact with the theta-modulation to affect learning in recurrent networks (Clopath et al., 2010; Graupner and Brunel, 2012), although at a single-synapse one can find qualitatively similar regimes for an array of plasticity rules in the presence of pre- and post-synaptic oscillations (Albers et al., 2013).

Our results clearly do not depend on the actual spike timing since our model neurons generate spikes as Poisson processes; rather, all lasting changes in the connectivity are due to time-varying modulations of the firing rates. In fact, recent work with a spiking neuron model suggests that such modulations in the firing rate, as opposed to exact spike timing, are sufficient to explain the effect of plasticity from STDP and more realistic calcium-based plasticity rules in general (Graupner et al., 2016). In any case the contribution of pairwise spike correlations to the evolution of the recurrent connectivity can be formally taken into account in (Equation 4), that is via its affect on the auto-correlation of place-cell activity.

The increase in SC during theta activity in our model is due to the emergence of theta sequences, Figure 2—figure supplement 5, a consequence of the strengthening of the recurrent excitatory connectivity during exploration. The emergence of theta sequences also gives rise to phase precession of the place cell activity. Specifically, the peak of the place cell firing rate shifts to earlier phases of the theta rhythm as the animal enters the place field. This mechanism was first studied in (Tsodyks et al., 1996), in which theta sequences arose through asymmetry in the recurrent connectivity. A very similar effect can be achieved through short-term synaptic depression even without the need for asymmetric connectivity (Romani and Tsodyks, 2015). The reason is that the motion of the animal, and hence the sequential activation of place cells, generates an asymmetric pattern of synaptic depression (upstream synapses are more depressed than downstream ones). This is the mechanism responsible for theta sequences and phase precession in our model. Phase precession therefore emerges over time in our model, and is not present during early exploration. We have also studied the effect of phase precession on learning numerically by generating theta sequences directly via the place-field input itself (not shown). Preliminary simulations suggest that phase precession can actually speed up learning, although the mechanism is non-trivial and requires additional study.

Recurrent network models for place-cell activity provide a parsimonious explanation for electro-physiological phenomena associated with exploratory behavior as well as the generation of sharp-wave bursts during awake quiescence and sleep (Romani and Tsodyks, 2015; Tsodyks et al., 1996; Shen and McNaughton, 1996; Cutsuridis and Hasselmo, 2011; Jahnke et al., 2015). In recent theoretical work (Romani and Tsodyks, 2015) sharp-wave-like bursts were generated spontaneously by virtue of the spatial modulation of the recurrent connectivity, which drives an instability to traveling waves in the absence of place-field input. The presence of short-term depression modulates the amplitude of the waves, leading to narrow bursts. This is the same mechanism we have used here. Alternatively, recent work with a biophysically detailed spiking network model focused on the role of nonlinear dendritic integration on the generation of replay during SWRs (Jahnke et al., 2015). In that work the authors found that the exploration of a virtual linear track in the presence of pairwise STDP led to highly asymmetric connectivity; this could generate replay activity given a sufficiently synchronous external input which recruited nonlinear dendritic events. In our work, we have sought to explain the replay as an emergent phenomenon which depends only on the network-wide organization of synaptic structure. In doing so we have considered a simple stochastic firing rate model which allowed us to fully characterize how interplay between the plasticity rule and the place-cell activity affects learning. Nonetheless, a detailed reproduction of the phenomenology of SWRs certainly requires mechanisms we have not included here. In particular, while our model generates sharp-wave-like bursts, it does not generate high-frequency ripples, which are likely generated by networks of inhibitory interneurons in CA1 (Buzsáki, 2006).

It seems reasonable that the learning of tasks which depend on spatial information require the formation of an internal representation of the relevant environment. This is the process we have studied here. While we have not modeled any particular cognitive task, we propose that the network-wide organization of synaptic structure, in order that it be in concordance with the place-field distribution of place cells, should be a necessary step in spatial learning tasks. Our results suggest that this process is dramatically sped up by modulating place-cell activity in the theta range, which is one possible role of this prominent rhythm. Interestingly, while theta modulation is prominent in land-going rodents, it is absent in the bat (Yartsev and Ulanovsky, 2013). Although we cannot purport to know why this is, we find that it actually may be consistent with the theoretical mechanism for learning rates we put forth here. Namely, to attain a high learning rate the AC of place-cell activity should be modulated on a time-scale commensurate with the window for plasticity. In the case of land mammals like rats and mice, the behavioral time-scale related to spatial navigation may be too slow, therefore necessitating oscillatory modulation in the form of theta. In the case of flying animals, like the bat, changes in sensory input alone are already likely to occur on the order of tens to hundreds of ms due to the high velocity of flight, thereby obviating the need for theta. This speculation is based on the modulatory effect of ongoing oscillations on learning in our computational model, and is therefore not yet supported by experimental evidence. Indeed, our theoretical prediction is that more generally, for learning to occur on behaviorally relevant time scales, neuronal activity should vary somehow on a time-scale commensurate with synaptic plasticity mechanisms. One means of achieving this is to have a broad window for synaptic plasticity, on the order of seconds, so-called behavioral time-scale plasticity (Bittner et al., 2017). Alternatively, internally generated rhythms may serve to modulate neuronal activity on non-behavioral time-scales, with a similar effect.

We simulate a network of $N$ place-cells. The firing rate of neuron $i$ is given by

where the time constant $\tau =10$ ms and $N=100$ for all simulations except for those in Figure 6 in which $N=200$. We take ${\displaystyle \varphi \left(I\right)=\alpha \mathrm{l}\mathrm{n}\left(1+e^{I/\alpha }\right)}$ with $\alpha =1$ Hz for all simulations. The input

where ${\theta }_i^{\mu }$ is the center of the place field for cell $i$ on track $\mu$. These place field positions are random and uncorrelated for each cell $i$ on different tracks. The position of the animal is given by ${\displaystyle X\left(t\right)}$. For the simulations in all figures except Figure 2—figure supplement 2 the position was calculated by assuming that the velocity of the animal was random with a characteristic time constant of 10 sec. Specifically the velocity is given by $v=v_0+v_1\left(t\right)$ where the mean velocity ${\displaystyle v_0=0.5}$ rad/sec and

where ${\xi }_{\nu }$ is a Gaussian white-noise process with mean zero and unit variance, with ${\tau }_{\nu }=10$ sec and ${\sigma }_{\nu }=2$ rad/sec. The position of the animal is then given by solving ${\displaystyle \dot{X}=v}$. For Figure 2—figure supplement 2 the velocity is taken to be a constant $v=1$ rad/sec and the position ${\displaystyle X=vt}$.

The short-term synaptic depression variable $x_i$ obeys

with ${\tau }_x=800$ ms and $U_0=0.0008$.

We first consider a simple case where the place cell activity is given by ${\displaystyle r\left(\theta ,t\right)=r_0+r_{PF}\cos \left(\theta -vt\right)}$, that is the rate dynamics trivially follow the motion of the virtual animal and are not ‘theta’-modulated. Also, we are taking a continuum limit in which the index of a neuron $i$ is replaced with the position of its place field along the track $\theta$. In (Equation 3) we have the product

where $f\left(t\right)$ is a time-periodic function. The first two terms in (Equation 18) will therefore determine the mean growth of the even and odd modes of the connectivity on a slow time-scale, while the remainder will lead to fluctuations about this mean on a faster time scale. (This time-scale separation will be conducted more formally in a later section) Averaging over the fast time eliminates these fluctuations. Therefore the rate of change of the connectivity can be expressed as

where the brackets on the l.h.s. indicate that we have taken a time-average to eliminate time-oscillating terms. Therefore, (Equation 3) reduces to

Performing these integrals gives

We see from (Equation 22) that the balance condition $A_{+}{\tau }_{+}=A_{-}{\tau }_{-}$ must hold to avoid the synapses from all reaching their maximum or minimum value. At this point we would also note that for any reasonable velocity the growth rate of the even mode is extremely small and far from the maximal possible growth rate. Specifically, the velocity has unit of radians. For a three meter track, a rat running at 50 cm/s would have a velocity of about 1 rad/sec. Typical time scale for STDP are 50 ms or 0.05 s. Therefore the non-dimensional quantity ${\tau }^2{v}^2=0.0025\ll 1$ and, given the balance condition, ${\displaystyle {\dot{w}}_{\text{even}}\sim 0}$.

Including the triplet interactions by using (Equation 17) gives

From (Equation 25) we see there are two additional balance conditions to avoid saturation or decay to zero of synaptic weights: ${\tau }_x={\tau }_y$ and $A_3^{+}{\tau }_{+}{\tau }_y=A_3^{-}{\tau }_{-}{\tau }_x$. Once we apply the balance conditions the growth rate of the mean connectivity $w_0$ is identically zero.

Here we consider the case where the place cell activity is modulated periodically with a frequency $f$, that is ${\displaystyle r\left(t\right)=r_0+r_{\text{PF}}\cos \left(\theta -vt\right)\left(1+r_{\text{theta}}\cos \left(2\pi ft\right)\right)}$. This can be rewritten as

where $v_{+}=v+2\pi f$ and $v_{-}=v-2\pi f$. Now when we plug (Equation 28) into (Equation 3) or (Equation 17) we find that each of the three cosine terms makes an independent contribution. The reason is that the cross terms (which arise due to the product of rates) lead to time periodic terms which vanish once we average the equation. Therefore we end up with equations of the form (Equation 22–24) or (Equations 25–27) but with three times the number of terms, that is for $v$, $v_{+}$ and $v_{-}$. These equations can, however, be simplified by considering the limit where the ratio $ϵ=\frac{v}{2\pi f}\ll 1$. For example, when $f=8$ Hz and $v=1$ rad/sec (50 cm/sec on a three meter track) then $\frac{v}{2\pi f}\sim 0.02$. Slower running speeds or faster oscillations will make the following approximation even more accurate. In this limit we have (for the pairwise rule)

Furthermore, given that $\tau v$ is also a small parameter, we can neglect the first term in (Equation 29), which means that the growth rate of the even mode is entirely dominated by the modulation frequency to leading order. Because $f$ could take on any value the ${\displaystyle {\dot{w}}_{\text{even}}}$ need not be close to zero. In fact, there is a critical value of f which maximizes the growth rate $f=\frac{1}{2\pi \sqrt{{\tau }_{+}{\tau }_{-}}}$. On the other hand the growth of the odd mode is not dependent on the modulation frequency to leading order, but only on the motion of the animal.

Now we return to the full problem, namely, (Equations 1–3) and we take $\varphi$ to be linear for simplicity. In principle (Equations 1–3) are extremely challenging to analyze by virtue of the nonlinear combinations of the rates and the connectivity, both of which evolve in time. However, if we assume that the connectivity evolves much more slowly than the rates, the problem simplifies considerably. Specifically, we assume that the changes to the synaptic weights due to each spike pair are small. Formally, we can introduce a small parameter $ϵ$ and write ${\displaystyle A_{+}=ϵ{\tilde{A}}_{+}}$ and ${\displaystyle A_{-}=ϵ{\tilde{A}}_{-}}$, and so ${\displaystyle K\left(T\right)=ϵ\tilde{K}\left(T\right)}$, where all of the ‘tilded’ quantities are order one. In the limit $ϵ\to 0$ we would then find that (Equation 1) is unchanged while (Equation 3) becomes ${\displaystyle \dot{w}\left(\theta -{\theta }^{{\mathrm{\prime }}^{}}\right)=0}$, that is the connectivity is a constant. This clearly is an approximation since the connectivity will change over time, just slowly. In fact it is precisely this slow evolution that we would like to describe. In order to do this we can define a new, slow time ${\tau }_s=ϵt$ and allow both the firing rate and the connectivity to evolve on both the fast and the slow time scales, which are now taken to be independent. This is known as a multi-scale approach.

Formally we make the following ansatz

where $t_s=ϵt$. Therefore, the temporal derivatives in (Equations 1–3) can be written ${\displaystyle \dot{r}={\mathrm{\partial }}_{t}r+ϵ{\mathrm{\partial }}_{t_s}r}$ by using the chain rule. This leads to the system of equations

To leading order we have the equations

(Equation 36) shows that, to leading order, the connectivity is independent of the fast time $t$, that is $w_0=w_0(\theta -{\theta }^{\text{ ′}},t_s)$. Given this, the term $w_0$ can be treated as a constant in (Equation 35), and the rate is given by $r_0=r_0(\theta ,t_s,t)$. At order $\vartheta \left(ϵ\right)$ the equations are

The solution to (Equation 37), $r_1$ gives the next order correction to the firing rates. On the other hand, (Equation 38) would seem to require solving for both $w_0$ and $w_1$ simultaneously. However, we note that the product of rates $r_0(\theta ,t_s,t)r_0({\theta }^{\text{ ′}},t_s,t+T)$ can actually be written as $g(\theta -{\theta }^{\text{ ′}},t_s;T)+f(\theta -{\theta }^{\text{ ′}},t_s,t;T)$. That is, there is a part which is independent of time, except for the dependence of the slow time through the connectivity, and a part which varies on a regular, non-slow time scale. This fast, time-varying part, has zero mean (it is periodic here) and hence can be eliminated by averaging on an appropriate intermediate time scale, which is still much faster than the slow time scale, that is ${\displaystyle {⟨f⟩}_t=0}$. Therefore, we can define $w_1$ as that component of the connectivity which is driven by $f$, and hence also will have zero mean and will vanish upon averaging. Finally, combining the order one equation for the rate and the averaged order epsilon equation for the connectivity yields

(Equations 39 and 40) are the system we wish to solve here. We first note that the connectivity can always be decomposed into a Fourier series. Anticipating the fact that the only terms which will not have zero mean in time will be the first cosine and sine terms, we write

Given that the forcing term can be written

then the firing rates $r_0$ will be

where, after some linear algebra one finds

Taking the product of the rates and keeping only those terms with non-zero mean yields the function $f$, which when plugged into (Equation 40) gives

In the limit $v/\left(2\pi f\right)\ll 1$ the leading order terms of (Equations 44 and 45) are given by (Equations 6 and 7).

The solution for the triplet rule can be similarly found by replacing $r_{\text{PF}}$ in (Equations 25–27) by $(r_{\text{even}}^2+r_{\text{odd}}^2)/2$ as well as the analogous terms corresponding to $v_{+}$ and $v_{-}$.

Up until now we have considered a constant velocity $v$ for ease in calculation. In reality, the motion of the animal will be much more erratic, with both positive and negative speeds during the exploration of the track. In the extreme case of a random walk, in which both clockwise (CW) and counterclockwise (CCW) motion is equally likely, the resulting connectivity is much more symmetric. This occurs because contributions from CW and CCW motion to the even mode simply add, while the same contributions to the odd mode cancel out. We can illustrate this more clearly by considering the continuum learning rule, (Equation 3), with a given distribution of velocities. Specifically, we take the simple case where the firing rates can be written ${\displaystyle r_0+r_{\text{PF}}\cos \left(\theta -vt\right)}$ (we have seen that the solution to the full model can be written in an analogous way to the solution of this simpler problem) and assume a distribution of firing rates $\rho \left(v\right)$. From (Equations 20 and 21) we have

If we choose $\rho \left(v\right)=\delta (v-v^{\ast })$, that is, a fixed velocity $v^{\ast }$, then we recover our results from previous sections. We can allow for a range of velocities by considering a uniform distribution of velocities, that is

For this choice of distribution we find ${\displaystyle {\dot{w}}_{\text{even}}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}dT\cdot K\left(T\right)\frac{\sin \left(v^{\ast }T\right)}{v^{\ast }T}}$, while ${\displaystyle {\dot{w}}_{\text{odd}}=0}$. The function ${\displaystyle \sin \left(x\right)/x}$ is even and decays more slowly than ${\displaystyle \cos \left(x\right)}$ around zero lag, meaning that the AC of the firing rate is qualitatively similar whether $v$ is a constant or distributed. On the other hand, the CC of neurons with disparate place field location will be zero on average if the animal samples both directions equally. Nonetheless, any asymmetry in the velocity distribution (which is expected for any finite length of exploration) will lead to a non-zero odd mode in the connectivity.

In the previous sections we have seen how the connectivity grows in time due to the non-stationarity of the firing rates. In particular, the AC of the place-cell activity drives growth of the even mode, and this growth depends on the STDP rule parameters and the modulation (theta) frequency. Here we show that once the amplitude of the even mode reaches a critical value there is an instability of the network activity in the absence of place field input, which gives rise to spontaneous traveling waves. This is the origin of replay activity in the model.

We consider the firing rate equation in the absence of place field input, that is spontaneous dynamics.

and a small perturbation about a steady state solution of the form ${\displaystyle r\left(\theta ,t\right)=r_0+\left(\delta r_0+\delta r_{\text{even}}\cos \left(\theta \right)+\delta r_{\text{odd}}\sin \left(\theta \right)\right)e^{\lambda t}}$. Plugging this into (Equation 49) leads to a characteristic equation for $\lambda$ the solutions of which are