1 Introduction

As a rat navigates in an environment, place cells fire sequentially during one theta cycle (∼100 ms) and form time-compressed 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. Non-local 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 out-of-field firing at reward locations (Sasaki et al., 2018). These remote activations of place cells reflect the underlying circuit connectivity rather than the actual locomotive state of 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 hard-wired 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 short-term 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 Tsodyks et al. (1996) model use a fixed asymmetrical weight matrix pre-designed to align with one movement trajectory. 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 short-term synaptic plasticity mechanisms, while the intrinsic sequences are played out by the CA3-DG recursive loops with fixed asymmetrical weights, as inspired by the experimental evidence that lesions in DG abolish non-local 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 pre-designed CA3-DG 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 sensory-motor experience is lacking and may reflect the mechanistic basis for offline replay.

2 Methods

Neuronal model

Neurons are 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) ← c and the adaptation variable by u ← 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

Spatial input

The place field centers 80 × 80 = 6400 excitatory CA3 cells equally tile the 80 by 80 cm square arena. Place cell firing rates are modelled direction-sensitive, with best directions semi-randomized among each2 × 2 tile of place cells by randomly rotating a set of 4 equally spaced direction angles by a uniformly distributed angle ξ, i.e.

The sensory input into the i-th neuron depends on the instantaneous position, I^MEC(t) if p(t) = x(t), y(t), and heading direction ψ(t) of the animal as

where A^pos is the amplitude of positional tuning and d(·) 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 5cm from the field center. Directional tuning is implemented as an additional amplitude gain A^dir to the positional current depending on the circular difference between the animal’s heading and the neuron’s best direction . The sensory input is assumed to be modulated bytheta oscillations from medial entorhinal cortex (MEC) I^MEC(t) with a phase shift of 70^° (Mizuseki et al., 2009).

The sensory input JS is subsequently transformed to the input current I^S via short-term facilitation (STF)

where the facilitation variable decays to with a time constant τ ^F = 500ms and increases to when the sensory input is present. Φ controls the strength of the STF. The facilitation variable is squared to include non-linear interactions in presynaptic calcium dynamics. As a result, facilitated sensory input increases over time and becomes stronger in the later part of the field.

Note that only the CA3 place cells receive the sensory input. 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 is conductance-based, and follows the equations:

The conductance of a post-synaptic cell i is increased by the spike arrivals at times from the pre-synaptic cell j, and decay with a time constant τ ^E = 12ms. N_J = 6, 400 is the number of presynaptic place cells, V ^E = 0mV is the reversal potential of the excitatory synapses and τ₀ = 2ms 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 angles, i.e.,

where B^pos and σ correspond to the maximum strength and width of the location-specific interaction, respectively. B^dir and K^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 2-D variant of Tsodyks et al. (1996) model in Figure 1 and otherwise turned off in the rest of our analysis.

Figure 1:

Furthermore, the recurrent synaptic conductances underwent short-term synaptic depression (STD), as was proposed in Romani and Tsodyks (2015) to serve as sequence generator in 2-D 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 which represents the available synaptic resourceD and follows the dynamics:

where recovers to 1 with a time constant τ ^D = 500 ms and is depleted by a fraction U^D every time a spike occurs. The STD only applies to synaptic connections when CA3 place cells are pre-synaptic. is fixed at 1 when the pre-synaptic cells are inhibitory interneurons or DG place cells.

DG layer

We simulated N_DG = 40 × 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 and directional tunings are determining synaptic strengths to and from the CA3 layer. The directional tuning is randomized as in 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 best angles between the CA3 and DG populations. The variable defines the path corresponding to the intrinsic sequence by choosing

where varies with the intrinsic path direction θ_DG as .

The excitatory synaptic strengths from DG to CA3 are chosen such that DG cells projects back to CA3 cells with place field centers shifted by a vector r = [4 cos(θ_DG), 4 sin(θ_DG)] of fixed length of 4 cm along the intrinsic path, i.e.,

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 (CA3-Inh-CA3 and DG-Inh-DG). The dynamics of their synaptic currents mirrors the excitatory synapses, i.e.,

with τ ^I = 10ms, V ^I = −80mV. CA3 and DG have all-to-all connections to their inhibitory populations with uniformly randomized strengths, i.e. , with ξ ∼ U (0, 1). is the maximum synaptic strength, and the notation X-Y corresponds to Inh-CA3 and Inh-DG connections. There is no synaptic connection between interneurons, i.e. W^Inh−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 all-to-all excitatory currents from their respective layer. Those currents are modelled according to eqs.(1-2). The synaptic weights are constant and denoted by and .

Parameters of the models

Model parameters that are adjusted in different analyses are listed in table 1.

Table 1:

Cross-correlation analysis

Cross-correlation represents the probability that a spike of one place cell would occur following a certain time lag from the spike of the another cell. Cross-correlation is always 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^°, 45^° and 90^° trajectory in each comparison group, and in Figure 2, where the time lag direction is from left to right cells.

Figure 2:

Correlation lag is derived by band-pass (4-12 Hz) filtering the cross-correlation 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 correlations in a pair of place fields following Yiu et al. (2022). We compare the cross-correlation signals of a field pair for a running direction along the DG-loop (θ_DG) and opposite to the loop (θ_DG + 180^°). Extrinsicity (Ex) is computed as the Pearson’s correlation (r) between two cross-correlation signals, and intrinsicity (In) between the signal of θ_DG and the horizontally flipped signal of θ_DG + 180^°. The Pearson’s correlation is then transformed (r^′ = (r + 1)/2) to be in the range of 0 and 1. An extrinsic pair would give an extrinsicity near 1, since the effect of DG loop is minimal and correlation signals are similar in both θ_DG and θ_DG + 180^° directions. An intrinsic pair would see correlation signal horizontally flipped in the θ_DG + 180^° 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. is a post-synaptic potential (PSP) kernel with decay and rising time constants of τ = 5ms and τ_r = 1.25ms respectively. V₀ is a factor which normalizes the PSP kernel to 1. A spike is said to occur if V (t) crosses the firing threshold V_Θ = 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 pre-synaptic place cells in a theta cycle (100ms).

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=+20cm at y=+20cm in direction θ_DG = 0^°, while the loop is absent in the lower half of the arena.

During training, we applied “non-moving” spatial inputs to the CA3 place cells at the with-loop (0 cm, 20 cm) and no-loop (0 cm, -20 cm) locations for 1 second, 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 20cm squared boxes centered at the two locations. Each population contains 400 pre-synaptic cells, forming the input space for the tempotron. The spikes from the with-loop population will train the first tempotron and those from the no-loop population will train the second tempotron. Prior to training, the input spikes are sub-divided to 10 patterns based on their theta cycles. Each pattern has a window of 100ms. We added noise to the patterns by jittering the spikes with Gaussian noise 𝒩∼ (0, (2ms)²) for 100 times. As a result, each tempotron receives 10 × 100 = 1000 training patterns from the activity evoked by the location. All training patterns are (+) patterns and there is no (−) pattern.

After training, trajectories (20cm long, 1s duration) with running directions from 0^° to 360^° with 15^° increment were simulated to cross each of the locations. The trajectories produce intrinsic sequences in the with-loop population and extrinsic sequences in the no-loop population. The patterns evoked by the running trajectories were separately applied to the tempotrons. The input spikes for testing were also subdivided in to 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.

3 Results

Theta sequence directionality in intrinsic and extrinsic models

Theta-scale 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. For extrinsic models, the spike correlation is determined by sensory-motor inputs. We first illustrate how these mechanisms work for two exemplary representatives of these two major model classes.

For intrinsic models, we refer to the original 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 one preferred running direction. Here we simulate a network of CA3 place cells with fixed asymmetrical connectivity as suggested in Tsodyks et al. (1996) model (see Methods for the implementation) and applied our model to behavioral running trajectories in a 2-D open space (Figure 1A). Phase precession and spike correlations (Figure 1B-D) are compared for opposite running directions. In our simulation, all place cells project excitatory synapses to their counterparts with rightward neighboring place fields (see CA3 recurrent connections in Methods). Phase precession therefore is determined by how closely the running direction matches the preferred direction imposed by the intrinsic connections. The closer the trajectory angle aligns with this preferred direction, the more negative is the slope of phase precession (Figure 1C). 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 movement direction (Figure 1D, see Cross-correlation analysis in Methods). 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 ratebased model by Romani and Tsodyks (2015), where phase precession was explained by symmetric recurrent connections that undergo running direction-dependent attenuation by short-term 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. We simulated our spiking variant of the Romani and Tsodyks (2015) model with the same trajectories as the intrinsic model (Figure 1A, E), and recovered direction-independent phase precession (Figure 1F). Since now, the theta sequences are played out in the same direction as the movement, theta spike correlations are symmetrically reversed (Figure 1G) 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 best directions (See Spatial input in Methods). The sensory input is assumed to arise from MEC, and hence, it is also theta-modulated and phase-shifted by 70^° with respected to the peak of theta cycle (Mizuseki et al., 2009). Furthermore, since the precession slope observed in Romani’s and Tsodyks’ (2015) model is limited (-1.13 radians per field size, see Figure 1F) as compared to the experimental reports (−4.44 radians (Yiu et al., 2022) and −266 degrees (Schmidt et al., 2009) per field size), we introduced short-term 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 lower phase tail and increasing the steepness of slopes. Finally, we designated the synaptic weights to be stronger between place cells with similar preferred 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 best direction of the field (red) and one for which the trajectory runs along the worst direction, which is opposite to the best direction of the field (Figure 2F). Consistent with experimental data (Yiu et al., 2022), phase precession has a lower onset and marginal spike phase along best direction than along the worst (Figure 2G-H), reflecting that larger depolarizations generally yield shorter latencies. Directionality of the input, however, 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 DG-CA3 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. 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, again through the DG loop. This scheme produces a rigid sequence whose activation order is independent of the movement direction. The connection pattern of DG-CA3 projections could be determined by pre-existing network structure or past experience through associative learning, or both.

Figure 3:

Figure 3, provides schematic illustrations, for a DG layer that either only projects CA3 activity to their rightward neighbours (θ_DG = 0^°, Figure 3A) or only to their left-ward neighbors (θ_DG = 180^°, Figure 3E). Simulations for both cases (θ_DG = 0^° and θ_DG = 180^°) 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. DG-loop connectivity is additionally modulated by firing rate directionality of the CA3 place fields. Fields with similar best direction are more strongly connected via the loop than those with opposite best directions (see DG layer in Methods).

We found that when the simulated animal is running in the same direction as the DG-CA3 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 along best-direction remains lower than along the worst direction (Figure 3D). When, however, the animal is running against the DG-CA3 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 (Figure 3F). The latter appears as higher phase at the tail of 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 DG-loop projection (θ_DG = 180^°), 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, we use the measures extrinsicity and intrinsicity (Yiu et al., 2022) (see Methods for details). Consistent with experimental data (Yiu et al., 2022), 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 (Figure 3J), since along the best direction, cells receive more sensory depolarization, and thus, the movement-dependent extrinsic sequences are more activated. The model also explains, why pairs with similar best place field directions may be less extrinsic than pairs with approximately opposite best direction (dissimilar pairs), since the DG loop preferentially connects CA3 place cells with similar best directions.

Thus, by introducing feedback excitation via the DG layer, intrinsic sequences are able to propagate in fixed directions on top of the movement-dependent 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.

Lesion in DG reduces pair correlation

One prediction of the DG-loop 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 short-term synaptic facilitation (Figure 4A).

Figure 4:

For a simulated rightward trajectory, we plotted the pair correlation lags versus the distance between the centers of two fields (Figure 4B-C). Theta compression (Dragoi and Buzsáki, 2006) measures how much theta phase encodes a certain interval in space and is quantified by magnitude of the linear-circular regression slopes (a, in radians per maximum pair distance) of the lag-distance plots. Theta compression is lower for the DG lesion case (a = 0.54 for both θ_DG = 0^° and θ_DG = 180^°), as compared to the control case (a = 2.23 for θ_DG = 0^° and a = −0.77 for θ_DG = 180^°) reproducing the finding (Ahmadi et al., 2022) that spatial encoding via theta sequences crucially depends on intact DG and predicting that loss of DG inputs is compensated for by the increase of release probability in the spared afferent synapses from the MEC.

Theta sequences in 2-D

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 2-D space by examining oblique trajectories which cross the DG-loop projection at certain angles.

We first arrange the DG-loop connections such that they cross a rightward trajectory at 45^° and 225^° (Figure 5A-F). Similar to the cases of θ_DG = 0^° and θ_DG = 180^° (Figure 3D and 3H), precession slopes are steeper and onsets higher when the trajectory aligns more with the orientation of the DG-loop, but with a smaller effect size for oblique crossings (Figure 5A and 5C) since DG-loop connectivity only covers part of the trajectory near the intersection. We further resolve the precession slope, onset and marginal phase for each place cell into 2-D maps (Figure 5B and 5D). Intrinsic sequences with a higher marginal spike phase can be clearly seen along the belt of DG-loop projections and are even extended to the outside of trajectory predicting “off-track” spikes at high phases. Depending on the movement alignment with the DG-loop orientation, the slope becomes either more negative (θ_DG = 45^°) or more positive (θ_DG = 225^°). 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 extrinsic pair center differences are mostly oriented horizontally (along the running direction) while intrinsic pair center differences are oriented along the DG-loop orientation θ_DG = 45^°, as by design (Fig 5F).

Figure 5:

The analysis above is repeated for the geometric configurations that DG-loop connectivity is minimally interacting with the movement direction, i.e., when they are perpendicular to each other (θ_DG = 90^° and θ_DG = 270^°, Figure 5G). Similar effects as in Figure 5B and 5D on precession slope, onset and marginal phases are also observed in the 2-D map, except 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 extrinsic-intrinsic ratios) due to the smaller overlapping area between DG-loop projection and the trajectory. Lastly, the pair center difference orientation confirms that extrinsic pairs follow the trajectory direction while intrinsic pairs are biased towards the DG-loop orientations (90^°).

The results demonstrate the distinct roles of extrinsic and intrinsic sequences in 2-D 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 the change in phase precession properties, most notably the higher spike phases from the DG-CA3 recurrent input, as well as increased intrinsicity of pair correlation and extended firing fields along the orientation of the DG-loop projections.

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 lesser expressed intrinsic sequence contributions is not straight forward. Simulations of trajectories in 2-D (Figure 5) suggest intrinsic activity may serve a role to identify certain location-direction 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. 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 Tempotron section in Methods). Two tempotrons were trained to recognize the spike patterns from the place cells, one taking input from a model with DG-loop connectivity at θ_DG = 0^°, and one without DG-loop connectivity separately (Figure 6A). Non-moving spatial inputs were applied to the CA3 place cells at the centers of with-loop and no-loop populations and their spike patterns in each theta cycle are used as training patterns, mimicking a situation in which network activity is evoked without sensory-motor input as, e.g., in a situation before the animal has seen the environment. The training patterns have only (+) labels, which the tempotrons are trained to recognize by firing a spike (Figure 6B). We then test the tempotrons with spike patterns induced by trajectories from 0^° to 360^° going through the two current injection sites (Figure 6C). All spikes in each training and testing pattern are individually jittered by adding a noise term σ ∼ N (0, 2ms²) 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. We found that for the population involving a DG loop, the tempotron is able to recognize the sequence patterns produced for all running directions, while for the no-loop population, the tempotron fails to identify the sequences in most of the trajectories (see accuracies in Figure 6D). Note that the tempotron performed its classification while the extrinsic sequences were not disabled. 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 6E), while spike patterns of the no-loop network can no longer be recognized (Figure 6F).

Figure 6:

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 as a mixture with the extrinsic sequences.

4 Discussion

We presented a model of hippocampal theta sequences in 2-d environments, suggesting that both extrinsic and intrinsic mechanisms are required to explain experimental reports that phase precession and spike timing correlations are non-homogeneous 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 recursive 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 spike-based temporal pattern decoder (tempotron), we showed that the trajectory-independent 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 1C, 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 (θ_DG = 180^°) but at a lower probability as indicated by the larger fraction of positive phase-position 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 (worst-angle) pairs and extrinsicity (intrinsicity). These predictions of our model are corroborated by past reports of higher spike phases in the non-preferred arm of a T-maze (Kay et al., 2020) as well as the association of rigid correlations with upward shifts in spike phases and an increase in worst-angle pairs (Yiu et al., 2022).

Since intrinsic sequences can also propagate outside the trajectory (Figure 5) and activate place cells non-locally, our model predicts direction-dependent expansion of place fields. 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 short-range 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 2-d 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).

The back-projection from CA3 to DG is a crucial anatomical prerequisite of our model, but was rarely explored compared to the feed-forward inputs via the perforant pathway. The proposed CA3-DG 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, 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 out-of-ensemble activity. Indeed, past studies showed that reliable excitatory effect could be observed when granule cells were depolarized (Scharfman, 1995) and when they received back-propagation 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. Moreover, lesions of the DG layer 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 of 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 2-d sensory features to cell ensembles in the loop (Leibold, 2020). The resulting topology can exhibit discontinuous long-range 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 head-direction 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 spike-based model, where stronger depolarization causes earlier spiking. This phase-rate 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 highand low-spiking runs, suggesting that the phase is not single-handedly determined by firing rate. In our model, firing rate is determined by both low-phase spiking from sensory input and high-phase spike arrivals of DG-CA3 loops, both producing opposing effects on the phase distribution. Thus, depending on the strength and geometry of the DG-CA3 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) with multiple fixed sequences that form a decodable spatial representation in a reinforcement learning paradigm.

We speculate that the functional roles of intrinsic sequences may not be limited to spatial memories. While, in the spatial domain, intrinsic sequences could be interpreted as planning of alternative trajectories during navigation or prospective planning of future pathways (Kay et al., 2020; Sasaki et al., 2018), virtual non-spatial trajectories could represent working memories contents (Jensen et al., 1996) available for general decision making processes.

Acknowledgements

This work was funded by the German Research Association (DFG) nder grant LE2250/13-1.