|Home | About | Journals | Submit | Contact Us | Français|
Hippocampal neurons can display reliable and long-lasting sequences of transient firing patterns, even in the absence of changing external stimuli. We suggest that time-keeping is an important function of these sequences, and propose a network mechanism for their generation. We show that sequences of neuronal assemblies recorded from rat hippocampal CA1 pyramidal cells can reliably predict elapsed time (15-20 sec) during wheel running with a precision of 0.5sec. In addition, we demonstrate the generation of multiple reliable, long-lasting sequences in a recurrent network model. These sequences are generated in the presence of noisy, unstructured inputs to the network, mimicking stationary sensory input. Identical initial conditions generate similar sequences, whereas different initial conditions give rise to distinct sequences. The key ingredients responsible for sequence generation in the model are threshold-adaptation and a Mexican-hat-like pattern of connectivity among pyramidal cells. This pattern may arise from recurrent systems such as the hippocampal CA3 region or the entorhinal cortex. We hypothesize that mechanisms that evolved for spatial navigation also support tracking of elapsed time in behaviorally relevant contexts.
Tracking time is of fundamental importance in a wide range of brain operations, including sensory perception and motor actions, learning, memory, planning, decision-making and language (Gibbon et al., 1997; Buonomano and Karmarkar, 2002; Ivry and Spencer, 2004; Mauk and Buonomano, 2004; Buhusi and Meck, 2005). Despite the central importance of temporal processing, its underlying neural mechanisms remain unknown. At the systems level, two competing ideas have been put forward: timing is generated by a central mechanism and distributed to various brain regions (Church, 1984) or each subsystem produces its own timing (Mauk & Buonomano 2004). With regard to timing and duration, a distinction is made between subsecond (perceptual-motor) and suprasecond (cognitively mediated) scales (Michon, 1985; Lewis et al., 2003). At the level of mechanisms, two models are typically distinguished; clocks and ramping time keepers, with neuronal substrates in the cerebellum, basal ganglia, prefrontal, motor and parietal cortical regions, cf. (Mauk and Buonomano, 2004; Buhusi et al., 2005). The hippocampus has also been implicated in timing (Clark and Isaacson, 1965; Thompson and Krupa, 1994; Young and McNaughton, 2000), although the mechanism has remained elusive.
We report here on a novel form of time-tracking mechanism, which is manifested by evolving transiently active cell assemblies and accurate for periods of tens of seconds. First, we examine the ability of evolving neuronal sequences to predict elapsed time in a memory task. Second, we propose a simple network model with Mexican-hat-type connectivity and adaptation of the membrane potential thresholds for action potential generation that is similar to what has been observed in rodent hippocampus (Henze and Buzsaki, 2001) and in fish (Chacron et al., 2007). The threshold-adaptation model reproduces the key properties of the observed sequences, suggesting that time-keeping in the hippocampus may arise from the same cellular and network mechanisms that support spatial navigation.
The experimental data used in this paper were adopted from (Pastalkova et al., 2008), where all relevant experimental methods and protocols are described. The animals were all male rats.
Two versions of time prediction models were fit from experimental data: a ‘rate-only’ model and a ‘phase-only’ model. These models used firing rates or theta phases of spikes to fit a probability distribution for population spiking activity in 0.5s time windows. Maximum-likelihood estimation was then used on single trials to predict the most likely time reflected by the population activity at each time bin. In all cases, the models were fit from complementary trials, never including the trial on which time was subsequently inferred. The time estimates corresponding to each time bin were obtained independently, so that only knowledge of the current population activity was needed to estimate elapsed time. See the Supplementary Text for a more detailed description of the time prediction models.
We modeled network dynamics using a standard firing rate model, with threshold nonlinearity [y]+ = max(y,0). At any point in time, the vector represents a population vector of firing rates for each of N neurons. A key ingredient of the model is the activity-dependent adaptation of the ‘spike’ thresholds of individual neurons, represented by the dynamic variables . The model can thus be described by a system of 2N equations:
where τm = 30ms and τa = 2s are membrane and threshold-adaptation time constants, respectively, J is the matrix of synaptic weights for the recurrent network, and I is a (time-dependent) vector of external inputs. The constant c = 0.5 controls the strength of the activity-dependent adaptation, whereas τa determines the timescale with which the thresholds recover in the absence of firing.
The connectivity matrix J is constant in time and can be written as a sum of two components, J = J0 + Jhet, where J0 has Mexican-hat-type connectivity on a two-dimensional lattice of neurons with periodic (torus) boundary conditions, and Jhet is a matrix of heterogeneous weights, sampled randomly and independently from a normal distribution with mean zero. In simulations, we generated three different instances of Jhet, resulting in three different connectivity matrices J1 J2 and J3.
In simulations, two types of behavioral conditions were distinguished: ‘task’ and ‘homecage’. In all cases initial conditions consisted of the same “bump” of firing rate activity in the center ofr the two-dimensional lattice of neurons, and differed in the adaptation variables only. Neurons with adapted thresholds were chosen to lie at the left, bottom, right or top of the initial bump of activity, resulting in initial conditions ‘A’, ‘B’, ‘C’ and ‘D’, respectively. For ‘task’ trials, the initial conditions were consistent from trial to trial, and the model was driven by temporally and spatially unstructured noise I(t); different instances of noise was thus the only difference between trials of the same initial condition type. In the ‘homecage’ trials, the initial conditions A, B, C and D were randomized across trials, and the model was driven by spatially unstructured noise that had temporal correlations on the order of 125ms. (See the Supplementary Text for further information.)
In order to investigate the ability of a downstream layer to ‘inherit’ the sequence generated by the threshold adaptation model (‘layer 1’), we simulated activity in a second layer connected to the first via sparse and random feedforward projections. The dynamics in this layer were governed by the same equations (1) and (2) in the previous layer, with two differences. First, the connectivity matrix J in layer 2 represented an overall global inhibition (to ensure sparse firing) and had no spatial structure. Second, the input vector I(t) had two components: temporally and spatially unstructured noise; and feedforward input derived from the activity in the previous layer via random and sparse connections (10% connection probability between layers). (See the Supplementary Text for further information.)
A reliability measure was used to compare population vectors from individual trials to the mean. First, we normalized firing rates so that each neuron had the same average firing rate, calculated over all time bins and all trials.
Let be the population vector corresponding to the i-th trial at time t. For each fixed time bin t0, we computed a “reliability score” Ri(t0) for each trial by computing the squared-distance between the population vector and the mean vector across trials , . This distance was then compared to the average squared-distance between and the mean population vectors at all other times, , where Nt is the total number of time bins. If a trial is reliable at time t0, then Wi(t0) < Bi(t0). The reliability was thus defined to be , with the denominator chosen so that –1 ≤ Ri(t0) ≤ 1. A reliability of 1 corresponds to the “best case” scenario of , and -1 corresponds to the “worst case” scenario of for all t ≠ t0, and . If is closer to the mean vector than it is to the mean vectors from other time bins , then Ri(t0) will be positive. For random (and thus completely unreliable) population vectors, Ri(t0) ≈ 0. The average reliability R(t) was computed as the population trial-average of the time series vectors Ri(t).
A reliable pattern of sequential activation of neuronal activity was observed in the CA1 region of hippocampus during the delay period of a memory task (Pastalkova et al., 2008) (c.f. (Gill et al., 2010; Macdonald et al., 2010)). Three rats were trained to run for approximately 20 seconds in a running wheel during the delay period prior to making a choice (left or right) for the next run through a T maze. Action potentials for pyramidal cells were recorded together with local field potentials. The pattern of sequential activation of simultaneously recorded neurons for a given trial type (Fig. 1A,B) was reliable across trials and lasted 10-20 sec without repeating itself. Therefore, we hypothesized that the population spiking activity of pyramidal neurons in CA1 at any point in time during a trial could be used to infer elapsed time.
To verify this hypothesis, we designed two probabilistic models for inferring elapsed time from instantaneous neural activity: one based on the cells’ firing rates, and the other using phases of spikes with respect to the theta oscillation (see Methods). Both models were good predictors of elapsed time on single trials (Fig. 1C). The average error of time estimation by the rate model was 1-2 seconds in rats 2 and 3 and 2-3 seconds in rat 1 (Fig. 1D). Similar errors were observed using the phase model (Supplementary Fig. 1). The accuracy of time-estimation from the models increased with the number of cells used in each animal (not shown). This observation suggests that by recording from a much larger fraction of hippocampal neurons, the accuracy of time estimation can be improved further. It also suggests that a greater amount of information is available to structures downstream from CA1 than what was available for our method of inference. Since we obtained timing precision on a behaviorally relevant scale, we hypothesize that the brain could use population activity to estimate elapsed time.
To investigate how behavioral relevance might influence timekeeping, we also recorded the activity of CA1 cells during running in a wheel placed adjacent to the homecage of the rat (control condition). The rat could enter the wheel and run at its leisure and was not required to keep track of elapsed time. While the patterns of neuronal activity on individual runs displayed some semblance of sequential activation near the beginning of the wheel run, the overall sequences were not consistent from trial to trial (Supplementary Fig. 2); as a result, the time prediction models did not yield any statistically significant prediction of time on the control data (not shown).
While the reliability R(t) of sequential activity during the memory task was well above 0 (the expected value for random, unreliable data; see Methods) (Fig. 1E), R(t) was not significantly positive for the control (homecage) data (Fig. 1F). This suggests that the patterns of neuronal activity during wheel runs reflect timing information only in behaviorally relevant contexts.
The experimentally observed sequences have several important features that make them particularly suitable for timekeeping. First, they are internally-generated; that is, the sequences are not brought about by changing, temporally-structured environmental or body-derived inputs. Second, the sequences are reliable from trial to trial, which allows for time inference on single trials. Third, the sequences are context-dependent, with ‘left’ and ‘right’ trials producing different sequences. Fourth, the sequences are long, on the order of 10-20 seconds, lasting for the entire duration of the required delay period. To explore how the entorhino-hippocampal circuit might be capable of generating these sequences, we considered a minimalistic model of a recurrent network that we call the Threshold Adaptation Model (TAM). This model captures all of the above properties.
We began with a firing rate model (Methods, equation (1)). Firing rate models provide fairly accurate descriptions of the dynamics of large recurrent networks when exact spike timing is not important. In our model, the recurrent network structure is a superposition of two components: J = J0 = Jhet, where J0 is the ‘correlated’ component and Jhet is a component uncorrelated across neurons. J0 has connection strengths that depend only on the distance between neurons arranged on a two-dimensional torus-like grid (Figure 2A) and follow a ‘Mexican-hat’ pattern of connectivity reflecting short-range excitation and long range effective inhibition. Jhet is a random matrix that adds heterogeneity to the pattern of connections. Although the synapses represented by J0 are chosen to be weak relative to Jhet (see Supplementary Text), they strongly influence network dynamics due to their highly correlated structure. The matrix Jhet disrupts the perfect symmetry of the Mexican-hat connectivity; this ensures that our results do not depend on the fine-tuned symmetry of J0, as this symmetry is unrealistic and may produce misleading results (Zhang, 1996; Seung et al., 2000; Renart et al., 2003).
Mexican-hat connectivity, as in J0, and the associated continuous attractor dynamics have been hypothesized as an underlying network mechanism of spatial working memory, spatial navigation and path integration (Samsonovich and McNaughton, 1997; Tsodyks, 1999; Constantinidis and Wang, 2004; McNaughton et al., 2006; Burak and Fiete, 2009). For a wide range of perturbed Mexican-hat connectivity matrices J = J0 + Jhet, the network activity will quickly converge to a “bump attractor” in the presence of constant input (Seung et al., 2000; Renart et al., 2003). If the input stays approximately constant, the bump will not move. Therefore, the dynamics of equation (1) alone cannot produce self-generated sequences, since the “bump” only moves in response to significant changes in external inputs.
To overcome this limitation, we added an activity-dependent adaptation of the spike thresholds (Methods, equation (2)). Threshold adaptation in hippocampal pyramidal neurons was observed experimentally (Henze and Buzsaki, 2001), and evolves on a relatively slow timescale (~1s). In our model, threshold adaptation has the important effect of destabilizing the “bump” attractor states of the network. In the fast-timescale dynamics (equation (1)), the system still evolves to a bump attractor, but as the firing rates of the neurons in the bump of activity increase, so do the corresponding thresholds, and this in turn decreases each neuron's ability to continue firing. The threshold adaptation thus forces the bump to move to a new location, at which point the process repeats itself, resulting in a continuously-moving bump that never stabilizes (Figure 2A, C; see also Supplementary Movie). The moving bump of activity is what produces sequential firing of the neurons.
To understand how context-dependence of sequences may arise in our model, we considered bump trajectories under different initial conditions, mimicking the ‘left’ versus ‘right’ trial types in the behavioral task. Similar to the behavioral experiments, we kept the firing rate initial conditions identical in all simulations, initiating them as a bump of activity in the center of the grid of neurons (Figure 2C, top left). Different initial conditions thus differed only in the initial values for the threshold adaptation variables, as these depend on the recent spiking history of the neurons. Initial conditions A, B, C, and D represents neurons with adapted thresholds to the left, bottom, right, and top of the activity bump, respectively (Figure 2B; see also Supplementary Text). In each case, the activity tends to move away from the neurons with adapted thresholds, as it is more difficult for these to fire. Different initial conditionsthus produce different center-of-bump trajectories. The precise contours of these trajectories are determined by the synaptic heterogeneities Jhet. Importantly, the trajectories for the same initial condition are reliable across trials despite being driven by noisy, temporally and spatially unstructured inputs (Figure 2D). Different instances of the heterogeneous connectivity matrixJhet also produce reliable, context-dependent bump trajectories, while the trajectories vary significantly between matrices (Supplementary Fig. 3).
To better assess the length and reliability of cell assembly sequences produced by the model, we simulated multiple ‘task’ trials for each of four initial conditions (A, B, C and D). Reliability was quantified using the same reliability measure as in Figure 1E,F (see Methods). Sequential activity on the order of 15-20s can be seen for different initial conditions (Figure 3A,B), and the sequences are quite different (Figure 3C). As expected from the reliability of bump trajectories (Figure 2D), sequential activity on single trials closely resembled that of the average (Figure 3D,E) and had a high degree of reliability (Figure 3G,H). Average and single-trial bump trajectories also showed reliability in the task conditions (Figure 3J,K). Note that the timescale of the adaptation controls the speed of the moving bump. A shorter timescale would result in a faster-moving trajectory, leading to shorter sequences (not shown). We also simulated a ‘homecage’ condition (see Methods) that resulted in large trial-to-trial variability, which bore little resemblance to the average activity across trials (Figure 3F,I). The lack of sequential structure in the homecage condition resulted from the lack of consistency in initial conditions, and the temporally correlated noise (‘task’ trials had uncorrelated noise). Finally, we investigated the reliability of sequences as a function of the strength of the input noise to each neuron. As expected, sequences generated in the presence of 5-fold and 10-fold increases in input noise were less reliable, and the reliability decreased more rapidly as a function of time (Figure 3L). Average and single-trial sequences in these higher noise conditions are shown in Supplementary Figure 6. In summary, the sequences for all ‘task’ (but not ‘homecage’) initial conditions were long-lasting and reliable, making them suitable for accurate time estimation to tens of seconds.
Although we initially introduced heterogeneities Jhet to our matrix of synaptic weights to ensure that our results did not depend on the fine-tuned symmetry of the Mexican-hat connectivity matrix J0, we found that synaptic heterogeneity had the unexpected benefit of lengthening the sequences (Supplementary Figure 4). This is because the heterogeneities ‘carved out’ irregular bump trajectories, allowing the bump to travel for a longer period of time without repeating itself. To verify that our results did not depend on the particular instance of heterogeneity we chose for the matrix J1 = J0 + Jhet, we repeated the analyses in Figure 3 and Supplementary Figure 4 using two more instances of Jhet to obtain matrices J2 and J3 (Supplementary Figure 5.)
We have shown that reliable and context-dependent sequences of neuronal activation similar to what we have observed in CA1 may arise from a recurrent network with torus-like architecture and a weakly correlated pattern of Mexican-hat connectivity. However, the architecture of the CA1 region, with its super-sparse recurrent excitation, does not fit with this pattern of connectivity. For this reason, we investigated whether or not reliable sequences generated in one layer can be ‘inherited’ by a downstream layer. In contrast to the first layer, the second layer we devised had no recurrent excitation and only a global, non-specific recurrent inhibition. Layer 2 was driven by both the output of the previous layer and noisy, temporally and spatially unstructured inputs. The feedforward connections between the first, torus-like layer and the second layer were random and sparse (see Methods). Figure 4 shows that despite the lack of structure in layer 2 the sequential activity from the first layer was perfectly ‘inherited’ by the second layer with a similar reliability profile. The reliability remained unchanged even when the magnitude of the noisy inputs to the second layer was increased 5- or 10-fold (see Supplementary Figure 7).
Neural correlates of elapsed time on a suprasecond scale have been documented in several cortical regions (Kojima and Goldman-Rakic, 1982; Fuster, 2001; Brody et al., 2003; Janssen and Shadlen, 2005; Lebedev et al., 2008; Mita et al., 2009). Surprisingly, hippocampal circuits have not been considered as timers, despite the critical role of the hippocampus in timing behavior (Clark and Isaacson, 1965; Young and McNaughton, 2000) and the key importance of temporal context in episodic memory (Tulving, 1972) and navigation (McNaughton et al., 1996).
Our experimental observations and modeling results suggest that in the hippocampus the same cell populations that keep information about past memories and planned travel directions of the animal (Pastalkova et al., 2008) also provide information about elapsed time. Elapsed time was reliably inferred from the population firing rate vector of the recorded neurons at any time point of the memory task. Although the time estimation error from the neuronal population increased over time, it did not increase proportionally to the duration of elapsed time, in contrast to Weber's law (cf. (Staddon, 2005)). In our network model (TAM), the sequences emerged as a natural byproduct of a network with a perturbed Mexican-hat connectivity pattern and adaptation of the spike thresholds. The TAM model does not involve learning of sequences, synfire chains, or a ‘hidden’ feedforward network structure (Abeles, 1991; Levy et al., 2005; Ganguli et al., 2008; Liu and Buonomano, 2009; Fiete et al., 2010). In principle, however, any network that allows for self-sustained, sequential activation of neurons can potentially be a substrate for time keeping. The adapting spike threshold mechanism was adopted because of its simplicity and because dependence of the spike threshold on prior spiking activity has been demonstrated experimentally (Henze and Buzsaki, 2001). We emphasize though that other mechanisms such as short-term synaptic depression may play a similar role (Abbott and Regehr, 2004). Although the mechanism for sequence generation in our model relies on connectivity patterns unlikely to be present in CA1, we have shown that sequences generated in one area with this kind of architecture (potentially in the entorhinal cortex or CA3) can be inherited by another area via sparse, random connections.
The model we have described here may be specific to the hippocampal system, where time keeping is needed on the scale of tens of seconds, and is different from the timing mechanisms in sensory and motor systems. Since evolving neuronal assemblies, or sequences, have been observed in other systems (Luczak et al., 2007; Fujisawa et al., 2008; Long and Fee, 2008; Johnson et al., 2010), it is possible that our modeling results apply to them as well. In general, our findings support the view that each neuronal system generates its own timing, providing temporal frames for its operations (Buonomano and Karmarkar, 2002). In summary, we found that a simple network mechanism can generate long-lasting, reliable sequences that may be used for timekeeping in the hippocampus.
V.I. was supported by the Swartz foundation and NSF DMS 0967377.
C.C. was supported by NSF DMS 0920845 and a Courant Instructorship.
E.P. was supported by the Patterson Trust and HHMI.
G.B. was supported by J.D. McDonnell Foundation, NSF SBE 0542013, NIH NS034994, and MH54671.