Abstract
Telling time is fundamental to many forms of learning and behavior, including the anticipation of rewarding events. Although the neural mechanisms underlying timing remain unknown, computational models have proposed that the brain represents time in the dynamics of neural networks. Consistent with this hypothesis, changing patterns of neural activity dynamically in a number of brain areas—including the striatum and cortex—has been shown to encode elapsed time. To date, however, no studies have explicitly quantified and contrasted how well different areas encode time by recording large numbers of units simultaneously from more than one area. Here, we performed large-scale extracellular recordings in the striatum and orbitofrontal cortex of mice that learned the temporal relationship between a stimulus and a reward and reported their response with anticipatory licking. We used a machine-learning algorithm to quantify how well populations of neurons encoded elapsed time from stimulus onset. Both the striatal and cortical networks encoded time, but the striatal network outperformed the orbitofrontal cortex, a finding replicated both in simultaneously and nonsimultaneously recorded corticostriatal datasets. The striatal network was also more reliable in predicting when the animals would lick up to ∼1 s before the actual lick occurred. Our results are consistent with the hypothesis that temporal information is encoded in a widely distributed manner throughout multiple brain areas, but that the striatum may have a privileged role in timing because it has a more accurate “clock” as it integrates information across multiple cortical areas.
SIGNIFICANCE STATEMENT The neural representation of time is thought to be distributed across multiple functionally specialized brain structures, including the striatum and cortex. However, until now, the neural code for time has not been compared quantitatively between these areas. Here, we performed large-scale recordings in the striatum and orbitofrontal cortex of mice trained on a stimulus–reward association task involving a delay period and used a machine-learning algorithm to quantify how well populations of simultaneously recorded neurons encoded elapsed time from stimulus onset. We found that, although both areas encoded time, the striatum consistently outperformed the orbitofrontal cortex. These results suggest that the striatum may refine the code for time by integrating information from multiple inputs.
Introduction
Anticipating events that will happen in the future is among the most important functions the brain performs. Indeed, it has been increasingly stressed that learning and memory are prospective brain functions; that is, they are only adaptive to the extent that they help animals anticipate and prepare for the future (Dudai and Carruthers, 2005; Schacter and Addis, 2007). To anticipate when events will happen, the brain has evolved mechanisms to tell time across a wide range of temporal scales (Buhusi and Meck, 2005; Buonomano, 2007).
Timing on the scale of hundreds of milliseconds to a few seconds is of particular importance in that it allows animals to predict and prepare for events unfolding within the immediate future. Within this range, animals discriminate the temporal features of sensory stimuli, such as those used for communication, and generate timed motor responses to prepare for external events such as expected rewards. The neural mechanisms underlying the brain's ability to tell time on the scale of seconds remains unknown (Mauk and Buonomano, 2004; Merchant et al., 2013a), but a rapidly growing literature has reported that dynamically changing patterns of neural activity encode information about the amount of time elapsed since a given stimulus. These patterns of activity, which have been referred to as population clocks (Buonomano and Karmarkar, 2002; Buonomano and Maass, 2009; Buonomano and Laje, 2010), have now been observed in a wide range of different brain areas, including the striatum (Matell et al., 2003; Chiba et al., 2008; Jin et al., 2009; Gouvêa et al., 2015; Mello et al., 2015; Bakhurin et al., 2016), prefrontal cortex (Brody et al., 2003; Oshio et al., 2008; Genovesio et al., 2009; Jin et al., 2009; Merchant et al., 2011; Kim et al., 2013; Carnevale et al., 2015), parietal cortex (Janssen and Shadlen, 2005; Crowe et al., 2010), and hippocampus (Pastalkova et al., 2008; Kraus et al., 2013), as well as in the bird song system (Hahnloser et al., 2002; Long et al., 2010). In addition, pharmacological, lesion, and neuroimaging work suggests a role of the basal ganglia (Meck, 1996; Coull et al., 2011) and prefrontal cortex (Dietrich and Allen, 1998; Kim et al., 2009; Xu et al., 2014) in timing.
The diversity of areas implicated in timing likely reflects the range of tasks and temporal scales examined, but it is also possible that, even within the same task, different areas track time in parallel (Matell et al., 2003; Jin et al., 2009). To date, however, no single study has quantified directly the degree to which two different circuits encode time through simultaneous multiple-region recordings. Here, we contrast directly the ability of two circuits, the striatum and orbitofrontal region (OFC) of the prefrontal cortex, to encode time.
We examined the neural representation of time during a Pavlovian conditioning task in which a food reward is presented at a specific interval after a conditioned stimulus (CS). Mice exhibited anticipatory licking during the fixed cue–reward delay period. Silicon microprobe recordings of dozens of units from either the striatum or OFC or both simultaneously revealed that population activity in both circuits encoded an internal representation of elapsed time. This code was quantified by feeding the trial-by-trial spike pattern into a pattern classifier and training it to read out elapsed time. The quality of the striatal population code for time was significantly better than that of the OFC. Our results support the hypothesis that many different brain areas encode time simultaneously, but the striatum may play a privileged role in timing relative to the OFC because it holds a more accurate clock. We hypothesize that, by sampling the changing patterns of activity unfolding throughout the cortex and other inputs continuously, the striatum implements a robust code for elapsed time via a temporal “winners-take-all” mechanism.
Materials and Methods
Animals and surgical procedures.
All procedures were approved by the University of California–Los Angeles Chancellor's Animal Research Committee. Singly housed male C57BL/6J mice (n = 11, 15–22 weeks old at the time of recording; The Jackson Laboratory) were used in the experiments. Animals underwent an initial head bar implantation surgery under isoflurane anesthesia in a stereotaxic apparatus to fix stainless steel head restraint bars bilaterally on the skull with dental cement. After training, animals underwent a second surgery under isoflurane anesthesia on the recording day to make craniotomies for acute microprobe recordings. An additional craniotomy was made over the posterior cerebellum for placement of an electrical reference wire. All behavioral training and recording sessions were performed in fully awake, head-restrained animals.
Behavioral task.
After a 1 week recovery period following the initial head bar implantation surgery, animals were food restricted and fed daily after each training session to maintain ∼90% of their baseline weight. Water access was ad libitum. During daily training sessions, animals were mounted on the head bar restraint bracket on the recording rig and stood on a polystyrene treadmill ball (200 mm diameter; Graham Sweet Studios) that rotated along a single axis during forward/backward ambulation. Animals were initially habituated to the head-fixed recording rig and trained to consume a liquid reward (5 μl, 10% sweetened condensed milk). The reward was delivered from a tube positioned between an infrared lick meter (Island Motion) by actuation of an audible solenoid valve (Neptune Research). During daily reward-only training sessions, animals consumed 100 rewards and were exposed to a constant stream of pure air through a tube positioned next to the nose [100 rewards per session, 13–21 s intertrial interval (ITI), sampled from a normal distribution, 1.5 L/min air flow]. Once animals could consume ≥90% of the rewards for 2 consecutive days, they began conditioning with olfactory cues using an olfactometer. Odorants were introduced by bubbling air (0.15 L/min) through aromatic odorants diluted 1:10 in mineral oil (Sigma-Aldrich) and merging this product with the 1.5 L/min stream of pure air. The constant flow of pure air into which odors are introduced decreased the possibility that animals used decaying concentrations of odorant as a temporal cue. During daily training sessions, animals received pseudorandom presentations of each odor stimulus (1 s duration, 17–29 s ITI, sampled from a normal distribution). Isoamyl acetate served as the CS+ odor because its offset was followed by a 1.5 s delay and a reward delivery. Citral served as the CS− odor because it was not followed by any explicit outcome. Animals received 100 presentations of each trial type in random order during each training session. The solenoid valves controlling the odors were sound isolated and thus inaudible to the mouse. Typically, during the first or second day, animals began predicting the delivery of the reward by licking in anticipation during the interval between the odor and the reward. Correct CS+ trials were defined as those trials during which licking was initiated before reward delivery (between 0.7 and 2.5 s after stimulus onset). Correct CS− trials were defined as those containing no licking activity for 5 s after stimulus onset. False alarm CS− trials were defined as those trials during which licking was initiated between 0.7 and 2.5 s after stimulus onset. Once animals demonstrated correct responding on ≥90% of trials, they underwent surgery for recording. During the recording session, animals received 100 CS+ trials with 85% reward probability and 100 CS− trials. Animals performed between 54 and 99 correct CS+ trials and between 1 and 56 false alarm CS− trials.
Electrophysiological recordings.
Procedures for developing and recording with silicon microprobes have been described previously (Shobe et al., 2015). One recording was performed per animal. Each area was targeted with a silicon microprobe containing a total of 256 electrodes that were divided across four or five prongs. The electrodes spanned between 0.825 to 1.05 mm of the distal tip of the prongs. Data in this study were aggregated from two groups of animals. In the first group (n = 5), recordings took place in the anterior striatum only (silicon prong tip positions: 1.2 mm anterior, 0.8 to 2.2 mm lateral, −3.4 to −5.7 mm ventral relative to bregma). In the second group (n = 6), we simultaneously recorded from the orbitofrontal region of the prefrontal cortex (2.2 mm anterior, 0.26 to 2.05 mm lateral, −3.6 mm ventral relative to bregma) and both the anterior and posterior regions of the striatum (anterior striatum: 1.2 mm anterior, 0.78 to 2.1 mm lateral, −5 mm ventral; posterior striatum: −0.5 mm anterior, 2.4 to 3.2 mm lateral, −4.3 mm ventral tip position relative to bregma) using multiple 256 electrode probes attached together (Shobe et al., 2015). Therefore, the striatal dataset analyzed in this study was composed of the anterior striatal recordings performed in the first group, combined with anterior and posterior striatal recordings performed in the second group. The OFC dataset was composed of orbitofrontal recordings performed in the second group. Because of the wide spatial distribution of recording sites above the prong tips, the anterior striatal dataset contained units sampled from both dorsal and ventral striatal areas. Positions of units included in analysis are illustrated in Figure 5A. Spike sorting was performed on the data using custom, semiautomated software written in MATLAB (The MathWorks). The placement of silicon probes was confirmed histologically at the end of each experiment by coating the prongs with a fluorescent dye (Di-D; Thermo Fisher) before implantation.
Delineation of anatomical subregions.
In each animal, the recordings in the anterior striatum consisted of predominantly ventral or dorsally positioned units, with one recording containing units evenly distributed in each area. The mean electrode position of −4.2 mm DV was used to divide the anterior striatal recordings into dorsal or ventral regions. To divide the OFC into medial and lateral subregions, the mean electrode position of 1.19 mm in the OFC was used.
Unit classification.
Analysis was performed on putative principal neuron populations; that is, pyramidal cells in the OFC and medium spiny neurons (MSNs) in the striatum. Spike waveform trough-to-peak duration was used to distinguish putative MSNs and pyramidal neurons from nonprincipal neurons. Putative fast-spiking interneurons (FSIs) were separated from principal cells in both the OFC and the striatum by their narrow waveform (maximum FSI ttr-pk = 0.475 ms, minimum principal neuron ttr-pk = 0.55 ms, and maximum MSN ttr-pk = 1.25 ms; Bakhurin et al., 2016). A measure of firing rate regularity (coefficient of variation, CV) was also used to exclude putative tonically active neurons from the striatal recordings (maximum CV = 1.5; Bennett and Wilson, 1999). A total of 690 putative MSNs of a total of 1115 striatal units and 505 putative pyramidal cells of a total of 654 cortical units were recorded.
Identification of lick-modulated units.
Licking modulated units were determined by correlating estimated firing rates with licking rate around lick episodes that occurred throughout the recording, including within and outside of trial periods. Licking episodes were defined as containing 2 licks that were separated by at most 250 ms (4 Hz). Licking episodes could not occur within 5 s of each other. To calculate the correlations, individual licks occurring within a 2 s window around each lick episode were binned into 50 ms time bins. For each unit, spikes occurring around each licking episode were binned within a 2 s window into 50 ms bins. The resulting episode vectors reflecting licking and spiking counts for each episode were concatenated into two vectors and convolved using a Gaussian function (SD = 100 ms) to obtain licking and spiking rate estimates across all lick episodes in the recording. A Pearson correlation was performed between the lick-rate vector and each spiking rate vector for each unit. A unit was considered to be lick-rate modulated if it demonstrated a positive correlation coefficient with a p-value <0.01.
Elapsed time prediction analysis.
All analyses were performed independently on data collected from each animal and each brain region using correctly performed CS+ trials or CS− trials with false alarm licking. All decoding models were generated using only simultaneously recorded cells from individual animals. For each trial, neural population activity was analyzed over the 2.5 s interval between cue onset and reward delivery. Over this interval, the activity of each neuron in the simultaneously recorded population was transformed into an analog rate code estimate by: (1) convolving its spike train with a decaying exponential function (τ = 100 ms) and (2) calculating its firing rate estimate as a binned average (100 ms time bins) of its convolved spike train. This procedure resulted in 25 population firing rate vectors, one per 100 ms time bin, in the trial.
Elapsed time was decoded from the population firing rates in each trial by requiring a classifier to label each rate vector in the trial as coming from one of the 25 time bins. The classification task was performed with a multiclass support vector machine (SVM) with a radial-basis function (RBF) kernel, as implemented in the LIBSVM library (version 3.20; Chang and Lin, 2011). This SVM uses a one-against-one multiclass approach to distinguish the population firing rates encoding a given time bin from those encoding each of the 24 other time bins (Knerr et al., 1990; Kreβel, 1999; Hsu and Lin, 2002). In the one-against-one multiclass approach, binary classifiers are trained to distinguish between the population codes for each pair of distinct time bins (i, j), for a total of 300 binary classifiers. SVM output is represented in 25 readout units, one per time bin. Given a test population rate vector, readout i generates a classification score indicating how closely this vector resembles the population code encoding bin i. It is calculated as an aggregate of the outputs of the 24 binary classifiers (i, 1) (i, 2), … (i, I − 1) (i, i + 1), … (i, 25). The SVM predicts that the test vector encodes time bin k whenever readout k produces the highest score of all 25 readout units (see Fig. 2).
Individual animals showed varying numbers of correctly performed trials. To ensure that the decoding performance across animals was compared under equivalent conditions, the predicted time bins in all figures were generated with a Monte Carlo cross-validation strategy. The rate vectors from each trial were tested on 30 independently trained SVMs and each SVM was trained on the rate vectors from M randomly sampled trials excluding the test trial. Because the minimum number of correct CS+ trials for an individual animal was 54, M was chosen to be 53.
The number of simultaneously recorded units used to train and test the models was controlled. The number of simultaneously recorded cells, N, used to generate each decoding model and the number of animals used for averaging is always indicated on the figure or in the figure caption. N varied from 29 to 55 because of subregion-specific limitations in the number of simultaneously recorded units. Furthermore, to test the effect of the population size on model performance (see Figs. 3, 4, and 9), random samples sizes of 5, 10, 15, 20, and 40 units taken from the entire striatum or entire OFC were compared. During each of the 30 repetitions of the Monte Carlo cross-validation, N distinct units from the population were randomized for training and testing. To maximize decoder performance, the RBF SVM regularization parameters were optimized for each brain region of each animal. Specifically, the misclassification cost parameter, C, and the data complexity parameter, γ, were optimized via a grid search with fivefold cross-validation. Across all datasets, the predominant value of C was 4 (range: 1–16) and of γ was 0.25 (range: 0.0156–0.25).
Comparing population coding between correct CS+ and false alarm CS− trials.
To determine the extent to which the CS+ code for time generalized to CS− trials, the classifier was trained in the same way as described above using 55 cells per area and using M = 53 trials per Monte Carlo cross-validation repetition. The models were then tested on the 25 rate vectors generated for each false alarm CS− trial available for each animal (identical binning and rate estimation procedure as done for CS+ trials). This procedure was repeated 30 times and random combinations of 55 units and 53 trials were used in training the model.
Lick onset prediction analysis.
For each trial, neural population activity was transformed into estimated population firing rate vectors using 100 ms bins, as in the elapsed time prediction analysis. This sequence, or trajectory, of neural population activity started 1 s before cue onset and ended 200 ms after the latest lick onset time of all correctly performed CS+ trials. As a result, the number of time bins (and population rate vectors) analyzed per trial varied between 31 and 37 across animals.
Lick onset bins were predicted from the population firing rates in each trial with an RBF SVM binary classifier. The SVM's output is represented by a single readout that scores how closely each population vector in the test population trajectory predicts lick onset. The predicted lick onset bin was the one in which the readout was at its highest value. Testing was performed with a Monte Carlo approach similar to the elapsed time prediction in which each trial was tested on 30 SVMs independently trained on M = 53 randomly sampled trials. The dataset contains disproportionately fewer lick onset bins than nonlick onset bins because only a single bin of the 31 to 37 bins per trial can be a lick onset bin. To avoid the resulting bias in the SVM model, the training set for each SVM was altered by randomly down-sampling the subset of nonlick onset bins by 75% and expanding the set of SVM target bins to include one bin immediately preceding and one bin immediately after the actual lick onset bin in each trial, for a total of three target bins per trial. The misclassification cost and data complexity regularization parameters for the RBF SVMs were optimized for each brain region of each animal similarly to the elapsed time prediction analysis. Across all datasets, the predominant value of C was either 2 or 8 (range: 2–128) and of γ was 0.125 (range: 10−7 to 0.5).
The binary SVMs were retrained for each prelick time to determine how far in advance the neural trajectory could predict lick onset (see Fig. 10D). At each prelick time, the SVMs were retrained to predict a new set of target bins that were appropriately shifted backward in time from the actual lick onset bin. During training, the data down-sampling procedure was accordingly altered to down-sample the subset of nontarget bins.
Trial shuffling.
Trial shuffling was used as a control for elapsed time prediction. This procedure disrupts correlations in simultaneously recorded population activity, but preserves the correct bin order for each unit. To create trial shuffled activity, each unit's firing rate estimate in each time bin of each trial was replaced with the same unit's firing rate estimate in the same time bin of a randomly selected trial. This control could not be performed with the lick onset prediction analysis because of the resulting dissociation between lick onset times and unit activity.
Bin shuffling.
Bin shuffling was used to generate population responses that were dissociated from their correct temporal order. To create bin shuffled activity, each unit's firing rate estimate in each time bin was replaced with the same unit's firing rate estimate in a randomly selected bin of the same trial.
Prediction analyses from trial and bin shuffled data involved training and then testing on the respectively transformed datasets. To allow for direct comparisons between observed, trial-shuffled, and bin-shuffled controls, care was taken to make sure that the same units and trials were subsampled for analysis.
Temporal warping of internal time representation.
Given our hypothesis that the population code for elapsed time and lick-onset time share a common internal timing representation, the population's encoding of the animal's internal representation of time should covary with the lick onset time. To measure this effect, correct CS+ trials were divided into three approximately equally sized sets corresponding to each of the terciles of the animal's lick onset distribution and SVMs were trained to classify elapsed time in the first (third) tercile trials and then tested on the second and third (first) tercile trials. Biases in the resulting error distribution would then reveal an underlying comodulation. A more direct measurement of this effect was performed by comparing the temporal relationship between trial-averaged trajectories of first and third tercile trials; if the third tercile trajectory was consistently slower than the first tercile trajectory, then this would indicate that the two timing variables comodulate one another. Population spike trains were convolved with a Gaussian function (mean = 0, SD = 100 ms) and then trial averaged separately over the first and third tercile trials to produce characteristic first (T1) and third (T3) tercile trajectories. These trajectories were compared by temporally aligning them as follows: (1) a Euclidean distance matrix was constructed by comparing the population code at each moment along T1 to the population code at each moment along T3, resulting in an NT × NT distance matrix (see Fig. 7C), where NT = 2500, given the time resolution of the spike trains (1 ms); (2) T3 and T1 were temporally aligned with a dynamic time-warping procedure that calculated the deviation of T3 from T1 over the course of time as the path along the distance matrix between the beginning and the end of T1 with minimum cumulative distance (see Fig. 7C, black trace). The relative speed (temporal warping) of T3 with respect to T1 was indicated by the difference between the respective times at which the two trajectories were temporally aligned (see Fig. 7D). When T3 ran slower than T1, this difference would be positive and monotonically increase and, when it ran faster, the difference would be negative and monotonically decrease.
Effective dimensionality.
The effective dimensionality of each recorded population was calculated from trial-averaged population firing rate estimates. To control for the difference in the number of units measured across different recordings, the effective dimensionality for each recording was calculated as a mean over 30 randomly sampled subpopulations of size 55. Performing principal component analysis (PCA) on the dynamics of a single such sample produced a list of 55 principal components (PCs) ordered by the percentage of variance in the population dynamics explained by each PC. The effective dimensionality was calculated as the minimum number of PCs required to explain 95% of the variance in the dynamics (Rajan et al., 2016).
Statistical analysis.
To determine the efficacy of the SVM models in elapsed time prediction, the correlation between the correct bin number and the predicted bin number was calculated. A single correlation coefficient was calculated from all the test data (i.e., 25 time bins per correct CS+ trial in the dataset, 30 repetitions each) for a given brain region of a given animal. For clarity, graphs display correlation coefficients, but statistical analysis was performed using Fisher's z-statistic for correlation coefficients (Fisher transformation). For the lick onset time prediction analysis, model accuracy was measured by the root mean squared error (RMSE) of the predicted lick onset bins. A single RMSE value was calculated from all the test data (i.e., 30 repetitions for the lick onset bin in each of the correct CS+ trials) for a given brain region of a given animal. During hypothesis testing, we assumed that the population size used in the analysis represented a repeated measure because units were sampled from the same population of units. Brain region (i.e., striatal vs OFC networks) was considered a repeated measure only when recorded in the same animal (see Figs. 5D, 6B,E, 9, 10). Two-way repeated-measures and mixed-model ANOVA analysis was performed using GraphPad Prism version 6.0 software. Two-sided paired and unpaired t tests were performed using standard functions in MATLAB.
Results
Behavior
We obtained large-scale recordings from the striatum and OFC in head-fixed mice (n = 11 mice) previously trained to perform an odor discrimination task (Shobe et al., 2015; Bakhurin et al., 2016). In this task, mice were presented for 1 s with one of two olfactory stimuli. One of the odors (CS+) was followed by a reward delivered 2.5 s from cue onset. The delivery of the reward was not contingent on any instrumental actions of the animal. The second odor (CS−) was followed by no specific outcome (Fig. 1A). After repeated presentations of the CS+ trials, animals learn to generate anticipatory licking behavior that preceded the reward delivery (Fig. 1B). Previous experiments from our group have demonstrated that animals time their anticipatory licking response depending on the cue–reward delay duration (data not shown), consistent with timed reward-guided behavior found in many other studies (Bermudez and Schultz, 2014). Our recordings were performed in animals that had experienced five to 10 training sessions and were performing above a criterion of at least 90% correctly performed trials (see Materials and Methods) before the recording day. The onset of anticipatory CS+ licking responses was concentrated during the cue–reward delay period for all animals studied (mean lick onset time = 1.8 s, SD = 0.25; Fig. 1C). We focused our analysis on correct CS+ trials because these displayed discrete behavioral evidence that animals timed their behavior to anticipate the reward.
Large-scale recording of OFC and striatal networks during reward-predictive behavior. A, Task schema. Mice received pseudorandomly ordered presentations of a CS+ odor that predicted reward delivery 2.5 s after odor onset and an unrewarded CS− odor. Rectangles represent odor-on time. Red triangle and vertical blue dashed line indicate reward delivery. B, Example of anticipatory licking behavior of one mouse during CS+ trials. Shaded blue rectangle represents odor presentation time. Black tick marks indicate individual licks and red ticks denote lick onset times that are used for subsequent analysis. Trials are sorted by descending latency to first lick. C, Cumulative distributions of lick onset times during CS+ trials for all mice included in the study (n = 11 mice). D, Distribution of the trough-to-peak width (milliseconds) recorded from striatal units (top) and OFC units (bottom). Vertical dotted lines depict the threshold margin (0.475 to 0.55 ms) for segregating putative FSIs (red histograms) from putative principal cells (striatal MSNs and OFC pyramidal cells, blue histograms). Gray bars reflect unclassified cells. E, Individual population-level recordings from the striatum (top) and the prefrontal cortex (bottom) during correctly performed CS+ trials. Each row in a matrix represents the mean normalized firing rate of one recorded putative projection neuron in the corresponding brain area. Units are sorted by their latency to maximum firing rate. Blue rectangles indicate CS+ odor presentation time and red triangles mark the time of reward delivery.
Large-scale striatal and orbitofrontal recordings
After animals reached criterion performance on the task, we used silicon microprobes (Shobe et al., 2015) to record population activity from the striatum, OFC, or simultaneously from both of these areas as the mice performed the task. We focused our analysis on putative principal cells in these brain regions: striatal MSNs and cortical pyramidal cells. If these brain areas contain a code for time, then principal cells would be the most likely to transmit that signal to downstream brain regions (Buonomano and Merzenich, 1995). To identify these populations, we measured the action potential duration of each unit and used a threshold margin to segregate putative principal cells from FSIs. In both the striatum and OFC, the distribution of spike widths across all cells was bimodal (Fig. 1D). Based on the separation of these distributions, we only included putative principal cells in our analysis. We analyzed data from animals containing at least 55 principal units per region (n = 9 striatal recordings, and 6 OFC recordings). Our datasets contained between 55 and 120 simultaneously recorded principal neurons. We found that, on average, the population of striatal and prefrontal neurons exhibited highly heterogeneous firing activity during the cue–reward interval (Fig. 1E). This observation is qualitatively similar to the sequential firing patterns reported from other cortical recordings (Crowe et al., 2010; Harvey et al., 2012; Stokes et al., 2013) and striatal recordings (Gage et al., 2010; Thorn and Graybiel, 2014; Mello et al., 2015; Rueda-Orozco and Robbe, 2015; Bakhurin et al., 2016). However, from the average firing rate representation, it was not evident whether the dynamics were robust at the single-trial level. We thus examined whether it was possible to decode elapsed time and lick onset time on a trial-by-trial basis during correctly performed CS+ trials.
Decoding time from network dynamics
To investigate the possibility that neural network activity could provide a mechanism for the stable representation of time, we used a SVM decoder to detect and measure the reoccurrence of dynamic population activity in striatal and OFC networks on a trial-by-trial basis. An SVM was trained to identify population activity in each of the 25 time bins (100 ms bin duration) between stimulus onset and reward delivery.
Each unit's firing rate for the 25 time bins of a given trial was estimated from its spike train over that trial (see Materials and Methods). Next, population firing rate dynamics across multiple trials were used to train the SVM classifier (Fig. 2A). During testing, population activity from time bins of novel trials were presented to the trained SVM. SVM output for the population activity in a given time bin was represented by a vector of values generated by 25 readout units, where each readout value i represented a prediction score that the input pattern was from time bin i. This resulted in a vector of 25 readout values per test time bin (Fig. 2B). For each test time bin, the SVM predicted its bin label as the index of the maximal readout (Fig. 2C). Testing was performed with a Monte Carlo cross-validation approach that controlled for the variance in the number of trials and size of the simultaneously recorded population across brain regions and animals (see Materials and Methods).
Schematic of the SVM decoding of elapsed time. A, Training the SVM. Single-trial spiking activity of each unit in a simultaneously recorded population (only 3 units represented) is transformed into a firing rate estimate for the unit during the 2.5 s interval after odor presentation onset (data not shown here). The rate estimates are binned (100 ms time bins) to construct 25 population activity patterns per trial. Using a one-against-one multiclass strategy, the SVM trains a set of binary classifiers to distinguish the population activity pattern in each time bin from every other time bin. SVM output is conceptualized as 25 readout units, one per target time bin, that learn to distinguish activity patterns in their respective target time bin from those in all other bins. B, The model is tested using a Monte Carlo cross-validation approach in which each activity pattern from novel trials (i.e., those excluded from the training set) is tested on trained SVM models. Illustrated is the testing of bin #2 of the test trial. C, Readout units score each test activity pattern for how closely it corresponds to their respective target bins. The target time bin of the readout with the maximal value is chosen as the predicted time in a winners-take-all manner (marked with a red vertical line). Actual readout values are depicted here.
Elapsed time encoding by striatal and cortical networks
We first investigated the ability of striatal MSN dynamics during single CS+ trials to be sorted into the correct temporal order by an SVM. Strong SVM performance would suggest that striatal neuron populations stably encode an internal representation of time elapsed from stimulus onset and may provide a mechanism by which downstream regions could read out temporal information from striatal activity. We found that the highest SVM readout values during testing generally fell along the diagonal line in single trial cross-temporal classification matrices (Fig. 3A).
Striatal networks encode elapsed time. A, Cross-temporal classification matrices visualizing SVM model performance on striatal network data recorded during individual correctly performed CS+ trials. Each column represents the normalized readout values normalized across SVM readout units for the activity pattern from the corresponding correct time bin (x-axis). Peaks in each column reflect the predicted time chosen by the model. The black dotted line lies along the diagonal. B, Top, Average of classification matrices generated across all correct CS+ trials for one striatal recording. Center, Average classification matrix across all correct CS+ trials after bin shuffling each unit's activity in the same recording. Bottom, Average classification matrix across all correct CS+ trials after trial shuffling each unit's activity in the same recording. C, Scatter plot of predicted versus correct time bins across 80 correctly performed CS+ trials for one striatal recording. Predicted bin numbers (y-axis) were jittered (Gaussian noise, mean = 0, SD = 0.2) to separate overlapping points. The blue solid line represents the regressed line describing the correlation between actual and predicted time. The red dotted line lies along the identity line. D, Mean correlation coefficients between predicted and correct time bins across all striatal recordings (55 units per animal, n = 9) for observed, bin-shuffled, and trial-shuffled data types. SVM classification of population activity was repeated 30 times (see Materials and Methods). SVM models trained on trial-shuffled activity performed better than when trained on observed (nonshuffled) activity patterns (p = 0.023, paired t test). Bin-shuffled models performed at chance level significantly worse than nonshuffled models (p < 0.0001, paired t test). E, Comparison of SVM performance using nonshuffled and trial-shuffled network activity as a function of the number of units used for training and testing. There was a significant effect of data type (F(1,8) = 7.9, p = 0.023) and number of units (F(6,48) = 109.7, p < 0.0001, two-way repeated-measures ANOVA). F, Bin-shuffled models performed worse than nonshuffled models for each population size used in the model (F(1,8) = 178.0, p < 0.0001, two-way repeated-measures ANOVA). Error bars indicate SEM.
Figure 3B, top, illustrates the average classification matrix over all trials in a single striatal recording and reveals the presence of a time code in the recorded dynamics. For each recording, we repeated the analysis on two different control patterns. First, to evaluate the temporal encoding efficacy of striatal population dynamics, we trained and tested an SVM on the dataset after scrambling its temporal dynamics by bin shuffling, wherein the sequence of firing activity for each unit within each trial was shuffled independently (see Materials and Methods). The control confirmed that bin shuffling completely eliminated the ability of the SVM to identify a code for elapsed time in the population activity (Fig. 3B, center). Because the above analysis was based on a set of 55 simultaneously recorded cells, we are able to determine the effect of noise correlations on the time code (Nirenberg et al., 2001; Schneidman et al., 2003; Averbeck and Lee, 2006; Averbeck et al., 2006); in other words, does decoding based on simultaneously recorded cells hamper or improve performance. To do this, we measured decoding performance after independently shuffling the firing activity of each unit across trials (see Materials and Methods). Although bin-shuffling population activity rendered time bin predictions entirely random, trial-shuffled controls performed very similarly to models trained on observed data (Fig. 3B, bottom).
To quantify performance and the effects of bin and trial shuffling on the quality of the time code, we calculated the Pearson correlation coefficient between the correct and predicted time bin values in each recording (Fig. 3C). Across all striatal datasets, population dynamics were highly predictive of elapsed time during the task (mean Pearson correlation coefficient = 0.85, SD = 0.069, n = 9). However, whereas bin shuffling reduced time prediction to chance levels, elapsed time decoding performance on the trial shuffled control was slightly, but significantly, better than on the observed data (mean Pearson correlation coefficient = 0.88, SD = 0.064, p = 0.024, paired t test on Fisher transformed coefficients; Fig. 3D). This shows that the neurons are not noise independent and that the noise correlations, the within-trial correlations between neurons, impair decoding.
Next, we examined whether performance was dependent on the size of the striatal population used in decoding time. A two-way, repeated-measures ANOVA detected a significant effect of population size on classification performance (F(6,48) = 109.7, p < 0.0001). The analysis again revealed a significant effect of trial shuffling (F(1,8) = 7.9, p = 0.023; Fig. 3E). These results show that, under physiological conditions, striatal noise correlations are detrimental for neural coding of elapsed time, in agreement with the detrimental role of correlations found in other studies (Averbeck and Lee, 2006; Averbeck et al., 2006; Cohen and Maunsell, 2009; Mitchell et al., 2009; Tremblay et al., 2015). A separate two-way, repeated-measures ANOVA comparing observed and bin-shuffled data at different population sizes also revealed a significant effect of bin shuffling (F(1,8) = 178.0, p < 0.0001; Fig. 3F).
We next applied these same analyses to OFC pyramidal cell dynamics using the same procedures and numbers of cells. We found that OFC network dynamics also encoded elapsed time during the task (Fig. 4A). Interestingly, in contrast to the striatal code, we found no significant difference in the encoding efficacy between observed (mean Pearson correlation coefficient = 0.7, SD = 0.104, n = 6) and trial-shuffled OFC network activity (mean Pearson correlation coefficient = 0.72, SD = 0.13, p = 0.21, paired t test; Fig. 4B). Although a two-way, repeated-measures ANOVA demonstrated that decoding performance using OFC population dynamics also depended on population size (F(6,30) = 49.6, p < 0.0001), the analysis did not show a significant difference between the encoding efficacy of observed and trial-shuffled data (F(1,5) = 2.4, p = 0.18; Fig. 4C). These findings suggest that temporal encoding in OFC dynamics is potentially less sensitive to noise correlations than in the striatum. Finally, as expected, bin shuffling the OFC dynamics eliminated the temporal code and resulted in chance level decoding performance (two-way, repeated-measures ANOVA, F(1,5) = 109.5, p < 0.0001; Fig. 4D).
Striatal networks encode elapsed time better than OFC networks. A, Average cross-temporal classification matrix across all correct CS+ trials for one OFC recording. Color scale is the same as in Figure 3B. B, Mean correlation coefficients across all OFC recordings (55 units per animal, n = 6) for observed, bin-shuffled and trial-shuffled data types. SVM classification of population activity was repeated 30 times (see Materials and Methods). SVM models trained on trial-shuffled activity were not significantly different from those trained on nonshuffled activity patterns (p = 0.21, paired t test). Bin-shuffled models performed at chance level and significantly worse than the nonshuffled models (p < 0.0001, paired t test). C, Comparison of SVM performance using nonshuffled and trial-shuffled network activity as a function of the number of units. There was no significant effect of data type (F(1,5) = 2.4, p = 0.18), but we observed a significant effect of the number of units (F(6,30) = 49.6, p < 0.0001, two-way repeated-measures ANOVA). D, Bin-shuffled models performed worse than nonshuffled models for each population size used in the model (F(1,5) = 109.5, p = 0.0001, two-way repeated-measures ANOVA). E, Comparison of SVM model performance between all striatal and OFC recordings (55 units per region, n = 9 striatal recordings and 6 OFC recordings) showed that the classification performance of models trained on striatal network data was significantly better (p = 0.0092, unpaired t test). F, Mean performance of SVM classification as a function of number of units used in training and testing for each brain region. A mixed-model ANOVA revealed a significant effect of number of units (F(5,65) = 191.9, p < 0.0001) and a significant effect of brain region (F(1,13) = 9.0, p = 0.01). The ANOVA excluded the “all units” column because it contained inconsistent numbers of cells between regions. Error bars indicate SEM.
Striatal networks outperform prefrontal networks in encoding elapsed time
Consistent with the striatal results above, other experimental studies have reported the presence of a time code in the striatum (Gouvêa et al., 2015; Mello et al., 2015). We found that OFC networks also encode time, suggesting that this information is distributed throughout multiple brain areas. An important and unaddressed question pertains to the relative quality of this neural code in the striatum and OFC. We thus compared the performance of OFC and striatal network dynamics in encoding elapsed time. SVM classification performance was significantly better when trained and tested on striatal activity than on OFC activity (p = 0.0092, unpaired t test; Fig. 4E). A two-way mixed-model ANOVA between brain region and population size revealed that this effect was consistent across a broad range of population sizes (F(1,13) = 9.5, p = 0.01; Fig. 4F). These results suggest that striatal networks show a significantly more robust representation of time compared with the OFC.
Dorsal and ventral striatum equally encode elapsed time
In the above analysis, we adopted an unbiased approach for quantifying temporal coding in the striatum in that we incorporated units from both anterior and posterior areas of this structure (Fig. 5A, left and center). Most of our recorded units were from the anterior striatum, but it is unclear to what extent this subregion by itself contained a better neural code than the OFC. We therefore repeated our comparative analysis after excluding posterior striatal MSNs (this reduced the minimum number of simultaneously recorded cells from 55 to 48). We found that the anterior striatum alone continued to have an improved code for time over the OFC (Pearson correlation coefficients: mean anterior striatum = 0.83, SD = 0.083, n = 9; mean OFC = 0.69, SD = 0.083, n = 6, p = 0.0083, unpaired t test; Fig. 5B). Next, we focused on differences between dorsal and ventral areas of the anterior striatum. Few studies have investigated whether the ventral striatum encodes time; however, the role of this area in reward prediction suggests that it may have a code for time. We took advantage of our widely distributed recording positions to compare the decoding performance of dorsal and ventral striatum MSNs (we used datasets with at least 35 simultaneously recorded MSNs). We found that SVM models trained on dorsal or ventral units performed as well as models trained with units taken randomly from either dorsal or ventral areas (F(2,16) = 0.02, p = 0.98, one-way ANOVA; Fig. 5C). Together, these results suggest that the quality of temporal coding appears to be evenly distributed across the striatum and that this area consistently outperforms the OFC.
Population encoding of elapsed time is distributed throughout striatum and OFC. A, Illustrations of recording positions of all principal units included in analysis from posterior striatum (left), anterior striatum (center), and OFC (right). Dotted red lines indicate boundaries used to separate units recorded in dorsal and ventral striatum (center) or those recorded in lateral and medial OFC (right). Scale bar, 1 mm. AP positions are distance from bregma. Section diagrams were adapted from Franklin and Paxinos (2008). B, Comparison of elapsed time decoding performance between models trained on recordings from OFC and anterior striatal neurons showed that anterior striatum performs better than OFC (p = 0.0083, unpaired t test). C, Recordings in the anterior striatum were grouped based on whether they included predominantly dorsal or ventrally recorded neurons (n ≥ 35 cells), with one recording being distributed into both subregions. Dorsal and ventral populations performed as well as populations containing 35 cells drawn uniformly at random from both areas (F(2,16) = 0.02, p = 0.98, one-way ANOVA). D, All recordings in the OFC were bisected into lateral and medial populations. Lateral and medial populations performed and populations containing 29 cells drawn uniformly at random from both areas (F(2,10) = 0.48, p = 0.64, one-way repeated-measures ANOVA). Error bars indicate SEM.
Medial and lateral OFC equally encode elapsed time
Our cortical recordings were mostly positioned within the OFC (Fig. 5A, right). However, this area is composed of several different anatomical subdivisions, raising the possibility that certain subregions encode time better than others. We therefore examined whether medial or lateral fields within our OFC recordings had a differential neural representation of time (we used datasets with at least 29 simultaneously recorded pyramidal cells). We found that models trained on medially or laterally positioned OFC units were just as effective at representing time as models using units taken randomly from either medial or lateral areas (F(2,10) = 0.48, p = 0.64, one-way, repeated-measures ANOVA; Fig. 5D). These findings suggest that the encoding of time via population dynamics is not localized to specific regions of the OFC.
Lick-related movement does not explain the striatum's improved encoding of time
Timing and movement are intimately related. Indeed, in the current study, task licking should have been driven in part by an internal representation of time, but it is possible that some of the code for time that we observed might reflect neurons encoding motor behaviors directly. If the encoding or planning of motor activity were the primary basis for the observed code for time during reward-anticipatory licking after CS+ cues, then we would predict that any licking episode would also encode time. We therefore examined whether population coding for time transferred to false alarm CS− trials in which animals errantly licked after CS− odor presentations. Licking onset time was conserved between CS+ trials and false alarm CS− trials (mean CS+ lick onset time = 1.8 s, SD = 0.25 s; mean CS− lick onset time = 1.8 s, SD = 0.30 s; p = 0.80, paired t test; Fig. 6A). To quantify the extent to which time-related coding could be detected during false alarm trials, we trained the SVM decoder on correct CS+ trials and tested the model on the 0 to 2.5 s interval during false alarm trials. The performance of these models tested on false alarm trials was severely attenuated compared with their performance when tested on correct CS+ trials (Fig. 6B). A two-way, mixed-model ANOVA revealed a significant effect of trial type used for testing (F(1,13) = 33.0, p < 0.0001). The ANOVA did not show a significant interaction between brain region and trial type, demonstrating that both striatum and OFC saw an equal decrement in model performance when tested on CS− trials (F(1,13) = 0.5, p = 0.48). These results suggest that temporal encoding is not only sensitive to licking, but also to other task variables such as the cue context.
Population coding of elapsed time is specific to CS+ trials and is not fully explained by licking behavior. A, Mice showed similar licking onset times during CS+ trials and CS− false alarm trials (p = 0.80, paired t test). B, Comparison of performance in decoding elapsed time for SVM models trained on correct CS+ trials and tested on either correct CS+ trials or on CS− false alarm trials (55 units per region, n = 9 striatal recordings and 6 OFC recordings). There was a significant effect of trial type (F(1,13) = 33.0, p < 0.0001, two-way, mixed-model ANOVA) and a significant effect of brain region (F(1,13) = 18.3, p = 0.00091), with no significant interaction (F(1,13) = 0.5, p = 0.48). C, Mean fraction of recorded principal cell populations showing significant activity modulation by licking in each brain region (p = 0.044, unpaired t test). D, Example licking-modulated principal cells recorded in each region (left, striatal MSN; right, OFC pyramidal). Shaded blue rectangle represents odor presentation time. Black tick marks indicate individual spikes, red ticks denote lick onset times, and blue dotted line shows reward delivery time. Trials are sorted by descending latency to first lick. E, Comparison of elapsed time decoding performance between models generated using all cells or all non-lick-modulated cells. Performance showed a significant decrease with the exclusion of lick-modulated cells (F(1,13) = 17.2, p = 0.0011, two-way, mixed-model ANOVA). The striatum maintained an improved code for time over the OFC after excluding lick-modulated cells (F(1,13) = 7.4, p = 0.017). We did not observe a significant interaction between region and population (F(1,13) = 0.9, p = 0.35). Error bars indicate SEM.
To further examine the contribution of lick-related movement to our data, we identified principal cells that were positively correlated with lick rate. Although we identified lick-rate-modulated cells in both areas, the striatum contained a significantly greater proportion of these cells than the OFC (mean striatal lick-modulated fraction = 0.35, SD = 0.127; mean OFC lick-modulated fraction = 0.226, SD = 0.049, p = 0.044, unpaired t test; Fig. 6C). Figure 6D depicts two examples of lick-rate-modulated neurons from the striatum (left) and the OFC (right). We retrained and tested the decoder after excluding these cells from the population. We found that removing lick-rate-modulated cells reduced decoder performance below what would be expected after removing the same number of randomly selected cells (F(1,13) = 17.2, p = 0.0011, two-way, mixed-model ANOVA; Fig. 6E). However, crucially, the decoder still performed significantly above chance levels, demonstrating that a code for time was still present without lick-rate-modulated cells. In addition, we found that the striatum still performed better at representing time over the OFC despite the exclusion of lick-modulated cells (F(1,13) = 7.4, p = 0.017, two-way, mixed-model ANOVA). The ANOVA did not reveal a significant interaction between brain region and the type of population used in analysis (F(1,13) = 0.936, p = 0.35). Together, the results in Figure 6 show that, although movement does indeed contribute to the observed code for elapsed time in both the striatum and OFC, it is not sufficient to fully explain the neural representation of time in these areas. Furthermore, we demonstrated that our main finding that striatal ensembles outperform OFC ensembles in terms of temporal coding is robust even after controlling for lick-rate-modulated cells.
Striatal population codes for elapsed time covaries with lick onset time
To this point, our decoding analysis had been performed on all correct CS+ trials regardless of the animal's actual lick onset time. However, because we found that lick-related movement partially contributed to the neural code for time, this implies that the neural code may vary on a trial-to-trial basis depending on the precise timing of lick onset. If the population dynamics are sensitive to lick onset time, then a prediction is that the encoding trajectories are respectively traversed faster (slower) when an animal licks earlier (later) than the mean. To test this prediction, we took advantage of the trial-to-trial variability in the time at which animals initiated licking during CS+ trials (Fig. 1C). We determined whether population dynamics in the striatum and OFC reflected this variable lick onset time. For each animal, we divided trials into three evenly sized groups representing early (first tercile), intermediate (second tercile), and late (third tercile) lick onset time trials (Fig. 7A). We then trained SVM models on trials in the first or third terciles and tested each separately on trials in the remaining terciles. In the striatal population, we found that, when testing the first tercile's model versus the third tercile's model on trials from the second tercile, the evaluations showed opposing classification error biases (p = 0.00034, paired t test; Fig. 7B). In other words, the model trained on the first tercile consistently classified time bins in the second tercile as having occurred later than they had. Conversely, the model trained on the third tercile consistently classified time bins in the second tercile as having occurred earlier than they actually had. Furthermore, when testing the first tercile's model on the third tercile's trials or testing the third tercile's model on the first tercile's trials, these evaluations also showed opposing classification error biases (p = 0.002, paired t test). Altogether, these results show that internal representation of time in the striatum appears to covary with the timing of lick onset, consistent with earlier work suggesting that the latency of the motor response was driven by the neural code for time (Gouvêa et al., 2015).
Striatal population coding of elapse time shows higher sensitivity to lick onset variability than OFC. A, Schematic illustrating the division of correct CS+ trials into three sets based on terciles of the lick onset distribution. B, Mean prediction biases of SVM decoders trained to predict elapsed time from striatal population data recorded in first tercile trials (orange) and tested on second and third tercile trials. Green bars show decoder biases when trained on third tercile trials and tested on data from first and second tercile trials. Training on first and third tercile trials and testing on second tercile trials produces opposing biases (p = 0.00034, paired t test), as does training on first tercile trials and testing on third tercile trials compared with training on third tercile trials and testing on first tercile trials (p = 0.002, paired t test). C, Mean prediction biases of SVM decoders trained to predict elapsed time from OFC data under similar conditions as in B. No significant difference in biases were observed when training on first and third tercile trials and testing on second tercile trials (p = 0.22, paired t test) or when training on first tercile trials and testing on third tercile trials compared with training on third tercile trials and testing on first tercile trials (p = 0.06, paired t test). D, Illustration of temporal alignment procedure on one striatal recording (88 cells). Distance matrix represents the Euclidean distance between all pairs of population activity patterns in the trial-averaged trajectories for the first and third tercile trials. Red line traces the minimum distance path along the distance matrix between the beginning and the end of the mean first tercile trajectory. A deviation (red arrows) of this path from the diagonal (dashed yellow line) measures the temporal warping of the mean third tercile trajectory relative to the mean first tercile trajectory. The upward shift observed here indicates that the mean third tercile trajectory is consistently slower. E, Mean temporal warping of striatal (black) and orbitofrontal (blue) third tercile trajectories relative to their respective first tercile trajectories. Error bars indicate SEM.
In contrast to the striatum, in the OFC, we did not find any significant effects of training classifiers on the first or third terciles and testing those models on the second tercile's trials (p = 0.22, paired t test; Fig. 7C). Testing first or third tercile classifiers on the third or first terciles' trials, respectively, also did not result in biased classification error deviations (p = 0.06, paired t test), although there was a trend. Therefore, in contrast to the striatum, the temporal code in the OFC may not covary as effectively with movement onset time.
It was possible that the decoded biases quantified above did not fully establish the extent of the underlying relationship between the internal representation of time and lick time due to potential artifacts imposed by binning and smaller training datasets after grouping by terciles. To better determine the extent of temporal covariation between licking and neural dynamics, we compared population trajectories averaged over the trials in the first tercile with population trajectories averaged over trials in the third tercile. A temporal alignment procedure applied to the two trial-averaged trajectories (see Materials and Methods) revealed that, whereas the two trajectories remained close to each other over the course of the trial interval, they were not uniformly aligned in time (Fig. 7D). Instead, the third tercile trajectory consistently lagged behind the first tercile trajectory, illustrated in Figure 7D as an upward shift of the minimum distance curve between the two trajectories away from the diagonal line. The magnitude of this shift is a measure of the temporal warping, or speed of progression, of the third tercile trajectory with respect to the first tercile trajectory. In the striatum, temporal warping emerges very early on in the trial relative to lick onset timing (Fig. 7E), which suggests that the striatal activity encoding an internal representation of time undergoes “subjective” fluctuations that may drive trial-to-trial variability in lick onset. In contrast, warping was less prevalent, particularly near the beginning of the trial, in the OFC. Together, these results suggest that the internal representation of time as encoded in the striatal dynamics are comodulated by the elapsed time and the lick onset time and these effects are less evident in the OFC.
Striatal ensembles predict movement onset time
Because striatal ensemble dynamics possessed a better code for time, we hypothesized that the lick onset time could also be predicted with better accuracy from patterns of striatal activity than OFC activity. Using the ensemble firing rate pattern in each 100 ms time bin of a trial, a binary SVM classifier was trained to discriminate the population activity in the first time bin when an animal licked (i.e., the lick onset bin) from the activity in all other time bins (Fig. 8A). SVM output for the population activity in a given bin is represented by a single readout unit with a value that captures the propensity of lick onset occurring in that bin. To establish how well network activity predicted lick onset times, we used a Monte Carlo cross-validation method to test trained SVM classifiers on population activity patterns in novel trials (Fig. 8B). The classifier generates one readout value for the activity pattern from each bin in a trial and the predicted lick onset bin for the trial is chosen as the one with the maximal readout value. Figure 8C illustrates the readout value distributions decoded from the striatal dataset of an animal and its observed lick onset bins (red ticks) for all correct CS+ trials. To quantify the classification performance, we measured the RMSE of the predicted lick onset times across all correct CS+ trials as generated by the Monte Carlo cross-validation approach (Fig. 8D). SVM models trained on the striatal network datasets (observed mean RMSE = 4.07, SD = 1.90; bin-shuffled mean RMSE = 15.32, SD = 1.48, n = 9, p < 0.0001, paired t test; Fig. 8E) and the OFC network datasets (observed mean RMSE = 6.50, SD = 1.96; bin-shuffled mean RMSE = 14.71, SD = 1.12, p = 0.0002, n = 6, paired t test; Fig. 8E) performed well above chance levels in predicting lick onset time. However, consistent with our hypothesis, the SVM models trained on striatal activity outperformed those trained on OFC activity in predicting lick onset times during the task (p = 0.032, unpaired t test; Fig. 8E).
Striatal networks outperform OFC networks at predicting lick onset time. A, Illustration of lick onset time prediction analysis. Raster plots show the same MSN population's activity during different correct CS+ trials. Top schematic shows odor on time (blue rectangle), reward delivery (red triangle), and actual lick times (red/black lines) that correspond to the recorded raster plots. Each correctly performed CS+ trial has a lick onset time indicated by a red line. As in the elapsed time prediction analysis, in each trial, spiking activity of each unit was transformed into corresponding firing rate estimates (data not shown) and the firing rates of simultaneously recorded units were binned (100 ms time bins) to construct population firing patterns for the trial. In each trial, the bin during which the first lick occurred is labeled as its lick onset bin (violet shading). A binary SVM classifier, represented here by a readout unit, was trained to distinguish between lick onset bins and non-lick onset bins (green shading). B, The model is tested using a Monte Carlo cross-validation approach. Population activity patterns for all time bins in a trial are presented to the classifier, which predicts the lick onset bin for the trial as the time bin with the maximal readout value. C, Heat plot showing normalized trial-averaged readout values generated by the SVM trained and tested on striatal network activity of one mouse. Trials are sorted by decreasing latency to lick onset time, indicated by a red tick mark. D, 2D density plot showing the joint distribution of actual lick onset times and those predicted by the SVM from striatal network activity, for one mouse. Prediction performance is measured as the RMSE. Lick onset bin classification was repeated 30 times for each trial (see Materials and Methods). Actual and predicted lick onset bins were jittered (Gaussian noise with 0 mean, 0.3 SD) to separate overlapping points. E, Comparison of mean predicted lick onset bin RMSEs across all striatal and OFC recordings (55 units per region, n = 9 striatal recordings and 6 OFC recordings) showed that models trained on striatal network data performed significantly better (p = 0.032, unpaired t test). Bin-shuffled models based on striatal recordings performed significantly worse than corresponding nonshuffled models (p < 0.0001, paired t test). Bin-shuffled models based on OFC recordings also performed worse than corresponding nonshuffled models (p = 0.0002, paired t test). Error bars indicate SEM.
Simultaneous OFC and striatal recordings exhibit a superior code for elapsed time in the striatum
The above analyses suggest that the network dynamics of the striatum constitute a better “clock” than the dynamics of the OFC. However, it is possible that these observations are partially influenced by differences in neural coding performance across animals. In a subset of our recordings (n = 4), we were able to measure at least 55 OFC pyramidal cells and 55 striatal MSNs simultaneously within the same animal and session (Shobe et al., 2015). Therefore, we investigated whether the observation that striatal dynamics contain a better code for elapsed time than the OFC was supported in these simultaneous dual region recordings. One particular advantage of this within-animal comparison is that the SVM models are trained and tested with network data in two brain regions that were recorded using identical behavioral conditions and trials. Therefore, the networks share the same stimulus inputs, interval durations, and lick onset times. A cursory comparison of the cross-temporal classification matrices for simultaneously recorded brain regions in a single animal indicated that its striatal population encoded elapsed time more robustly than its OFC population (Fig. 9A). An accuracy comparison of the elapsed time decoded from population activity in the two brain regions across all simultaneously recorded network activity datasets reasserted that the striatal networks' temporal encoding efficacy was consistently better (mean striatal correlation coefficient = 0.90, SD = 0.041; mean OFC correlation coefficient = 0.71, SD = 0.13, n = 4, p = 0.013, paired t test; Fig. 9B). We also measured whether this effect was consistent at different sizes of the decoded neural population. A two-way, repeated-measures ANOVA between brain region and the decoded population size showed a significant effect of brain region (F(1,3) = 58.1, p = 0.0047) and population size (F(5,15) = 73.4, p < 0.0001; Fig. 9C).
Simultaneous multiregion recordings indicate that striatum encodes elapsed time better than OFC. A, Left, Average cross-temporal classification matrix showing mean performance of the elapsed time classifier across all correct CS+ trials for one striatal recording that occurred in parallel with a OFC recording in the same mouse. The classification matrix for the corresponding OFC recording is shown at right. B, Mean correlation coefficient across simultaneous striatal and OFC recordings (55 units per region, n = 4) for each brain region. SVM classification of population activity was repeated 30 times (see Materials and Methods). SVM models trained on striatal activity performed better than when trained on OFC activity patterns (p = 0.013, paired t test). C, Performance comparison of SVM models trained and tested on striatal and OFC network activity from simultaneous recordings as a function of number of units. There was a significant effect of brain region (F(1,3) = 58.1, p = 0.0047) and number of units (F(5,15) = 73.4, p < 0.0001, two-way repeated-measures ANOVA). Error bars indicate SEM.
Simultaneous OFC and striatal recordings exhibit a superior lick onset time prediction in the striatum
We also explored whether lick onset time prediction was significantly better using striatal population activity within the simultaneously recorded striatal and OFC datasets. We again observed that SVM models trained to identify population activity encoding lick onset time appeared to be more precise when decoding from striatal population activity than from OFC population activity (Fig. 10A). When comparing lick onset prediction performance between SVM models trained on simultaneously recorded networks as a function of the brain region and the size of the decoded population, we found that the striatal networks encode lick onset time with a significantly higher efficacy than OFC networks (F(1,3) = 18.9, p = 0.022, two-way, repeated-measures ANOVA). There was also a highly significant effect of population size on lick onset prediction performance across the two brain regions (F(5,15) = 178.4, p < 0.0001, two-way, repeated-measures ANOVA; Fig. 10B).
Simultaneous multiregion recordings show distinct prelick dynamics across striatal and OFC networks. A, Heat plots showing normalized readout values generated by SVM models trained to detect lick onset times. Heat plots reflect trial-averaged readout values of SVM models trained and tested on striatal (left) and OFC (right) network activity from simultaneous recordings from the same mouse (55 units per region, n = 4). Trials are sorted by decreasing latency to lick onset time, indicated by a red tick mark. B, Mean performance of lick onset bin prediction as a function of number of units included in training and testing the SVM models for each simultaneously recorded brain region (55 units per region, n = 4). A two-way, repeated-measures ANOVA revealed a significant effect of number of units (F(5,15) = 178.4, p < 0.0001) and a significant effect of brain region (F(1,3) = 18.9, p = 0.022). The ANOVA excluded the “all units” column because it contained inconsistent numbers of cells between simultaneously recorded regions. C, Heat plots showing normalized readout values generated by SVM models trained to detect time bins occurring 500 ms before actual lick onset times. Heat plots reflect trial-averaged readout values of SVM models trained and tested on striatal (left) and OFC (right) network activity from simultaneous recordings from the same mouse. Trials are sorted by decreasing latency to actual lick onset time, indicated by a red tick mark. Magenta tick marks indicate 500 ms before lick onset (55 units per region, n = 4). D, Mean RMSE values across all simultaneous striatal and OFC recordings (55 units per region, n = 4) quantifying performance of SVM models trained and tested to predict time bins that occurred in advance of actual lick onset times. A two-way, repeated-measures ANOVA revealed a significant effect of time bin (F(15,45) = 8.8, p < 0.0001) and brain region (F(1,3) = 16.0, p = 0.028). Error bars indicate SEM.
Last, we investigated how far in advance the neural activity could predict lick onset time. This allowed us to further investigate whether lick onset prediction was not simply a product of neural activity directly driving motor responses. For this analysis, we trained separate SVM classifiers on increasingly earlier target bins, moving the target bin backward in time with respect to the actual lick onset bin (see Materials and Methods). In each brain region, a comparison of the readout value distributions when the classifier is trained on the actual lick onset bin (Fig. 10A) versus on the bin occurring 500 ms ahead of the actual lick onset bin (Fig. 10C) indicates that either time bin can be decoded with similar reliability (Fig. 10D). A two-way, repeated-measures ANOVA revealed a significant effect of brain region (F(1,3) = 16.0, p = 0.028). SVM models trained to decode prelick onset time bins from OFC population activity were consistently worse than those based on striatal population activity. There was also a significant effect of time (F(15,45) = 8.8, p < 0.0001) in which increasing the “look-ahead” time resulted in a decay in classification performance at approximately 1 s. This effect was more pronounced in the striatum, which maintained a relatively constant RMSE until ∼1 s before actual lick onset. Conversely, the OFC showed higher error at actual lick onset compared with the striatum. Interestingly, the OFC showed a slight improvement of classification further away in time from the actual lick onset. These results suggest that striatal and OFC dynamics are causally related to behavior and that the activity patterns preceding the lick by up to 1 s encode when the animal will lick. Again, the results show that striatal dynamics provide more information about lick onset time.
An important question related to both the nature of the time code and the mechanisms underlying the neural dynamics in the striatum and OFC relates to the “complexity” of the dynamics. One way to measure the complexity of patterns of neural activity is through its effective dimensionality, a PCA-based measure calculated as the number of PCs required to explain 95% of the variance in the trial-averaged population dynamics (Rajan et al., 2016; see Materials and Methods). The calculations show that the OFC dynamics (mean = 15.75, SD = 2.06) are of a significantly higher dimension than MSN dynamics (mean = 11.25, SD = 1.5, p < 0.0001, paired t test) and suggest that the OFC may encode other variables that are not immediately relevant to the task.
Discussion
This study took advantage of the high single-unit recording throughput of silicon microprobes (Shobe et al., 2015) to examine the neural coding properties of large prefrontal and striatal ensembles during a conditioning paradigm in which mice learned to anticipate the timing of reward. We used a machine-learning algorithm to quantify the ability of dynamically changing patterns of network activity to encode time at the single trial level within individual animals. Using activity patterns in these two brain regions, we decoded time elapsed from the onset of a reward-predictive cue. We also could predict the time of anticipatory licking onset. Our results show that the striatum consistently outperformed the OFC in terms of the ability to encode time. This superior time representation of striatal ensembles was confirmed via simultaneous recordings in the OFC and striatum.
The large scale of the recordings (at least 55 simultaneously measured units per animal) enabled quantitative comparisons of decoding performance without the need for pooling units across sessions or subjects. Although a few studies have recorded simultaneously in the cortex and striatum (Matell et al., 2003; Jin et al., 2009), to our knowledge, this is the first study to compare quantitatively the quality of the time code between the striatum and one area of the cortex, the OFC, while recording simultaneously from large populations of neurons. Therefore, we anticipate that the approaches outlined in this work offer new opportunities for understanding the coding properties of neural ensembles across multiple brain areas during behavior (Brown et al., 2004; Buzsáki, 2004).
Computational models (Medina et al., 2000; Buonomano and Laje, 2010; Laje and Buonomano, 2013) and recent experimental work (Stokes et al., 2013; Crowe et al., 2014; Carnevale et al., 2015) have suggested that motor timing may be encoded in dynamic patterns of neural activity: a “population clock”. In vivo, population clocks in the form of either simple sequential patterns of activity or complex high-dimensional patterns have now been observed in many different brain areas, including the striatum (Matell et al., 2003; Chiba et al., 2008; Jin et al., 2009; Gouvêa et al., 2015; Mello et al., 2015; Bakhurin et al., 2016) and multiple areas of the prefrontal cortex (Dietrich and Allen, 1998; Fuster, 2001; Brody et al., 2003; Oshio et al., 2008; Genovesio et al., 2009; Jin et al., 2009; Merchant et al., 2011; Kim et al., 2013; Xu et al., 2014; Carnevale et al., 2015). Our data are consistent with these previous reports of dynamic time-varying coding properties of neural ensembles and support population clock models for the coding of time. This study provides evidence that different brain regions, here the striatum and the OFC, may be part of a distributed but regionally specialized network for encoding time.
The brain's code for time and the timing of movements are highly interdependent; animals use a timing mechanism to determine when to generate actions, but those actions in turn cause changes in brain activity, potentially influencing any observable code for time. We have addressed such interdependences in our data and show that, whereas licking-related activity can contribute to the population clock, there is a robust code in a number of control analyses aimed at removing potential motor influences. Because we trained mice to learn to time a single interval, we exploited the fact that licking behavior demonstrated variable onset timing, similar to the timing variability of lever pressing during fixed-interval tasks (Matell et al., 2003; Mello et al., 2015). Our study benefited from this natural variability in that we could demonstrate that population codes reflected early or late onset times within the single interval. Interestingly, we found evidence that licking behavior that occurred particularly late or early strongly covaries with the speed at which temporal codes evolved along the entire duration of the trial. This was particularly apparent in the striatum, which supported results in earlier work on population coding in the dorsal striatum (Gouvêa et al., 2015). Overall, these analyses support the notion of high-dimensional multiplexed representations within the striatum and OFC (Rigotti et al., 2013; Fusi et al., 2016).
Our study focused on the OFC region of the prefrontal cortex, an area that has not received extensive attention with respect to encoding of time. Parts of the prefrontal cortex, including the dorsolateral, medial, and premotor areas, have been shown previously to encode time (Onoe et al., 2001; Kim et al., 2009; Merchant et al., 2013b; Crowe et al., 2014; Merchant et al., 2015). It is therefore possible that these other areas exhibit a better code for time than both the regions that we investigated in the OFC and the striatum. Despite not having been studied extensively in the context of coding time per se, neurons in the OFC has been shown previously to be sensitive to time during reward expectation and other related behaviors such as temporal discounting, in which rewards received earlier in time are preferred to those associated with a greater delay period (Roesch et al., 2006; Moorman and Aston-Jones, 2014, but see Jo et al., 2013). In this study, we now show that the OFC is capable of encoding time using population clocks. We also add to a growing literature that the OFC contains movement-related coding (Feierstein et al., 2006; Furuyashiki et al., 2008; Simon et al., 2015). A further implication of our data is that the OFC may be better suited for representing higher-dimensional information about the behavioral task. This was reflected in the effective dimensionality analysis and is consistent with the role of the OFC in complex cognitive processing (Fuster, 2001), which has been hypothesized to involve recurrent circuit activity (Mante et al., 2013; Rigotti et al., 2013).
A major finding of this study is that ensembles in the striatum outperform the OFC in terms of time encoding, even after controlling for motor effects. A fundamental question thus pertains to how these differences in coding arise. Is time encoding generated within the striatum, through a corticostriatal loop (Merchant et al., 2013a), or is the striatum instantiating a readout of the dynamics generated within neocortical areas? These questions cannot yet be answered, but we suggest that the most consistent interpretation given the superior performance of the striatum observed here is that the striatum is serving as a readout of the dynamics generated in the cortex. The striatum is innervated by a plethora of areas, including prefrontal regions and motor areas (McGeorge and Faull, 1989; Voorn et al., 2004; Hintiryan et al., 2016). As a result of its functionally diverse inputs, the striatum as a whole integrates information that spans the cue–reward delay period, including sensory stimuli, reward prediction, and action initiation (Tremblay et al., 1998; Jog et al., 1999; Nicola et al., 2004; Roitman et al., 2005; Jin and Costa, 2010; Rueda-Orozco and Robbe, 2015). Our results therefore suggest an important role for sensorimotor integration in the striatum (Reig and Silberberg, 2014), which may lead to a more refined representation of time in this structure than what was found in an upstream frontal cortical area.
We subdivided our dataset into dorsal and ventral striatal ensembles and showed that each subregion separately performed as well in encoding time as when both subregions were pooled together. Most studies of timing in the striatum have focused on the dorsal striatum (Chiba et al., 2008; Jin et al., 2009; Bartolo et al., 2014; Mello et al., 2015), with the exception of a study that investigated how dopaminergic signaling in dorsal and ventral areas contributed to timing behavior (Meck, 2006). Our results show that the ventral striatum encodes time as effectively as the dorsal striatum. We do not suggest that the ventral and dorsal areas are encoding the same kinds of information (Bakhurin et al., 2016), nor would they need to do so to represent time in their distinct patterns of dynamic activity. It is also possible that the synchronization of a temporal code could be attributed to local striatal microcircuitry (Bakhurin et al., 2016; Barbera et al., 2016) or to basal ganglia feedback loops (Haber et al., 2000).
Several studies using fMRI approaches in humans have shown that dorsal areas of the striatum, including both the caudate and the putamen, are involved in interval timing tasks (Harrington et al., 2004; Wiener et al., 2010; Coull et al., 2011). It is important to point out that, in contrast to the current results, the ventral striatal subregions were not often modulated significantly in these studies. This could be explained by differences between the timing tasks performed by subjects across these studies. For example, successful performance in our task relies on animals anticipating the delivery of a reward, thus explaining why the ventral striatum, a region commonly implicated in reward processing (Day and Carelli, 2007), may be recruited. Tasks developed for humans may allow for more explicit control of the brain's capacity to use time, but also may rely less on simple Pavlovian associations that we used to train mice. In addition, fundamental differences between single-unit spike and BOLD signal measurements could make direct comparison of our data with fMRI experiments challenging.
The neuronal architecture of striatum differs dramatically from that of the neocortex: striatal circuits are characterized by recurrent inhibition (Tepper et al., 2004), whereas neocortical circuits contain recurrent excitation. Theoretical studies have established that, in contrast to inhibitory circuits, excitatory recurrent circuits are ideally suited to generate self-sustaining time-varying patterns of activity (Jaeger and Haas, 2004; Sussillo and Abbott, 2009; Laje and Buonomano, 2013), although it is possible for such patterns to emerge from circuits that exhibit recurrent (feedback) inhibition (Mauk and Donegan, 1997; Medina et al., 2000). Because of the lateral inhibition interactions among MSNs (Taverna et al., 2008) and the influence of local interneurons (Tepper et al., 2010), the striatal microcircuit may be well suited to refine those signals into an improved time code through a temporal, winners-take-all mechanism by ensuring that the time-varying patterns of activity within cortical areas only activate a subpopulation of MSN cells at a time (Humphries et al., 2009; Carrillo-Reid et al., 2011; Ponzi and Wickens, 2012).
Footnotes
This work was supported by the National Institutes of Health (Ruth Kirschstein National Research Service Award NIH T32-NS058280 to K.I.B.; Grants MH060163 and NS100050 to D.V.B.; and Grants DA005010, DA034178, and NS100050 to S.C.M.), and the National Science Foundation (Grant IIS-1420897 to D.V.B. and Grant CBET 1263785 to S.C.M.). S.C.M. acknowledges support from a McKnight Technical Innovations in Neuroscience Award.
The authors declare no competing financial interests.
- Correspondence should be addressed to Sotiris C. Masmanidis, Department of Neurobiology, University of California, 650 Charles E. Young Drive South, Los Angeles, CA 90095. smasmanidis{at}ucla.edu
