Intrinsic timescales across the basal ganglia

Recent studies have shown that temporal stability of the neuronal activity over time can be estimated by the structure of the spike-count autocorrelation of neuronal populations. This estimation, called the intrinsic timescale, has been computed for several cortical areas and can be used to propose a cortical hierarchy reflecting a scale of temporal receptive windows between areas. In this study, we performed an autocorrelation analysis on neuronal populations of three basal ganglia (BG) nuclei, including the striatum and the subthalamic nucleus (STN), the input structures of the BG, and the external globus pallidus (GPe). The analysis was performed during the baseline period of a motivational visuomotor task in which monkeys had to apply different amounts of force to receive different amounts of reward. We found that the striatum and the STN have longer intrinsic timescales than the GPe. Moreover, our results allow for the placement of these subcortical structures within the already-defined scale of cortical temporal receptive windows. Estimates of intrinsic timescales are important in adding further constraints in the development of computational models of the complex dynamics among these nuclei and throughout cortico-BG-thalamo-cortical loops.


Intrinsic timescales of STN, PANs, and GPe.
To assess the spike-count autocorrelation values as a function of time lags, we required a non-zero mean activity for each neuron in each 50 ms bin during the baseline period. In particular, 77/78 neurons in the STN, 103/158 PANs in the striatum, and 92/92 neurons in the GPe fulfilled this first requirement (see "Material and methods"). We further reduced the neuronal sample by requiring each single timescale larger than zero, an R 2 obtained from the fit procedure larger than 50%, and by an accurate eye inspection as described in "Material and methods". After applying this selection criteria, we reduced our neural population to 14/78, 39/158, and 71/92 neurons in STN, striatum, and GPe, respectively. The electrophysiological characteristics of each neuronal population are summarized in Table 1. Figure 1 (left) shows the autocorrelation values as a function of time lags averaged across neurons for each brain structure, with the exponential fit superimposed along with the estimated timescale τ and its interval at 95% confidence level. In particular, the GPe showed a shorter timescale (118 ms [113][114][115][116][117][118][119][120][121][122][123][124]) than both the striatum (250 ms ) and the STN (201 ms ). We noticed that the timescale estimate of STN showed a broad confidence interval, and this is likely the effect of the small size of the neuronal sample (14 neurons). To overcome this issue and to assess whether GPe timescale was significantly smaller than the striatum and the STN timescales, we compared the sample mean of the single timescale distributions (see "Material and methods"). The mean (± standard error of the mean, sem) of the single timescale distributions is 212 ± 46 ms, 223 ± 28 ms, and 131 ± 6 ms for STN, striatum, and GPe, respectively (see Figs. 1B, D, F). We computed the absolute value of the difference between GPe-striatum timescales (91 ± 29 ms), and GPe-STN timescales (80 ± 47 ms), and both differences are not compatible with zero within the error, meaning that GPe showed a significantly shorter timescale than both the striatum and STN (see "Material and methods" for quantification of the difference between the mean of the timescales' distributions based on the general propagation of the error rule). Moreover, the difference between STN-striatum timescales (11 ± 54 ms) is compatible with zero within the error, meaning that the intrinsic timescale of the STN is compatible with the timescale of the striatum.
Comparing the timescales of the BG structures, it is worth noticing that the standard error of the GPe mean timescale is smaller than the other two BG structures. This result might be explained by the lower heterogeneity in the autocorrelation structures of the GPe single neurons compared to that of the neurons of the other two BG structures. We found a significant difference in the heterogeneity of the timescale values of GPe neurons and PANs (two-sample F test for equal variances, p = 9e−18), and GPe and STN neurons (two-sample F-test for equal variances, p = 3e−12). No significant difference in heterogeneity was found between STN neurons and PANs timescale distributions (two-sample F test for equal variances, p = 0.9). In addition, to investigate the degree of heterogeneity in the single timescale distributions for each BG structure, we computed the coefficient of variation (CV) of the timescale distribution using formula (7) in "Material and methods". We found a higher degree of heterogeneity in STN neurons (CV = 15%) and PANs (CV = 12%) τ distributions than in the GPe neurons τ distribution (CV = 7%), in line with the previous individual timescale distribution parameters, i.e. mean and standard error of the mean (sem), that were smaller for GPe neurons than for STN neurons and PANs.

Discussion
To our knowledge, our study is the first to report estimations of intrinsic timescales of neuronal populations at the subcortical level, and reveals timescale differences between the input structures of the BG, the striatum and the STN, and the GPe. Earlier studies have used the autocorrelation function to understand the firing pattern properties of single cells within the BG [13][14][15] and midbrain dopaminergic neurons 16 . Specifically, autocorrelograms of single cells have been used to classify these cells into different subpopulations, for example into GPe neurons 17 or different types of dopaminergic neurons 16 , and to assess firing rate rhythmicity of single neurons from BG nuclei in healthy and diseased conditions 13,14,18,19 . In this study, we used the autocorrelation function to characterize the properties of neuronal populations in the BG and place them in the context of already-known intrinsic timescales throughout the cortex. Table 1. Electrophysiological characteristics of the selected neuronal populations. Each row represents a structure of the basal ganglia (STN for Subthalamic Nucleus, PANs for the Phasically Active Neurons of the Striatum and GPe for the external segment of the globus pallidus). The columns show the number of neurons selected for the analysis, the number of trials during which they were recorded and the mean and standard deviation of their electrophysiological characteristics. The coefficient of variation of the inter-spike interval (ISI) was computed by dividing the standard deviation of the ISI by the mean of the ISI.  The striatum and STN intrinsic timescales reported here place these structures on a comparable level as the timescales assigned to prefrontal areas. Importantly, in both of these BG structures recordings were mainly from their associative parts, known for receiving their major inputs from the prefrontal cortex [22][23][24][25] . In contrast, the GPe exhibited an intrinsic timescale of only 118 ms, which places it at the same intermediate level as the PMd, LIP, and S2. One possibility is that this lower timescale could reflect lower temporal information integration ability of the GPe among BG structures. Our study shows a gradient of receptive time windows within BG circuitry (Fig. 2, red circles). As suggested for areas of the visual system 5 , the differences found among the BG nuclei could reflect a broader window of temporal information storage in the striatum and the STN compared to in the GPe. Indeed, intrinsic timescales could be described as "a temporal counterpart of the spatial hierarchy" 7 . Previous studies 4,26 shed light on the existence of a parallel between the anatomical hierarchy of cortical areas and their intrinsic timescales. Drawing a parallel with the cortex, the interpretation of a functional hierarchy based on timescales described for cortical areas might apply to the BG. This would suggest a greater need for information accumulation within the BG input nuclei rather than in the GPe. The convergence of information from multiple cortical areas could explain the necessity for longer timescales in both input structures because of the need of striatal projection neurons and subthalamic neurons to gate and maintain information from prefrontal neurons 27,28 .
A large body of work in primate neurophysiology has shown that, at the neural level, choosing corresponds to crossing a firing rate threshold in the cortex. Depending on the environment and on the decision being made, the threshold level can be regulated. A computational study 29 implemented a biophysically-based network model of decision thresholds of the cortical-BG-superior colliculus (SC) pathway. After comparing the different nodes of these networks, the authors concluded that the all-or-none activity of SC neurons is triggered by a threshold crossed by cortical neurons that can be optimally tuned by the strength of cortico-striatal synapses. Indeed, through this pathway, the output structures of the BG inhibit cortex and SC activity to preclude inappropriate motor outputs 30 . Some authors 31 have hypothesized that "the BG may convert cortical representations of sensory evidence into evaluative quantities, " allowing generation and adjustment of decisions. They suggest that the BG can modify the decision rules by modifying the decision threshold, but also modify the "value of a developing decision variable." Different BG models hypothesize that the main role of the BG nuclei is to act as a central  www.nature.com/scientificreports/ selection device 32 that examines each action requested based on its urgency and salience 33 and that, with their unique anatomical properties, allows the allocation of motor resources to the appropriate actions. The striatum and the STN have distinct roles in these processes. According to Frank and colleagues, the former has a crucial role in gating sensory input for updating working memory in the prefrontal cortex and then in maintaining it, preventing the influence of distracting information 27 , especially when adaptive gating is necessary for the processing of multiple goal demands. The STN is supposed, particularly during high-conflict decisions, to reduce premature responses and to refine the selection process that takes place via cortico-striatal pathways 28,34,35 . This special role could be conferred by the direct afferences it receives from prefrontal areas 24,25 and the top-down signal it sends to the striatum through the GPe. A computational study suggests that during a countermanding task, after the appearance of an infrequent stop signal, STN acts as a brake to cancel forthcoming actions and keeps accumulating the evidence of the stop signal over time along the decision process 36 . As a consequence, the STN is indirectly involved in higher cognitive functions that require long timescales as for prefrontal areas 26,37 .
On the other hand, models of BG action selection assign another role to GPe neurons. Gurney and colleagues 33,38 have proposed that the GPe, mainly through its massive projections into the STN, automatically limits the activity of BG output structures and allows the network to make a selection. The control that the GPe exerts on other BG nuclei is supported by in vitro electrophysiological studies showing that the firing rate of GPe neurons can be approximated by a linear function of the injected current 39 . In contrast, the firing of the STN neurons fits better with an exponential function of their inputs 33 . This understanding is also consistent with in vivo electrophysiological studies in primates that report independent encoding of task variables by GPe neurons 10,40 , suggesting more a parallel processing of information by GPe neurons rather than an integration of different variables 10 .
Taken together, these results support a lack of integration of information at GPe level and are in line with the intermediate intrinsic timescale exhibited by the GPe in the present study. However, a recent computational study 26 suggests caution around assigning timescales to brain areas too rigidly, because processing different sensory inputs may lead to different timescales based on their model. Moreover, the three BG nuclei studied here exhibited different degrees of heterogeneity in their single unit timescales (Fig. 1, left panel). In particular, the input structures displayed a higher degree of heterogeneity than the GPe. Some studies have shown that within each cortical area, the individual intrinsic timescale computed during a baseline period predicted the strength of response modulation during following task periods in the LIP 41 , PMd 6 , and dlPFC 21 , suggesting that neurons with longer timescales are more involved in the encoding of task-related information. We could hypothesize that the high degree of heterogeneity found in the input structures of the BG could serve to support the heterogeneity of information that these structures have to process with different temporal integration requirements, although further study is necessary to reach conclusions about these functions. Our study is the first to quantify intrinsic timescales of BG nuclei, and some limitations should be considered. First, our datasets are relatively small compared to others used for the cortex and have high variability at the population level. This is mainly true for the input structures, which showed a higher degree of heterogeneity at the single-cell timescale level. Moreover, our selection criterions excluded a large number of cells, especially in the input structures, suggesting that the decay constant of only a subpopulation of PANs and STN neurons can be well estimated with an exponential fit. Second, the BG nuclei are known to be partially specialized in sensorimotor, associative, and limbic territories, and our datasets cover mainly the associative and the limbic parts of these structures. For future studies, it remains to be investigated how our BG nuclei timescale estimates could be generalized at the whole-structure level, and whether our estimates are consistent with the timescales of other BG nuclei not studied here. Third, we performed our analysis during the baseline period as in 4 . It offered the possibility to compare the timescales values of subcortical and cortical structure but did not fully represent the neuronal activity at rest. It is also important to investigate the relationships within each structure at the single-cell level between timescales and show persistent representations of task-relevant signals, as has been done in the cortex 6,20,21,39 . We believe that notwithstanding these limitations, the timescales reported here could be useful as a first approximation for the validation of computational models of action selection based on evidence accumulated through cortico-striatal/subthalamic synapses and architecture of cortico-BG-cortical loops, as has been done for the cortex 26 .

Material and methods
Resource availability. Further information and requests for resources and reagents should be directed and will be fulfilled by Simon Nougaret (simon.nougaret@univ-amu.fr).
The datasets and code supporting the current study have not been deposited in a public repository but are available from the corresponding author on request.
Use of experimental animals. All experimental procedures were conducted in conformity with the French laws on animal experimentation, with the European directive 2010/63/UE on the protection of animals used for scientific purpose and with the ARRIVE guidelines. All experimental protocols were reviewed and approved by the local ethical committee of the Institut de Neurosciences de la Timone.
Experimental model and subject details. Two male rhesus monkeys (Macaca mulatta) were used in this study. Relevant details about the animals have already been reported in 9,10 . Dataset. The experimental details of the datasets used in the current study have already been reported 9,10 .
The single-unit activity of two male rhesus monkeys (Macaca mulatta), recorded from three populations of neurons within the (BG), were analyzed during a foreperiod, considered as a baseline period in which no cognitive process was engaged. During this 1 s period, both monkeys needed to maintain a basal pressing force on a lever www.nature.com/scientificreports/ and wait for the presentation of a pair of visual stimuli indicating the amount of force needed and the amount of reward to be expected at the end of the trial. The data set consisted of 1) a population of 158 putative projection neurons recorded in the striatum, also called phasically active neurons (PANs); 2) a population of 78 neurons from the subthalamic nucleus (STN); and 3) a population of 92 irregular neurons from the external globus pallidus (GPe).
Quantification and statistical analysis. Neuronal sample selection and population timescale analysis. We analyzed the recorded activity of neurons in the STN, GPe, and striatum (PANs) from 100 ms after the beginning of the trial until the end of the baseline period for a total of 900 ms of baseline activity. The same baseline period had already been chosen 9 for electrophysiological analysis of neurons recorded for the same task. We included only correct trials in the following analyses. For each structure, we selected neurons satisfying the following criterium: We performed all analyses with MatLab (The MathWorks, Inc., Natick, MA, USA). For all the following fit procedures, we used the Trust-Region algorithm for Nonlinear Minimization problem implemented in the Optimization Toolbox in Matlab.
To assess the single neuron spike-count autocorrelation structure, we calculated the spike count during the baseline period in 50 ms time bins. It is worth noting that the results did not change within a difference of 20% of the bin length. Given a neuron, the spike-count autocorrelation value across trials between time bins k and j (k, j as integer numbers) at a time lag equal to |k-j|x Δ (Δ = 50 ms), the Pearson's correlation coefficient r is defined as follows 4 : where N(k) and N(j) are the spike counts computed in the k and j time bins, respectively, and N(k) and N j are the spike counts averaged across trials in k and j time bins, respectively. The covariance (Cov), the variance (Var), and the autocorrelation value r were computed for each possible combination of pair-bins (k, j). We calculated the autocorrelation values as a function of the time lags for each neuron satisfying the criterion in (1). We subsequently performed an exponential fit to each neuron's autocorrelation as defined below 4 : where nΔ indicates the time lag between the time bins k and j, with n =|k-j | (n = 1,2, …, 18); r is the autocorrelation value at time lag nΔ; A is the amplitude; τ is the decay constant of the exponential function, called intrinsic timescale; and B is the offset that mirrors the value of r in the limit of time lag nΔ → ∞ (i.e., time lag values much larger than our 900 ms baseline length). Throughout this paper, we refer to intrinsic timescale as simply timescale or τ.
We further reduced the neuronal sample by selecting those neurons satisfying both the criterion in (1) and the following requirements: where R 2 is the coefficient of determination obtained by the fit. The first requirement was introduced because a negative or 0 ms τ value is meaningless; the second requirement of an R 2 larger than 50% was a trade-off between the need to keep as many neurons as possible and the importance of having a good fit. We also excluded outliers, defined as neurons having an intrinsic timescale below the 5th percentile and above the 95th percentile of the τ distribution, from the sample. This last requirement was established due to the heterogeneity of timescale values within each brain structure and to avoid having an estimate of the mean of the intrinsic timescales biased towards the outlier values.
Moreover, we decided to discard those neurons poorly fit by the exponential function by eye inspection 20 . In more detail, we removed a neuron whether the autocorrelation values had an oscillatory behavior with increasing time lags, or the autocorrelation values showed a flat shape with increasing time lags.
In order to assess the timescale of the whole neuronal population, we first averaged the coefficient r defined in formula (2), at a fixed time lag, across neurons satisfying all the previous selection criteria. We subsequently fit the averaged autocorrelation values at increasing time lags with the exponential function defined in (3). The estimate of the decay constant τ was the timescale assigned to the neuronal population for each brain structure.
Single-neuron intrinsic timescale analysis. We moved further by characterizing the single neuron timescale distribution. For each brain structure, we computed the mean of the single timescale distribution. To compare the sample mean of the timescale distributions across brain structures, we applied the following procedure. Given 2 sample means τ1 and τ2, with 2 standard errors of the mean, sem1 and sem2, respectively, we computed the difference ∆: (1) 1. each 50 -ms time bin in the baseline period with non -zero mean activity across trials.
(2) r = Cov(N(k), N(j)) Var(N(k)) × Var(N(j)) = < (N(k) − N(k))(N(j) − N(j)) > Var(N(k)) × Var(N(j)) 2. τ > 0 ms; 3. R 2 > 50%, www.nature.com/scientificreports/ In order to assess whether the difference ∆ is compatible with zero (i.e. no significant difference between the 2 mean timescales of the distributions), we quantified the error associated to ∆ as follows: that stems from the general propagation of the error rule. If the value of ∆ is compatible with 0, within the error range, we would say that there is no significant difference between the mean of the timescales distributions.
We further investigated the degree of heterogeneity of the single timescale distribution for each brain structure. To quantify this, we used the coefficient of variation (CV), defined as follows: where sigma and mean are the standard deviation and mean of the timescale distribution, respectively.