Introduction

Time estimation is an essential function in the brain^1–3, since many crucial cognitive functions implicitly require a record of time, such as motor control^4,5 or memory^6,7. The neural substrate for time estimation in the brain has been studied for decades, and several modes of time coding emerged from these studies: 1) Time can be represented by the gradual change of activity levels in certain cells (ramping)^8–10; 2) Individual cells show selective activation at a specific time point (sequential)^11–15; and 3) Population coding that showed complex patterns but stable dynamics in latent space^16–18. Interestingly, by training the animal to learn different lengths of waiting periods, several groups found that these “time codes” exhibit scaling properties so that the number of cells for coding different lengths of time remains constant, while the activity can be compressed or stretched according to the duration of target time^11,19–21. While these findings provide strong evidence for a neural mechanism of time coding in the brain, true causal evidence at single-cell resolution remains beyond reach due to technical limitations. Although inhibiting certain brain regions (such as medial prefrontal cortex, mPFC,²²) led to disruption in the performance of the timing task, it is difficult to attribute the effect specifically to the ramping or sequential activity patterns seen in those regions as other processes may be involved.

Lacking direct experimental evidence, one potential way of testing the causal involvement of “time codes” in time estimation function is to examine their correlation at a finer resolution. However, a limitation in the experimental protocols of previous studies is that the animal learns a fixed time target, so that the scaling phenomenon is observed at the group level. Thus, there is a lack of evidence regarding whether the scaling happens rapidly at single-trial level to support time estimation.

In this study, we utilized a novel timing task in rat that allowed the animal to control the waiting period on its own will, so that we can observe the correlation of “time code” scaling with behavioral waiting responses at individual trial resolution. We found robust sequential activities in the mPFC when the rats performed the task. To the best of our knowledge, we provide the first piece of evidence that the scaling effect is dynamic to account for the variation of waiting periods at individual trial level. And the rats were capable of subjectively adjusting the scaling factor for accurate time estimation. Intriguingly, we found that cells coding for the start or end of the waiting period undergo cross-session shifts, while the sequential time code remains stable over weeks. This surprising stability of the time code suggests an underlying mechanism that is different from apparently similar sequential activities seen in place coding or time coding in hippocampus CA1^23,24, which showed dynamic shifts across sessions. Altogether, our study provides strong evidence for the online utilization of sequential time code in rats mPFC during a timing task. The unique rapid scaling and cell-identity stability of these sequential time code suggest a designated cell population for coding time in this region.

Results

Calcium imaging in mPFC during rats perform the timing task

We trained rats to perform a modified version of the timing task used in previous studies²¹. In order to get the water rewards, water-deprived rats must poke their nose in a designated hole and maintain position (Figure 1A). The rat was free to start and end the nose poke at its own will, but only when the duration of the nose poke was above the minimal threshold would the rat receive a water reward. The rat did not receive additional punishment for nose poke duration below the minimal threshold, and can start a new trial when it was ready. Importantly, once the nose poke duration was over the minimum threshold, the amount of water reward was proportional to the total length of the nose poke duration. We trained the rat to perform this task in two phases. After the initial shaping for rats to associate nose poke with water rewards, the rats were trained on a short duration phase in which the minimal threshold was set to be 300 ms. When the rats’ correct performance reached above 70%, they advanced into a long duration phase and the threshold was set to be 1500 ms (Figure 1A).

Figure 1.

Rats that underwent the training paradigm described above showed gradual increase in their nose poke durations (Figure 1B). And the median durations in each phase were slightly above the respective threshold in that phase (Figure 1C), indicating that the rats were capable of learning the different minimal threshold. Notably, there was a wide distribution of the nose poking durations in both phases of the task, which allowed us to study the trial-to-trial variation of the timing behavior.

We recorded neuronal activities from 4 different sessions (1 from short phase and 3 from long phase, Figure 1D). To obtain neuronal activity, we injected adeno-associated viruses expressing calcium sensor GCaMP6S in the medial prefrontal cortex (mPFC) of the rat and installed a 1.8 mm diameter GRIN lens after aspiring the cortical tissue on top (Figure 1E). Using a miniaturized microscope, we were able to record calcium activities from pyramidal neurons in the mPFC (Figure 1F), while the rats can freely move and perform the timing task. Our preparation yielded stable imaging data with good signal-to-noise ratio and cell counts (ranging from 50 to 150 cells per session).

Scaling of the mPFC time sequences at single-trial resolution

We started our analysis by examining the neural activities correlated with different features of the nose-poking behavior. To this end, we aligned calcium traces by nose-poking events for each cell, and calculated the Pearson’s correlation between the averaged activity trace to the binarized event trace (start, during and end). To determine the statistical significance of the correlation, we performed random shuffles of the calcium traces and repeated the calculations above, generating a distribution of correlation coefficients for each cell. If the actual correlation coefficient reached beyond 95% confidence interval of this distribution, we then assigned this cell to the group that codes for the specific feature of the nose-poking event. With this method, we identified neurons that showed selective activation at the beginning, the middle, or the end of the nose poke events (Figure 2A), which we designated as “start cell”, “duration cell” and “end cell”, respectively. A small portion of these coding cells (∼10%) showed significant correlation between two features of the task (“start” and “duration”; “end” and “duration”). In later analysis, we excluded these cells from the duration cell group. Notably, we observed gradual decays of GCaMP signal in start cells and sometimes gradual rise in end cells. However, we cannot definitively separate these “ramping-like” activities from the potential artifact due to slow kinetics of calcium sensor. These three types of cells constituted more than half of the observed neurons (Figure 2B). By training a support-vector machine classifier using activities from the start cells and end cells, we were able to predict nose-poking events with above 90% accuracy, indicating that their activities were highly specific to the nose-poking events (Figure 2C, D).

Figure 2.

We then examined the activity patterns of the duration cells. We found that many duration cells showed activation during a certain proportion of the nose-poking events, and the timing of the activation seemed to be modulated by the total length of the events (Figure 2E, Supp. Figure 1A). When we normalized the actual time to the total length of each nose-poking event and transformed the calcium traces accordingly, we found that many duration cells showed selective activation at a fixed point in the normalized time scale (Figure 2E), consistent with the view that their actual activity was scaled by the total event duration ^11,18. The ensemble activities from duration cells tiled across the normalized nose-poking time, showing a sequential activation pattern (Figure 2F). The duration cell activity exhibited robust sequentiality as measured by peak entropy ^25,26, which was not seen in activities from start or end cells (Figure 2G). This activity sequence was stable across all trials within a session, regardless of whether rat reached minimal reward threshold (Figure 2H, Supp. Figure 1B-C). Furthermore, we were able to train a Gaussian process regression model to predict the progress of each trial with high accuracy (Figure 2I-J). While activities from start and end cells can also decode time, this might be due to slow calcium dynamics arising from either ramping activity or GCaMP kinetics as peak-extracted traces showed no decoding power (Figure 2K). On the other hand, peak-extracted activities from duration cells maintained high decoding power, indicating a more reliable sequential activity pattern. Together, these results indicated that distinct cell groups existed in the mPFC to code for different features of the nose-poking timing task. Duration cells showed selective activation towards specific timepoints in the normalized time scale, and their actual trial-by-trial activity was scaled by total event duration.

mPFC sequential time code remains stable over weeks

Previous studies reported that sequential activity patterns in the hippocampus can represent temporally ordered events²⁷ and time points²⁴. Interestingly, these sequences were highly dynamic, and the ordered pattern showed significant session-to-session variations^28,29. Therefore, we then examined whether the sequential time code we observed in the mPFC showed similar instability across different sessions. To achieve this, we aligned imaging data from multiple sessions based on matching unique blood vessel patterns in the field of view (Supp. Figure 2). This allowed us to examine the coding properties in the aligned cells. We were able to see cells that maintained their coding features in different sessions (Figure 3A-C). And particularly for duration cells, we found that cells can represent a fixed timepoint in the duration across sessions, or shifted to another timepoint while still being a duration cell (Figure 3C, Supp Fig. 2).

Figure 3.

We next formally quantified these dynamics. Below 30% of the start cells or end cells continued to maintain the same type of coding in the next session (Figure 3D). When we used the activities from the aligned start and end cells to classify the nose poking state, we found chance level performance for cross-session decoding despite of high within-session accuracy (Figure 3E). These results indicated that cells code for nose poking start and end are highly dynamic and they did not form stable representation over days. In contrast, around 70% of the duration cells kept their identity in the next session. And within these stable duration cells, more than 70% of them coded for a fixed timepoint in the normalized time scale (Figure 3F). Using activities from aligned duration cells, we were able to make GPR models that show good cross-session decoding of the nose poking time (Figure 3G), suggesting a stable code in the duration cells regarding the sequential activity pattern and time. We found largely the same results from three pairs of sessions: two sessions in the long phase with 2-day interval, two sessions in the long phase with 14-day interval and one short phase session and one long phase session with 16-day interval. Thus, the dynamism and stability we observed were not influenced by elapsed time or task structure, but may be a reflection of intrinsic properties of those cells. Given that previous studies in the hippocampus²⁹and cortex³⁰ all showed unstable sequential pattern across days, to the best of our knowledge our data for the first time demonstrated a stable time code in the brain.

Active scaling of the mPFC sequences represents subjective time estimation

Having found that mPFC showed time-associated stable sequential activities in the duration cells, we next investigated whether this activity could serve as a neural substrate for time estimation in rats performing our task. Previous studies showed that cooling of the mPFC in rats disrupted their performance in tasks that require timing ^21,22. However, in that study it was difficult to separate the influence of other cognitive functions affected by cooling the mPFC. Lacking the technology of single-neuron manipulation in freely-moving rats, we argue that further evidence on the necessity of mPFC sequential activities for time estimation could come from the analysis of the trials in that rats made errors. If the sequential time code we observed is causally linked with time estimation function, we should see reflection of behavior errors in the neural activity in some form. We speculated that there might be three types of coding errors in the mPFC: type I, absence of sequential activities that could be due to inattention or disengagement of the rats; type II, disordered sequences that may lead to errors in time estimation; and type III, scaling error with intact sequential code but wrong target time.

We then examined the above error types in detail. We performed a decomposition of the neural activity by partial least square regression of the neural activities in individual trials to normalized time, and found that almost all trials showed similar trajectories, indicating a lack of type I error described above (Figure 4A). In addition, we found that duration cells were more important for this trajectory pattern compared to other cells that we observed (Figure 4B). The constructed trajectories could explain ∼52% of the variance in the normalized time, which was primarily contributed by the activities from duration cells (Figure 4C). When we examined the raw activities from duration cells in individual error trials, we also found only a small fraction showed the absence of sequential activities (Supp. Figure 3). Furthermore, we can find duration cells that showed stable coding timepoints regardless of correct or incorrect trials, and some duration cells with their activity modulated by trial types (Figure 4D), suggesting that precise time estimation from these sequential activities may be affected. We then quantified the time estimation errors between the predicted time from duration cell neural activity and actual time using our previously established GPR models (Figure 4E). Although we observed some fluctuations in the accuracy, in general we did not find significant difference between correct and incorrect trials, and we were able to train a common model that showed accurate time estimation for both trial types (Figure 4F), indicating a lack of type II error.

Figure 4.

To estimate the scaling errors associated with the type III errors, we calculated the cumulative errors between GPR model predicted time based on neural activities and actual trial time. When we plotted these data against actual individual trial durations, we found that most trials were scattered around zero, indicating relatively high accuracy in the time estimation. Interestingly, incorrect trials showed a cluster centered below zero, and correct trials showed one above zero (Figure 4G), indicating that during the majority of the incorrect trials the rats timed too short. Surprisingly, the duration time at which the polarity of the cumulated error reversed was 1500 ms, exactly the minimal time threshold required for the rats to receive reward. This result strongly suggests that the rat was capable of perceiving the minimal time threshold and adjusting trial-by-trial scaling factor of the mPFC sequential activities during the task, and that errors in the online scaling accounted for most of the incorrect trials that we observed.

Discussion

The neural basis of time representation is a fundamental question in neuroscience. As timing is closely intertwined with various cognitive functions, such as short-term memory, decision-making, etc., it is perhaps not surprising that neuron activities correlated with time were found throughout the brain. Previous studies established that sequential or ramping activity can be feasible neural code for time representation. Indeed, time-related sequential or ramping patterns of neuron activity have been reported in prefrontal cortex ^31,32, motor cortex ^33,34, sensory cortex^35–37, striatum^11,38–40, thalamus⁴¹, hippocampus^12,42,43, and entorhinal cortex^30,44. A fascinating feature of time coding is the capability of scaling, in which the same sequential or ramping patterns were compressed or stretched to represent different durations of time^2,45. While these previous findings set forth substantial insights into the neural substrate of time in the brain, two major limitations exist in the field.

One aspect is the lack of direct causal evidence regarding the neuron activity and the perception of time. This would not be an easy task as the causal evidence between similar sequential patterns of place cell firing and spatial perception was not fully demonstrated either, although the phenomenon was described about 50 years ago. Recent advances in the technology indicated a possibility of single-cell optogenetic control for such experiments⁴⁶, yet the low throughput and the difficulty for inhibitory control still greatly limit its application. While lacking direct causal evidence, results from our current study provides a finer correlation between sequential firing and timing function. As previous studies generally trained the animal to learn distinct categories of durations^21,47, the scaling effect of sequential activity was usually described as a group effect. We discovered that scaling is correlated with durations at individual trial level, and that failure in reaching the minimal timing threshold for reward can be primarily attributed to scaling errors, in which the sequential pattern prematurely reached the end. Furthermore, the probability of high scaling error was disproportionally increased when the duration of the trial is close to the timing threshold for reward, suggesting that the animal provides trial-by-trial adjustment of the scaling factor that reflected its estimation of the timing threshold. The real-time scaling of sequential activity, and the enrichment of scaling errors around the timing threshold strongly suggest that the animal actively uses this mechanism, and thus the sequential activity we observed is likely to be the neural code for time estimation. Future advances in single-cell manipulations may provide more definitive experimental evidence.

The second limitation in our understanding for the neural basis of time perception is the lack of biophysical mechanism for the emergence of sequential patterns and the scaling effect. In this regard, several theoretical models have been proposed^18,25,48, yet testing these models in biological systems is currently not possible. Our data demonstrated a unique stability of the duration cells and their sequential activity patterns for coding time. Large proportions of duration cells’ coding type and their sequential order remained unchanged over weeks, even when the target time threshold was changed. In contrast, within the same field of view, cells that code for the start or end phase of the nose poke showed dynamic changes, although we could not separate the start and end-coding cells from time-coding cells with ramping patterns and could not assess the long-term stability of time coding by ramping activities. Previous studies of sequential activity did not report such stability, and activities related to event trajectories are generally considered unstable^28,29. It is possible that this difference reflects different mechanisms underlying sequence emergence in different brain regions. This prolonged stability of duration cells suggests that in the mPFC, participation in the sequential pattern was unlikely a result of input-driven flexible coding, but a reflection of certain intrinsic properties. In other words, our data implied the existence of predetermined factors such as a genetic program for time coding neural population in the mPFC. Such feature could the advantageous so that time coding can be a relatively independent module and can be inert to constant plastic changes happened in the brain. Future studies isolating molecular signatures of the time-coding population could provide a powerful tool in dissecting the neural mechanism of time representation in the brain.

Supplementary Figure 1:

Supplementary Figure 2:

Supplementary Figure 3:

Author contributions

CL.M. and BM.L. designed the study. YT.L., WQ.L. and JW.L. performed the experiments. P.Y., SL.Z. and YT.L. performed data analysis. YT.L. drafted the manuscript. P.Y., CL.M. and BM.L. revised the manuscript.

Competing interest statement

The authors claim no competing interest.

Acknowledgements

This study was supported by Jiangxi Natural Science Foundation 20171ACB20002 to BM.L., and National Natural Science Foundation of China 31960171 to CL.M.. P.Y. was additionally supported by Shanghai Pilot Program for Basic Research – FuDan University 21TQ1400100 (22TQ019), the Lingang Laboratory (grant no. LG-QS-202203-09) and National Natural Science Foundation of China (32371036).

Materials and Methods

Subjects

Our study used Sprague Dawley rats. Equal number of male and female rats were used in the study. The environment temperature was controlled at 25-degree Celsius with 12/12 light-dark cycle. All the rats were raised in single-housing after surgery with ad libitum access to food and water. During behavior training, water access was restricted according to an experimental procedure, in which the body weight loss was monitored and kept under 20%. Animal procedures were conducted in accordance with the Guiding Principles for the Care and Use of Laboratory Animals. All experiments were approved and monitored by the Ethical Committee of Animal Experiments at the Institute of Life Science, Nanchang University.

Apparatus

Behavioral training and testing were conducted in a 40×20×30 cm plastic chamber. One side of the chamber had a semicircular hole (2.5 cm diameter) which was 13 cm above the floor and an infrared detector was placed to detect the nose poking, this hole was called as the operant hole. Under the operant hole is a reward hole where rats can get a certain volume of water, and water delivery is controlled by a solenoid valve. If the detector had detected a nose poking longer than the minimal time threshold from a rat, the solenoid valve will deliver a certain amount of water. We had written scripts to control those processes and record the timestamp of the nose poking.

Virus injection

The rat was anesthetized with isoflurane and placed into a stereotactic apparatus. Rats were injected with carprofen solution (5 mg/kg) and dexamethasone solution (0.2 mg/kg) to minimize inflammation. We then shaved the head area and disinfected 3 times with alternating iodophors and alcohol. An opening on the skull was made with a dental drill. Then we unilaterally microinjected 500 nl of rAAV-CaMKIIa-GCaMP6S-WPRE-PA virus (1.37e⁺¹² vg/ml, BrainVTA) at 80 nl/min into the prelimbic cortex with the following stereotactic coordinates: 3.7 mm anterior to bregma, 0.5 mm lateral to the midline, and -3.2 mm ventral to the skull surface. After the virus injection, we waited 10 minutes to remove the microelectrode from the prelimbic cortex, and then sutured the scalp and disinfected. The rat was placed back in the homecage after waking up from anesthesia.

GRIN lens implantation

Our surgical method is based on the surgical procedure reported in mice⁴⁹. Two weeks after the virus injection, the GRIN lens was implanted above the previous injection site. For the procedure, after disinfecting the head area, we removed the scalp with scissors and cleaned up the skull with hydrogen peroxide and saline, and then inserted three skull screws on the skull. The skull screws were placed in a way that would not block the lens. Then, we used saline to clean up the fragment, and drilled a hole with 1.9 mm diameter on top of the virus injection site. We slowly aspired brain tissue above the prelimbic cortex by connecting a syringe needle with the vacuum pump. The aspiration stopped at 2.5mm below brain surface. We then inserted a GRIN lens (1.8 mm diameter, N.A: 0.54., Edmund) to the target position with the micromanipulators. The placement was checked using a miniscope (Labmaker) and adjusted to maximize the fluorescence signal.

We then removed excess liquid with the vacuum pump, and glued the GRIN lens with skull and skull screw using cyanoacrylate glue. After the glue was fully cured, we then applied black dental cement and covered the top of the lens with Kwik-Sil (WPI) after the dental cement is dry. The rat was placed back to their homecage and amoxicillin was added to their drinking water for 7 days.

After GRIN lens implantation, we waited 4 to 6 weeks to install baseplate. The rat was anesthetized with isoflurane again, and then we used forceps to remove the Kwik-Sil from the top of the lens, cleaned the lens with wet lens paper. We then installed a miniscope with a micro-manipulator holder and found the best view with clear vasculature and cells. At this position, we applied glue around the baseplate. After the glue was set, we continued to cover the larger area with dental cement. When the dental cement is dry, we removed the miniscope. Finally, we covered the plastic cap on the top of the baseplate and fixed it with a screw.

The timing task

To motivate the rats to learn the task, their drinking water was limited during the training period, so that rats could only obtain water through task operation. Specifically, the rats need to probe their nose into the operant hole and stay longer than the minimal threshold time. After rats withdrew their nose from the hole, they can get reward from the reward hole. If the rats performed this action and acquired water reward, we called this as a correct trial. Conversely, an incorrect trial is when the rats started a nose poke but did not reach designated time for reward.

The training was separated into three phases. In the pre-training phase, once rats produced the nose poking action, they could immediately acquire water reward at a small amount (60 μl). The rats moved on to the short phase after they mastered the nose poking action. In this phase, rats would acquire the reward only if they kept the action of nose poking for at least 300 ms. Rats received 50 μl water with 300 ms duration. The rats graduated from this phase when their performance reached above 70% correct. In the long phase, the nose poking duration must reach at least 1500 ms. Rats received 100 μl water with 1500 ms duration. In these two phases, with each additional 500 ms nose poke time in a trial, the reward was increased by 60 μl.

Calcium image recording

We did not anesthetize the rat for installing miniscope. Instead, we hand-held the rat to maintain its position and installed the miniscope and fixed it with screws. We monitored the calcium signal so that if the recording was not stable, the experiment would be interrupted. The calcium imaging of animals was recorded using a UCLA miniscope (V3), and the animal behavior was recorded using a webcam. Data acquisition software recorded the timestamps of the behavior camera and miniscope at the same time in order for subsequent alignment. Video streams were recorded at 20 frames per second. And excitation was set at ∼4.5 mW. In our experiment, we used an MCX adapters to link the coaxial cable to ensure the animals can move freely and avoid the winding of cables.

Extracting calcium activity

All analyses of the calcium imaging data were performed using custom software written in MATLAB. First, the raw imaging data was motion corrected using a previously published the NoRMCorre algorithm⁵⁰. The corrected movies have been processed by median filter and down-sampled prior to extract the calcium traces. Here, to accurately extract the calcium signals of individual cells, we applied a supervised and robust the EXTRACT algorithm^51,52. Extracted units were checked manually to exclude artifacts.

The identification of cell types

To identify cells that encode different periods of the nose-poking event, we calculated the Pearson’s correlation coefficients between calcium traces and binarized event traces. In our task, we focused on three types of events, the start of nose pocking was defined as the 5 frames prior to the detection of the nose-poking event start. The end of nose probing was defined as the 5 frames after the rat began to withdraw from the operant hole. The duration of the nose poke was defined as the interrupted time of the infrared sensor in each trial. In order to calculate the statistics of the correlation coefficients, we generated a null hypothesis for each cell by shuffling the calcium traces with random delays for one hundred times, and calculated the resulted correlation coefficients. A cell was deemed to respond to certain events when the observed coefficient reached above the 95th percentiles of the corresponding coefficient distribution from the shuffled dataset. Those neurons responding specifically the three events were called respectively the start cells, the end cells and the duration cells.

The alignment of calcium imaging cross sessions

We performed cross-session alignment by manually select unique blood-vessel or neuron patterns as fixed control points, and then calculated the transformation between sessions by linear intrapolation. We then applied the transformation to the centroid coordinates of each extracted cell to obtain the aligned cell positions in the target session. By calculating the pair-wise distances between the aligned centroids and that from the extracted cells in that session, we selected the minimal distance for each cell as the most likely candidate matches. And if this distance is below 5 pixels, which is about the radius of a cell, we recorded these two cells as a matched pair. In order to maximize the number of aligned cells, in all cross-session analysis, we used matched paired between two sessions.

Support vector machine model for classification

We trained a linear kernel support vector machine (SVM) model using the activities of the start and end cells to predict the binarized nose poking state for each timepoint. In order to avoid biases caused by the imbalance of the samples, we chose random subsamples of recorded timepoints with equal number of nose-poking and non-nose-poking states. From this subsample, we then selected 80% of the timepoints for model training, and used the remaining 20% to evaluate the performance of the model. For performance measurement, we tabulated the differences between the model’s predicted states from the actual states and calculated percentages of accurate predictions. In order to avoid variations introduced by sampling, we repeated our random sampling process for 100 times for each dataset and recorded the average performance. For generating a null hypothesis, the actual states were randomly shuffled during the training process, and the resulted model was then tested on real testing dataset. For testing cross-session performance, the within-session training and testing procedures were similar as described above, but using only cells that can be matched to the target session as model input. After training, the model was tested on the timepoints of the target session using activities from the activities of the matched neurons, and accuracy calculation was the same as that described above.

Gaussian process regression model

Neuron activities from each nose poking trial were isolated by taking the frames of each nose-poking bouts detected by the sensors and the amplitude of each cell’s activity was normalized into a 0-1 scale within the trial. We then used linear interpolation algorithm (MATLAB “interp1” function) to make each trial into the same number of “normalized” timepoints. In our case we used 30 frames as the normalized trial length. Before model fitting, we selected subsamples of the same number of correct and incorrect trials, and used 80% of the trials for model training and saved 20% for testing. We then fitted a Gaussian process regression (GPR) model with a principal analysis pre-processing with a rational quadratic kernel function to describe the relationship between neuron activities and the progression of trial in normalized timepoints.

Rational Quadratic Kernel:

where σ is the characteristic length scale, α is a positive-valued scale-mixture parameter, r is the Euclidean distance between x_i and x_j.

We feed the testing set data to the fitted model to calculate the predicted normalized time. The performance was defined as the root mean square error (RMSE) between the actual normalized time and the predicted time across all tested trials. Similar to the SVM model, we repeated this process for 100 times to eliminate the sampling variations and shuffled the normalized time in the training set to generate a null hypothesis. It is notable that in this case, the null model output was generally around 15, which is at the middle of our normalized trial length. While this number gave no prediction power, the RMSE (∼9) was lower than a random number set between 1 and 30 (∼12).

Partial least square regression

The VIP score was computed to authenticate the importance of neuronal ensembles on time estimation. VIP score is a measure used in PLS regression to identify the most important variables that predict a response variable. VIP score values range from 0 to 1, where higher values indicate greater importance of the variable in predicting the response. To calculate VIP score using the Partial Least Squares Regression method, several steps are involved. First, the dataset is divided into a training set and a test set. The training set is used to build the PLS regression model, while the test set is used to evaluate its performance. Next, the PLS regression model is trained using the training set. During this step, the model identifies the most relevant variables that explain the variance in the response variable. The VIP score for each variable is then calculated based on its contribution to explaining the variance in the response variable.

The VIP score is importance in predicting the response variable.

XS is an orthonormal matrix including n time points by m components. Each row of XS corresponds to one time point, and each column corresponds to one component. YL is an n-by-one matrix, where n is the number of response variables and one is the number of PLS components. Each row of YL contains coefficients that define a linear combination of PLS components approximating the original response variables. W is n-by-m components matrix of PLS weights.

Statistics analysis

GraphPad Prism version 9.00 was used for statistical analyses. All data are presented as mean ± standard error (SE). Statistical significance was assessed by two-tailed and non-parametric Student’s t-tests. p < 0.05 was considered statistically significant. ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001