A theory of hippocampal theta correlations accounting for extrinsic and intrinsic sequences
Abstract
Hippocampal place cell sequences have been hypothesized to serve as diverse purposes as the induction of synaptic plasticity, formation and consolidation of longterm memories, or navigation and planning. During spatial behaviors of rodents, sequential firing of place cells at the theta timescale (known as theta sequences) encodes running trajectories, which can be considered as onedimensional behavioral sequences of traversed locations. In a twodimensional space, however, each single location can be visited along arbitrary onedimensional running trajectories. Thus, a place cell will generally take part in multiple different theta sequences, raising questions about how this twodimensional topology can be reconciled with the idea of hippocampal sequences underlying memory of (onedimensional) episodes. Here, we propose a computational model of cornu ammonis 3 (CA3) and dentate gyrus (DG), where sensorimotor input drives the directiondependent (extrinsic) theta sequences within CA3 reflecting the twodimensional spatial topology, whereas the intrahippocampal CA3DG projections concurrently produce intrinsic sequences that are independent of the specific running trajectory. Consistent with experimental data, intrinsic theta sequences are less prominent, but can nevertheless be detected during theta activity, thereby serving as runningdirection independent landmark cues. We hypothesize that the intrinsic sequences largely reflect replay and preplay activity during nontheta states.
eLife assessment
This important work presents an interesting perspective for the generation and interpretation of phase precession in the hippocampal formation. Through numerical simulations and comparison to experiments, the study provides a convincing theoretical framework explaining the segregation of sequences reflecting navigation and sequences reflecting internal dynamics in the DGCA3 loop. This study will be of interest for researchers in the spatial navigation and computational neuroscience fields.
https://doi.org/10.7554/eLife.86837.4.sa0Introduction
As a rat navigates in an environment, place cells fire sequentially during one theta cycle (∼100 ms) and form timecompressed representations of behavioral experiences (Skaggs et al., 1996), called theta sequences. Theta sequences were proposed to be driven by extrinsic (extrahippocampal) sensorimotor input (Foster and Wilson, 2007; Huxter et al., 2008; Romani and Tsodyks, 2015; Yiu et al., 2022), since they are played out in the direction of travel during locomotion and, hence, represent current behavioral trajectories. In contrast, various types of hippocampal sequences have also been proposed to arise from intrinsic hippocampal connectivity. Nonlocal activation of place sequences during immobile periods was observed in replay of past locations after the space has been explored (Skaggs and McNaughton, 1996; Lee and Wilson, 2002) as well as in preplay (Dragoi and Tonegawa, 2011) of prospective locations before the animal explores a novel environment. In addition, some CA3 place cells exhibit outoffield firing at reward locations (Sasaki et al., 2018). These remote activations of place cells reflect the underlying circuit connectivity rather than the actual location and movement of the animal. Furthermore, a subset of CA3 cell pairs shows rigid theta correlations with peak lags that are independent of the traversal order of their place fields (Yiu et al., 2022), suggesting the existence of hardwired sequences even when sensorimotor drive is present. Such intrinsic sequences that are driven by intrahippocampal connectivity (Tsodyks et al., 1996), although less predominantly observed during theta (Yiu et al., 2022), are generally interpreted as reflecting spatial memories or planning (Kay et al., 2020).
Existing models for theta sequences are, however, either fully extrinsic or intrinsic. The former often employ shortterm plasticity (Romani and Tsodyks, 2015; Thurley et al., 2008), which creates synaptic couplings that are temporally stronger along the instantaneous forward direction. In contrast, intrinsic models, such as the Tsodyks et al., 1996 model, use a fixed asymmetrical weight matrix predesigned to align with one movement trajectory (for review see Maurer and McNaughton, 2007; Jaramillo and Kempter, 2017). Neither of these models alone can explain the simultaneous presence of rigid and flexible correlations in theta sequences. Here we present a network model that accounts for both types of correlations by separating their generation into two anatomically distinct layers: CA3 and dentate gyrus (DG). Extrinsic sequences are generated in the CA3 layer by shortterm synaptic plasticity mechanisms, while the intrinsic sequences are evoked by the CA3DG feedback loop with fixed asymmetrical weights, as inspired by the experimental evidence that lesions of DG abolish nonlocal activation of CA3 place cells (Sasaki et al., 2018) and CA3 theta correlations (Ahmadi et al., 2022).
In this paper, we present a model for theta correlations that unifies both extrinsic and intrinsic mechanisms. Extrinsic and intrinsic sequences can propagate simultaneously in separate directions, along the movement trajectory and the predesigned CA3DG feedback loops, respectively. As a result, spike correlations display directionality as the two sequences cross each other at various angles: The more parallel they are, the stronger the correlation. Our simulations are in quantitative agreement with directionality properties found in experimental data (Yiu et al., 2022) and propose that rigid correlation structure can serve as a stable temporal pattern, which is recognizable across multiple movement directions. This temporal ‘landmark’ pattern allows spatial encoding even if sensorymotor experience is lacking and may reflect the mechanistic basis for replay in nontheta states.
Results
Dependence of theta sequences on heading directions: Extrinsic and intrinsic sequences
Thetascale correlations of place cells have been explained by previous models using two different types of network mechanisms, intrinsic and extrinsic ones. For intrinsic models spike correlations are explained by only the recurrent connectivity of the neuronal network (Figure 1A). For extrinsic models, the spike correlation is defined by sensorymotor parameters such as movements (Figure 1B). We first illustrate how these mechanisms work for two exemplary representatives of these two major model classes.
For intrinsic models, we refer to the Tsodyks et al., 1996 model where phase precession is generated by the fixed asymmetrical connectivity between place cells. Spike phases of the place cells ahead of the animal decrease as the excitatory drive is gradually increasing, but only along the direction in which the connection strength is asymmetrically stronger (e.g. rightward in Figure 1A), called the asymmetry direction. Here we simulate a network of CA3 place cells with fixed asymmetrical connectivity (see Methods section: CA3 recurrent connections) as suggested in the Tsodyks et al., 1996 model and applied our model to behavioral running trajectories in a 2d open space (Figure 1C). Phase precession and spike correlations (Figure 1D–E) are compared for opposite running directions. In our simulations, all place cells project excitatory synapses to their counterparts with rightward neighboring place fields. Phase precession therefore is determined by how closely the running direction matches the asymmetry direction imposed by the intrinsic connections. The closer the trajectory angle aligns with this asymmetry direction, the more negative is the slope of phase precession (Figure 1E). Since in this case, the theta sequence only propagates rightwards as place cells are sequentially activated from left to right, the signs of spike correlations between cell pairs remain invariant to the movement direction (Figure 1F, see Methods section: Crosscorrelation analysis). Intrinsic models thus cannot explain experimentally observed directional independence of phase precession and directional dependence of theta spike correlations (Huxter et al., 2008).
Our example of an extrinsic model is based on our spiking simulations of the originally ratebased model by Romani and Tsodyks, 2015, where phase precession was explained by symmetric recurrent connections that undergo running directiondependent attenuation by shortterm synaptic depression (STD): place fields with centers behind the current animal position on the trajectory thereby received largely reduced recurrent input resulting in recurrently driven theta sequences to play out only in forward direction (see Methods section: CA3 recurrent connections). We simulated our spiking variant of the Romani and Tsodyks, 2015 model with the same trajectories as the intrinsic model (Figure 1G–I), and recovered directionindependent phase precession (Figure 1H). Since now, the theta sequences are played out in the same direction as the movement, theta spike correlations are symmetrically reversed (Figure 1I) as shown experimentally in CA1 neurons (Huxter et al., 2008; Yiu et al., 2022).
In area CA3, however, theta spike correlations are neither solely extrinsic (Yiu et al., 2022; Kay et al., 2020), since phase precession properties change in relation to running directions, nor are they solely intrinsic since reversal of correlation is still observed in most of the sequences (Huxter et al., 2008; Yiu et al., 2022). We therefore propose a new theory of phase precession for CA3 incorporating both intrinsic and extrinsic factors.
Directional sensory input
To, however, fully explain directional properties of theta phase precession and theta spike correlations by a model, also directional modulations of firing rates (Yiu et al., 2022) need to be taken into account.
We therefore included both directional and positional modulation of the sensory input to the model place cells (Figure 2A–B) with randomized preferred heading directions (see Methods section: Spatial input). The sensory input is assumed to arise from MEC, and hence, it is also thetamodulated and phaseshifted by ${70}^{\circ}$ with respected to the peak of theta cycle (Mizuseki et al., 2009). Furthermore, since the precession slope observed in the Romani and Tsodyks, 2015 model is limited (1.13 radians per field size, see Figure 1H) as compared to the experimental reports (4.44 radians (Yiu et al., 2022) and about 2.0 radians (Harris et al., 2002) per field size), we introduced shortterm synaptic facilitation (STF) to the sensory input (Berretta and Jones, 1996; Thurley et al., 2008) generating temporally asymmetric depolarization as suggested by intracellular recordings in vivo (Harvey et al., 2009; Figure 2B). STF amplifies the sensory current at the later part of the field, thus creating phase precession with steeper slopes thereby extending the phase range (see Methods section: Spatial input). Finally, we designated the synaptic weights to be stronger between place cells with similar preferred heading directions (Figure 2C) as has been proposed (Brunel and Trullier, 1998) as a result of Hebbian plasticity applied to directional firing fields.
A simulation of the place cell network was performed for a rightward trajectory through the arena based on our variant of the extrinsic Romani and Tsodyks, 2015 model (Figure 2D). We focus on two sets of place cells, one for which the trajectory aligns with the preferred heading direction of the field (red, denoted as best direction) and one for which the trajectory runs opposite the preferred heading direction (denoted as worst direction; Figure 2F). Phase precession has a lower onset and marginal spike phase along best direction than along the worst (Figure 2G–H), which is consistent with experimental data (Yiu et al., 2022 report mean spike phases ± SEM for best and worst direction of 1.61 ± 0.02 and 2.22 ± 0.03 in radians respectively), reflects that larger depolarizations generally yield shorter latencies. Directionality of the input, although it yields lower spike phases through higher depolarization, does not affect spike pair correlations, which remains solely extrinsic (Figure 2I). Thus, even though rate directionality and directional bias in recurrent connectivity can render phase precession directionally dependent, they are not sufficient to account for intrinsic sequences.
Generation of intrinsic sequences by the DGCA3 recurrent network
To explain the expression of intrinsic sequences in CA3, we propose them to be generated by the interaction of two networks, CA3 and DG (Figure 3A). DG is a good candidate region to be involved in phase precession, since lesions of it were shown to reduce prospective spiking (Sasaki et al., 2018) and to lower the onset phase of phase precession (Ahmadi et al., 2022). In our model, the neurons in DG receive excitatory synaptic inputs from CA3 place cells (putatively via hilar mossy cells) and project back to the CA3 cells with place field centers at a different location (Equation 3) to induce propagation of intrinsic sequences along a specific spatial direction. The CA3 cells at the target location of DG input are then activated and evoke higher depolarization in cells with place fields at the next DG target locations through the feedback. This scheme produces a rigid sequence whose activation order is independent of the movement direction. The connection pattern of DGCA3 projections (for brevity, we also refer to it as ‘DG loop’ in the subsequent text) could be determined by preexisting network structure or past experience through associative learning, or both.
Figure 3, provides schematic illustrations, for a DG layer that either only projects CA3 activity to their rightward neighbours (${\theta}_{\mathrm{DG}}={0}^{\circ}$, Figure 3A) or only to their leftward neighbors (${\theta}_{\mathrm{DG}}={180}^{\circ}$, Figure 3E). Simulations for both cases (${\theta}_{\mathrm{DG}}={0}^{\circ}$ and ${\theta}_{\mathrm{DG}}={180}^{\circ}$) assume a rightward trajectory. Apart from the addition of the DG layer, the model architecture and parameters of CA3 layer are the same as in Figure 2 (including best and worst direction in place field firing rate), which only generates extrinsic sequences through STD in the CA3 recurrent synaptic connections. DGloop connectivity is additionally modulated by firing rate directionality of the CA3 place fields. Fields with similar preferred heading directions are more strongly connected via the loop than those with opposite preferred heading directions (see Methods section: DG layer).
We found that, when the simulated animal is running in the same direction as the DGCA3 projection, phase precession starts from a higher phase (Figure 3C–D) due to the forward activation of place cells through DG layer (recovering the effect of asymmetric connectivity in the original Tsodyks et al., 1996 model), as compared to the model without DG layer (Figure 2G–H). Spike phases in best direction remain lower than along the worst direction (Figure 3D). When, however, the animal is running against the DGCA3 projection (Figure 3E), extrinsic sequences are still present in forward direction, evoked by the movement of the animal, but the intrinsic sequences are played out backward as determined by the direction of fixed recurrence (see cyan and yellow shaded regions in Figure 3F). The latter is reflected by the higher phase at the end of the phase position plots (Figure 3G) which leads to flatter precession slopes and decreases the fraction of phase precession (slope <0) of all traversal trials (Figure 3H). A closer look into pair correlation reveals that for trajectories opposite to the DGloop projection (${\theta}_{\mathrm{DG}}={180}^{\circ}$), spike probability is added to positive time lags (Figure 3I). Therefore, introducing fixed recurrence through DG loops elicits both extrinsic and intrinsic sequences and qualitatively changes theta sequences.
To quantify the degrees of extrinsic and intrinsic sequence firing in a way allowing comparison to experimental reports, we use the measures extrinsicity and intrinsicity (Yiu et al., 2022) that are based on pairs of place cells with overlapping place fields (see Methods section: Extrinsicity and intrinsicity, and Discussion). In our simulation, extrinsically and intrinsically driven cell pairs are both present in the population (Figure 3J), indicating a coexistence of extrinsic and intrinsic sequences. Our model reproduces a greater extrinsicity for cell pair activity when running direction aligns with both best place field directions as compared to when it aligns to both worst field directions, since along the best direction, cells receive more sensory depolarization, and thus, the movementdependent extrinsic sequences are more activated. The model also explains, why pairs with similar preferred heading directions may be less extrinsic than pairs with approximately opposite preferred heading direction (dissimilar pairs), since the DG loop preferentially connects CA3 place cells with similar preferred heading directions. Both of the results follow the same trend as found in experimental data (Yiu et al., 2022 report ratios of extrinsically to intrinsically correlated CA3 field pairs of 1.57 for bothbest directions, 0.41 for bothworst directions, 0.87 for similar pairs and 2.43 for dissimilar pairs).
Thus, by introducing feedback excitation via the DG layer, intrinsic sequences are able to propagate in fixed directions on top of the movementdependent extrinsic sequences. Theta sequence directionality is reflected through the change in spike correlation, which varies as a function of the difference between the direction of DG feedback and movement direction. The combination of extrinsic and intrinsic theta sequence activity is robust regarding changes of the parameters of shortterm synaptic plasticity (Figure 3—figure supplement 1), as well as running speed (Figure 3—figure supplement 2), as long as place fields are wide enough to allow spikes in sufficiently many theta cycles.
Lesion of DG reduces theta compression and phase precession range
One prediction of the DGloop model, consistent with findings from DG lesion experiments (Ahmadi et al., 2022), is that DG would contribute to the temporal organization of spike sequences in CA3. To verify this hypothesis also in the model, we implemented a lesion of DG by disabling activity in the DG layer. To compensate for reduced excitatory drive caused by the lesion, we then increased probability of release of the sensory inputs thereby increasing the initial input amplitudes but removing shortterm synaptic facilitation (Figure 4A).
We found that a DG lesion would reduce theta compression in sequence activity. Theta compression (Dragoi and Buzsáki, 2006) refers to the compression of secondslong behavioural experience of placefield crossing into a neural representation of spike sequences at a shorter (theta) timescale. To quantify the strength of theta compression, we plotted the pair correlation lags versus the distance between the centers of two fields (abbreviated as ‘lags’ and ‘pair distance’ respectively), after simulating a rightward trajectory (Figure 4B–C). The magnitude of the linearcircular regression slope $a$ measures how much theta phase encodes a certain interval in space, and therefore, the strength of theta compression. As a result, theta compression is reduced for the DGlesioned case ($a=0.053$ radians/cm for both ${\theta}_{\mathrm{DG}}={0}^{\circ}$ and ${\theta}_{\mathrm{DG}}={180}^{\circ}$), as compared to the control case ($a=0.183$ radians/cm for ${\theta}_{\mathrm{DG}}={0}^{\circ}$ and $a=0.059$ radians/cm for ${\theta}_{\mathrm{DG}}={180}^{\circ}$) reproducing the finding (Ahmadi et al., 2022) that spatial encoding via theta sequences crucially depends on intact DG and suggesting that the loss of DG inputs could be compensated for by the increase of release probability in the spared afferent synapses from the MEC. The DG lesion also reduces the spike phase and phase range of phase precession (Figure 4D), which indicates the participation of DG loops in highphase spiking. Both weaker phase precession and theta compression stress the important role of DG in temporal organization of CA3 sequences.
Theta sequences in 2d and outoffield firing
So far, the model was only evaluated on bidirectional linear tracks, where running directions completely overlapped with the orientation of the DG loop connectivity. Now, we extend our analysis to 2d space by examining oblique trajectories which cross the orientation of DGloop projection at certain angles.
We first arrange the DGloop connections such that the DGloop orientation crosses a rightward trajectory at 45° and 225° (Figure 5A–F). Similar to the cases of ${\theta}_{\mathrm{DG}}={0}^{\circ}$ and ${\theta}_{\mathrm{DG}}={180}^{\circ}$ (Figure 3D and H), precession slopes are steeper and onsets higher when the trajectory direction aligns more with the orientation of the DGloop, but with a smaller effect size for oblique crossings (Figure 5A and C) since DGloop connectivity area only overlaps with part of the trajectory near the intersection. We further resolve the precession slope, onset and marginal phase for each place cell into 2d maps (Figure 5B and D). Intrinsic sequences with a higher marginal spike phase can be clearly seen along the belt of DGloop projections and are even extended to the outside of trajectory predicting ‘offtrack’ spikes at high phases. Depending on the alignment between movement direction and DGloop orientation, the slope becomes either more negative (${\theta}_{\mathrm{DG}}={45}^{\circ}$) or more positive (${\theta}_{\mathrm{DG}}={225}^{\circ}$). Analysis of extrinsicity and intrinsicity was conducted for all field pairs and confirmed the same trend as in Figure 3 that best and dissimilar pairs are more extrinsic than worst and similar pairs, respectively (Figure 5E). As a quantitative prediction, we computed the angle differences between field centers of cell pairs for the extrinsic and intrinsic populations, and observe that place field center differences in extrinsically correlated field pairs are mostly oriented horizontally (along the running direction) while place field center differences from intrinsically correlated field pairs are oriented along the DGloop orientation ${\theta}_{\mathrm{DG}}={45}^{\circ}$, as by design (Figure 5F).
The analysis above is repeated for the geometric configurations that DGloop connectivity is minimally interacting with the place cell activity induced by movement, that is when DGloop orientation and the movement direction are perpendicular to each other (${\theta}_{\mathrm{DG}}={90}^{\circ}$ and ${\theta}_{\mathrm{DG}}={270}^{\circ}$, Figure 5G). Similar effects as in Figure 5B and D on precession slope, onset and marginal phases are also observed in the 2d map, except that the effects are further restricted to the intersection area in the middle. Also, the whole population has become more extrinsic as compared to the 45° and 225° cases (Figure 5H, see the numbers for extrinsicintrinsic ratios) due to the smaller overlapping area between DGloop projection and the trajectory. Lastly, the pair center difference orientation confirms that field pairs with extrinsic correlations follow the trajectory direction while those with intrinsic correlation are biased towards the DGloop orientations (90°).
The results demonstrate the distinct roles of extrinsic and intrinsic sequences in 2d spatial encoding. The former represents trajectory direction while the latter the associative memory towards specific locations. They can be played out at the same time separately in different directions and only interact with each other when they overlap. The interaction is reflected in directional dependence of phase precession properties, most notably the higher spike phases from the DGCA3 recurrent input, as well as increased intrinsicity of pair correlation and extended firing fields along the orientation of the DGloop projections. Intrinsic sequences also triggered outoffield firing (Figure 5B, D and G) at late theta phases. In this case, the DGloop connects to cells with remote place fields. These cells could even display multiple separated place fields, with high spike phases indicating the target location of the intrinsic sequence.
Topologyfree mechanisms of extrinsic phase precession
A wellknown problem of phase precession models based on recurrent connectivity that applies to both, the original intrinsic Tsodyks et al., 1996 and Romani and Tsodyks, 2015 model, is that they do not explain how the topological connectivity matrix (in our case ${W}^{\mathrm{CA3}}$) is generated (Lisman and Redish, 2009; Jaramillo and Kempter, 2017). Extrinsic theta sequences in a first exposure to a novel environment should therefore be missing. Although Feng et al., 2015 find that theta sequences on a first exposure of a novel linear track are indeed much weaker (maybe only reflecting intrinsic sequences), their results nevertheless indicate a very fast learning time scale that is hard to reconcile with recurrent learning of a full spatial topology (particularly the generalization to 2d). Also their result might be hampered by place field plasticity that biases the decoder towards backwardshifted place maps of later trials (ParraBarrero and Cheng, 2023). We therefore explored, whether extrinsic 2d sequences could also be generated by a model that is not relying on 2d topology in the recurrent weights. To this end, we disabled the CA3 recurrence and compensated the missing level of excitation by an increased strength of the spatial input and the DG loops (see Methods section: Parameters of the models). Our simulations show that extrinsic sequences can still be generated by spatial input alone (Figure 6A), relying only on the shortterm facilitation mechanism. Simulating CA3 activity with lesioned DG similarly abolishes the temporal organization of theta sequence and reduces the phase range (Figure 6B). The results demonstrate that the temporal order of extrinsic sequences could be coordinated solely by sensorimotor drive and does not necessarily require CA3 recurrence.
Functional role of intrinsic sequences
While the function of extrinsic theta sequences in encoding the actual trajectory of an animal (connecting the recent past, present and near future locations) is obvious, the potential role of the less readily apparent intrinsic sequences is not straight forward. Simulations of trajectories in 2d (Figure 5) suggest intrinsic activity may serve a role to identify certain locationdirection pairs independent of the current trajectory. Here, we follow this idea by evaluating the hypothesis that the intrinsic sequences signal a stable ‘landmark’ (location/direction pair) cue by a temporal code that is invariant to different directions of approach.
To test our hypothesis, we constructed a downstream readout neuron that would reliably identify the presence of the intrinsic sequence independently of the animal’s running direction, whereas it would not be able to do so for only extrinsic sequences. To this end, we trained the synaptic weights using the tempotron learning rule (Gütig and Sompolinsky, 2006), which is able to implement binary classification based on temporal relations of input spike patterns (see Methods section: Tempotron). Two tempotrons were trained to recognize the spike patterns from the place cells, one taking input from a model with DGloop connectivity at ${\theta}_{\mathrm{DG}}={0}^{\circ}$, and one without DGloop connectivity to serve as a control only having access to extrinsic sequences (Figure 7A). Nonmoving spatial inputs were applied to the CA3 place cells at the centers of withloop and noloop populations and their spike patterns in subsequent theta cycle were used as training patterns, mimicking a situation in which network activity is evoked without sensorymotor input as, for example in a offline situation before the animal walks or maybe even has seen the environment. The training patterns have only (+) labels, which the tempotrons are trained to recognize by firing a spike (Figure 7B). We then test the tempotrons with spike patterns induced by the animal running on different trajectories through the trained location with running directions varying between 0° to 360° (Figure 7C). All spikes in each training and testing pattern are individually jittered by adding a noise term $\sigma \sim \mathcal{N}(0,(2\text{}\mathrm{m}\mathrm{s}{)}^{2})$ for 100 times, producing 100 samples for each pattern. The tempotron is said to successfully recognize the sequence of a trajectory direction if any of the theta cycles throughout the trajectory elicits a spike, that is while running, the readout cell would evaluate the place cell sequence in every theta cycle for information on the trained landmark.
We found that the tempotron trained on the intrinsic sequence from the DG loop is able to recognize the sequence patterns produced for all running directions, while the tempotron trained without a DGloop fails to identify the extrinsic sequences most of the time (see accuracies in Figure 7D). The reason is that the spike patterns induced by intrinsic sequences remain similar to the training pattern despite being approached in other directions (see sequential contributions in Figure 7E), while spike patterns for the noloop condition are different between training and testing (Figure 7F). The distinction is further illustrated in Figure 8, where 2d maps of spike time gradients in one theta cycle are plotted with respect to running direction/training condition. Sequences always contain components that propagate along the projection direction of DG loop, while, without such a loop, they only propagate along the running direction. Moreover, during training, the noloop condition evokes concentric waves, reflecting the 2d topology of the recurrent weights. Similar results were also achieved without CA3CA3 recurrence (Figure 7—figure supplement 1 and Figure 8C).
Our results show that intrinsic sequences can provide a stable correlation signal which allows reliable decoding of locations through temporal correlations. The intrinsic temporal code remains detectable even when mixed with extrinsic sequences.
Discussion
We presented a model of hippocampal theta sequences in 2d environments, suggesting that both extrinsic and intrinsic mechanisms are required to explain experimental reports that phase precession and spike timing correlations are nonhomogeneous across running directions. Although phase precession already becomes directional by including directiondependent sensory input into a purely extrinsic model, directionality of spike timing correlations cannot be explained by such a model. We, however, demonstrated that the correlation preference could be implemented by fixed recurrent loops via a model DG layer. We further supported the model assumptions by showing that DG lesions plus compensatory sensory drive can abolish the theta compression effect in CA3 spiking activity (Ahmadi et al., 2022). By employing a spikebased temporal pattern decoder (tempotron), we showed that the intrinsic sequences could function as stable signatures that act as anchors of the spatial code.
Early intrinsic models (Tsodyks et al., 1996) were challenged owing to their inability to generate phase precession in backward travel (Figure 1E, also see Cei et al., 2014), as well as the predominantly extrinsic correlations observed in CA1 (Huxter et al., 2008). In our hybrid model, phase precession still occurs during backward travel (${\theta}_{DG}={180}^{\circ}$) but at a lower probability as indicated by the larger fraction of positive phaseposition slopes (Figure 3H). Also, extrinsic sequences still dominate over intrinsic sequences as indicated by the majority of field pairs being extrinsic (Figure 3J). Both the reduced expression of phase precession in reverse runs and the dominance of extrinsic sequences are in accordance with the experimental data (Yiu et al., 2022).
The mixture of extrinsic and intrinsic mechanisms in our theory, naturally gives rise to the directionality of spike correlations and phase distributions. As the trajectory aligns itself with the DG loops, the ratio of intrinsic to extrinsic sequences increases. As a result, spike correlations become more rigid and the phase distribution is shifted upward due to the accumulated synaptic transmission delay from the reverberating activity between CA3 and DG populations. Adding directional sensory input activates extrinsic sequences in the best direction more strongly, and hence, leads to an association between bestangle (worstangle) pairs and extrinsicity (intrinsicity). These predictions of our model are corroborated by past reports of higher spike phases in the nonpreferred arm of a Tmaze (Kay et al., 2020) as well as the association of rigid correlations with upward shifts in spike phases and an increase in worstangle pairs (Yiu et al., 2022). The experimental distinction between extrinsic and intrinsic components in theta sequences has so far only been achieved in pairs of place cells, owing to the limited number of simultaneously recorded place cells with overlapping fields. Our model predicts that similar distinctions should also be observable in higherorder statistics, obtained from overlapping fields of a larger number of cells. Instead of correlation lags, we suggest to use temporal pattern detection methods (e.g. Chenani et al., 2019) to unveil the respective sequence contributions.
Since intrinsic sequences can also propagate outside the trajectory (outoffield firing in Figure 5) and activate place cells nonlocally, our model predicts directiondependent expansion of place fields, or even multiple place fields, with the intrinsic sequence’s target location exhibiting late spike phases and higher phase precession onsets. Remote activation during locomotion has already been observed in a previous study (Sasaki et al., 2018) where CA3 place cells preferentially firing at one arm of the maze were also activated at reward locations at other arms. In our model, only shortrange intrinsic connectivity was considered, thus, place field boundaries expand locally but in a skewed manner matching the sequence direction. Skewness of place fields has been reported by a number of studies (Mehta et al., 1997; Shen et al., 1997; Mehta et al., 2000; Ekstrom et al., 2001; Lee et al., 2004; Burke et al., 2008; Cei et al., 2014; Roth et al., 2012; Dong et al., 2021) showing place fields to be asymmetrically expanded opposite to the direction of travel. This effect was connected to plasticity as it develops after repeated traversal, and due to its dependence on NMDA receptor activation (Ekstrom et al., 2001; Burke et al., 2008; Shen et al., 1997). These plasticity studies show that the hippocampal place code is shaped by intrinsic synaptic computations including temporal activation patterns in theta sequences (Feng et al., 2015). Apart from being conducted on linear tracks and not 2d environments, most of this work focused on CA1 and associated Schaffer collateral plasticity. Yet some prior studies (Lee et al., 2004; Roth et al., 2012) did show that place fields in CA3 were more skewed than in CA1, which our model would explain by CA3 expressing more intrinsic sequences than CA1 consistent with prior experimental observations (Yiu et al., 2022 reported ratios of extrinsically to intrinsically driven cell pairs of 1.44 in CA1 and 1.23 in CA3).
A further prediction of hardwired DG loops is that the resulting activity patterns (intrinsic sequences) should not remap under conditions of global or partial remapping (Leutgeb et al., 2004). Instead the same intrinsic sequence components should be observable in multiple environments, however, they might only be seen in a small fraction and thus this prediction is potentially hard to test.
The backprojection from CA3 to DG is a crucial anatomical prerequisite of our model, but was rarely explored compared to the feedforward inputs via the perforant pathway. The proposed CA3DG recurrent structure of this model, albeit simplified, is consistent with the anatomical evidence. Pyramidal cells in CA3 innervate the mossy cells at the DG hilus (Scharfman, 1994; Scharfman, 2016), which then project to granule cells through both excitatory and inhibitory pathways (Hsu et al., 2016; Scharfman, 1995; Larimer and Strowbridge, 2008; Soriano and Frotscher, 1994), and subsequently back to CA3 pyramidal cells. An optogenetic study (Hsu et al., 2016) showed that the net effect of mossy cells on granule cells was predominantly inhibitory, suggesting that the DG ensembles excited by mossy cell synaptic drive are sparsified by suppressing unwanted outofensemble activity. Indeed, past studies showed that reliable excitatory effect could be observed when granule cells were depolarized (Scharfman, 1995) and when they received backpropagation of sharp wave bursts from CA3 population (Penttonen et al., 1998). This indicates that the excitatory recurrent pathway from CA3 via DG exists and might allow activity reverberation between two layers. While our model, owing to its simplicity and generality does not require any DG specific pathways and would work equally well with any other anatomical interpretation of the CA3 feedback, we hypothesize the intrinsic feedback connectivity to arise via the DG, particularly because DG lesions were shown to eliminate the coordinated temporal structure of CA3 activity and to be instrumental to sequence organization (Figure 4 and Ahmadi et al., 2022).
Our model assumed a connectivity pattern in the DG loops, in which neurons activate the neighbours along a specific direction, as inspired by Hebb’s phase sequences (Hebb, 1949) and, hence, replay of the loop would activate a spatially plausible virtual trajectory. The loop connectivity could either arise from previous learning, or might be present already beforehand (Dragoi and Tonegawa, 2013), with spatial topology inherited by associating 2d sensory features to cell ensembles in the loop (Leibold, 2020). The resulting topology can exhibit discontinuous longrange jumps to other locations (Sasaki et al., 2018) or consist of a discrete set of (behaviorally relevant) locations (Pfeiffer, 2022).
Different from other phase precession models, we also included heading direction as part of the sensory input, as inspired by past literature that CA1 (Markus et al., 1995; Acharya et al., 2016; Stefanini et al., 2020), CA3 (Mankin et al., 2019), and DG place cells (Stefanini et al., 2020) exhibit directional selectivity in firing rates, potentially inherited from the upstream headdirection cells in the medial entorhinal cortex (Giocomo et al., 2014) and postsubiculum (Taube et al., 1990). As a result, the directional drive immediately translates to phase directionality in theta sequences, partly contributing to the upward shift of the phase distribution in the worst angles. Such phase directionality arises naturally from the intracellular dynamics of a spikebased model, where stronger depolarization causes earlier spiking. This phaserate dependence has already been used in previous models (Harris et al., 2002; Mehta et al., 2002; Thurley et al., 2008), where the increasing depolarization within place fields directly relates to decreasing spike phases. The causal effect of firing rate on spike phases, however, was disputed by Huxter et al., 2003 as they showed that precession slopes and spike phases remained the same between high and lowspiking runs, suggesting that the phase is not singlehandedly determined by firing rate. In our model, firing rate is determined by both lowphase spiking from sensory input and highphase spike arrivals of DGCA3 loops, both producing opposing effects on the phase distribution. Thus, depending on the strength and geometry of the DGCA3 connectivity, spike phases are not fully determined by firing rate.
By using a tempotron to decode the spike patterns, we show that the spike patterns of intrinsic sequences can serve as a stable landmark which remains decodable across multiple running directions. The invariant temporal patterns could serve as anchors of spatial memories in a novel environment, since place fields only stabilize after the animal becomes familiar with the environment (Wilson and McNaughton, 1993). The preexisting sequence motifs, even at times when the spikes of the neurons are not spatially tuned to a location, can still encode the position based on their temporal relations alone. The idea has previously been spelled out (Cheng, 2013) and numerically verified (Leibold, 2020; ParraBarrero and Cheng, 2023) with multiple fixed sequences that form a decodable spatial representation.
Intrinsic sequences may thus act as a scaffold around which a new spatial code can be built for new but similar behavioral contexts, where similarity could for example be identified by a salient feature. Once the behavioral context of a situation changes, new intrinsic sequences would be observable. These intrinsic landmarks need to be stable across time, as shown for some dentate gyrus representations (Hainmueller and Bartos, 2018). We speculate that offline sequences observed during replay and preplay (for review see Buhry et al., 2011; Dragoi and Tonegawa, 2014), would correspond to the intrinsic activity patterns and indicate the context expectation of an animal (which can be detected by a tempotron). The functional roles of intrinsic sequences may thus not be limited to spatial memories. While, in the spatial domain, intrinsic sequences could be interpreted as spatial trajectories (Kay et al., 2020; Sasaki et al., 2018), virtual nonspatial trajectories could represent working memories contents (Jensen et al., 1996) available for general decision making processes.
Methods
Neuronal model
Generation of neuronal action potentials is modelled according to Izhikevich, 2003. The soma potential $v$ and the adaptation variable $u$ of unit $i$ at time $t$ (in ms) follows the equations:
Any time $v(t)$ crosses the threshold 30 mV from below, we register a spike for the neuron and reset the soma potential by $v(t)\leftarrow c$ and the adaptation variable by $u(t)\leftarrow u(t)+d$. For the excitatory pyramidal place cells, we use parameters $a=0.035$, $b=0.2$, $c=60$ mV, $d=8$, which provides the neuron with burst firing characteristics. For the inhibitory interneurons, the parameters were $a=0.02$, $b=0.25$, $c=65$ mV, and $d=2$, which corresponds to fast spiking patterns. $I(t)$ is the total sum of recurrent ${I}^{R}(t)$, sensory ${I}^{S}(t)$ and oscillatory theta input
We chose to use the phenomenological spike generation model of Izhikevich, 2003, since it allows to adjust burst firing properties with only few parameters that efficiently emulate the bifurcation structure of spike generation. Synaptic integration below threshold is not affected by the spike generation model and will thus be treated by conventional synaptic models.
Spatial input
The place field centers ${\mathbf{p}}_{i}^{CA3}=[{x}_{i}^{CA3}(t),{y}_{i}^{CA3}(t)]$ of $80\times 80=6400$ excitatory CA3 cells equally tile the 80 by 80 cm square arena. Place cell firing rates are modelled directionsensitive, with preferred heading directions ${\psi}_{i}^{CA3}$ semirandomized among each 2×2 tile of place cells by randomly rotating a set of four equally spaced direction angles by a uniformly distributed angle $\xi $, that is
The sensory input ${J}_{i}^{S}(t)$ into the $i$th neuron depends on the instantaneous position, $\mathbf{p}(t)=[x(t),y(t)]$, and heading direction $\psi (t)$ of the animal as
where ${\mathrm{A}}^{\mathrm{pos}}$ is the amplitude of positional tuning and $d(\cdot )$ computes the Euclidean distance between two positions. The positional tuning curve is implemented as a rectangular box function, where the place cell only receives sensory input if the animal is within 5 cm from the field center. Directional tuning is implemented as an additional amplitude gain ${\mathrm{A}}^{\mathrm{dir}}$ to the positional current depending on the circular difference between the animal’s heading and the neuron’s preferred heading direction ${\psi}_{i}^{\mathrm{CA3}}$. The sensory input is assumed to be modulated by theta oscillations from medial entorhinal cortex (MEC) ${I}^{\mathrm{MEC}}(t)$ with a phase shift of 70° (Mizuseki et al., 2009).
The sensory input ${J}^{S}$ is subsequently transformed to the input current ${I}^{S}$ via shortterm facilitation (STF)
where the facilitation variable ${s}_{i}^{F}$ decays to ${\mathrm{S}}_{0}^{F}$ with a time constant $\tau}^{F}=500\text{}\mathrm{m}\mathrm{s$ and increases to ${\mathrm{S}}_{1}^{F}$ when the sensory input ${J}_{i}^{S}$ is present. The time constant ${\tau}^{F}$ of facilitation of neocortical synapses was in the range suggested by Tsodyks et al., 1998 following previous experimental reports (Mejías and Torres, 2008; Zucker and Regehr, 2002). ${\mathrm{\Phi}}^{F}$ controls the strength of the STF. The facilitation variable is squared to include nonlinear interactions in presynaptic calcium dynamics. As a result, facilitated sensory input ${I}_{i}^{S}$ increases over time and becomes stronger in the later part of the field, thus effectively generating a spatially graded input strength.
Note that only the CA3 place cells receive the sensory input. ${I}_{i}^{S}(t)$ is not applied to the place cells in DG and all of the inhibitory interneurons.
CA3 recurrent connections
Place cells in CA3 connect with each other by excitatory synapses. The excitatory synaptic current ${I}_{i}^{E}(t)$ is conductancebased, and follows the equations:
The conductance ${g}_{i}^{E}$ of a postsynaptic cell $i$ is increased by the spike arrivals at times ${t}_{j}^{(f)}$ from the presynaptic cell $j$, and decay with a time constant $\tau}^{E}=12\text{}\mathrm{m}\mathrm{s$. ${N}_{J}=6,400$ is the number of presynaptic place cells, $V}^{E}=0\text{}\mathrm{m}\mathrm{V$is the reversal potential of the excitatory synapses and $\tau}_{0}=2\text{}\mathrm{m}\mathrm{s$ is the synaptic transmission delay.
The synaptic weights ${W}_{ij}$ from cell $j$ to cell i depend on the distance between place cell centers and on the similarity of their preferred heading angles, i.e.,
where ${B}^{\mathrm{pos}}$ and $\sigma =2\text{}\mathrm{c}\mathrm{m}$ correspond to the maximum strength and width of the locationspecific interaction, respectively. ${B}^{\mathrm{dir}}$ and ${K}^{\mathrm{CA3}}$ control the maximum strength and the concentration of the directional dependence, respectively. ${J}_{ij}$ models the rightward asymmetry of the cell connections, which was only turned on when we simulated the 2d variant of Tsodyks et al., 1996 model in Figure 1C–F and otherwise turned off in the rest of our analysis.
Furthermore, the recurrent synaptic conductances underwent shortterm synaptic depression (STD), as was proposed in Romani and Tsodyks, 2015 to serve as sequence generator in 2d space. The mechanism penalizes the recurrent input into the place cells behind the animal. As a result, the differential recurrence strengths translate to a gradient of spike phases and produces extrinsic sequences in the direction of travel. We model the STD by the variable ${s}_{i}^{D}(t)$ which represents the available synaptic resource and follows the dynamics:
where ${s}_{i}^{D}$ recovers to 1 with a time constant ${\tau}^{D}=500$ ms and is depleted by a fraction ${U}^{D}$ every time a spike occurs. The recovery time constant is comparable to experimentally obtained values of cortical neurons (200–800 ms in Tsodyks and Markram, 1997; Markram et al., 1998; Abbott et al., 1997; Zucker and Regehr, 2002) and previous modelling work (450–800 ms in Romani and Tsodyks, 2015; Haga and Fukai, 2018; Tsodyks and Markram, 1997; Tsodyks et al., 1998). The STD only applies to synaptic connections when presynaptic cells are CA3 place cells. ${s}_{i}^{D}(t)$ is fixed at 1 when the presynaptic cells are inhibitory interneurons or DG place cells.
DG layer
We simulated ${N}_{\mathrm{DG}}=40\times 40=1600$ place cells in the DG layer, with place field centers equally tiling the environment. The DG cells do not receive sensory input. Their positional (${x}_{i}^{\mathrm{DG}},{y}_{i}^{\mathrm{DG}}$) and directional (${\psi}_{i}^{\mathrm{DG}}$) tunings are determining synaptic strengths to and from the CA3 layer. The directional tuning is semirandomized as described for CA3. The synaptic current dynamics follow Equations (1) and (2). Excitatory synaptic weights from CA3 place cells to DG place cells are defined as
which are dependent on the differences in the place field centers and preferred heading angles between the CA3 and DG populations. The variable ${C}_{j}^{\mathrm{CA3}}$ strengthens outgoing connections from CA3 place cells on the path corresponding to the intrinsic sequence by choosing
where $\mathbf{p}}_{k}^{C$ varies with the intrinsic path direction ${\theta}_{DG}$ as ${\mathbf{p}}_{k}^{C}=[2k\phantom{\rule{thinmathspace}{0ex}}\mathrm{cos}({\theta}_{\mathrm{D}\mathrm{G}}),2k\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}({\theta}_{\mathrm{D}\mathrm{G}})]$.
The excitatory synaptic strengths from DG to CA3 are chosen such that DG cells project back to CA3 cells with place field centers shifted by a vector $\mathbf{r}=[4\mathrm{cos}({\theta}_{\mathrm{D}\mathrm{G}}),4\mathrm{sin}({\theta}_{\mathrm{D}\mathrm{G}})]$ of fixed length of 4 cm along the intrinsic path, that is
The model has no synaptic connections between DG excitatory neurons.
Inhibitory synapses
The model additionally contains ${N}_{I}=250$ inhibitory interneurons (denoted as Inh) each for the CA3 and the DG layer. They provide inhibitory feedback separately to the excitatory cells within each layer (CA3InhCA3 and DGInhDG). The dynamics of their synaptic currents mirrors the excitatory synapses, that is
with $\tau}^{I}=10\text{}\mathrm{m}\mathrm{s$, $V}^{I}=80\text{}\mathrm{m}\mathrm{V$. CA3 and DG have alltoall connections to their inhibitory populations with uniformly randomized strengths, i.e. ${W}_{ij}^{XY}={W}_{0}^{XY}\xi $ , with $\xi \sim \mathcal{U}(0,1)$ is the maximum synaptic strength, and the notation XY corresponds to InhCA3 and InhDG connections. There is no synaptic connection between interneurons, that is ${W}^{\mathrm{Inh}\mathrm{Inh}}=0$.
The total recurrent current entering each excitatory neuron is thus the sum of the excitatory and inhibitory current:
Excitatory synapses to interneurons
Interneurons only receive alltoall excitatory currents from their respective layer. Those currents are modelled according to Equations 1; 2. The synaptic weights are constant and denoted by ${W}_{0}^{\mathrm{CA3}\mathrm{Inh}}$ and ${W}_{0}^{\mathrm{DG}\mathrm{Inh}}$.
Parameters of the models
Model parameters that are adjusted in different analyses are listed in Table 1. The values of the synaptic weights and spatial input were chosen to allow for a large range of phase precession and stability of the network activity. For the analyses including the DG layer, weights are adjusted to allow coexistence of extrinsic and intrinsic sequences.
Crosscorrelation analysis
Crosscorrelation represents the probability that a spike of one place cell would occur following a certain time lag from the spike of the another cell. Crosscorrelation is always empirically computed as a histogram of time lags between spike pairs with a resolution of 5ms in a window of 200ms. Throughout the present study, the direction of a time lag is designated as the lag of the first encountered cell relative to the next cell along the trajectory, except in Figure 1, where the direction of time lag follows the cell order along the ${0}^{\circ}$, ${45}^{\circ}$, and ${90}^{\circ}$ trajectory in each comparison group, and in Figure 2, where the time lag direction is from left to right cells.
Correlation lag is derived by bandpass (4–12 Hz) filtering the crosscorrelation histogram and applying a Hilbert transform on the filtered signal. The phase of the analytic signal at time lag 0 is the correlation lag.
Extrinsicity and intrinsicity
We apply quantitative measures for the extrinsic or intrinsic nature of crosscorrelations in a pair of place fields following Yiu et al., 2022. We compare the crosscorrelation histograms of a field pair for a running direction along the DG loop (${\theta}_{DG}$) and opposite to the loop (${\theta}_{DG}+{180}^{\circ}$). Extrinsicity (Ex) is computed as the Pearson’s correlation ($r$) between two crosscorrelation histograms, and intrinsicity (In) between the histogram of ${\theta}_{DG}$ and the horizontally flipped histogram of ${\theta}_{DG}+{180}^{\circ}$. The Pearson’s correlation is then transformed (${r}^{\prime}=(r+1)/2$) to be in the range of 0 and 1. An extrinsic correlation would give an extrinsicity near 1, since the effect of DG loop is minimal and correlation histograms are similar in both ${\theta}_{DG}$ and ${\theta}_{DG}+{180}^{\circ}$ directions. A pair of place fields with intrinsic correlation would see crosscorrelation horizontally flipped in the ${\theta}_{DG}+{180}^{\circ}$ condition due to the large effect of DG loop, and hence, give an intrinsicity near 1. We classify a pair as extrinsic if its extrinsicity exceeds intrinsicity, and vice versa.
Tempotron
A tempotron is a neuronally inspired classifier (readout neuron) whose dendritic synaptic weights can be adapted to recognize temporal patterns of spikes arriving at the afferents (for details, see Gütig and Sompolinsky, 2006). Briefly, the soma potential of the tempotron follows the equations
where w_{i} is the adaptable weight of the afferent fiber conveying spikes from place cell $i$ to the tempotron. $K(t{t}_{i}^{(f)})$ is a postsynaptic potential (PSP) kernel with decay and rising time constants of $\tau =5\text{}\mathrm{m}\mathrm{s}$ and $\tau}_{r}=1.25\text{}\mathrm{m}\mathrm{s$ respectively. V_{0} is a factor which normalizes the PSP kernel to 1. A spike is said to occur if $V(t)$ crosses the firing threshold ${V}_{\mathrm{\Theta}}=2$ from below. After threshold crossing, the afferents will be shunted and spike arrivals will not evoke more PSPs for the rest of the pattern. A pattern is defined as the set of spike times of all the presynaptic place cells in a theta cycle (100 ms).
The weight w_{i} follows the update rule
where t_{max} is the time at the peak of the soma potential $V(t)$. The learning rule assigns credit to the afferents based on spike timing. Spike times closer to the peak are considered to have higher contribution to the tempotron firing, hence their afferents are incremented by a larger step. After training, spike times with similar temporal correlations as the (+) patterns would be able to evoke enough PSP in the tempotron’s soma and elicit a spike as a positive response of binary classification, while those similar to () patterns would not elicit a spike from the tempotron.
We trained the tempotrons to identify the spike patterns of place cells at locations with and without intrinsic connectivity separately. To this end, we modified our network such that DG loops are present at the upper half of the arena, spanning the space from x=20cm to x=+20 cm at y=+20 cm in direction ${\theta}_{DG}={0}^{\circ}$, while the loop is absent in the lower half of the arena.
During training, we applied ‘nonmoving’ spatial inputs to the CA3 place cells at the withloop (0 cm, 20 cm) and noloop (0 cm, –20 cm) locations for 1 s, as if the animal were standing still at the locations, evoking the activities representing the two location cues. For computational efficiency, we restricted our analysis to the populations of CA3 place cells within the 20 cm squared boxes centered at the two locations. Each population contains 400 presynaptic cells, forming the input space for the tempotron. The spikes from the withloop population will train the first tempotron and those from the noloop population will train the second tempotron. Prior to training, the input spikes are subdivided to 10 patterns based on their theta cycles. Each pattern has a window of 100 ms. We added noise to the patterns by jittering the spikes with Gaussian noise $\mathcal{N}\sim (0,(2\text{}\mathrm{m}\mathrm{s}{)}^{2})$ for 100 times. As a result, each tempotron receives 10×100 = 1,000 training patterns from the activity evoked by the location. All training patterns are (+) patterns and there is no () pattern.
After training, trajectories (20 cm long, 1 s duration) with running directions from 0° to 360° with 15° increment were simulated to cross each of the locations. The trajectories produce a mix of extrinsic and intrinsic sequences in the withloop population and only extrinsic sequences in the noloop population. The patterns evoked by the running trajectories were separately applied to the tempotrons. The input spikes for testing were also subdivided into theta cycles and jittered in the same manner as during training, forming 1000 testing patterns for each running direction. A sequence is said to be correctly identified if the tempotron fires at at least 1 out of 10 theta cycles along the trajectory. The accuracy rate for each running direction of trajectory is computed across the 100 jittered realizations.
Code availability
We used Python 3 for simulations and visualization. The codes are available from a github repository (https://github.com/yyhhoi/directionalnet, copy archived at Yiu, 2023).
Data availability
The current manuscript is a computational study, so no data have been generated for this manuscript. We used Python 3 for simulations and visualization. The codes are available from https://github.com/yyhhoi/directionalnet (copy archived at Yiu, 2023).
References

Reactivation, replay, and preplay: how it might all fit togetherNeural Plasticity 2011:203462.https://doi.org/10.1155/2011/203462

Reversed theta sequences of hippocampal cell assemblies during backward travelNature Neuroscience 17:719–724.https://doi.org/10.1038/nn.3698

The CRISP theory of hippocampal function in episodic memoryFrontiers in Neural Circuits 7:88.https://doi.org/10.3389/fncir.2013.00088

Selection of preconfigured cell assemblies for representation of novel spatial experiencesPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 369:20120522.https://doi.org/10.1098/rstb.2012.0522

Topography of head direction cells in medial entorhinal cortexCurrent Biology 24:252–262.https://doi.org/10.1016/j.cub.2013.12.002

The tempotron: a neuron that learns spike timingbased decisionsNature Neuroscience 9:420–428.https://doi.org/10.1038/nn1643

Simple model of spiking neuronsIEEE Transactions on Neural Networks 14:1569–1572.https://doi.org/10.1109/TNN.2003.820440

Phase precession: a neural code underlying episodic memory?Current Opinion in Neurobiology 43:130–138.https://doi.org/10.1016/j.conb.2017.02.006

Quantifying circularlinear associations: hippocampal phase precessionJournal of Neuroscience Methods 207:113–124.https://doi.org/10.1016/j.jneumeth.2012.03.007

Nonrandom local circuits in the dentate gyrusThe Journal of Neuroscience 28:12212–12223.https://doi.org/10.1523/JNEUROSCI.361208.2008

Prediction, sequences and the hippocampusPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 364:1193–1201.https://doi.org/10.1098/rstb.2008.0316

The hippocampal code for space in Mongolian gerbilsHippocampus 29:787–801.https://doi.org/10.1002/hipo.23075

Interactions between location and task affect the spatial and directional firing of hippocampal neuronsThe Journal of Neuroscience 15:7079–7094.https://doi.org/10.1523/JNEUROSCI.151107079.1995

The role of synaptic facilitation in spike coincidence detectionJournal of Computational Neuroscience 24:222–234.https://doi.org/10.1007/s1082700700528

Learning to predict future locations with internally generated theta sequencesPLOS Computational Biology 19:e1011101.https://doi.org/10.1371/journal.pcbi.1011101

The enigmatic mossy cell of the dentate gyrusNature Reviews. Neuroscience 17:562–575.https://doi.org/10.1038/nrn.2016.87

The effect of aging on experiencedependent plasticity of hippocampal place cellsThe Journal of Neuroscience 17:6769–6782.https://doi.org/10.1523/JNEUROSCI.171706769.1997

Phase precession through synaptic facilitationNeural Computation 20:1285–1324.https://doi.org/10.1162/neco.2008.0706292

Neural networks with dynamic synapsesNeural Computation 10:821–835.https://doi.org/10.1162/089976698300017502

Directional tuning of phase precession properties in the HippocampusThe Journal of Neuroscience 42:2282–2297.https://doi.org/10.1523/JNEUROSCI.156921.2021

SoftwareDirectionalnet, version swh:1:rev:7d87577ddb437b0c5941a29fae1b2f1bed909b33Software Heritage.

Shortterm synaptic plasticityAnnual Review of Physiology 64:355–405.https://doi.org/10.1146/annurev.physiol.64.092501.114547
Peer review
Reviewer #1 (Public Review):
In the manuscript entitled "A theory of hippocampal theta correlations", the authors propose a new mechanism for phase precession and thetatime scale generation, as well as their interpretation in terms of navigation and neural coding. The authors propose the existence of extrinsic and intrinsic sequences during exploration, which may have complementary functions. These two types of sequences depend on external input and network interactions, but differ on the extent to which they depend on movement direction. Moreover, the authors propose a novel interpretation for intrinsic sequences, namely to signal a landmark cue that is independent of direction of traversal. Finally, a readout neuron can be trained to distinguish extrinsic from intrinsic sequences.
The study puts forward novel computational ideas related to neural coding, partly based on previous work from the authors, including published (Leibold, 2020, Yiu et al., 2022) and unpublished (Ahmedi et al., 2022. bioRxiv) work. The manuscript will contribute to the understanding of the mechanisms behind phase precession, as well as to how we interpret hippocampal temporal coding for navigation and memory.
https://doi.org/10.7554/eLife.86837.4.sa1Reviewer #2 (Public Review):
Place cells fire sequentially during hippocampal theta oscillations, forming a spatial representation of behavioral experiences in a temporallycompressed manner. The firing sequences during theta cycles are widely considered as essential assemblies for learning, memory, and planning. Many theoretical studies have investigated the mechanism of hippocampal theta firing sequences; however, they are either entirely extrinsic or intrinsic. In other words, they attribute the theta sequences to external sensorimotor drives or focus exclusively on the inherent firing patterns facilitated by the recurrent network architectures. Both types of theories are inadequate for explaining the complexity of the phenomena, particularly considering the observations in a previous paper by the authors: theta sequences independent of animal movement trajectories may occur simultaneously with sensorimotor inputs (Yiu et al., 2022).
In this manuscript, the authors concentrate on the CA3 area of the hippocampus and develop a model that accounts for both mechanisms. Specifically, the model generates extrinsic sequences through the shortterm facilitation of CA3 cell activities, and intrinsic sequences via recurrent projections from the dentate gyrus. The model demonstrates how the phase precession of place cells in theta sequences is modulated by running direction and the recurrent DGCA3 network architecture. To evaluate the extent to which firing sequences are induced by sensorimotor inputs and recurrent network architecture, the authors use the Pearson correlation coefficient to measure the "intrinsicity" and "extrinsicity" of spike pairs in their simulations.
I find this research topic to be both important and interesting, and I appreciate the clarity of the paper. The idea of combining intrinsic and extrinsic mechanisms for theta sequences is novel, and the model effectively incorporates two crucial phenomena: phase precession and directionality of theta sequences. I particularly commend the authors' efforts to integrate previous theories into their model and conduct a systematic comparison. This is exactly what our community needs: not only the development of new models, but also understanding the critical relationships between different models.
https://doi.org/10.7554/eLife.86837.4.sa2Author response
The following is the authors’ response to the previous review
In response to the additional concerns voiced by Reviewer# 2, we have conducted control simulations. The outcomes are summarized in the new supplements to Figure 3. They show that the model is robust under changes of shortterm plasticity parameters and running speed.
Below, we give a pointbypoint response to the remaining comments of the editors and reviewers.
Editorial Assessment: This important work presents an interesting perspective for the generation and interpretation of phase precession in the hippocampal formation. Through numerical simulations and comparison to experiments, the study provides a convincing theoretical framework explaining the segregation of sequences reflecting navigation and sequences reflecting internal dynamics in the DGCA3 loop. This study will be of interest for researchers in the spatial navigation and computational neuroscience fields.
We would like to thank the Editors very much for this positive assessment of our work!
Reviewer #1
In the manuscript entitled ”A theory of hippocampal theta correlations”, the authors propose a new mechanism for phase precession and thetatime scale generation, as well as their interpretation in terms of navigation and neural coding. The authors propose the existence of extrinsic and intrinsic sequences during exploration, which may have complementary functions. These two types of sequences depend on external input and network interactions, but differ on the extent to which they depend on movement direction. Moreover, the authors propose a novel interpretation for intrinsic sequences, namely to signal a landmark cue that is independent of direction of traversal. Finally, a readout neuron can be trained to distinguish extrinsic from intrinsic sequences.
The study puts forward novel computational ideas related to neural coding, partly based on previous work from the authors, including published (Leibold, 2020, Yiu et al., 2022) and unpublished (Ahmedi et al., 2022. bioRxiv) work. The manuscript will contribute to the understanding of the mechanisms behind phase precession, as well as to how we interpret hippocampal temporal coding for navigation and memory.
I am very pleased to have seen major improvements in the manuscript regarding (i) a clarification of the concepts of extrinsic and intrinsic mechanisms, and (ii) overall arrangement of Figures but also (iii) expanding on some important concepts such as the role of experience in determining the asymmetric connectivity that is necessary for intrinsic models of sequence generation.
We are delighted to have been able to amend the Reviewer’s concerns voiced after the initial submission. We are very grateful for their many good suggestions that allowed us to make important additions to the revised manuscript.
Reviewer #2
Place cells fire sequentially during hippocampal theta oscillations, forming a spatial representation of behavioral experiences in a temporallycompressed manner. The firing sequences during theta cycles are widely considered as essential assemblies for learning, memory, and planning. Many theoretical studies have investigated the mechanism of hippocampal theta firing sequences; however, they are either entirely extrinsic or intrinsic. In other words, they attribute the theta sequences to external sensorimotor drives or focus exclusively on the inherent firing patterns facilitated by the recurrent network architectures. Both types of theories are inadequate for explaining the complexity of the phenomena, particularly considering the observations in a previous paper by the authors: theta sequences independent of animal movement trajectories may occur simultaneously with sensorimotor inputs (Yiu et al., 2022).
In this manuscript, the authors concentrate on the CA3 area of the hippocampus and develop a model that accounts for both mechanisms. Specifically, the model generates extrinsic sequences through the shortterm facilitation of CA3 cell activities, and intrinsic sequences via recurrent projections from the dentate gyrus. The model demonstrates how the phase precession of place cells in theta sequences is modulated by running direction and the recurrent DGCA3 network architecture. To evaluate the extent to which firing sequences are induced by sensorimotor inputs and recurrent network architecture, the authors use the Pearson correlation coefficient to measure the ”intrinsicity” and ”extrinsicity” of spike pairs in their simulations.
I find this research topic to be both important and interesting, and I appreciate the clarity of the paper. The idea of combining intrinsic and extrinsic mechanisms for theta sequences is novel, and the model effectively incorporates two crucial phenomena: phase precession and directionality of theta sequences. I particularly commend the authors’ efforts to integrate previous theories into their model and conduct a systematic comparison. This is exactly what our community needs: not only the development of new models, but also understanding the critical relationships between different models.
We also would like to express our gratitude to Reviewer 2 for their numerous constructive criticisms that led to a very much improved revised manuscript!
Reviewer #2
1. The choice of timescale parameters for input facilitation and synaptic depression is still not fully justified in my opinion. The authors themselves mention that previous experiments suggest wide ranges for both timescales. Given that the generation of intrinsic and extrinsic sequences in their model is primarily driven by these two mechanisms, their chosen timescales should significantly impact the simulation results. I urge the authors to discuss the potential effects of selecting different sets of timescales and the possible limitations of the current selection of 500ms for both.
For instance, the authors state in the caption of Fig 1 that all simulated rat trajectories were set at a speed of 20cm/s, which is a rat’s walking speed. However, the running speed of rats can exceed 3m/s. In this case, none of the CA3 cells in the model would produce any extrinsic sequences since the animal would traverse the place fields much more rapidly, preventing the sensorimotor input from increasing as it does in the model.
The reviewer raised the valid point that our simulations may be sensitive to the shortterm plasticity time constants and running speeds. We therefore conducted new simulations illustrated in Figure 3—figure supplements 1 and 2.
In agreement with the reviewer’s assertion, using the current model parameters, a higher running speed would not elicit extrinsic sequences due to the lack of depolarization from spatial input (Figure 3—figure supplement 2A). However, an increase of running speed also requires sensory inputs to be available on a larger spatial scale (width of the spatial input box in our case). ParraBarrero et al., eLife 10:e70296 and ParraBarrero & Cheng 2023, PLOS Comput. Bio. 19:e1011101, e.g., showed that place field sizes become larger under higher running speeds and consequently lengthen the theta sequences. With such modification, along with a longer DG projection length (r), we were able to recover the theta sequence at a higher speed (100 cm/s), using the same STD and STF time constants (Figure 3—figure supplement 2B). Furthermore, it has been shown that theta frequency increases with running speed (e.g., Rivas et al., 1996, Exp Brain Res 108:1138). In our analysis, a higher theta frequency (12Hz instead of 10Hz) is also able to counteract the effect of running speed and leads to controllevel like phase precession (Figure 3—figure supplement 2C).
Consistent with this finding, the original study of Romani & Tsodyks 2015, Hippocampus 25:94105, found a fourfold increase of speed (from 0.05 to 0.2 fraction of the track per second) to not affect phaseposition relations (with UD = 0.8 and 800ms STD time constant), likely due to the large place field sizes covering 1/3 of the track. Thus, phase precession may only be affected by high speeds in narrow place fields in which activity would only be present for few theta cycles thus naturally having limited capacity for phase coding.
We further refrain from increasing the running speed beyond 1m/s (e.g. 3m/s as suggested by the reviewer), as the typical running speed of a rat in an 80cm square environment is between 2040cm/s (Mankin et al. 2012, PNAS, 109:1946219467). Even on linear tracks, reported running speeds hardly exceed 120 cm/s (e.g. Ahmed and Mehta, 2018, J Neurosci 32:73737383; Schmidt et al., 2009 J Neurosci 29:1323213241). To our knowledge phase precession for speeds above 1.2 m/s has not been reported so far at all, certainly also owing to experimental challenges. We, however, would speculate that beyond 120 cm/s phase precession could be meaningful in large environmental enclosures with wide place cells. Thus a version of our input model with very large place field sizes should generally be able to also cover very high running speeds.
To conclude, STD and STF time constants do not need to be in a precise range to accommodate the behavioural time scales if the sensory input changes on accordingly larger spatial scales.
Following up on the reviewer’s additional concern, we also checked the effect of time constants on the theta sequences (while keeping the running speed unchanged). Decreasing the time constant of STF (τF) to 100ms would degrade the theta sequence due to a lack of depolarization, as sensory input reverts to its resting value (=0) too fast, but at 250ms, the temporal correlation of theta sequences is largely maintained (Figure 3—figure supplement 1A). However, such effects can be compensated for by an increase in sensory input which promotes input facilitation (Figure 3—figure supplement 1B). Further increasing τF does not significantly affect theta sequences as the sensory input amplitude have asymptotically reached their target values (Figure 3—figure supplement 1A bottom). The temporal correlation of theta sequences is not sensitive to the change in the time constant of STD (τD) (Figure 3—figure supplement 1C), possibly because the synaptic resource of the place cells behind the animal is reliably depleted by strong depolarization despite a fast recovery time (τD=100ms).
Since the relation between running speed and theta sequences has been thoroughly studied in ParraBarrero et al. 2021 and ParraBarrero & Cheng 2023, and the precise range of STD and STF time scales does not play a critical role in the temporal structure of theta sequence, we refrain from substantially revising the manuscript and only briefly add these points after Figure 3.
1. This is a point I overlooked in the initial review. The synaptic depression fraction UD is set at 0.9 or 0.7, implying that the synaptic coupling weight between CA3 excitatory cells (and CA3 to DG) is almost entirely depleted within a few hundred milliseconds. To my limited neuroscience knowledge, I am not aware of any experimental results that corroborate this potentially bold setup, and I urge the authors to provide relevant experimental and theoretical references if they exist.
Most crucially, I find this setup biased towards supporting the authors’ theory for intrinsic sequences because it essentially eliminates the possibility of any CA3 cell producing an effective output to other neurons after it fires. Hence, I question whether the simulation results would be much less clean if a more moderate depression factor UD were utilized.
We thank the reviewer for giving us the opportunity to further clarify. (1) Probabilties of synaptic release (here called UD for consistency with the original work by Romani and Tsodyks), can attain a very wide range and indeed achieve values up to 0.9 (for review see e.g. Dobrunz LE, Stevens CF, 1997, Neuron, 18: 9951008). (2) Contrary to the reviewer’s impression, a higher UD (0.70.9 in our case) would bias the simulation towards even more extrinsicity. Larger UD produces steeper phase precession in extrinsic sequences, because it (temporarily) generates an even stronger asymmetrical connectivity. (3) The extreme value of 0.9 was only used in Figure 1 to best illustrate the original Romani and Tsodyks 2015 idea. (4) Our simulations without recurrent synaptic connections (Figure 6) do not even require shortterm synaptic depression. In view of these arguments we refrained from making further additions to the paper and refer the critical reader to this comment.
I have a few final suggestions for the authors in the hopes of further improving the manuscript for the neuroscience community:
line 62: sensorimotor input is present or ABSENT?
Intrinsic activity signatures are found ”EVEN when sensorimotor feedback is present”, as one may assume that this input may be able to completely override the intrinsic patterns.
line 76: played out. colloquial, consider rewriting/explaining
We use ”evoked” now.
line 104: second part of motivation for Izhikevichtype model is wordy, and grammatically incorrect.
We have shortened the sentence.
on potential limitations of the model lines 116120: is the use of a box an important assumption, as opposed to a more graded function, exponential or gaussian?
Usingspike based input, it is not straightforward how to implement a graded input. One way would be to employ a stochastic point process with graded firing probability. We, however, chose to use a nonlinear facilitation function (see below).
line 124 (equation) and 129130: How crucial is the nonlinearity in the synaptic variable for the results? This is a strong assumption, as the nonlinearity is the dominant effect (as opposed to a correction/perturbation). Are there any other contributions for this ramp of activity due to sensory input?
We found results to fit best with a nonlinear facilitation function (see above), and, as argued in the manuscript, facilitation indeed acts nonlinearly owing to the calciumdependence of synaptic release. We have added a comment to the Methods section explaining that we use facilitation to generate a graded spatial input.
line 187: ’...neglecting gamma activity in the model.’ I suggest removing this part of the sentence, unless you motivate why gamma would be relevant and the conditions for its generation.
We have followed the reviewer’s suggestion.
https://doi.org/10.7554/eLife.86837.4.sa3Article and author information
Author details
Funding
Deutsche Forschungsgemeinschaft (LE2250/131)
 Christian Leibold
Deutsche Forschungsgemeinschaft (INST 39/9631 FUGG)
 Christian Leibold
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Senior Editor
 Panayiota Poirazi, Foundation for Research and Technology Hellas, Greece
Reviewing Editor
 Brice Bathellier, CNRS, France
Version history
 Preprint posted: February 5, 2023 (view preprint)
 Sent for peer review: March 2, 2023
 Preprint posted: May 17, 2023 (view preprint)
 Preprint posted: August 14, 2023 (view preprint)
 Preprint posted: September 18, 2023 (view preprint)
 Version of Record published: October 4, 2023 (version 1)
Cite all versions
You can cite all versions using the DOI https://doi.org/10.7554/eLife.86837. This DOI represents all versions, and will always resolve to the latest one.
Copyright
© 2023, Yiu and Leibold
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 260
 Page views

 47
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
The joint storage and reciprocal retrieval of learnt associated signals are presumably encoded by associative memory cells. In the accumulation and enrichment of memory contents in lifespan, a signal often becomes a core signal associatively shared for other signals. One specific group of associative memory neurons that encode this core signal likely interconnects multiple groups of associative memory neurons that encode these other signals for their joint storage and reciprocal retrieval. We have examined this hypothesis in a mouse model of associative learning by pairing the whisker tactile signal sequentially with the olfactory signal, the gustatory signal, and the tailheating signal. Mice experienced this associative learning show the whisker fluctuation induced by olfactory, gustatory, and tailheating signals, or the other way around, that is, memories to multimodal associated signals featured by their reciprocal retrievals. Barrel cortical neurons in these mice become able to encode olfactory, gustatory, and tailheating signals alongside the whisker signal. Barrel cortical neurons interconnect piriform, S1Tr, and gustatory cortical neurons. With the barrel cortex as the hub, the indirect activation occurs among piriform, gustatory, and S1Tr cortices for the secondorder associative memory. These associative memory neurons recruited to encode multimodal signals in the barrel cortex for associative memory are downregulated by neuroligin3 knockdown. Thus, associative memory neurons can be recruited as the core cellular substrate to memorize multiple associated signals for the firstorder and the secondorder of associative memories by neuroligin3mediated synapse formation, which constitutes neuronal substrates of cognitive activities in the field of memoriology.

 Chromosomes and Gene Expression
 Neuroscience
Mathys et al. conducted the first singlenucleus RNAseq (snRNAseq) study of Alzheimer’s disease (AD) (Mathys et al., 2019). With bulk RNAseq, changes in gene expression across cell types can be lost, potentially masking the differentially expressed genes (DEGs) across different cell types. Through the use of singlecell techniques, the authors benefitted from increased resolution with the potential to uncover cell typespecific DEGs in AD for the first time. However, there were limitations in both their data processing and quality control and their differential expression analysis. Here, we correct these issues and use bestpractice approaches to snRNAseq differential expression, resulting in 549 times fewer DEGs at a false discovery rate of 0.05. Thus, this study highlights the impact of quality control and differential analysis methods on the discovery of diseaseassociated genes and aims to refocus the AD research field away from spuriously identified genes.