Inhibitory control of frontal metastability sets the temporal signature of cognition
Abstract
Cortical dynamics are organized over multiple anatomical and temporal scales. The mechanistic origin of the temporal organization and its contribution to cognition remain unknown. Here, we demonstrate the cause of this organization by studying a specific temporal signature (time constant and latency) of neural activity. In monkey frontal areas, recorded during flexible decisions, temporal signatures display specific area-dependent ranges, as well as anatomical and cell-type distributions. Moreover, temporal signatures are functionally adapted to behaviourally relevant timescales. Fine-grained biophysical network models, constrained to account for experimentally observed temporal signatures, reveal that after-hyperpolarization potassium and inhibitory GABA-B conductances critically determine areas’ specificity. They mechanistically account for temporal signatures by organizing activity into metastable states, with inhibition controlling state stability and transitions. As predicted by models, state durations non-linearly scale with temporal signatures in monkey, matching behavioural timescales. Thus, local inhibitory-controlled metastability constitutes the dynamical core specifying the temporal organization of cognitive functions in frontal areas.
Editor's evaluation
The paper investigates the temporal signatures of single-neuron activity (the autocorrelation timescale and latency) in two frontal areas, MCC and LPFC. These signatures differ between the two areas and cell classes, and form an anatomical gradient in MCC and, moreover, the intrinsic timescales of single neurons correspond with their coding of behaviorally relevant information on different timescales. The authors develop a detailed biophysical network model which suggests that after-hyperpolarization potassium and inhibitory GABA-B conductances may underpin the potential biophysical mechanism that explains diverse temporal signatures observed in the data. The proposed relationship between the intrinsic timescales, coding of behavioral timescales, and anatomical properties (e.g., the amount of local inhibition) in the two frontal areas is novel, the use of the biophysically detailed model is creative and interesting and the claims are convincingly supported by the data.
https://doi.org/10.7554/eLife.63795.sa0Introduction
Large-scale cortical networks are anatomically organized in hierarchies of interconnected areas, following a core-periphery structure (Markov et al., 2013). Within this large-scale organization, the dynamical intrinsic properties of cortical areas seem to also form a hierarchy in the temporal domain (Chaudhuri et al., 2014; Murray et al., 2014). The temporal hierarchy arises from increasing timescales of spiking activity from posterior sensory areas to more integrative areas including notably the lateral prefrontal and midcingulate cortex (MCC). Intrinsic areal spiking timescales are defined from single unit activity autocorrelation (Murray et al., 2014). Long spiking timescales potentially allow integration over longer durations, which seems crucial in the context of higher cognitive functions, learning, and reward-based decision-making (Bernacchia et al., 2011). Recent studies uncovered links between single unit working memory and decision-related activity and spiking timescales in the lateral prefrontal cortex (LPFC) (Cavanagh et al., 2018; Wasmuht et al., 2018). However, the mechanisms that causally determine the timescale of cortical neuron firings and their role in the functional specificity of areas remain to be described.
To address this question, we recorded in the MCC and LPFC, because these two frontal areas both display particularly long spiking timescales and are functionally implicated in cognitive processes operating over extended timescales. These interconnected regions collaborate in monitoring performance and in integrating the history of outcomes for flexible decisions (Kennerley et al., 2006; Khamassi et al., 2015; Kolling et al., 2018; Medalla and Barbas, 2009; Rothé et al., 2011; Seo and Lee, 2007; Womelsdorf et al., 2014a). Recent anatomical and physiological investigations revealed that the cingulate region has relatively higher levels of synaptic inhibition on pyramidal neurons than LPFC, with higher frequency and longer duration of inhibitory synaptic currents (Medalla et al., 2017), suggesting that excitatory and inhibitory cell types differentially contribute to the specific dynamics of distinct frontal areas. Moreover, MCC also seems to have a longer spiking timescale than the LPFC (Cavanagh et al., 2018; Murray et al., 2014).
In this context, we sought to understand the relationship between temporal features of spiking activity, local neural network dynamics, and the computations implemented by frontal neural networks. We focused on whether and how different temporal features play distinct roles in different frontal areas. To this aim, we addressed the following questions: what are the exact differences in the temporal organization of spiking in the LPFC and MCC? How do they relate to the distinct roles of excitation and inhibition? Do they reflect cognitive operations, and can they be adjusted to current task demands? Can they be accounted for by local biophysical circuit specificities? If so, do distinct collective network neurodynamics emerge from such areal biophysical characteristics and what are their functional implications?
We examined the contribution of single unit temporal signatures to dynamical differences between LPFC and MCC in monkeys. After clustering units based on spike shape (putative fast spiking [FS] and regular spiking [RS] units), we computed spike autocorrelograms and their temporal signatures (time constant and latency). We discovered that LPFC and MCC differed not only in average time constant, but also specifically in the autocorrelogram latency of their RS units.
Regular and FS MCC neurons showed different temporal signatures. Remarkably, through these signatures, neurons contributed to encoding information at different timescales, that is, information relevant between trials or across multiple trials. Exploring constrained biophysical recurrent network models, we identified the ionic after-hyperpolarization potassium (AHP) and inhibitory GABA-B receptor conductances as critical determinants mechanistically accounting for the difference in spiking temporal signatures between LPFC and MCC. The models predicted how differences in temporal signature amounts to the ability of networks to undergo metastable states with different properties. Indeed, we found, in monkey data, long-lasting states in primate MCC activity but not in the LPFC.
Critically, we show that by controlling states stability and transitions, local inhibition – rather than synaptic excitation (Chaudhuri et al., 2015) – is the major factor setting temporal signatures. Moreover, inhibitory-mediated temporal signatures did not require specific disinhibition between molecularly identified subnetworks of interneurons but naturally emerged from inhibitory weight variability (Wang, 2020).
Results
We analysed 570 units recorded in MCC and LPFC (298 and 272 units, respectively). Using the autocorrelogram of binned spike counts (see Materials and methods), we were able to extract population spiking timescales for a subset of this population (140 and 159 units, respectively, for MCC and LPFC) and observed population autocorrelograms similar to those obtained with other datasets (Cavanagh et al., 2018; Murray et al., 2014; Wasmuht et al., 2018; Figure 1a). At the population level, the characteristic timescale of spiking fluctuation over time, scTAU (the time constant from the exponential fit of the spike count autocorrelogram), was longer for MCC than for LPFC (MCC = 519 ± 168 ms, LPFC = 195 ± 17 ms). In addition, MCC single units exhibited longer individual scTAU than LPFC units (medians, MCC = 553 ms, LPFC = 293 ms; two-sided Wilcoxon signed-rank test on log(scTAU), W = 15,192, p<10–8; Figure 1b), as in previous datasets (Figure 1c in Cavanagh et al., 2018). Aside from being characterized by a slow decay (long scTAU), the MCC population autocorrelation displayed a distinctive feature: a positive slope at the shortest time lags equivalent to a latency in the autocorrelogram, that can be observed in previous publications (see Figure 1c in Murray et al., 2014, Figure 1d in Cavanagh et al., 2018). However, the method we employed above (derived from Murray et al., 2014) cannot resolve the fine dynamics of neuronal activity at short time lags. To improve upon this approach, we developed a method based on the spike autocorrelogram of individual units from all spike times (named spike autocorrelogram below), which provides high temporal precision in parameter estimation and is computed using the spike time series in the entire or subset segments of recordings (see Materials and methods).

Midcingulate cortex (MCC) and lateral prefrontal cortex (LPFC) spike count autocorrelograms.
(a) Population exponential fit. Autocorrelograms were computed for each unit and the fit was performed on all the units of the MCC (dark grey) and the LPFC (blue) to extract the decay parameter scTAU (as in Murray et al., 2014). (b) Single unit fits were used to capture individual spiking timescales and produce the distribution of scTAU values for each region. Dotted lines represent the median of scTAU. (c) Clustering of spike shape. After extracting the spike width and amplitude from each unit average waveform, we performed a hierarchical clustering revealing the presence of three groups of units (coloured groups RS1, RS2, FS; see Materials and methods). Fitting Gaussian mixed model on the population (lines) confirmed the presence of the three clusters. In the paper, units with narrow spike width were termed as fast spiking (FS), whereas units with broader waveform were marked as regular spiking (RS: RS1 + RS2).
One basic assumption to explain local dynamical properties is that interactions between cell types (e.g. pyramidal cells and interneurons) might induce specific dynamics in different areas (Medalla et al., 2017; Wang, 2020; Womelsdorf et al., 2014b). To separate putative cell populations in extracellular recordings, we clustered them using single unit waveform characteristics (Nowak et al., 2003). Although associating spike shapes to cell types is not a fully reliable methods for cell-type identification (Vigneswaran et al., 2011), several studies have shown that on population data different cell types and coding properties can be clustered in this way (Krimer et al., 2005; Trainito et al., 2019). Clustering our dataset discriminated three populations, with short, large, and very large spikes (Figure 1c). The results below were obtained using two clusters (small, and large + very large), as detailed analyses showed no clear difference between large and very large spike populations (see Figure 2—figure supplement 1). We classified units as FS (short spikes; nMCC = 41, nLPFC = 57 units) or RS (long spikes; nMCC = 257, nLPFC = 215 units) which, in previous studies, were associated to putative interneurons and pyramidal cells, respectively. In the rest of the paper, and especially for the purpose of modelling, we thus assume simplistically an equivalence between FS vs. RS and interneurons vs. excitatory neurons.
MCC temporal signatures differ for RS units
From spike autocorrelograms we extracted multiple metrics, namely the peak latency (LAT; the time lag of the peak of the autocorrelogram) and time constant (TAU) (see Materials and methods). Together, TAU and LAT constituted the temporal signature of single neurons spiking dynamic. The success rate of fitting an exponential function on spike autocorrelograms using the whole recordings was 91.4% and largely outperformed the alternative method (see Materials and methods). Figure 2a shows comparative examples. All subsequent analyses of this study were performed on this pool of units (nMCC-FS=39, nMCC-RS=225, nLPFC-FS=55, nLPFC-RS=202). Note that because of the methodological criteria on spike numbers required for good fitting, the sample size of units can change depending on the analysis, especially when restricting recordings to specific time periods. Note also that in the pool of neurons where TAU was successfully extracted using both methods (n=280, see Materials and methods for criteria), we found a correlation between the two measures (Murray methods – scTAU – vs. spike autocorrelograms; Spearman’s correlation: rho(282) = 0.46, p<10–15) although scTAU were overall larger, as observed by another recent study using a different method (Spitmaan et al., 2020). Importantly, TAU was not correlated with firing rate across units (Figure 2b, Figure 2—figure supplement 2).

Spike autocorrelogram and temporal signatures in midcingulate cortex (MCC) and lateral prefrontal cortex (LPFC).
(a) Three single examples of spike count (purple, scTAU) vs. normalized spike autocorrelograms (green) contrasting the outcome of the two methods. The measured time constant (TAU) is indicated for both when possible. Numbers of spikes used for each method is also indicated. (b) TAU values extracted from each methods are significantly correlated (n=280, Spearman’s rho(282) = 0.46, p<10–15). (c) Distributions of TAUs (upper histograms) and peak latencies (LAT – lower histogram) for fast spiking (FS) (left) and regular spiking (RS) (right) units. ‘n’ indicates the number of units. Vertical dashed lines indicate medians of respective populations. Boxplots on the right show the respective population data. TAU values were longer in MCC (dark grey) than in LPFC (blue) for both FS and RS (linear model fit on BLOM transformed TAU for normality, TAU = region * unit type, region: t=−4.68, p<10–6, unit type: ns, interaction: ns). Peak latencies significantly differed between MCC and LPFC for RS but not for FS units (medians: MCC FS = 48.5 ms, RS = 102.0 ms, LPFC FS = 48.5 ms, RS = 51.8 ms; linear model fit on BLOM transformed latency for normality, latency = region * unit type, interaction: t-value=−3.57, p<10–3).
TAU was higher on average in MCC than in LPFC for both RS and FS cells (medians ± sd: MCC FS = 284.7 ± 132 ms, RS = 319.5 ± 199 ms, LPFC FS = 175.1 ± 67 ms, RS = 191.6 ± 116 ms; linear model fit on Blom transformed TAU for normality, TAU = area * unit type, area: F(1,520)=18.36, p<10–4, unit type: F(1,520)=2.72, p=0.12, interaction: F(1,520)=0.19, p=0.79) (Figure 2c; individual monkey data in Figure 2—figure supplement 3).
Additionally, our new approach allowed us to extract LAT, which captures other aspects of neurons temporal dynamics. Importantly, it differed significantly between MCC and LPFC for RS but not for FS units, with MCC RS units having particularly long latencies (median ± sd: MCC FS = 48.5 ± 30 ms, RS = 108.7 ± 64 ms, LPFC FS = 48.5 ± 35 ms, RS = 51.9 ± 46 ms; linear model fit on Blom transformed LAT for normality, LAT = area * unit type, interaction: F(1,520) = 11.81, p<0.005) (Figure 2c).
TAU and LAT both reflect temporal dynamics, but those measures were significantly correlated only in LPFC RS units (Spearman’s correlations with Bonferroni correction, only significant in LPFC RS: rho(203) = 0.29, p<10–3). The absence of correlation suggested TAU and LAT likely reflect different properties of cortical dynamics. Moreover, the data also suggested that the different temporal signatures of RS units could reflect differences in the physiology and/or local circuitry determining the intrinsic dynamical properties of MCC and LPFC.
MCC temporal signatures are modulated by current behavioural state
A wide range of temporal signatures might reflect a basic feature of distributed neural processing (Bernacchia et al., 2011). But do different temporal signatures play distinct roles in terms of neural processing in different areas? And, are these signatures implicated differentially, depending on task demands? These questions are unresolved, although recent studies suggest a lack of relationship between individual neuron timescale and selectivity to task-relevant signals (Spitmaan et al., 2020). As single units were recorded while monkeys performed a decision-making task (described in Stoll et al., 2016; Figure 3a), we extracted each unit’s temporal signature separately for periods in which monkeys were either engaged in the cognitive task or were pausing from performing the task (Figure 3b). TAU extracted during engage and pause periods were significantly correlated across neural populations (MCC FS n=19, LPFC FS n=21, MCC RS n=80, LPFC RS n=97, Pearson correlation: r(215)=0.20, p=3.0e-3), indicating that TAU reflects stable temporal properties across conditions (corrected from time-on-task, see Materials and methods). The MCC RS population exhibited a significant modulation of TAU, expressing longer TAU during engage periods compared to pause periods, suggesting that engagement in cognitive performance was accompanied by a lengthening of temporal dynamics for RS neurons in MCC (Figure 3c top) (Wilcoxon signed-rank test with Bonferroni correction, only significant for MCC RS: MD = 1.06, V=2467, p=3.9e-7). To control for a time-on-task effect on such timescale modulation, we contrasted pause periods with engaged periods that occurred at similar times within sessions (i.e. considering only engaged periods occurring after the first pause – see limits in Figure 3b, red marks) (nMCCFS = 19, nLPFCFS = 21, nMCCRS = 80, nLPFCRS = 97, Wilcoxon signed-rank test with Bonferroni correction, only significant for MCC RS: MD = 1.06, V=2467, p=3.9e-7).

Behavioural engagement in task and spiking timescale changes.
(a) Schematic representation of the task. At the start of each trial, animals can either initiate a delayed response task (WORK option) which can lead to one reward delivery, or use the CHECK option to check the current size of the gauge (or collect the bonus reward). Each reward in the task contributes to increase the gauge size and bring the bonus availability closer. The graph (bottom) schematized the speed of increase of the gauge size which varies between blocks (fast or slow blocks). (b) Distribution of pauses in sessions. Each line represents the time course of behaviour for one session for monkey H and monkey A. Grey zones represent engagement in the task, black zones represent pauses in work. Red marks indicate the start limit of the first pause in session which defined the beginning of the period taken for control analyses. (c) Boxplots of indices for each unit type and region calculated to estimate potential changes in TAU between engage and pause (top), and between empty and full gauge (bottom). TAUs increased in engage vs. pause only for midcingulate cortex (MCC) regular spiking (RS) units.
We observed no significant variation of LAT with task demands.
Temporal signatures are linked to cognitive processing
Contrary to MCC, LPFC temporal signatures were not modulated by engagement in the task. Multiple cognitive models propose a functional dissociation between MCC and LPFC and indeed empirical data reveal their relative contribution to feedback processing, shifting, and decision-making (Khamassi et al., 2015; Kolling et al., 2018; Stoll et al., 2016). One important question is thus whether temporal signatures observed for a given area and/or cell type contribute to selected aspects of cognitive processing. For example, temporal signatures might be adjusted to the current functional context and timescale required to perform a task. In our experiment monkeys gained rewards by performing trials correctly in a categorization task while each success also brought them closer to obtaining a bonus reward (Figure 3a, right panel, see Materials and methods for task description). By touching a specific lever at trial start, animals could either enter a categorization trial or check the status of a visual gauge indicating the proximity of the bonus reward availability. The number of rewards (i.e. correct categorization trials) needed to get the bonus, and thus the speed of the gauge increase, varied across blocks (i.e. either fast or slow). Previous analyses revealed that feedback influenced the likelihood of checking in the following trial (Stoll et al., 2016). Thus, feedback can be considered as information used on a short timescale (within the intertrial period). The animals also built an estimation of the gauge size that was updated upon checking in order to regulate the frequency of checks during blocks, allowing animals to seek and collect the bonus in a cost-efficient manner (Stoll et al., 2016). Gauge size can thus be considered as information used and carried over long timescales.
We first hypothesized that blocks of different speeds and/or gauge encoding could engage neurons and modulate their spiking timescale. This was not the case. TAU values were not significantly modulated depending on the state of the gauge (less vs. more than half full, Figure 3c, bottom), nor related to different speeds (Wilcoxon signed-rank test, median = 1) with Bonferroni correction, for gauge state and gauge speed, all (p>0.6).
Conversely, we assessed whether temporal signatures observed for certain cell types contributed to code-specific aspects of the task. We used mixed effect models on groups of single units to test the contribution of population activity to encoding task-relevant information: feedback in categorization trials (i.e. reward vs. no-reward), and gauge size. The rationale was that feedback information was relevant within the intertrial period, whereas gauge information was relevant across trials between two successive checks. Previous analyses had revealed that both MCC and LPFC units encode such information, although MCC units showed greater contributions (Stoll et al., 2016). We used data from the whole recordings (all periods) and classified both FS and RS units as either short or long TAU units using a median split. We used a time-resolved generalized mixed linear models (glmm; Figure 4a) to reveal notable dissociations between these populations that we complemented using a more classical approach at the single unit level, using Poisson glm and weighted proportion of variance explained (wPEV; Figure 4—figure supplement 1).

Encoding of feedback and gauge size for different unit types and spiking timescales and rostro-caudal distribution.
(a) Regression weights (β-coefficients) for the midcingulate cortex (MCC) (grey) and lateral prefrontal cortex (LPFC) (blue) unit populations obtained from time-resolved glmm for feedback (reward vs. no reward; top graphs) and gauge size (bottom) (see group analyses using glmm’ in Materials and methods). Regression weights are obtained at successive time points covering the entire intertrial period between feedback onset and the lever onset in the following trial. Significant effects are indicated by a red triangle (p<0.05 corrected) when more than two successive bins are concerned, shadings indicate standard deviations. Positive values depict a population activity bias towards negative feedback (top) and positive slope of linear coding for gauge size (bottom). Data are presented for fast spiking (FS) and regular spiking (RS) units (left and right respectively for each panel) and have been obtained on subpopulations with short or long TAU values (determined by a median split). Short and long TAU populations are represented by light and dark colour intensity, respectively. Thick bars above the x-axes indicate significance of the coloured corresponding data compared to a null distribution generated through permutations of median split unit identity. Note in particular the dissociation for RS MCC units with short and long TAU respectively coding for feedback and gauge size.
Early phases of feedback encoding recruited MCC long TAU populations for both FS and RS units (Figure 4a, upper graphs). This discrepancy was confirmed by a difference between early feedback coding in short and long TAU population at the single unit level (Figure 4—figure supplement 1). Interestingly FS units in the MCC were mostly engaged in the first second after feedback onset, with a strong bias towards encoding negative feedback (Figure 4a, upper left, positive estimates). Effects were more transient and involved short TAU units in the LPFC (Figure 4a).
During the intertrial interval, feedback valence was represented in different directions between short and long TAU RS populations in MCC, coding being positive for short TAU populations (higher activity for incorrect feedback) and becoming negative for long TAU populations (higher activity for correct feedback). Conversely, solely the population of MCC long TAU RS units coded for the gauge during the intertrial period (Figure 4a, lower graphs). Single unit analyses confirmed the higher contribution of the long TAU population to gauge encoding (Figure 4—figure supplement 1).
Spiking timescales are anatomically organized in MCC
Spiking timescales measured in MCC and LPFC covered several orders of magnitudes (10–1000 ms; Figure 2c). Because single unit recordings spanned large regions, such wide range could reflect anatomical organization of segregated populations with distinct homogeneous intrinsic properties. Such organization has been observed in MCC with human fMRI (Meder et al., 2017). We indeed found that average TAU values in MCC were higher in more posterior parts, in particular for RS units (ANOVA on Blom transformed TAU: MCC, monkey A: F(5,112)=2.8, p=0.041, monkey H: F(5,54)=3.09, p=0.033; linear regression on Blom transformed TAU: MCC, monkey A: t(1,112)=8.99, p=0.0067, monkey H: t(1,54)=2.22, p=0.28; all p-values are FDR corrected for n=2 comparison per monkey) (Figure 4—figure supplement 2a). This suggests an antero-posterior gradient or heterogeneity of spiking timescales. No such effect was observed in our LPFC data (ANOVA on Blom transformed TAU: LPFC, monkey A: F(6,110)=0.34, p=1, monkey H: F(6,64)=2.49, p=0.066; linear regression on Blom transformed TAU: LPFC, monkey A: t(1,110)=1.09, p=0.60, monkey H: t(1,64)=0.25, p=1; all p-values are FDR corrected for n=2 comparison per monkey). Note that the so-called LPFC data covered several subparts of posterior LPFC (see Stoll et al., 2016). Similar analyses for LAT revealed no consistent heterogeneity within MCC or LPFC (Figure 4—figure supplement 2b).
The consequence of such an organization, knowing the respective functional involvement of units with long and short TAU (Figure 4a), should be an antero-posterior functional gradient. We tested this by separating MCC cells in posterior vs. anterior subgroups and tested their contribution to feedback and gauge encoding (Figure 4—figure supplement 2c). Indeed, posterior RS units activity contributed to positive encoding of gauge size, preceded in time by encoding of positive feedback (negative estimates) (Figure 4—figure supplement 2c, lower and upper right), while anterior RS units showed primarily a contribution to feedback encoding (upper right). Finally, anterior FS units were primarily (in time and in strength) contributing to encoding negative feedback. This remarkable contribution of FS to feedback encoding is studied and discussed further below.
In summary, MCC RS units with relatively short or long TAU contributed to the encoding of task elements relevant over short and long terms, respectively. The spiking timescales seemed to be organized along the rostro-caudal axis in MCC. This suggests a correspondence between cell type, temporal signatures, and their functional involvement in processing specific aspects of cognitive information in different functional subdivisions of cortical regions. The crucial questions thus remain of the mechanistic origin of temporal signatures and of how they relate to cognitive functions.
Biophysical determinants of temporal signatures in frontal network models
To uncover the source and consequences of distinct temporal spiking signatures in the LPFC and MCC, we designed a fine-grained model of local recurrent frontal networks. This model is unique in combining (1) highly detailed biophysical constraints on multiple ionic channels, synaptic receptors, and architectural frontal specificities, and (2) the cardinal realistic features of mammals cortical neurodynamics including the excitation/inhibition balance, high-conductance state of neuronal activity, and asynchronous irregular regime characterizing the awake state (Brunel, 2000; Destexhe et al., 2003; Hennequin et al., 2017). Our specific goal was to evaluate whether biophysical circuit specificities could mechanistically account for differences in LPFC and MCC temporal signatures. We also assessed whether these specificities induce distinct collective network neurodynamics and functional implications, possibly explaining the empirical relationships between temporal signatures, cell type, and information processing. Note that for modelling purposes we equate FS units to GABAergic interneurons and RS units to excitatory neurons while acknowledging that it is a crude simplification.
We first explored, using Hodgkin-Huxley cellular models (see Materials and methods), whether specific frontal temporal signatures may arise from ionic or synaptic properties of individual neurons. Extensive explorations of these models identified the maximal cationic non-specific conductance (gCAN) and potassium after-hyperpolarization conductance (gAHP) as the sole couple affecting both LAT and TAU (Figure 5—figure supplement 1a-b). By contrast, conductance couples setting spiking adaptation, post-inhibitory rebound, and slow synaptic transmission were ineffective in changing LAT and TAU (Figure 5—figure supplement 2). However, we could not find any region of the gCAN and gAHP parameter space that yielded reasonable values for both LAT and TAU (Figure 5—figure supplement 1b). Therefore, the temporal signature of the monkey dataset (Figure 5—figure supplement 1c) was poorly reproduced by the cellular model (Figure 5—figure supplement 1d). Thus, we then assessed whether collective dynamics at the level of recurrent networks models could better account for frontal temporal signatures (Figure 5a–b, see Materials and methods). One-dimensional explorations of the large parameter space failed to identify single biophysical determinants accounting, alone, for differences between monkey LPFC and MCC (RS and FS) temporal signatures (Figure 5—figure supplement 3; Figure 5—source data 1). However, these explorations targeted four parameters of interest regulating either LAT or TAU confirming those already revealed in cellular models (gCAN and gAHP) and uncovering, in addition, NMDA and GABA-B maximal conductance (gNMDA and gGABA-B) whose slow time constants strongly affected network dynamics.

Temporal signature of LPFCm and MCCm recurrent network biophysical models.
(a) Scheme of the frontal recurrent networks modelled, with 80% excitatory (green) and 20% inhibitory (red) neurons and sparsity of synaptic connections. (b) Membrane potential in the 484 excitatory (lower part) and 121 inhibitory (upper part) neurons of LPFC and MCC example network models (respectively LPFCm and MCCm ; 'm' for model) with parameter set to approximate LPFC dynamics (gCAN = 0.025 mS·cm–2, gAHP = 0.022 mS·cm–2, gGABA-B=0.0035 mS·cm–2; see text and legend of Figure 6b for the choice of LPFCm and MCCm standard gAHP and gGABA-B maximal conductances) and MCC dynamics (gCAN = 0.025 mS·cm–2, gAHP = 0.087 mS·cm–2, gGABA-B=0.0143 mS·cm–2). (c) (Upper left) Membrane potential of an example excitatory neuron of LPFCm. Scaling bars 1 s and 10 mV (spikes truncated). (Lower left) Autocorrelogram of this LPFCm example excitatory neuron (black) and its exponential fit (red, see Materials and methods). (Right) Bivariate probability density distribution of autocorrelogram parameters in LPFCm excitatory neurons. Contour lines at 50%, 75%, and 90% of the maximum of the bivariate probability density distribution in LPFC monkey regular spiking (RS) units. (d) Same as (c) for LPFCm inhibitory neurons, with contour lines from the bivariate probability density distribution in LPFC monkey fast spiking (FS) units. (e,f) Same as (c,d), for the MCCm and MCC.
-
Figure 5—source data 1
Summary of the effects of the main parameters determining TAU and LAT in the network model.
- https://cdn.elifesciences.org/articles/63795/elife-63795-fig5-data1-v2.docx
Two-dimensional explorations using these key parameters (Figure 5 and Figure 6—figure supplement 1) identified a single specific setup which demonstrated network dynamics that reproduced the shift from the LPFC-like temporal signature to that resembling the MCC with striking precision. An increase of both gAHP and gGABA-B, in the presence of gCAN, drove the model from an LPFC-like temporal signature (LPFCm) (Figure 5c–d; map and contours: bivariate probability density model and monkeys’ distributions, respectively) towards that of the MCC (MCCm, Figure 5e–f). Specifically, gAHP increased LAT and decreased TAU in excitatory (possibly equivalent to RS) neurons (Figure 6a left) and had no effect in inhibitory (putatively FS) neurons (Figure 6a, right). Besides, gGABA-B decreased LAT in both excitatory and inhibitory neurons (Figure 6a, top) and increased TAU in an intermediate range (Figure 6a, bottom). A bivariate similarity measure of probability density (see Materials and methods) revealed that monkey temporal signatures were robustly reproduced by the model in two large contiguous regions in the (gAHP, gGABA-B) space (from which best fits were drawn), with both conductances increased in the MCCm compared to LPFCm (Figure 6b).

Similarity to monkey lateral prefrontal cortex (LPFC) and midcingulate cortex (MCC) temporal signatures critically depends on AHP and GABAB conductance in the network model.
(a) Mean population LAT (top) and TAU (bottom) in Exc (left) and Inh (right) neurons, as a function of AHP and GABA-B maximal conductances. Blue and grey dots indicate the (gAHP, gGABA-B) parameter values of the best fits for LPFCm and MCCm, respectively. (b) Similarity of the temporal signature between the network model and monkey data in the LPFC (left) and MCC (right), as a function of AHP and GABA-B maximal conductances (see Materials and methods). In (a) and (b), the value for each (gAHP, gGABA-B) is averaged over five simulations. Contour line at 80% of maximum similarity. LPFCm and MCCm (gAHP, gGABA-B) parameter values calculated as coordinates of the contour delimited area’s weighted average. (c) Bivariate probability density distribution of the autocorrelogram LAT and first-order latency (the latency of the inter-spike interval [ISI] distribution) in regular spiking (RS) units in monkey LPFC (left) and excitatory neurons in the example LPFCm (right). The model accounts for two distinct neuronal subsets in RS neurons, where LAT is determined by first-order latency solely (due to gAHP-mediated refractoriness; diagonal band), or in conjunction with other factors (gGABA-B slow dynamics-mediated burstiness and recurrent synaptic weight variability; horizontal band). (d) Single excitatory neuron frequency/intensity relationship in LPFCm (blue) and MCCm (grey) in response to a constant injected current.
Several lines of evidence further indicated the model’s relevance. First, spiking statistics were similar to those of monkeys (Figure 7—source data 1). Then, the model properly accounted for the larger LAT variability in monkey RS vs. FS units (Figure 5). Moreover, it reproduced the complex relations between LAT and first-order latency (inter-spike interval [ISI] distribution latency) remarkably well, and in all populations (Figure 6c and Figure 6—figure supplement 2). Furthermore, both the firing frequency and input-output gain were lower in MCCm excitatory neurons (Figure 6d), because of its higher gAHP (Naudé et al., 2012), as found experimentally (Medalla et al., 2017).
Metastable states underlie LPFC and MCC temporal signatures
The asynchronous irregular (presumably chaotic) dynamics of network models was highly structured in time (Figure 5b). Hidden Markov models (HMMs) revealed that it organized through collective transitions between the so-called metastable (quasi-stationary) states in model neural populations (Figure 7a) or pseudo-populations (Figure 7—figure supplement 1; see Materials and methods) in the LPFCm and MCCm, as found in frontal areas (Abeles et al., 1995; Seidemann et al., 1996; Xydas et al., 2011). Moreover, while LPFCm states maximally lasted a few hundred milliseconds (Figure 7b, left, blue), MCCm states persisted up to several seconds (Figure 7b, grey). This suggested that such a difference in metastability may also parallel the difference of temporal signature in monkey LPFC and MCC areas. Applying HMM to neural pseudo-populations built from experimental data revealed that, as predicted by the model, neural activity was organized as metastable states at slower timescales in the MCC (vs. the LPFC, Figure 7b, right). State durations were globally shorter in models (compared to monkeys), as they contained neither temporal task structure nor learning (see Discussion) and were not optimized to fit data.

Properties of metastable states in the lateral prefrontal cortex (LPFC) and midcingulate cortex (MCC).
(a) LPFCm and MCCm spiking raster plots (black dots), with Hidden Markov model states (HMM, coloured bands). (b) State duration distributions: probability distributions of being in states of given durations in LPFCm (blue), MCCm (grey), MCCm with LPFCm gAHP (MCCmLPFC AHP pink), and MCCm with LPFCm gGABA-B (MCCmLPFC GABA-B, orange) models (left) and monkey LPFC (blue) and MCC (grey) areas (right). Each model was simulated 100 times and analysed via HMM, while monkey data was analysed via HMM with 100 different initiation parameter states. Periods above 300 s were excluded. (c, d) Regulation of state duration and short states: median state duration (c) and Kolmogorov-Smirnov one-sample test statistic or maximal distance of state duration probability distributions to log-normality, as a measure of the over-representation of short states (d), as a function of gAHP and gGABA-B maximal conductances. Coloured disks indicate parameter values of models LPFCm, MCCm, MCCmLPFC AHP, and MCCmLPFC GABA-B, respectively. Each point is the average of five simulations. (e) Separation between states: average distances between HMM states (averaged pairwise distance between neural centred standardized frequency centroids [temporal averages] of HMM states), as a function of median state durations. Distances calculated over 100 simulations in models and once for monkey LPFC and MCC data. (f) State segregation: projection of neural activity on the principal components of the principal component analysis (PCA) space of example model simulations and of monkey data. State colours as in (a). (g) Frontal processes and state regulation: schematic attractor landscapes in the LPFC and MCC. Horizontal and vertical arrows indicate possible regulations of AHP and GABAB conductance levels respectively by intrinsic/synaptic plastic processes or neuromodulation in the LPFC and MCC. Likely functional processes operating in these landscapes are indicated in blue for the LPFC and grey for the MCC. (h) Inhibitory control of state transitions: probability to escape an ongoing state (left) and to reach a target state (right), when the ongoing state is perturbed by substituting a given proportion of its excitatory (vs. inhibitory) neurons’ activity by that of the same neurons in the (perturbing) target state (see Materials and methods). Average (full line), ± s.e.m. (shaded areas, almost imperceptible).
-
Figure 7—source data 1
Spiking statistics comparison between monkey and model data.
Mean firing rates, coefficient of variations (CV), CV2, Lv, and Fano factors in monkey and model data for individual neurons, and averaged them across all neurons in each case. For the model, reported values are averaged over the network and across 100 simulations.
- https://cdn.elifesciences.org/articles/63795/elife-63795-fig7-data1-v2.docx
-
Figure 7—source data 2
Analysing the causal relationship between neural frequency drift and Hidden Markov model (HMM) state durations in monkey spike data.
Neurons are divided into two halves – most or least drifting – according to how much their frequency drifts across time to then analyse the spiking activity of each group via HMM (on data from 0 to 600 s). The neural frequency drift averaged across neurons is ~6.6× higher in most drifting vs. least drifting neurons across areas. The average HMM state duration increased by ~1.7× in most vs. least drifting neurons across areas, whereas the ratio of MCC vs. LPFC average state duration across groups was ~5.75×. Thus, neural frequency drift causally increases HMM state durations, but not enough to cause the difference between lateral prefrontal cortex (LPFC) and midcingulate cortex (MCC) average state durations.
- https://cdn.elifesciences.org/articles/63795/elife-63795-fig7-data2-v2.docx
Long states essentially required high gGABA-B in the MCCm, as they disappeared when gGABA-B was lowered to its LPFCm value (MCCmLPFC GABA-B model, Figure 7b left, orange curve). In contrast, they only marginally depended on gAHP. MCCm and an MCCm with the gAHP derived from that of LPFCm (MCCmLPFC AHP) showed state duration distributions that were essentially similar, although there was a small increase in the probability of short states at lower gAHP (pink vs. grey curves). In the (gAHP, gGABA-B) space, gGABA-B systematically proved to be essential in increasing the duration of states, with a border region that clearly separated short states (<0.1 s) from longer states (>1 s) (Figure 7c) At this intermediate border, lower gAHP increased the probability of short states (grey vs. pink dots; distributions were even bimodal at lowest gAHP values, not shown), as witnessed by departure from log-normality (Figure 7c). As such, the temporal structure of states in the LPFCm was dominated by short and unimodal state duration distributions (Figure 7c and d, blue dots), as in monkeys (Figure 7b, right) and previous studies (Abeles et al., 1995; Seidemann et al., 1996). In the MCCm, by contrast, the distribution displayed large durations and a slight departure from log-normality (Figure 7c and d, grey dots), resulting in a majority of long states (>1 s) coexisting with short states, as found in data (Figure 7b).
State duration, that is, stability, scaled with spatial separation in the neural space of activity (Figure 7e, see Materials and methods). Indeed, the shorter states of network models with lower gGABA-B (LPFCm and MCCmLPFC GABAB, blue and orange dots) were less distant, compared to those of networks models with higher gGABA-B (MCCm and MCCmLPFC AHP, grey and pink dots). While states were largely intermingled in the LPFCm and MCCmLPFC GABAB (Figure 7f, upper and middle left), they clearly segregated in the MCC and MCCmLPFC AHP (Figure 7f, upper and middle right). As predicted by the model, segregation between states was indeed higher in the monkey MCC (Figure 7e, large grey triangle, and Figure 7f, lower right), compared to the LPFC (Figure 7e, large blue triangle, and Figure 7f, lower left). This suggests that the higher stability of states in monkey MCC arose from a larger segregation of representations in the space of neural activity.
Altogether, these results suggested that itinerancy between metastable states constitutes a core neurodynamical principle underlying the diversity of computational processes and functions operated in primate frontal areas (Figure 7g, see Discussion). From this perspective, the conditions governing transitions between states is critical. We thus evaluated how perturbations of selective neuronal populations would escape ongoing states and reach specified target states (Figure 7h). In the MCCm, we substituted the membrane potentials and synaptic opening probabilities of a fraction of excitatory (vs. inhibitory) neurons of the ongoing HMM state by those of a target state. This could mimic the effect of internal chaotic fluctuations or external inputs aimed at reaching that target state. Surprisingly, escaping the ongoing state or reaching the target state remained quite unlikely when substituting excitatory neurons, whatever the fraction (Figure 7h, left). By contrast, both probabilities of escaping and reaching scaled with the fraction of substituted inhibitory neurons, with high maximal probabilities (mean: 0.89 and 0.59 for escaping and reaching, respectively – Figure 7h, right panel). Interestingly, the probability of escaping a state could attained 0.24 even with as few as 2% of substituted inhibitory neurons, indicating the significant impact of single inhibitory neurons on state itineracy.
Thus, inhibition is a major factor controlling targeted transitions between metastable states in the MCC network model and is also crucial in determining their stability. Excitation had no such role. This result is remarkable, especially considering that MCC FS neurons encoded negative outcomes immediately after feedback onset that triggered behavioural adaptive responses (Figure 4). This could reflect the involvement of MCC FS neurons in inducing state changes on feedback associated to behavioural flexibility.
Discussion
We showed LPFC and MCC displayed long population spiking timescales (TAU), with larger values in MCC, consistent with previous observations (Chaudhuri et al., 2015; Murray et al., 2014). In fact, LPFC and MCC express distinctive and complex temporal organizations of their activity, which cannot be solely captured by the population spiking timescale. The spiking timescale has been used as a measure characterizing intrinsic areal properties and an inter-area temporal hierarchy. However, the spiking timescale of single units varied over two orders of magnitude within each area (Cavanagh et al., 2018; Murray et al., 2014; Wasmuht et al., 2018). The latency of autocorrelogram also demonstrate informative variability, which suggest important underlying functional richness. Our study demonstrates that the temporal signature (TAU and LAT) of single units, measured through spike autocorrelogram metrics and cell-type segregation, can highlight specific local ionic and synaptic mechanisms. Differences in temporal signatures, for instance between LAT (the time lag of the peak of the autocorrelogram) of FS and RS in MCC, and within regions, provide important information on the functional properties of the underlying neural network.
Unravelling the multidimensional nature of LPFC and MCC temporal signatures at the level of individual neurons enabled us to constrain refined biophysical recurrent network models and reveal the local biophysical determinants mechanistically accounting for their specific temporal organization. Moreover, we showed that these determinants control neurodynamical features that constitute core computational foundations for the executive cognitive processes operated by these frontal areas.
Functional spatio-temporal organization of temporal signatures in frontal areas
The relationship between temporal signatures and behaviour suggests how such biophysical properties could contribute to functional specificities. Such functional relations are still debated. Spiking timescales distributions have been related to persistent activity, choice value, and reward history in the LPFC and MCC (Bernacchia et al., 2011; Cavanagh et al., 2018; Meder et al., 2017; Wasmuht et al., 2018), but in a recent study no correlation was observed between task-dependent and intrinsic timescales at the unit level (Spitmaan et al., 2020). In all those studies, however, cell types were not considered. Here, we could not estimate task-relevant timescales to correlate with TAU, but we found that the spiking timescales of MCC RS units increased on average during periods of engagement in cognitive performance, likely reflecting the global implication of neural processes in task performance at long timescales. MCC units with different temporal signatures differentially contributed to cognitive processes known to engage MCC, namely feedback/outcome processing and outcome history representations (Kennerley et al., 2009; Quilodran et al., 2008; Seo and Lee, 2007). Outcome processing generally enables rapid – trial by trial – adaptation of control and decisions, while outcome history representations contribute to the long-term – across trials – establishment of values guiding strategy adaptation (Behrens et al., 2007; Karlsson et al., 2012). Here, population analyses suggested that short spiking timescale units contributed to feedback processing, whereas long spiking timescale units and especially RS units contributed to encode gauge size, which linearly increase with the accumulation of rewards across trials. In MCC, this temporal dissociation coincided with a spatial organization along the antero-posterior axis: anterior units mainly encoded feedback valence, more strongly and earlier than posterior units, whilst posterior units mostly encoded the long-term information related to gauge size. This antero-posterior gradient strikingly resembles that observed in humans (Meder et al., 2017).
Local molecular basis of frontal temporal signatures
Through extensive parameter exploration of constrained biophysical frontal network models, we identified two conductances that precisely reproduced all monkey temporal signatures. In the model, higher TAU (i.e. MCC vs. LPFC, posterior vs. anterior MCC) was accounted for by stronger synaptic GABA-B levels, consistent with reported higher GABA-B receptor densities (Zilles and Palomero-Gallagher, 2017), stronger and slower inhibitory currents in the MCC (vs. LPFC) (Medalla et al., 2017), and stronger GABA-B receptor densities in the posterior (vs. anterior) MCC (Palomero-Gallagher et al., 2009). Excitatory synaptic transmission has been proposed to be a crucial determinant of longer spiking timescales in the temporal cortical hierarchy (Chaudhuri et al., 2015). We found that while stronger excitatory transmission increases TAU (possibly accounting for longer MCC TAUs), it also decreases LAT. LAT, however, was longer in the monkey MCC. This inability to reproduce the temporal signature pattern of frontal areas suggests that GABA-B inhibitory – rather than excitatory – transmission is likely the principal causal determinant of longer spiking timescales, at least in the LPFC and MCC. Noticeably, long timescales do not require strong inhibitory-to-inhibitory connections (Kim and Sejnowski, 2021) nor specific disinhibition between molecularly identified subnetworks of interneurons (Wang, 2020), but of strong slow inhibition to both excitatory and inhibitory neurons. Note also that long timescales naturally emerge from weights variability (see below) and does not require synaptic learning as found elsewhere (Kim and Sejnowski, 2021). The model also predicts that higher LAT in the MCC originate from increased refractoriness through higher AHP conductances in RS units (which increases first-order latency). Higher AHP implies lower input-output gains in MCC RS units, compared to the LPFC (Naudé et al., 2012), as found empirically (Medalla et al., 2017). Finally, reproducing appropriate temporal signatures required the cationic non-specific (CAN) conductance in the areas’ RS units. This was observed in RS of rodent medial frontal areas (Haj-Dahmane and Andrade, 1997; Ratté et al., 2018), where it regulates, together with AHP, cellular bistability and memory, network persistent activity, and computational flexibility (Compte et al., 2003; Papoutsi et al., 2013; Rodriguez et al., 2018; Thuault et al., 2013). Our conclusions do not preclude the contribution of other factors to temporal signatures such as different positions in the anatomical hierarchy, different proportion of excitatory to inhibitory neurons, large-scale hierarchical gradients of other neurotransmitter receptor or receptor subunit expression (Chaudhuri et al., 2015), distinct neuromodulations (see below), different extra-regional inputs, or inputs with different spectral contents to LPFC and MCC.
Frontal temporal signatures uncover metastable dynamics
The LPFC and MCC activity, both in models and in monkeys, was metastable, that is, organized in sequences of discrete, quasi-stationary states in which activity fluctuates around fixed-point attractors (Abeles et al., 1995; La Camera et al., 2019; Rich and Wallis, 2016; Seidemann et al., 1996). Such states were robustly found, whether analysing populations or pseudo-populations of neurons (see Materials and methods). As a general rule, the duration of states increases with the stability of their attractor (i.e. the depth/width of their basin of attraction) and decreases with spiking fluctuations. Fluctuations originate from stochastic inputs or chaotic noise (as in our model), and they trigger state transitions. Here, activity was always present as consecutive states occurred, that is, with no interruption, and therefore departed from UP/DOWN dynamics in which the network was either active or silent (Jercog et al., 2017).
States were longer in monkeys, likely because extensive training induced attractors that were more stable, whereas models displayed less stable attractors that simply resulted from just random connectivity without learning. Thus metastability genuinely emerged from synaptic heterogeneity and did not require strong network clustering (La Camera et al., 2019). We showed that high GABA-B levels are crucial to stabilize states because they amplify the heterogeneity of inhibition and widens attractors, as reflected by higher state separation in the MCCm. In addition, GABA-B’s long time constant naturally promotes burstiness, that is, stable discharge episodes. Finally, higher AHP levels, required for higher LAT in MCC RS units, limited the occurrence of the shortest states, limiting frequent transitions between states. AHP conductances have been implied in other computational functions such as in the linearization of neuronal input-output function (Wang, 1998), network decorrelation (Renart et al., 2010), or the complexity of network dynamics (Cartling, 1993). This diversity may emerge from differences in AHP gating dynamics considered, for example, fast (here) vs. slow (Jercog et al., 2017) AHP currents.
In monkeys and biophysical models, temporal signatures, which correlate with state stability, actually reflect the underlying temporal organization of neurodynamics into metastable states. Interestingly, state durations (up to >10 s) were longer than spiking timescales (<0.5 s), reconciling the apparent discrepancy between typical spiking timescales in frontal areas (<1 s) and the functional timescales at which those areas operate (up to tens of seconds, Bernacchia et al., 2011).
Functional significance of controlled metastable states in frontal areas
Metastable states can be linked to specific representations in the brain at a variety of levels of abstraction, from stimuli to mental states (Engel et al., 2016; La Camera et al., 2019; Mazzucato et al., 2015; Mazzucato et al., 2019; Rich and Wallis, 2016; Taghia et al., 2018). In general, state transitions contain appreciable randomness, with high transition rates signing internal deliberation, whilst more stable states predicting forthcoming decisions (La Camera et al., 2019). We suggest that controlling itinerancy among metastable states constitutes a core neurodynamical process supporting executive functions in frontal areas, which allows to scan choices and strategies, generate deliberation, and solve ongoing tasks.
Specifically, in the MCC (Figure 7g, grey landscape) GABA-B-mediated long metastable states underlying long spiking timescales may contribute to the maintenance of ongoing strategies (Durstewitz et al., 2010; Enel et al., 2016; Stoll et al., 2016) and to the integration of outcome history (Kennerley et al., 2006; Meder et al., 2017; Seo and Lee, 2007; Tervo et al., 2014). At shorter timescales, short states might instantiate dynamic coding, flexible computations, and rapid decision-making in the LPFC (Figure 7g, blue landscape) (Rich and Wallis, 2016; Rigotti et al., 2013; Stokes, 2015). Short states may be lengthened in the LPFC when AHP is increased (Figure 7g, orange landscape), favouring longer timescales and a global stabilization of, for instance, working memory processes (Cavanagh et al., 2018; Durstewitz and Seamans, 2008). Conversely, decreasing GABA-B destabilizes all long states in the MCC model, globally favouring fast transitions (Figure 7g, orange landscape). This mechanism might contribute to abandon prior beliefs and to rapid search for adapted representations, for example, in uncertain environments (Karlsson et al., 2012; Quilodran et al., 2008; Stoll et al., 2016). In the LPFC model with increased GABA-B or in the MCC model with decreased AHP, activity destabilizes certain long states, favouring transitions to remaining long states (Figure 7g, pink landscape). Such a configuration might be relevant for flexible behaviours, directed exploration, and switching (Durstewitz et al., 2010; Pasupathy and Miller, 2005; Russo et al., 2021; Stoll et al., 2016). Regulating GABA-B and AHP to dynamically adapt computations and temporal signatures could be achieved through neuromodulatory or fast plastic processes (Froemke, 2015; Satake et al., 2008).
Macroscopic gradients of inhibitions and excitations appear as important determinants of the large scale organization of cortical dynamics (Wang, 2020; Womelsdorf et al., 2014b). Our results indicate a complementary fundamental dual role of local inhibition in regulating state durations and stability on one hand, and setting the timing and direction of state transitions, on the other. Moreover, transitions can be easily triggered using very few inhibitory neurons. Our study suggests that interneurons and inhibition might be causal in error-driven state transitions in the MCC. Such transitions, initiated by FS neurons immediately after feedback onset (Figure 4), would allow escaping currently unsuccessful states, reaching alternatives or exploring new states.
In conclusion, we showed that local ionic and synaptic determinants specify the scale of temporal organization of activity in frontal cortical areas. These determinants might produce the particularly long states observed in monkey MCC dynamics and could explain its contribution to functions operating over extended behavioural periods. More generally, our results suggest that the diversity of spiking timescales observed across the cortical hierarchy reflects the local excitability- and synaptic inhibition-mediated regulation of metastability, which sets the temporal organization of computational processes.
Materials and methods
Subjects and materials
Request a detailed protocolThis project was conducted with two male rhesus monkeys (Macaca mulatta), monkeys A and H. All procedures followed the European Community Council Directive (2010) (Ministère de l’Agriculture et de la Forêt, Commission nationale de l’expérimentation animale) and were approved by the local ethical committee (Comité d’Ethique Lyonnais pour les Neurosciences Expérimentales, CELYNE, C2EA #42). Electrophysiological data were recorded using an Alpha-Omega multichannel system (AlphaOmega Engineering, Israel).
Recording sites
Request a detailed protocolRecording chambers (Gray Matter Research, Bozeman, MT) were centred on antero-posterior coordinates of +34.4 and+33.6 relative to ear bars (for monkeys A and H, respectively) (Stoll et al., 2016). MCC recording sites covered an area extending over 10 mm (anterior to posterior), and at depths superior to 4 mm from cortical surface (corresponding to the anatomically defined aMCC or functionally defined dACC). Recording sites in LPFC were located between the principalis and arcuate sulcus and just dorsal to the arcuate (areas 6DR, 8B, 8A, and 9/46) and at depths inferior to 2 mm from cortical surface (see supplemental figures in Stoll et al., 2016). Reconstructions of cortical surface, of MRI sections perpendicular to recording grids and of microelectrode tracks were performed using neuronavigation. Locations were confirmed with MRI reconstructions and stereotaxic measurements by keeping track of electrophysiological activity during lowering of electrodes.
Single unit activity
Request a detailed protocolElectrophysiological activity was recorded using epoxy-coated tungsten electrodes (1–2 MOhm at 1 kHz; FHC Inc, Bowdoin, ME) independently lowered using Microdrive guidance (AlphaOmega Engineering). Neuronal activity was sampled at 22 kHz resolution. Single units were sorted offline using a specific toolbox (UltraMegaSort2000, Matlab toolbox, Kleinfeld Lab [Hill et al., 2011], University of California, San Diego, CA). Metrics served to verify the completeness and purity of single unit activity. Each single unit activity was selected, recorded, and included in analyses on the basis of the quality of isolation only. We obtained 298 MCC units and 272 LPFC units while monkeys performed a checking task (Stoll et al., 2016). A subset of these data has been used in a previous publication (Stoll et al., 2016).
Spike shape clustering
Request a detailed protocolSpike shapes can be clustered in different groups that might correspond to different putative cell populations. For each single unit, we computed the average spike shape on which we extracted the spike width, represented by the time between the peak and the trough (maximal and minimal value, respectively), and the spike amplitude defined by the ratio between the minimum value of the waveform following the peak and the peak. To cluster units we first computed the spike width vs. spike amplitude distance matrix (dist function in R). The partitioning led to three clusters, one with narrow spike shapes, one with wide spikes, and one with very wide spikes. Narrow and wide spikes were considered FS and RS, respectively. Although clustering revealed three clusters, no differences were found between the two wide ones, both considered RS neurons (see Figure 2—figure supplement 1). To get a statistical confirmation of the numbers of retained clusters, we then fitted the distribution of spike widths using Gaussian Mixture Models (Mclust function from the package MClust in R, which uses the expectation-maximization algorithm). This method was previously applied for spike clustering (Torres-Gomez et al., 2020). We tested the presence of up to three mixture components with variable/unequal variance, comparing the different models using Bayesian information criterion (BIC). In this context BIC values are an approximation to integrated likelihood and should be maximized (Banfield and Raftery, 1993; Scrucca et al., 2016).
The model which fitted best the spike widths distribution was composed of three Gaussians (BIC unimodal: 255, BIC bimodal: 477, BIC trimodal: 588).
We decided not to use the firing rate for the clustering because we did not have a clear justification for choosing a specific period for firing rate calculation. Yet we found that the so-called FS population we extracted had a higher firing rate than the RS population (Figure 2—figure supplement 2a). This difference is in adequation with the literature and supported our decision to cluster units solely based on spike shape/duration. This difference of firing rate is reported when computed from the whole recording. Firing rate computed from different periods of the recordings (when monkeys are engaged in the task, taking a break or around key task events, etc.) are correlated but variable. Actually, in our recordings the correlation is lowest when considering firing rates between pauses and the fore period of the task, two of the periods which could have been logical candidates for a firing rate of reference.
Spike count autocorrelogram and timescale
Request a detailed protocolThe primary analysis of timescales was based on Murray et al., 2014. Spike counts (sc) were measured in 14 successive bins of 50 ms from the pre-cue period (700 ms) of each trial, when the monkey is in a controlled, attentive state awaiting stimulus onset. We first calculated the cross-trial bin cross-correlations. Each vector of spike counts from the 50 ms bin t was correlated with vectors of spike counts at subsequent bins (t+1, t+2, etc.) generating an autocorrelation matrix. Autocorrelograms were computed for negative and positive lags, producing a histogram symmetric along the zero axis. Timescales were computed using the autocorrelation defined over positive time lags. The autocorrelogram data was then fitted using non-linear least square (nls function in R) to a function of the form:
where R is the correlation coefficient and t the bin time. scTAU, representing the decay of the exponential function and thus the timescale, and A, a scaling constant, were obtained from the fit. We computed scTAU both at the population level, by using a global fit on all recorded units from a given area (as in Murray et al., 2014), and at the single unit level.
However, the above method cannot resolve the fine dynamics of neuronal activity at short time lags because it is based on counts pooled across trials and from coarse-grained time bins (50 ms). Moreover, the large variability of unit discharge resulted in a high variability of autocorrelograms, which could not be fitted in many cases (47.5% failures), as in other studies (52.1% and 48.4% failures in Wasmuht et al., 2018 and Cavanagh et al., 2018, respectively). Finally, tracking the causal determinants of LPFC and MCC temporal signatures in terms of local cellular and/or network dynamics requires a high temporal precision, because they rely on intrinsic and synaptic time constants, which often lie below the coarse time bin of the spike count method. To prevent these shortcomings, we directly computed the autocorrelogram of individual neurons from spike times, allowing for high temporal precision in parameter estimation. For this we leveraged all the data recorded for each neuron to reduce the large noise present at the level of individual neurons.
Spike autocorrelogram analysis
Request a detailed protocolTo capture the dynamics of neuronal activity, we computed autocorrelograms from individual unit spike time series and extracted their latencies (LAT; the time lag of the peak of the autocorrelogram) and time constants (TAU). The same method was applied to units from in vivo recordings and neurons from network models. To do so, we computed the lagged differences between spike times up to the 100th order, that is, the time differences between any spike and its successors (up to ) at the unit level. The lagged differences were then sorted in 3.33 ms bins from 0 to 1000 ms. The resulting counts allowed to build the probability density function of the autocorrelogram (AC) that we multiplied by the inverse of the time bin width so that it peaked at 1 and is graphically more understandable (as in Figure 2). We then smoothed the AC by local non-linear regression (loess method, with span 0.1; to filter high-frequency noise and correctly detect the peaks, see below) after removing its first 10 ms, to eliminate potential source data contaminations, such as ISIs shorter than the absolute refractory period. We defined the peak of the AC as its maximum, except when the maximum was the very first bin, in which case the peak was defined as the first local maximum after the first bin. The latency of that peak, LAT, was considered in further analyses as a structural parameter of the AC characterizing the temporal signature of the neuron/unit spiking set. For each AC, a global mono-exponential fit (GLOBAL fit) was then performed on the part of the AC situated after the peak using the Levenberg-Marquardt algorithm (nlsLM function in R) for monkey data or von-Neumann-Karmarkar interior-point algorithm (fmincon in Matlab) for network models (we checked that either algorithm on the same spiking sets gave similar results), as follows:
TAU, the time constant of the AC fit characterized the temporal signature of the neuron. , the amplitude of the exponential, and , the offset, are positive constants. Note that this mono-exponential fitting equation is strictly equivalent to that of Murray et al., 2014, here corresponding to in the Murray method. Choosing one or the other equation did not affect the resulting fit and we kept the present form as it is easier to interpret. Fits on each AC were performed 50 times, with random initial guesses in the range for , for , and [0, 1000] ms for TAU, from which the best fit was kept.
In a minority of cases (<3% of neurons), the AC following the peak (as defined above and denoted below the 1st peak) could present a shape that diverged from a simple exponential decay, because the first peak was followed by: a fast and large dip, then a second peak (local maximum), then the slower, final exponential decay. In this case, we developed a pipeline aiming at consistently choosing the peak (i.e. 1st or 2nd) from which the fit started. We defined the AC as having a dip if the first local minimum in the 100 ms after the 1st peak was below 75% of the global range of the AC, (to avoid modelling local troughs due to noise as dips). When two peaks were detected, the second peak was defined as the maximum of the AC after the dip. Two additional mono-exponential fits of the AC were then performed, one from the first peak to the dip (FAST fit) and a second one from the second peak to the end of the AC (SLOW fit).
To be valid, any individual fit (whether of the GLOBAL, FAST or SLOW type) had to display positive , , and TAU values. When neurons had a valid GLOBAL fit, two possibilities were considered. First, the valid GLOBAL fit was kept when (1) at least one of the FAST and SLOW fits were not valid or when (2) it was the best (i.e. its root-mean-square error was inferior to that of the sum of the valid FAST and SLOW fits). Neurons that did not have a valid GLOBAL fit were excluded from further analysis. Thus, while FAST and SLOW fits were de facto systematically excluded from further analysis, they were only used to ensure the quality of GLOBAL exponential fits. Note again that excluding <3% of neurons, this complex procedure was very conservative and designed for the sake of fitting performance.
All codes are freely available (Fontanier et al., 2020).
Statistical analysis
Request a detailed protocolAll analyses were performed using R (version 3.6.1) with the RStudio environment (R core team, 2014).
BLOM transformation. As some timescale measures are non-normally distributed, analyses required a robust non-parametric test. We opted for the BLOM transformation which is a subcase of rank-based inverse normal transformations (Beasley et al., 2009). Basically, the data is ranked and then back-transformed to approximate the expected normal scores of the normal distribution according to the formula:
where ri is the ordinary rank and Yi the BLOM transformed value of the ith case among the N observations. Φ−1 is the standard normal quantile (or probit) function and c a constant set to 3/8 according to Blom, 1958. Regular parametric analyses can then be performed on the transformed data. Since z-scores of the transformed data are normally distributed and differences are expressed in standard errors, main effects and interactions can easily and robustly be interpreted. As sanity checks we also ran more classical non-parametric tests (Wilcoxon test) on non-normally distributed data leading to the same conclusions.
Behaviour and context-dependent modulations
Request a detailed protocolBehavioural task. Monkeys were trained to perform a dual task involving rule-based and internally driven decisions (Stoll et al., 2016). Monkeys performed the task using a touch screen. In each trial they could freely choose whether to perform a rewarded categorization task or to check their progress towards a large bonus juice reward (Figure 3a). Upon checking (selection of a disk-shaped lever) progress was indicated by the onset of a visual ‘gauge’ (an evolving disk inside a fixed circle). Choosing the categorization task (selection of an inverted triangle lever) started a delayed response task in which an oriented white bar (cue) was briefly presented, followed by a delay at the end of which two bars oriented 45° leftward and rightward where presented. Selecting the bar matching the cue orientation led to a juice reward. An incorrect response led to no reward delivery. The gauge increased based on correct performance in the categorization task following seven steps to reach the maximum size. If the animal checked while the gauge was full, the bonus reward was delivered, and the gauge reset to step 1. The full gauge was reached after either 14, 21, 28, or 35 correct trials (=number of trials to complete the seven steps, pseudo-randomly chosen in each block). Thus, the gauge could increase at one of four different speeds.
Pause vs. engage periods. As each trial was self-initiated by the animal, monkeys could decide to take a break in their work. We defined pauses as periods of at least 60 s without trial initialization. Monkeys made on average 3.4±2.57 pauses per session (mean ± sd, monkey A: 3.44±2.55, monkey H: 3.34±2.63; see Figure 3b). We extracted spike times during the defined pause and engage time segments for each unit. To control for a time-on-task confound on timescale modulation in this analysis, we contrasted pauses with engaged periods that occurred at the same time of the session (after the first pause). Because engage periods were as long as pause periods for one monkey (monkey H, 53 sessions, MDengage = 392 s, MDpause = 396 s, Wilcoxon-paired test: V=790, p=0.51) and roughly twice as long for the other (monkey A, 96 sessions, MDengage = 638 s, MDpause = 372 s, Wilcoxon-paired test: V=406, p=2.19e-12), we decided not to further segment the data to avoid resampling biases. This analysis was conducted on units for which TAU could be extracted for both periods (nMCC-FS=19, nMCC-RS=80, nLPFC-FS=21, nLPFC-RS=97).
Fast vs. slow-paced blocks. We defined 14 and 21 correct trials blocks to be fast blocks and 28 and 35 correct trials blocks as slow blocks (Figure 3a, bottom). We considered neuronal activity from the first-time monkeys checked in a block until the end of the block. We excluded pause periods from this analysis. We extracted spike timing from the segments and computed timescales as previously, keeping only units with successful timescale extraction for both periods (nMCC-FS=33, nMCC-RS=165, nLPFC-FS=46, nLPFC-RS=165).
Emptier vs. fuller gauge size seen. In each block, monkeys used the gauge size observed upon checking to regulate their future decisions to check. The checking frequency increased with gauge size with a marked increase at steps >4. We thus compared neuronal activity in periods in which monkeys saw gauges of size <4, with periods in which they saw gauges >4, excluding the very beginning of blocks when monkeys have not seen the gauge yet, and pauses periods. We performed this analysis on 430 units (nMCC-FS=30, nMCC-RS=178, nLPFC-FS=47, nLPFC-RS=175).
To test whether current block speed had an influence on TAU at the unit level, we computed a modulation index for each unit: log(TAUslow)/log(TAUfast). Similarly, to test whether gauge filling state had an influence on TAU at the unit level, we computed a modulation index for each unit: log(TAUempty)/log(TAUfull), where TAUfull corresponds to TAU calculated on the spike data recorded during the time in blocks where the gauge was superior of equal to the 4th level.
Task-related analyses
Request a detailed protocolSingle unit activity. Each unit’s spikes were counted in sliding bins of 200 ms overlapping by 50 ms from feedback onset to 800 ms post-feedback and during the intertrial interval from 400 ms before the end of trial signal onset to 2000 ms after its onset.
Group analyses using a glmm. We used a glmm using a Poisson family. p-Values were corrected for multi-comparison with the false discovery rate algorithm with the number of comparisons being the number of timebins (p.adjust function in R).
The mixed models used were of the form:
where γ·Z is the random term, and CheckWork, Gauge, and PreviousFeedback are the fixed effects describing the Check vs. Work decision (0/1), the gauge size (1–7) and the feedback in the previous trial (0/1) with their respective parameters (β). In the glmm, the single unit identity was used as a random factor.
A persistent problem with Poisson models in biology is that they often exhibit overdispersion. Not accounting for overdispersion can lead to biased parameter estimates. To deal with overdispersion we used observation-level random effects, which model the extra variation in the response variable using a random effect with a unique level for every data point.
Median splits. To test the hypothesis that units with different timescales may encode feedback differently, we divided the units into two groups based on the median of the timescale metric. We computed the median of the metric (e.g. peak latency or TAU) in all the units of a given cell type. Then we put units with a metric value below the median into the ‘short’ group and units with a metric value above the median into the ‘long’ group. These splits led to the following population of units: LPFC: FS short: 37, FS long: 18, RS short: 148, FS long: 54 – MCC: FS short: 10, FS long: 29, RS short: 67, RS long: 161.
To assess differences between short and long TAU population coding for a given area and cell type, we have constructed null distributions of coding (z-values) by permuting TAU group allocation of units. Such permutations allowed us to retain differences in sampling (e.g. the population of MCC RS with long TAU is larger than the short TAU one). This procedure was performed 100 times for each area and cell type. We then compared the position of the true data relative to the cumulative distribution of the permutations and set a statistical threshold at α=0.05. Outcomes are shown as raster above x-axes in panel in Figure 4.
Single unit approach. To investigate whether each single unit activity encoded the different key variables of the task, we analysed variations of spike counts measured in each trial using a glm (using the libraries MASS and ggplot2 for graphics under R software) (see also Stoll et al., 2016). Spike counts were measured on successive bins of 200 ms moved smoothly by 50 ms around key event times in each trial. Because of the statistical properties of count data, the glm were applied using a Poisson regression (Poisson error structure). We checked for overdispersion by dividing residual deviance by the degree of freedom. In case of overdispersion, we applied a negative binomial regression using the glm.nb() function in R. To validate this choice for each set of data, we statistically compared the two models (Poisson and negative binomial) fitted for each set (likelihood ratio test, w2-test).
Proportions of significant single units were extracted from the sliding glm if they significantly discriminated the factor of interest for four consecutive bins (covering a time period of 350 ms).
We used also used the wPEV (computed with the function anova_stats from the sjstats package in R) as a statistical measure to quantify the extent to which the variability in neural firing rate was determined by feedback valence and gauge state. We then quantified the time-resolved proportion of cells coding for the task variables (p<0.05 for at least four successive bins) and the wPEV in populations of units with short or long TAU (median split by cell type). To assess differences between short and long TAU populations, we built null distributions by permuting 1000 times the TAU group allocation of units. We then compared the position of the true data relative to the cumulative distribution of the permutations and set a statistical threshold at α=0.05.
Timescale and coding variations along the antero-posterior axis
Request a detailed protocolWe considered the genu of the arcuate sulcus as an anatomical landmark from which we computed distances of recording location along the anterior-posterior axis from MRI reconstructions.
We questioned TAU antero-posterior variability keeping recording locations covering the same range in both monkeys. We ordered locations from the most posterior site for each area. We excluded FS units from statistical analysis due to their disparateness (RS units, monkey A: nMCC = 112, nLPFC = 110; monkey H: nMCC = 54, nLPFC = 64). This analysis was conducted separately between monkeys to account for inter-subject anatomical variability.
To test variation in population coding along the antero-posterior axis, we divided single units into a posterior and anterior group based on the range of locations of each area (posterior MCC from 4.5 to 7 mm, nMCCRSpost = 84, nMCCFSpost = 14; anterior MCC from 7 to 9.5 mm, nMCCRSant = 82, nMCCFSpost = 16; posterior LPFC from 2.5 to 6 mm, nLPFCRSpost = 77, nLPFCFSpost = 19; anterior LPFC from 6 to 8.5 mm, nLPFCRSant = 97, nLPFCFSant = 19). Population coding analysis is described in task-related analyses.
Cellular model of pyramidal neurons in frontal areas
Request a detailed protocolWe built a generic biophysical Hodgkin-Huxley model of the detailed dynamics of membrane potential and of ionic and synaptic currents of individual pyramidal neurons in frontal areas. The model was generic, being endowed with a large set of ionic voltage- and calcium-dependent conductances, to encompass the wide possible repertoire of spiking discharge patterning encountered in vivo. In the model, the membrane potential followed
where C is the specific membrane capacity and the membrane ionic current writes
in which the leak current is
and action potential (AP) currents ( , ) are taken from a previous model we devised to reproduce spike currents of frontal pyramidal regular-spiking neurons (Naudé et al., 2012). The high-threshold calcium current was
where the activation followed first-order kinetics
with a voltage-dependent time constant
where and were fitted from in vitro data (Helton et al., 2005). The infinite activation followed
where and respectively denote the half-activation potential and e-fold slope of the Boltzmann activation voltage dependence, estimated from in vitro data (Helton et al., 2005). The cationic non-selective () current and the medium after-hyperpolarization () current, responsible for frequency adaptation in pyramidal neurons were taken as in Rodriguez et al., 2018, with
and
The activation of both currents, () followed,
with
and
where and respectively denote activation and deactivation kinetic constants consistent with experimental data in layer 5 PFC pyramidal neurons (Faber and Sah, 2007; Haj-Dahmane and Andrade, 1997; Villalobos et al., 2004). The low-threshold calcium () and hyperpolarization-activated () currents were from reference Ritter-Makinson et al., 2019. To account for autocorrelogram parameters, we employed different versions of the model that contained distinct subsets of ionic currents, which have been implicated in adaptation and bursting (, ), rebound (, ), and regenerative and bistable discharge (, , ) in cortical pyramidal neurons. Calcium concentration dynamics resulted from the inward influx due to and and first-order buffering or extrusion (Rodriguez et al., 2018) through:
where F is the Faraday constant, is the basal intracellular calcium concentration, is the buffering time constant, and
is the surface area to volume ratio of an idealized intracellular shell compartment of thickness situated beneath the surface of a spherical neuron soma of radius .
The synaptic current () mimicked in vivo conditions encountered by neurons in the asynchronous irregular regime, summing random synaptic excitatory inputs, through AMPA and NMDA receptors, and inhibitory inputs, through GABAA and GABAB receptors. Thus,
For AMPA, GABAA, and GABAB,
where is the opening probability of channel receptors and the reversal potential of the current. The NMDA current followed
incorporating the magnesium block voltage dependence modelled (Jahr and Stevens, 1990) as
To simulate fluctuations encountered in vivo, all opening probabilities followed Ornstein-Uhlenbeck processes (Destexhe and Paré, 1999)
where is a Gaussian stochastic process with zero mean and unit standard deviation and and are the mean and standard deviation of the opening probabilities. For AMPA and GABAA, the mean was taken as the steady-state value of first-order synaptic dynamics described in the network model (see below):
with pre-synaptic neurons firing at a frequency (with , depending on the type of current considered), an instantaneous increase of opening probability upon each pre-synaptic spike and first-order decay dynamics with time constant between spikes. For NMDA and GABAB, the mean was taken as the steady-state value of second-order synaptic dynamics described in the network model (see below):
For all currents, standard deviations were taken as . Feed-forward excitatory and inhibitory currents were balanced (Xue et al., 2014), according to the driving forces and the excitation/inhibition ratio, through
Model similarity to monkey data
Request a detailed protocolThe bivariate probability density distribution of neuronal TAU and LAT autocorrelogram parameters was estimated in RS and FS units in monkey in both the LPFC and MCC, using bivariate normal kernel density functions. For cellular models, similarity maps to monkey data was determined as following: for each model parameter couple of the map, the similarity to the considered cortex (LPFC or MCC) was defined as the probability density of that cortex to display the TAU and LAT parameters produced by the model. Cellular models with mean firing frequency superior to 20 Hz were considered to discharge in an unrealistic fashion, compared to data, and were discarded. In network models, for each parameter value (one-dimensional explorations) or model parameter couples of the map (two-dimensional explorations), the similarity (S) was defined as the normalized Frobenius inner product between the bivariate probability density distributions of units in monkeys (U) and that of neurons in the network model (N), following
In order to account for the TAU and LAT autocorrelogram parameters for both RS and FS populations, the similarity was calculated separately as RS with Exc and FS with Inh. Seeing as excitatory neurons represent of the neurons in cortex (Beaulieu et al., 1992), the overall similarity was then calculated as
Parametric explorations in the pyramidal neuron model
Request a detailed protocolBiophysical properties of neurons can affect autocorrelation parameters in several ways. In principle, increasing the refractory period (through increased hyperpolarizing ionic conductance) shifts the distribution of first-order lags (ISIs), thus increasing LAT (Figure 5—figure supplement 1a i). Increasing burstiness of the spike discharge (through increased depolarizing conductance-mediated positive feedback) also increases the latency, because higher-order lag distributions are more peaked. Moreover, conductances with slow time constants (including many bursting-mediating conductances) increase that of the autocorrelogram itself. Finally, all these factors may interact in complex ways in vivo to set the spiking pattern that shapes autocorrelations.
We first explored these alternatives with a detailed biophysical Hodgkin-Huxley model of a generic frontal pyramidal cortical neuron, simulated in in vivo conditions. Pyramidal neurons display a huge electrophysiological diversity set by ionic channels, which, together with synaptic inputs, influences spiking patterns. Two conductances, that is, CAN and AHP, were the sole couple able to affect both the LAT and TAU of the autocorrelation (compare Figure 5—figure supplement 1a ii and iii). Interestingly, these conductances are prominent in monkey LPFC and MCC, as well as rodent prefrontal pyramidal neurons where they control regenerative discharge, bistability, and burstiness (Haj-Dahmane and Andrade, 1997; Medalla et al., 2017; Ratté et al., 2018; Rodriguez et al., 2018; Yang et al., 1996). Within physiological ranges, (1) the autocorrelogram LAT essentially increased with the maximal gAHP conductance, while (2) TAU increased in an intermediate range of gAHP and increased with gCAN (Figure 5—figure supplement 1b), possibly accounting for differences between LPFC and MCC in monkeys. Remarkably, the low-threshold calcium (CaT), high-threshold calcium (CaL), and hyperpolarization-activated H conductances, which are ubiquitous and govern spiking patterns through spiking adaptation and rebound, as well as NMDA and GABA-B synaptic input conductances, which display long time constants, were all ineffective in adequately modulating autocorrelation parameters (Figure 5—figure supplement 2).
Computing an estimation of the bivariate probability density distribution of neuronal autocorrelogram parameters for LPFC and MCC RS units (Figure 5—figure supplement 1c) allowed to build a map of the similarity of the cellular model to RS units temporal signatures in monkey LPFC and MCC, defined as the bivariate probability density observed for the LAT and TAU yielded by the cellular model, given a (gCAN, gAHP) couple of parameters (see Model similarity to monkey data). We found that the model displayed large (i.e. sub-maximal) similarity to the LPFC in a substantial region of (gCAN, gAHP) parameters (Figure 5—figure supplement 1d). By contrast, this was not true for the MCC (Figure 5—figure supplement 1d), because the model was unable to generate LAT in the 100–150 ms range that characterizes the MCC (Figure 5—figure supplement 1b), even when exploring large ranges of CAN and AHP conductance kinetic parameters.
Model of local recurrent neural networks in frontal areas
Request a detailed protocolWe built a biophysical model of a generic local frontal recurrent neural network, endowed with detailed biological properties of its neurons and connections. The network model contained neurons that were either excitatory (E) or inhibitory (I) (neurons projecting only glutamate or GABA), respectively (Dale, 1935), with probabilities and , respectively, and (Beaulieu et al., 1992). Connectivity was sparse (i.e. only a fraction of all possible connections exists Thomson et al., 2002) with no autapses (self-connections) and EE connections (from E to E neurons) drawn to insure the over-representation of bidirectional connections in cortical networks (four times more than randomly drawn according to a Bernoulli scheme; Song et al., 2005). The synaptic weights of existent connections were drawn identically and independently from a log-normal distribution of parameters and (Song et al., 2005). To cope with simulation times required for the massive explorations ran in the parameter space, neurons were modelled as leaky integrate-and-fire neurons, that is, the AP mechanism was simplified, compared to the cellular model (see above). Moreover, leveraging simulations at the cellular level, we only considered the and amongst the ionic currents of the cellular model (see above). Thus, the membrane potential followed
where repolarization occurred after a refractory period . The ionic current followed
with parameters and gating dynamics of ionic currents identical to the cellular model. The intra-somatic calcium concentration Ca evolved according to discrete spike-induced increments and first-order exponential decay:
where is the time of the th spike in the spike train of neuron , the Dirac delta function, the time constant of calcium extrusion, the basal calcium, and a spike-induced increment of calcium concentration. The recurrent synaptic current on post-synaptic neuron , from – either excitatory or inhibitory – pre-synaptic neurons (indexed by ), was
The delay for synaptic conduction and transmission, , was considered uniform across the network (Brunel and Wang, 2001). Synaptic recurrent currents followed
where is the synaptic weight, the opening probability of channel receptors, and Vx the reversal potential of the current. The NMDA current followed
with the magnesium block voltage dependence (see cellular model). AMPA and GABAA rise times were approximated as instantaneous (Brunel and Wang, 2001) and bounded, with first-order decay
To take into account the longer NMDA (Wang et al., 2008) and GABA-B (Destexhe et al., 1998) rise times, opening probabilities followed second-order dynamics (Brunel and Wang, 2001)
Recurrent excitatory and inhibitory currents were balanced in each post-synaptic neuron (Xue et al., 2014), according to driving forces and excitation/inhibition weight ratio, through
with approximating the average membrane potential.
The feed-forward synaptic current (putatively arising from cortical and sub-cortical inputs) consisted of an AMPA component
with a constant opening probability .
Numerical integration and parameters of the models
Request a detailed protocolModels were simulated and explored using custom developed code under MATLAB and were numerically integrated using the forward Euler method with time-steps in cellular models and in network models. Computational models are freely accessible on the Zenodo repository (https://doi.org/10.5281/zenodo.5707884) (Fontanier et al., 2020).
Unless indicated in figure legends, standard cellular parameter values were as follows. Concerning ionic currents, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , . Concerning synaptic currents, , , , , , , , , , , , , , , , , , , , , , .
Unless indicated in figure legends, standard parameter values in network models were identical to cellular model parameters, except for the following. Concerning the network, neurons, , so that and . Concerning the weight matrix, , , , . Concerning integrate-and-fire neuron properties and intrinsic currents, , , , , , . Concerning synaptic currents, , , , , , , a.u.
Parametric explorations in the network model
Request a detailed protocolWe first assessed whether variations of a single biophysical parameter could explain the differences in the temporal signature of the MCC, compared to the LPFC: an increased TAU for RS and FS units and an increased LAT for RS (but not for FS) units. To do so, we tested many biophysical parameters determining the architectural, synaptic, and intrinsic properties of the network, but none were able to account for these differences between frontal areas (not shown).
However, these explorations unraveled four model parameters of interest that were able, when varied within their physiological range (i.e. realistic regimes of network activity), to (1) affect either LAT in Exc (but not Inh) neurons, TAU in Exc neurons or TAU in Inh neurons, and (2) do so in a gradual fashion, that is, allowing some possible form of (developmental, homeostatic, or plastic) regulatory control. Indeed, several other parameters could vary TAU or LAT, but they did so abruptly, because their effects occurred at the vicinity of network bifurcations, where network activity dramatically saturated or was silenced, that is, in non-physiological regimes.
The parameters of interest were the maximal conductances, on the one hand, of two membrane ionic currents setting the intrinsic excitability and spiking pattern of cortical pyramidal neurons (gCAN and gAHP) and, on the other hand, of two neurotransmitter-gated channels that set slow synaptic neurotransmission in cortical networks (gNMDA and gGABA-B).
Firstly, decreasing gCAN, the maximal conductance of the spike-triggered calcium-dependent cationic current, which controls regenerative discharge and spiking bistability in pyramidal neurons (Haj-Dahmane and Andrade, 1997), gradually increased LAT in Exc neurons, although this occurred at low values where the network was nearby silence and displayed very low firing frequency (Figure 5—figure supplement 3a, upper left). The CAN current is absent in Inh neurons (Gorelova et al., 2002), so that changing gCAN (i.e. only in Exc neurons) left LAT constant in Inh neurons (Figure 5—figure supplement 3a, lower left). Besides, increasing the CAN maximal conductance gradually increased TAU in Exc neurons (Figure 5—figure supplement 3a, upper right). This arose because gCAN increases burstiness, a factor that can increase the autocorrelogram time constant (Bar-Gad et al., 2001). However, it had no effect on TAU in Inh neurons, where it is absent (Figure 5—figure supplement 3a, lower right). Thus, while gCAN possibly accounted for LAT in frontal areas (increased LAT in MCC Exc neurons, no change in Inh neurons), as well as for the increased TAU in MCC Exc neurons, it could not explain the TAU difference in Inh neurons (Figure 5—figure supplement 3). Moreover, accounting for LAT and TAU in the MCC required incompatible gCAN ranges of values and this was the same for LPFC (Figure 5—figure supplement 3).
Secondly, the maximal conductance of the medium AHP current (gAHP), a spike-triggered calcium-dependent potassium current, which balances the CAN current in the patterning of spiking in pyramidal neurons, increased LAT in Exc neurons (Figure 5—figure supplement 3b, upper left); this arose because AHP increases the refractory period, which can increase LAT (Bar-Gad et al., 2001). Similarly to CAN, the AHP current is absent in Inh neurons (Gorelova et al., 2002), so changing gAHP (i.e. in Exc neurons) left LAT constant in Inh neurons (Figure 5—figure supplement 3b, lower left). Besides, although gAHP is largely known for its effect on firing frequency adaptation, an important determinant of discharge temporal patterning, it displayed an extremely weak effect on TAU (Figure 5—figure supplement 3b, right). Thus, while gAHP possibly accounted for LAT in frontal areas (increased LAT in MCC Exc neurons, no difference in Inh neurons), it could not explain differences in TAU (Figure 5—figure supplement 3).
Together, these effects of intrinsic conductances at the network scale shared important trends with those in the cellular model, inasmuch as gAHP increased LAT and gCAN increased TAU in Exc neurons.
Thirdly, the NMDA receptor maximal conductance, gNMDA, displayed no effect on LAT (Figure 5—figure supplement 3c, left), but it increased TAU (Figure 5—figure supplement 3c, right) in Exc and, to a lesser extent, in Inh neurons, because of its slow synaptic action (decay time constant, 75 ms) on both neuronal types. However, these effects on TAU occurred at high gNMDA values where the network was nearby saturation and displayed unrealistic high-frequency activity. Thus, while gNMDA possibly accounted for TAU in frontal areas and for the absence of change in LAT in Inh neurons, it could not explain the difference in LAT in Exc neurons between LPFC and MCC (Figure 5—figure supplement 3).
Fourthly, the GABAB receptor maximal conductance, gGABA-B, as for the NMDA current, displayed no effect on LAT (Figure 5—figure supplement 3d, left), but it increased TAU (Figure 5—figure supplement 3d, right) both in Exc and Inh neurons, because of its slow synaptic action (rise and decay time constants, 90 and 160 ms, respectively) on both neuronal types. Thus, while gGABA-B possibly accounted for TAU differences in frontal areas and for the absence of change in LAT in Inh neurons, it could not explain the difference in LAT in Exc neurons between LPFC and MCC (Figure 5—figure supplement 3).
Interestingly, both NMDA and GABA-B currents, which had no effect at the individual level, were essential at the network scale, suggesting that the influence of slow synaptic transmission on recurrent collective network dynamics are central in determining the time constant TAU in frontal areas.
In summary, one-dimensional network explorations showed that: (1) gAHP was the sole biophysical parameter that changed LAT in Exc but not in Inh, while keeping network activity within the physiological regime (gCAN also changed LAT in Exc, but at the border of network silencing); (2) gGABA-B was the sole biophysical parameter that changed TAU in both Exc and Inh (gCAN and gNMDA also changed TAU, mainly in Exc, but at the border of unrealistic regimes, that is, network silence and saturation).
Together, these results pointed to gAHP and gGABA-B as major candidates, with the idea that their combined effect in the (gAHP, gGABA-B) space could account for the differences in temporal signature between the LPFC and the MCC. However, because ionic and synaptic conductances typically display strong non-linear interactions whereby some forms of counter-intuitive compensatory or amplificatory effects can emerge, we nevertheless conducted two-dimensional explorations in the (gAHP, gCAN), (gAHP, gNMDA), and (gAHP, gGABA-B) spaces, with the idea that the relative balance of gAHP, which affects LAT, on the one hand, and of either gCAN, gNMDA, or gGABA-B, which affect TAU, on the other hand, could synergistically account for both the larger LAT and larger TAU observed in the MCC, compared to the LPFC (Figure 5—figure supplement 1 and Figure 6—figure supplement 1).
Exploring the (gAHP, gCAN) space, we found that combined increases of gAHP and gCAN could both (1) increase LAT in Exc neurons but not in Inh neurons (Figure 6—figure supplement 1a, top) and (2) increase TAU in Exc neurons (Figure 6—figure supplement 1a, bottom), which translated, quantitatively, as two domains of smaller and larger (gAHP, gCAN) parameter values that displayed higher similarity to LPFC and MCC data, respectively (Figure 6—figure supplement 1b). However, increasing gCAN and gAHP only very weakly varied TAU in Inh neurons (as in one-dimensional explorations), one of the three major changes observed in FS units in the MCC (together with higher LAT and TAU in RS units). While this incapacity marginally reflected in the similarity measure (which integrates similarity in Exc (RS) and FS (Inh) neurons proportional to their relative abundance, that is, 0.2 for Inh neurons), the model therefore revealed qualitatively insufficient to account for the differences in LPFC and MCC temporal signatures.
The exploration of the (gAHP, gNMDA) space indicated a situation where combined increases of gAHP and gNMDA could increase LAT in Exc neurons but not in Inh neurons (Figure 6—figure supplement 1c, top) but could hardly reproduce TAU in Exc and Inh neurons (Figure 6—figure supplement 1c, bottom), so the qualitative agreement was weak. As a result, quantitatively, the domain of largest MCC similarity displayed modest similarities (Figure 6—figure supplement 1d). Moreover, as in one-dimensional explorations, gNMDA increased TAU mostly at high values (gNMDA ~ 1), where the network model was near saturation (unrealistic high frequency), while intracellular recordings show no difference in Exc post-synaptic current amplitudes between MCC and LPFC in monkeys (Medalla et al., 2017).
Contrarily to explorations in (gAHP, gCAN) and (gAHP, gNMDA) spaces, exploration in the (gAHP, gGABA-B) provided two domains of high similarity to LPFC and MCC that indicated a strong quantitative agreement of the model to monkey data (see main text). Moreover, these domains were large (relative to the mean values of gAHP and gGABA-B in said domains) indicating robustness to the inherent biological variability present in frontal cortical structures. Finally, qualitative agreement was present, in addition to the quantitative agreement revealed by the similarity measure, in the sense that all three main qualitative differences between LPFC and MCC (with higher LAT and TAU in RS units and a higher TAU in Inh neurons) were well reproduced in this setup.
HMM analysis
Request a detailed protocolWe used HMM to map the spiking set of neural network models and unit populations in monkeys onto discrete states of collective activity, based on previously established methods (Abeles et al., 1995; Seidemann et al., 1996). Monkey recordings occurred in separate sessions from which samples of typically three to four (maximum five) neurons were recorded simultaneously. HMMs were impractical on such small samples, so it was necessary to create pseudo-populations that could substitute for large simultaneous recordings of the population, to improve the amount of data leveraged and maximize the information extracted in terms of neural states. However, creating pseudo-populations ignores correlated activity between neurons, globally degrading the accuracy of state detection. In principle this degradation affects absolute quantitative state estimations in the LPFC and MCC equally. Therefore, it should not affect the conclusions drawn here, as they are based on relative comparisons between order of magnitudes of state estimates in the LPFC and MCC. Thus, applying HMM to large simultaneous populations would most likely not change the conclusions drawn from experimental data. This point is supported by the fact that we qualitatively obtained similar HMM results in populations and pseudo-populations in the model (see Results), as well as by the general observation that the amount of information that can be time-decoded in associative areas (such as in the LPFC) is globally the same, when using pseudo-populations or simultaneously recorded neurons (Anderson et al., 2007; Meyers et al., 2008). Extended results presented in supplementary materials were obtained using LPFC and MCC pseudo-populations built of groups of neurons drawn from simulations of distinct random neural networks (i.e. synaptic matrices), with group statistics (numbers of groups and numbers of excitatory and inhibitory neurons within each group) similar to those of monkey LPFC and MCC data.
HMM methods allow to determine the probability of the network to be in state at time . Typically, we found that, as previously shown in frontal areas, population activity organized into periods that lasted in the range , that is, transition probabilities were small and states were quasi-stationary. When all probabilities of being in a state , the network was considered to be in the null state , signifying that the network was not in any of the states. Periods in the state were typically short (mean: LPFCm =16 ms, MCCm =36 ms, not shown). Thus, when immediately preceded and followed by two periods in the same state , periods in were attributed the state . For each network spiking set assessed, we pooled the durations of all periods in all the states of the HMM, to build the overall probability distribution of period durations . We then used this probability distribution to compute
that is, the proportion of time spent in state periods of duration , that is, the probability, at any given instant in time, of being in state periods of duration . We could not find any generic univocal suitable method of stably determining the number of states . However, as a low number of states is more parsimonious in terms of data interpretation (Pohle et al., 2017) in general and because the task structure contains a low number of possible states in terms of actions (four), reward on the last trial (incorrect trial, first correct trial, correct trial after previous correct trials) and behavioural states (exploration, exploitation), we arbitrarily fixed . We checked that this choice was sound and evaluated Akaike information criterion (AIC) and BIC values for different numbers of states in both cortical areas and in both monkey and model data (Figure 7—figure supplement 4). We observed that AIC and BIC did not substantially change when varying the number of states in the range 2–10. We therefore kept the chosen value .
Each HMM analysis was conducted on a spiking set lasting 600 s, both in neural network models and unit populations in monkeys. For each monkey area, the activity of all neurons was pooled, regardless of their recording session. This was mandatory because the number of neurons simultaneously recorded in each session was typically inferior to 5, so that HMMs were inefficient in detecting states. Pooling all neurons allowed the detection of global states that corresponded to the combination of collective dynamics recorded during distinct sessions, that is, that were not time-locked together (phase information lost across sessions) and causally independent. Although chimeric, these HMM states were nevertheless able to indirectly capture the underlying temporal structure of collective spiking discharges in frontal areas in a similar way and thus allowed comparing LPFC and MCC collective temporal structure. In control HMMs, both the timing and neuron assignment of all spikes were randomly shuffled. The initial estimation of the average state duration across all periods in a given state was taken at a high value (300 ms), which was suggested to give better log-likelihood scores and converge to similar states across repetitions of the HMM (Seidemann et al., 1996). The time bin was . We checked that state durations were not caused by eventual drifts in data firing frequency (Figure 7—figure supplement 2; Figure 7—source data 2). We also checked that state durations were distributed exponentially, lending credence to the metastable nature of HMM states (Figure 7—figure supplement 3).
Principal component analysis
Request a detailed protocolThe principal component analysis (PCA) of LPFC and MCC of monkeys’ units and neural network models’ neurons spiking activity was computed from firing frequencies, in order to better visualize and characterize collective dynamics. PCA was achieved on the set of the spiking frequency vectors of all units/neurons in each case. Spiking frequency was estimated through convolution of spiking activity with a normalized Gaussian kernel with standard deviation , as average frequencies were typically in both areas. For each neuron, frequencies were then centred and standardized for optimal PCA. Cells with average frequencies <0.5 Hz were removed for the experimental data and for the model data, to avoid abnormal standardized frequencies when the neuron’s average frequency was too low (at most six cells per area).
Perturbation protocol for state transitions
Request a detailed protocolWe assessed the contribution of excitatory and inhibitory neural populations to the stability of HMM states. To do so, we estimated the probability to stay in a given ongoing (or perturbed, see below) HMM state or to switch towards a distinct target (or perturbing) state in response to specified perturbations. The perturbation was achieved by substituting the value of neural variables (membrane potential, spiking state, calcium concentration, downstream channel opening probabilities) of a random subset of excitatory (respectively inhibitory) neurons of the ongoing state by those of the same neurons taken from the (distinct) target state. Specifically, starting from an initial (unperturbed) 600 s simulation, perturbations were achieved by substituting state variables 50 ms after the onset of a randomly chosen period of a specified perturbed state by those taken 50 ms after the onset of a randomly chosen period of a distinct perturbing state and the resulting network states used as initial conditions for further ‘perturbation simulations’. For each perturbation simulation, the network was simulated from the perturbation time to the end of the period when the network was not perturbed and the HMM state was determined as the posterior state probability based on HMM transition and emission matrices obtained from the entire initial unperturbed simulation. The probability to escape the ongoing state (Figure 7h, left) and to reach the target state (Figure 7h, right) were then computed as the proportion of time spent, during the ongoing period, in a HMM state different from the ongoing perturbed state (escape ongoing state probability), and in the target perturbing state (reach target state probability), respectively. The effects of perturbations were tested by replacing either excitatory or inhibitory populations, where proportions of replaced neurons systematically varied in the range 0–1. For each neuron type and proportion tested, the perturbation protocol was applied and results averaged for 50 random combinations of periods (with period durations >100 ms), for each of the 12 possible pairs of the four HMM states (excluding pairs of repeated states), over 20 different randomly initialized MCCs. Probabilities were offset and normalized to remove the basal probability of escaping the ongoing (0.09) and reaching the target (0.01) states when no perturbation was applied (such transitions were due to random selection of simultaneous spikes when initiating the HMM analysis).
Data availability
All spike time series from monkey recordings, scripts for temporal signatures extraction and scripts of computational models are freely accessible on the Zenodo repository (https://doi.org/10.5281/zenodo.5707884).
-
ZenodoCode for spike analyses and biophysical modelling for 'Inhibitory control of frontal metastability sets the temporal signature of cognition".https://doi.org/10.5281/zenodo.5707884
References
-
Involvement of prefrontal cortex in visual searchExperimental Brain Research 180:289–302.https://doi.org/10.1007/s00221-007-0860-0
-
The neuronal refractory period causes a short-term peak in the autocorrelation functionJournal of Neuroscience Methods 104:155–163.https://doi.org/10.1016/s0165-0270(00)00335-6
-
Learning the value of information in an uncertain worldNature Neuroscience 10:1214–1221.https://doi.org/10.1038/nn1954
-
A reservoir of time constants for memory traces in cortical neuronsNature Neuroscience 14:366–372.https://doi.org/10.1038/nn.2752
-
Dynamics of sparsely connected networks of excitatory and inhibitory spiking neuronsJournal of Computational Neuroscience 8:183–208.https://doi.org/10.1023/a:1008925309027
-
Effects of neuromodulation in a cortical network model of object working memory dominated by recurrent inhibitionJournal of Computational Neuroscience 11:63–85.https://doi.org/10.1023/a:1011204814320
-
Control of the complexity of associative memory dynamics by neuronal adaptationInternational Journal of Neural Systems 4:129–141.https://doi.org/10.1142/s0129065793000122
-
Temporally irregular mnemonic persistent activity in prefrontal neurons of monkeys during a delayed response taskJournal of Neurophysiology 90:3441–3454.https://doi.org/10.1152/jn.00949.2002
-
Pharmacology and Nerve-endings (Walter Ernest Dixon Memorial Lecture): (Section of Therapeutics and Pharmacology)Proceedings of the Royal Society of Medicine 28:319–332.https://doi.org/10.1177/003591573502800330
-
Impact of network activity on the integrative properties of neocortical pyramidal neurons in vivoJournal of Neurophysiology 81:1531–1547.https://doi.org/10.1152/jn.1999.81.4.1531
-
The high-conductance state of neocortical neurons in vivoNature Reviews. Neuroscience 4:739–751.https://doi.org/10.1038/nrn1198
-
Reservoir Computing Properties of Neural Dynamics in Prefrontal CortexPLOS Computational Biology 12:e1004967.https://doi.org/10.1371/journal.pcbi.1004967
-
Selective modulation of cortical state during spatial attentionScience (New York, N.Y.) 354:1140–1144.https://doi.org/10.1126/science.aag1420
-
Functions of SK channels in central neuronsClinical and Experimental Pharmacology & Physiology 34:1077–1083.https://doi.org/10.1111/j.1440-1681.2007.04725.x
-
Plasticity of cortical excitatory-inhibitory balanceAnnual Review of Neuroscience 38:195–219.https://doi.org/10.1146/annurev-neuro-071714-034002
-
Mechanisms of dopamine activation of fast-spiking interneurons that exert inhibition in rat prefrontal cortexJournal of Neurophysiology 88:3150–3166.https://doi.org/10.1152/jn.00335.2002
-
Neuronal L-type calcium channels open quickly and are inhibited slowlyThe Journal of Neuroscience 25:10247–10251.https://doi.org/10.1523/JNEUROSCI.1089-05.2005
-
Inhibitory Plasticity: Balance, Control, and CodependenceAnnual Review of Neuroscience 40:557–579.https://doi.org/10.1146/annurev-neuro-072116-031005
-
Quality metrics to accompany spike sorting of extracellular signalsThe Journal of Neuroscience 31:8699–8705.https://doi.org/10.1523/JNEUROSCI.0971-11.2011
-
Voltage dependence of NMDA-activated macroscopic conductances predicted by single-channel kineticsThe Journal of Neuroscience 10:3178–3182.
-
Network resets in medial prefrontal cortex mark the onset of behavioral uncertaintyScience (New York, N.Y.) 338:135–139.https://doi.org/10.1126/science.1226518
-
Optimal decision making and the anterior cingulate cortexNature Neuroscience 9:940–947.https://doi.org/10.1038/nn1724
-
Neurons in the frontal lobe encode the value of multiple decision variablesJournal of Cognitive Neuroscience 21:1162–1178.https://doi.org/10.1162/jocn.2009.21100
-
Cortical computations via metastable activityCurrent Opinion in Neurobiology 58:37–45.https://doi.org/10.1016/j.conb.2019.06.007
-
Cortical high-density counterstream architecturesScience (New York, N.Y.) 342:1238406.https://doi.org/10.1126/science.1238406
-
Dynamics of multistable states during ongoing and evoked cortical activityThe Journal of Neuroscience 35:8214–8231.https://doi.org/10.1523/JNEUROSCI.4819-14.2015
-
Dynamic population coding of category information in inferior temporal and prefrontal cortexJournal of Neurophysiology 100:1407–1419.https://doi.org/10.1152/jn.90248.2008
-
A hierarchy of intrinsic timescales across primate cortexNature Neuroscience 17:1661–1663.https://doi.org/10.1038/nn.3862
-
A theory of rate coding control by intrinsic plasticity effectsPLOS Computational Biology 8:e1002349.https://doi.org/10.1371/journal.pcbi.1002349
-
Electrophysiological classes of cat primary visual cortical neurons in vivo as revealed by quantitative analysesJournal of Neurophysiology 89:1541–1566.https://doi.org/10.1152/jn.00580.2002
-
Induction and modulation of persistent activity in a layer V PFC microcircuit modelFrontiers in Neural Circuits 7:161.https://doi.org/10.3389/fncir.2013.00161
-
SoftwareR Foundation for Statistical ComputingThe R Project for Statistical Computing, Vienna, Austria.
-
The asynchronous state in cortical circuitsScience (New York, N.Y.) 327:587–590.https://doi.org/10.1126/science.1179850
-
Decoding subjective decisions from orbitofrontal cortexNature Neuroscience 19:973–980.https://doi.org/10.1038/nn.4320
-
Conditional Bistability, a Generic Cellular Mnemonic Mechanism for Robust and Flexible Working Memory ComputationsThe Journal of Neuroscience 38:5209–5219.https://doi.org/10.1523/JNEUROSCI.1992-17.2017
-
Coordination of high gamma activity in anterior cingulate and lateral prefrontal cortical areas during adaptationThe Journal of Neuroscience 31:11110–11117.https://doi.org/10.1523/JNEUROSCI.1016-11.2011
-
Coordinated Prefrontal State Transition Leads Extinction of Reward-Seeking BehaviorsThe Journal of Neuroscience 41:2406–2419.https://doi.org/10.1523/JNEUROSCI.2588-20.2021
-
Temporal filtering of reward signals in the dorsal anterior cingulate cortex during a mixed-strategy gameThe Journal of Neuroscience 27:8366–8377.https://doi.org/10.1523/JNEUROSCI.2369-07.2007
-
“Activity-silent” working memory in prefrontal cortex: A dynamic coding frameworkTrends in Cognitive Sciences 19:394–405.https://doi.org/10.1016/j.tics.2015.05.004
-
Specific frontal neural dynamics contribute to decisions to checkNature Communications 7:11990.https://doi.org/10.1038/ncomms11990
-
Prefrontal cortex HCN1 channels enable intrinsic persistent neural firing and executive memory functionThe Journal of Neuroscience 33:13583–13599.https://doi.org/10.1523/JNEUROSCI.2427-12.2013
-
SKCa channels mediate the medium but not the slow calcium-activated afterhyperpolarization in cortical neuronsThe Journal of Neuroscience 24:3537–3542.https://doi.org/10.1523/JNEUROSCI.0380-04.2004
-
Calcium coding and adaptive temporal computation in cortical pyramidal neuronsJournal of Neurophysiology 79:1549–1566.https://doi.org/10.1152/jn.1998.79.3.1549
-
Macroscopic gradients of synaptic excitation and inhibition in the neocortexNature Reviews. Neuroscience 21:169–178.https://doi.org/10.1038/s41583-020-0262-x
-
Dynamic circuit motifs underlying rhythmic gain control, gating and integrationNature Neuroscience 17:1031–1039.https://doi.org/10.1038/nn.3764
-
Revealing ensemble state transition patterns in multi-electrode neuronal recordings using hidden Markov modelsIEEE Transactions on Neural Systems and Rehabilitation Engineering 19:345–355.https://doi.org/10.1109/TNSRE.2011.2157360
-
Electrophysiological and morphological properties of layers V-VI principal pyramidal cells in rat prefrontal cortex in vitroThe Journal of Neuroscience 16:1904–1921.
-
Multiple Transmitter Receptors in Regions and Layers of the Human Cerebral CortexFrontiers in Neuroanatomy 11:78.https://doi.org/10.3389/fnana.2017.00078
Decision letter
-
Timothy E BehrensSenior and Reviewing Editor; University of Oxford, United Kingdom
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Decision letter after peer review:
Thank you for submitting your article "Inhibitory control of frontal metastability sets the temporal signature of cognition" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Timothy Behrens as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
As the editors have judged that your manuscript is of interest, but as described below
both reviewers found this study's findings to be interesting and novel, and they appreciated the integration of intrinsic timescale analysis, coding of behavioral signals, and exploration of mechanisms in a biophysical circuit model. However, they both raised serious concerns about methodological and interpretational aspects of the analyses which would need to be satisfactorily addressed in subsequent review. In consultation, both reviewers were in agreement with all points raised in the other's review. I view the two reviewers' requests and suggestions as appropriate and complementary, and all of them needed response. Given the quality of the reviewers' comments and their consensus about what needs to be addressed, I am appending both reviewers' full reviews at the end of this letter, to guide the authors' revision and response.
I highlight here two of the concerns raised by the reviewers (but all raised by the reviewers merit addressing). First is the need for intrinsic timescale analysis to not be confounded by slowly varying changes in firing rate coding task variables (Point 1.2 of Reviewer 2). This is important for interpretation on intrinsic timescale and its correlation with task variable coding as a major result. Second is the interpretation and implementation of Hidden Markov Models to non-simultaneously recorded neurons (Points 6 of Reviewer 1 and 3 of Reviewer 2). Here again the HMM states may be driven by task variable coding which is not corrected for, which could confound interpretation of results as in terms meta-stable states and its link to the circuit model. Furthermore, HMM analysis of circuit model does not match its methodology for experimental data, but could be through non-simultaneous model spike trains, which the manuscript does not justify. I will add my own suggestion that perhaps both of these methodological data analysis concerns could potentially be addressed through comparison to null model of a non-homogeneous Poisson process with firing rate given by the variable-coding PSTH. I do not consider it necessary to study a network model that performs the task, as suggested for consideration by Reviewer 1.
Reviewer #1:
This is an interesting manuscript which covers an important topic in the field of computational neuroscience – the 'temporal signatures' of individual neurons. The authors set out to address several important questions using a single-neuron electrophysiology dataset, recorded from monkeys, which has previously been published. The behavioural paradigm is well designed, and particularly well suited to investigating the functional importance of different temporal signatures – as it simultaneously requires the subjects to monitor feedback across a short timescale, as well as integrate multiple outcomes across a longer timescale. The neural data are of high quality, and include recordings from lateral prefrontal cortex (LPFC) and mid-cingulate cortex (MCC). First, the authors modify an existing method to quantify the temporal signatures of individual neurons. This modification appears helpful, and an improvement on similar previously published methods, as the authors are able to capture the temporal signatures of the vast majority of neurons they recorded from. The temporal signatures differ across brain regions, and according to the neurons' spike width. The authors argue that the temporal signatures of a subset of neurons are modified by the subjects' degree of task engagement, and that neurons with different temporal signatures play dissociable roles in task-related encoding. However, I have several concerns about these conclusions which I will outline below. The authors then present a biophysical network model, and show that by varying certain parameters in their model (AHP and GABAB conductances) the temporal signatures of the monkey data can be reproduced. Although I cannot comment on the technical specifics of their models, this seems to be an important advance. Finally, they perform a Hidden-Markov Model analysis to investigate the metastability of activity in MCC, LPFC, and their network model. However, there are a few important differences between the model and experimental data (e.g. neurons recorded asynchronously, and the network model not performing a task) that limit the interpretation of these analyses. Overall, I found the manuscript interesting – and the insights from the biophysical modelling are exciting. However, in its current form, the conclusions drawn from the experimental data are not supported by sufficient evidence.
1. The authors use a hierarchical clustering algorithm to divide neurons into separate groups according to their spike width and amplitude (Figure 1C). There are three groups: FS, RS1, and RS2. The authors ultimately pool RS1 and RS2 groups to form a single 'RS' category. They then go on to suggest that RS neurons may correspond to pyramidal neurons, and FS neurons to interneurons. I have a few concerns about this. Firstly, the suggestion that spike width determined from extracellular recordings in macaques can be used as an indicator of cell type is controversial. A few studies have presented evidence against this idea (e.g. Vigneswaran et al. 2011 JNeuro; Casale et al. 2015 JNeuro). The authors should at least acknowledge the limitation of the inference they are making in the Discussion section. Secondly, visualising the data alone in Figure 1C, it is far from clear that there are three (or two) relatively distinct clusters of neurons to warrant treating them differently in subsequent analyses. In the methods section, the authors mention some analyses they performed to justify the cluster boundaries. However, this data is not presented. A recent study approached this problem by fitting one gaussian to the spike waveform distribution, then performing a model comparison to a 2-gaussian model (Torres-Gomez et al. 2020 Cer Cortex). Including an analysis such as this would provide a stronger justification for their decision to divide cells based on spike waveform.
2. The authors conclude that the results in Figure 3 show that MCC temporal signatures are modulated by current behavioural state. However, this conclusion seems a bit of a stretch from the data currently presented. I can understand why the authors used the 'pause' periods as a proxy for a different behavioural state, but the experiment clearly was not designed for this purpose. As the authors acknowledge, there is only a very limited amount (e.g. a few minutes) of 'pause' data available for the fitting process compared with 'engage' data. Do the authors observe the same results if they constrain the amount of included 'engage' data to match the length of the 'pause' data? Also, presumably the subjects are more likely to 'pause' later on in the behavioural session once they are tired/sated. Could this difference between 'pause' and 'engage' data be responsible for the difference in taus? For instance, there may have been more across-session drift in the electrode position by the time the 'pause' data is acquired, and this could possibly account for the difference with the 'engage' data. Is the firing rate different between 'pause' and 'engage' periods – if so, this should be controlled for as a covariate in the analyses. Finally, it is not really clear to me, or more importantly addressed by the authors, as to why they would expect/explain this effect only being present in MCC RS neurons (but not FS or LPFC neurons).
3. At many points in the manuscript, the authors seem to be suggesting that the results of Figure 4 demonstrate that neurons with longer (shorter) timescales are more involved in encoding the task information which is used across longer (shorter) behavioural timescales (e.g. "long TAU were mostly involved in encoding gauge information", and "population of MCC RS units with short TAU was mostly involved in encoding feedback information"). However, I disagree that this conclusion can be reached based on the way the analysis has currently been performed. A high coefficient simply indicates that the population is biased to be more responsive depending on a particular trial type / condition – i.e. the valence of encoding. This does not necessarily tell us how much information the population of neurons is encoding, as the authors suggest. For instance, every neuron in the population could be extremely selective to a particular parameter (i.e. positive feedback), but if half the neurons encode this attribute by increasing their firing but the other half of neurons encode it by decreasing their firing, the effects will be lost in the authors' regression model (i.e. the β coefficient would equal 0). I would suggest that the authors consider using an alternative analysis method (e.g. a percentage of explained variance or coefficient of partial determination statistic for each neuron) to quantify coding strength – then compare this metric between the high and low tau neurons.
4. Similarly, in Figure 4 the authors suggest that the information is coded differently in the short and long tau neurons. However, they do not perform any statistical test to directly compare these two populations. One option would be to perform a permutation test, where the neurons are randomly allocated into the 'High TAU' or 'Low TAU' group. A similar comment applies to the different groups of neurons qualitatively compared in panel Figure 4C.
5. The authors make an interesting and well-supported case for why changing the AHP and GABA-B parameters in their model may be one mechanism which is sufficient to explain the differences in temporal signatures they observed between MCC and LPFC experimentally. However, I think in places the conclusions they draw from this are overstated (e.g. "This suggests that GABA-B inhibitory – rather than excitatory – transmission is the causal determinant of longer spiking timescales, at least in the LPFC and MCC."). There are many other biophysical differences between different cortical regions – which are not explored in the authors' modelling – which could also account for the differences in their temporal signatures. These could include differences in extra-regional input, the position of the region in an anatomical hierarchy, proportion of excitatory to inhibitory neurons, neurotransmitter receptor/receptor subunit expression, connectivity architecture etc. I think the authors should tone down the conclusions a little, and address some more of these possibilities in more detail in their discussion.
6. For the Hidden Markov Model, I think there are a couple of really important limitations that the authors only touch upon very briefly. Firstly, the authors are performing a population-level analysis on neurons which were not simultaneously recorded during the experiment (only mentioned in the methods). This really affects the interpretation of their results, as presumably the number of states and their duration is greatly influenced by the overall pattern of population activity which the authors are not able to capture. At this stage of the study, I am not sure how the authors can address this point. Secondly, the experimental data is compared to the network model which is not performing any specific task (i.e. without temporal structure). The authors suggest this may be the reason why their predictions for the state durations (Fig7b) are roughly an order of magnitude out. Presumably, the authors could consider designing a network model which could perform the same task (or a simplified version with a similar temporal structure) as the subjects perform. This would be very helpful in helping to relate the experimental data to the model, and may also provide a better understanding of the functional importance of the metastability they have identified in behaviour.
7. It is not clear to me how many neurons the authors included in their dataset, as there appear to be inconsistencies throughout the manuscript (Line 73, Figure 1A-b: MCC = 140, LPFC = 159; line 97-98: MCC = 294, LPFC=276; Figure 2: MCC = 266, LPFC = 258; Methods section line 734-735 and Fig2S2: MCC = 298, LPFC = 272). While this is likely a combination of typos and excluding some neurons from certain analyses, this will need to be resolved. It will be important for the authors to check their analyses, and also add a bit more clarity in the text as to which neurons are being included/excluded in each analysis, and justify this.
Reviewer #2:
The paper investigates the temporal signatures of single-neuron activity (the autocorrelation timescale and latency) in two frontal areas, MCC and LPFC. These signatures differ between the two areas and cell classes, and form an anatomical gradient in MCC. Moreover, the intrinsic timescales of single neurons correspond with their coding of behaviorally relevant information on different timescales. The authors develop a detailed biophysical network model which suggests that after-hyperpolarization potassium and inhibitory GABA-B conductances may underpin the potential biophysical mechanism that explains diverse temporal signatures observed in the data. The results appear exciting, as the proposed relationship between the intrinsic timescales, coding of behavioral timescales, and anatomical properties (e.g., the amount of local inhibition) in the two frontal areas is novel. The use of the biophysically detailed model is creative and interesting. However, there are serious methodological concerns undermining the key conclusions of this study, which need to be addressed before the results can be credited.
1. One of the key findings is the correspondence between the intrinsic timescales of single neurons and their coding of information on different behavioral timescales (Figure 4). However, the method for estimating the intrinsic timescales has serious problems which can undermine the finding.
1.1. The authors developed a new method for estimating autocorrelograms from spike data but the details of this method are not specified. It is stated that the method computes the distribution of inter-spike-intervals (ISIs) up to order 100, which was "normalized", but how it was normalized is not described. The correct normalization is crucial, as it converts the counts of spike coincidences (ISI distribution) into autocorrelogram (where the coincidence counts expected by chance are subtracted) and can produce artifacts if not performed correctly.
1.2. The new method, described as superior to the previous method by Murray et al., 2014, appears to have access to more spikes than the Murray's method (Figure 2). Where is this additional data coming from? While Murray's method was applied to the pre-cue period, the time epoch used for the analysis with the new method is not stated clearly. It seems that the new method was applied to the data through the entire trial duration and across all trials, hence more spikes were available. If so, then changes in firing rates related to behavioral events contribute to the autocorrelation, if not appropriately removed. For example, the Murray's method subtracts trail-averaged activity (PSTH) from spike-counts, similar to shuffle-correction methods. If a similar correction was not part of the new method, then changes in firing rates due to coding of task variables will appear in the autocorrelogram and estimated timescales. This is a serious confound for interpretation of the results in Figure 4. For example, if the firing rate of a neuron varies slowly coding for the gauge size across trials, this will appear as a slow timescale if the autocorrelogram was not corrected to remove these rate changes. In this case, the timescale and GLM are just different metrics for the same rate changes, and the correspondence between them is expected. Before results in Figure 4 can be interpreted, details of the method need to be provided to make sure that the method measures intrinsic timescales, and not timescales of rate changes triggered by the task events. This is an important concern also because recent work showed that there is no correlation between task dependent and intrinsic timescales of single neurons, including in cingulate cortex and PFC (Spitmaan et al., PNAS, 2020).
2. The balanced network model with a variety of biophysical currents is interesting and it is impressive that the model reproduces the autocorrelation signatures in the data. However, we need to better understand the network mechanism by which the model operates.
2.1. The classical balanced network (without biophysical currents such as after-hyperpolarization potassium) generates asynchronous activity without temporal correlations (Renart et al., Science, 2010). The balanced networks with slow adaptation currents can generate persistent Up and Down states that produce correlations on slow timescales (Jercog et al., eLife, 2017). Since slow after-polarization potassium current was identified as a key ingredient, is the mechanism in the model similar to the one generating Up and Down states, or is it different? Although the biophysical ingredients necessary to match the data were identified, the network mechanism has not been studied. Describing this network mechanism and presenting the model in the context of existing literature is necessary, otherwise the results are difficult to interpret for the reader.
2.2. Does the model operate in a physiologically relevant regime where the firing rates, Fano factor etc. are similar to the data? It is hard to judge from Figure 5b and needs to be quantified.
2.3. The latency of autocorrelation is an interesting feature in the data. Since the model replicates this feature (which is not intuitive), it is important to know what mechanism in the model generates autocorrelation latency.
3. HMM analysis is used to demonstrate metastability in the model and data, but there are some technical concerns that can undermine these conclusions.
3.1. HMM with 4 states was fitted to the data and model. The ability to fit a four-state HMM to the data does not prove the existence of metastable states. HMM assumes a constant firing rate in each "state", and any deviation from this assumption is modeled as state transitions. For example, if some neurons gradually increase/decrease their firing rates over time, then HMM would generate a sequence of states with progressively higher/lower firing rates to capture this ramping activity. In addition, metastability implies exponential distributions of state durations, which was not verified. No model selection was performed to determine the necessary number of states. Therefore, the claims of metastable dynamics are not supported by the presented analysis.
3.2. HMM was fit to a continuous segment of data lasting 600s, and the data was pooled across different recording sessions. However, different sessions have potentially different trial sequences due to the flexibility of the task. How were different trial types matched across the sessions? If trial-types were not matched/aligned in time, then the states inferred by the HMM may trivially reflect a concatenation of different trial types in different sessions. For example, the same time point can correspond to the gauge onset in one session and to the work trial in another session, and vice versa at a different time. If some neurons respond to the gauge and others to the work, then the HMM would need different states to capture firing patterns arising solely from concatenating the neural responses in this way. This confound needs to be addressed before the results can be interpreted.
https://doi.org/10.7554/eLife.63795.sa1Author response
I highlight here two of the concerns raised by the reviewers (but all raised by the reviewers merit addressing).
First is the need for intrinsic timescale analysis to not be confounded by slowly varying changes in firing rate coding task variables (Point 1.2 of Reviewer 2). This is important for interpretation on intrinsic timescale and its correlation with task variable coding as a major result.
We agree with the reviewers that the relationship between firing rate and timescale variation is not trivial and could constitute an important confound. We tackled this concern in the initial submission (Figure 2 —figure supplement 2) by showing that there is no linear correlation between the firing rate and TAUs. However, and as pointed out by the reviewers, this relationship depends on how the firing rate is calculated. In the analysis reported in the initial manuscript, we aimed to extract the firing rate when the units are the most active as we thought it would be a good measure of the physiological firing rate and this would avoid average issues. In response to the reviewers, we decided to go further into this analysis and we described TAUs and firing rate relationships between different recording periods, when monkeys are engaged, taking breaks or around key task events. We found little correlation around key task events (feedback and decision). These however were not reliable as segments used for data were short (2s for feedback, 4s for decision) resulting in noisy ISI computations leading to less spike-autocorrelogram being fitted and more variability in TAUs.
Conversely there were no limitations in firing rate computation (number of spikes given a duration) for larger segments (data during pauses or engage periods, or the whole data set). Using such large segments, we showed that there was no linear correlation between the firing rates computed at feedback or during decisions and TAUs extracted from the whole recording:
Consistent with the idea that TAUs measures are stable and descriptive of a unit’s properties, we found that TAUs extracted from different periods were quite correlated with each other, when data allowed for such correlations. Also, as pointed out in the manuscript, TAUs extracted from the spike autocorrelogram and from the spike-count method were correlated across units for which both measures were available.
The spike-count method to measure TAUs requires a repetition of task trials from which an autocorrelation can be computed between each time bin. In their supplementary materials, Murray and colleagues (Murray et al. 2014) demonstrated that the autocorrelation subtracts the mean spike-count of each time bin, thus controlling for firing rate fluctuation. Nevertheless, this does not remove within trials fluctuations of firing rate. In fact in our dataset we found a slight correlation between scTAUs and firing rate extracted from the precue period (p-value = 0.031, rho=0.14):
Reviewer #2 suggested to remove the average firing rates (across trials) by normalization when extracting TAUs in the spike autocorrelogram. This however might not be relevant as the spike autocorrelogram method is not trial-based but based on the ISIs distribution. Some confusion might arise from the term autocorrelation which refers to a true mathematical autocorrelation in the spikecount method, whereas it refers to the construction of spike autocorrelograms from the ISIs distribution in our method. We clarify the description of the method in the Result section.In conclusion, firing rate in our analyses appears to not be a confound for the analyses performed in this paper. We provide other more specific responses to reviewers below.
Second is the interpretation and implementation of Hidden Markov Models to non-simultaneously recorded neurons (Points 6 of Reviewer 1 and 3 of Reviewer 2). Here again the HMM states may be driven by task variable coding which is not corrected for, which could confound interpretation of results as in terms meta-stable states and its link to the circuit model. Furthermore, HMM analysis of circuit model does not match its methodology for experimental data, but could be through non-simultaneous model spike trains, which the manuscript does not justify. I will add my own suggestion that perhaps both of these methodological data analysis concerns could potentially be addressed through comparison to null model of a non-homogeneous Poisson process with firing rate given by the variable-coding PSTH. I do not consider it necessary to study a network model that performs the task, as suggested for consideration by Reviewer 1.
We have followed both of your propositions, regarding non-simultaneous model spike trains (i.e. built and run HMM analyses of pseudo-populations on model data) and task variable-coding (i.e. using non-homogeneous Poisson with firing rate given by the task variable-coding PSTH) to disambiguate these possible issues. These additional analyses clearly show that (1) pseudo-populations in experimental data vs populations in the model is not an issue in drawing conclusions on the fact that neural states display larger durations in the MCC, compared to the LPFC in both monkey and model activity data, and (2) task variable coding is not responsible for longer meta-stables states in the MCC (compared to the LPFC) in monkey data (see Specific responses to reviewers).
Reviewer #1:
This is an interesting manuscript which covers an important topic in the field of computational neuroscience – the 'temporal signatures' of individual neurons. The authors set out to address several important questions using a single-neuron electrophysiology dataset, recorded from monkeys, which has previously been published. The behavioural paradigm is well designed, and particularly well suited to investigating the functional importance of different temporal signatures – as it simultaneously requires the subjects to monitor feedback across a short timescale, as well as integrate multiple outcomes across a longer timescale. The neural data are of high quality, and include recordings from lateral prefrontal cortex (LPFC) and mid-cingulate cortex (MCC). First, the authors modify an existing method to quantify the temporal signatures of individual neurons. This modification appears helpful, and an improvement on similar previously published methods, as the authors are able to capture the temporal signatures of the vast majority of neurons they recorded from. The temporal signatures differ across brain regions, and according to the neurons' spike width. The authors argue that the temporal signatures of a subset of neurons are modified by the subjects' degree of task engagement, and that neurons with different temporal signatures play dissociable roles in task-related encoding. However, I have several concerns about these conclusions which I will outline below. The authors then present a biophysical network model, and show that by varying certain parameters in their model (AHP and GABAB conductances) the temporal signatures of the monkey data can be reproduced. Although I cannot comment on the technical specifics of their models, this seems to be an important advance. Finally, they perform a Hidden-Markov Model analysis to investigate the metastability of activity in MCC, LPFC, and their network model. However, there are a few important differences between the model and experimental data (e.g. neurons recorded asynchronously, and the network model not performing a task) that limit the interpretation of these analyses. Overall, I found the manuscript interesting – and the insights from the biophysical modelling are exciting. However, in its current form, the conclusions drawn from the experimental data are not supported by sufficient evidence.
We thank the reviewer for their clear appreciation of the gains provided by our study. We would like first to note that the dataset presented in the current study is an extension of the previous one published in Stoll et al. 2016, containing additional recordings from the same animals.
1. The authors use a hierarchical clustering algorithm to divide neurons into separate groups according to their spike width and amplitude (Figure 1C). There are three groups: FS, RS1, and RS2. The authors ultimately pool RS1 and RS2 groups to form a single 'RS' category. They then go on to suggest that RS neurons may correspond to pyramidal neurons, and FS neurons to interneurons. I have a few concerns about this. Firstly, the suggestion that spike width determined from extracellular recordings in macaques can be used as an indicator of cell type is controversial. A few studies have presented evidence against this idea (e.g. Vigneswaran et al. 2011 JNeuro; Casale et al. 2015 JNeuro). The authors should at least acknowledge the limitation of the inference they are making in the Discussion section.
We thank the reviewer for this comment which indeed raised an important point. We acknowledge that there is no one to one relationship between spike shape and cell type and that although several studies presented data in that sense, others, sometimes in different structures, showed that in some instances Fast Spiking (FS) units could be related to pyramidal cells. Also, several excitatory and inhibitory cell types can be found in the primate prefrontal cortex, with possibly different functional roles (Wang Nat Neursc Rev 2020), making the picture more complicated than just a binary dissociation. Yet, Trainito and colleagues found that (from a database of 2500 waveforms) several waveform clusters (4) could be dissociated and linked to different coding properties (Trainito et al. 2019 Curr Biol.) confirming that, using one of the only accessible cell specific features from extracellular recordings, one might extract some basic biological characteristics. Our more limited dataset revealed 3 clusters. Their data and the literature suggest that ‘spike waveform is a sufficiently sensitive and specific marker to dissociate more than two cell classes from extracellular recordings’ (Trainito et al. 2019). Other past studies have succeeded in separating FS and RS units using waveforms (in particular spike width) and observing clustered physiological properties (Nowak et al. 2003 J Neurophysiol, Krimer et al. 2005 J Neurophysiol). Combining spike width to other physiological measures obtained in vitro or with intracellular recordings are without doubt much better for isolating FS and RS. Krimer et al. (2005) found that 100% of pyramidal neurons were RS and 100% of chandelier cells (gabaergic interneurons) were FS, but some other intermediate interneuron types were found to be either RS or FS. Unfortunately, one can’t yet access such measures in behaving monkeys. Although spike waveform is a quite good parameter to separate unit types, it is not perfect. In our dataset a few elements suggest this dissociation might be relevant. First the proportion of FS and RS units observed in our database is 20/80 % which fits relatively well with the estimated proportions of interneurons vs excitatory neurons. Also, as mentioned in the litterature, the firing rate is slightly higher on average for interneurons than for pyramidal neurons, as in our models, and as observed for FS and RS units.
We now acknowledge the inherent pitfalls but also the potential of spike waveform clustering in the text and provide references according to the reviewer’s point (see Results, page 4). We also reformulated the text to clarify that we only hypothesized such equivalence for modeling purposes knowing this is a simplification, as in our spike segregation and modeling only 2 cell types were considered.
Secondly, visualising the data alone in Figure 1C, it is far from clear that there are three (or two) relatively distinct clusters of neurons to warrant treating them differently in subsequent analyses. In the methods section, the authors mention some analyses they performed to justify the cluster boundaries. However, this data is not presented. A recent study approached this problem by fitting one gaussian to the spike waveform distribution, then performing a model comparison to a 2-gaussian model (Torres-Gomez et al. 2020 Cer Cortex). Including an analysis such as this would provide a stronger justification for their decision to divide cells based on spike waveform.
We agree with the reviewers on the confusion about the clustering methods and the lack of justification for the number of clusters to be formed. For the sake of clarity, we decided to remove from the methods section the mention of analyses (PCA procedure) that were not shown as they were redundant with the hierarchical clustering. We now use the method from Torres-Gomez et al. 2020 Cereb Cortex suggested by the reviewer in order to obtain an objective statistical number of clusters formed from the spike widths distribution. We discuss this further below and have modified the methods section accordingly (see Spike Shape Clustering).
Furthermore, we decided not to use the firing rate for the clustering because we did not have a clear justification for choosing a specific period for firing rate calculation. Yet we found that the so-called FS population we extracted had a higher firing rate than the RS population (Figure 2 —figure supplement 2a). This difference is in adequation with the literature and supported our decision to cluster units solely based on spike shape/duration. This difference of firing rate is reported when computed from the whole recording. Firing rate computed from different periods of the recordings (when monkeys are engaged in the task, taking a break or around key task events…) are correlated but variable. Actually, in our recordings the correlation is lowest when considering firing rates between pauses and the foreperiod of the task, two of the periods which could have been logical candidates for a firing rate of reference.
We finally clustered units according to their spike width and V2/PiK. We first computed the spike width vs. V2/PiK Euclidean distance matrix (dist function in R). Then we performed hierarchical clustering using Ward’s method (hclust function in R). The number of retained clusters was statistically confirmed by fitting one-dimensional gaussian mixture models from 1 to 3 models with variable/unequal variance on the spike widths distribution using the function Mclust from the package Mclust in R which uses the expectation-maximization algorithm. The model that best fitted the spike widths distribution was composed of 3 gaussians (as shown in figure 1c).
2. The authors conclude that the results in Figure 3 show that MCC temporal signatures are modulated by current behavioural state. However, this conclusion seems a bit of a stretch from the data currently presented. I can understand why the authors used the 'pause' periods as a proxy for a different behavioural state, but the experiment clearly was not designed for this purpose. As the authors acknowledge, there is only a very limited amount (e.g. a few minutes) of 'pause' data available for the fitting process compared with 'engage' data. Do the authors observe the same results if they constrain the amount of included 'engage' data to match the length of the 'pause' data? Also, presumably the subjects are more likely to 'pause' later on in the behavioural session once they are tired/sated. Could this difference between 'pause' and 'engage' data be responsible for the difference in taus? For instance, there may have been more across-session drift in the electrode position by the time the 'pause' data is acquired, and this could possibly account for the difference with the 'engage' data. Is the firing rate different between 'pause' and 'engage' periods – if so, this should be controlled for as a covariate in the analyses. Finally, it is not really clear to me, or more importantly addressed by the authors, as to why they would expect/explain this effect only being present in MCC RS neurons (but not FS or LPFC neurons).
Indeed, pauses are a natural phenomena that reflects the animal's voluntary decision to engage or not. The distribution of pauses is indeed biased toward the second half of a session. We now provide a figure for both monkeys (Figure 3b). The across session drift is possible but unlikely to explain everything as we checked for each session the stability of spike shape and the approximate firing rate for each unit. A time on task effect is also possible and this is indeed a very interesting suggestion. Other data from our lab revealed time-on-task changes in frontal neural activity at the level of β oscillation power (Stoll et al. 2016) and at the level of LFP and unit activity (Goussi et al., in preparation).
To control for a time-on-task effect on timescale modulation, we now contrast pause periods with engaged periods that occurred at similar times within sessions (i.e. considering only engaged periods occuring after the first pause, so excluding the first segment of data where no pause occurred – See limits in figure 3b, yellow marks; see changes in Methods – Behaviour and context-dependant modulations).
Because engage periods were as long as pause periods for one monkey (Monkey H, 53 sessions, MDengage=392s, MDpause=396s, Wilcoxon-paired test: V=790, p=0.51) and roughly twice as long for the other (Monkey A, 96 sessions, MDengage=638s, MDpause=372s, Wilcoxon-paired test: V=406, p=2.19e-12), we decided not to further segment the data to avoid resampling biases. The contrast between TAU engage and pause led to results similar to those already reported in the article, with increased TAUs in engage relative to pause periods only in MCC RS population, thus discarding a confounding effect of time on task (nMCCFS=19, nLPFCFS=21, nMCCRS=80, nLPFCRS=97, Wilcoxon signed-ranks test with Bonferroni correction, only significant for MCC RS: MD=1.06, V=2467, p=3.9e-7). Because this analysis is based on a subpart of recordings, the number of units with correctly extracted timescale for those periods is not sufficient to conduct appropriate statistical tests separately for the two monkeys (several groups have n<10). We have however verified that the distribution of TAUs extracted from the whole recordings are similar for the two monkeys, with longer TAUs in the MCC (Linear model on BLOM transformed TAU, region x unit type, Monkey A (332 units), interaction: ns, unit type: t=1,89, p=0,059, region: t=-3,25, p=0,0013 ; Monkey H (189 units), interaction: ns, unit type: ns, region: t=-2,60, p=0,010).
The reviewer finally asked whether one would have an interpretation on why such effect engage/pause is observed for RS units in MCC and not for FS or for LPFC units. We do not have a definitive answer but the differentiation might be related to multiple factors including the large range of timescales available for RS units in MCC (compared to FS or LPFC RS units) which might be regulated by specific mechanisms including neuromodulation by noradrenaline or dopamine. The differential density of neuromodulatory fibers in MCC versus LPFC might be an influential factor. Finally, the different data for LPFC compared to MCC might also relate to the areas covered by recordings which include the most posterior parts of the LFPC (described in Stoll et al. 2016).
3. At many points in the manuscript, the authors seem to be suggesting that the results of Figure 4 demonstrate that neurons with longer (shorter) timescales are more involved in encoding the task information which is used across longer (shorter) behavioural timescales (e.g. "long TAU were mostly involved in encoding gauge information", and "population of MCC RS units with short TAU was mostly involved in encoding feedback information"). However, I disagree that this conclusion can be reached based on the way the analysis has currently been performed. A high coefficient simply indicates that the population is biased to be more responsive depending on a particular trial type / condition – i.e. the valence of encoding. This does not necessarily tell us how much information the population of neurons is encoding, as the authors suggest. For instance, every neuron in the population could be extremely selective to a particular parameter (i.e. positive feedback), but if half the neurons encode this attribute by increasing their firing but the other half of neurons encode it by decreasing their firing, the effects will be lost in the authors' regression model (i.e. the β coefficient would equal 0). I would suggest that the authors consider using an alternative analysis method (e.g. a percentage of explained variance or coefficient of partial determination statistic for each neuron) to quantify coding strength – then compare this metric between the high and low tau neurons.
We now supplement the population coding analysis with a more common single-unit coding analysis as reported in Stoll et al. 2016. The detail is in the Materials and methods section ‘Task-related analyses – Spike analyses’ and results are in figure 4 supplement 1. As suggested by the reviewer, we used the weighted proportion of variance explained (wPEV) as a statistical measure to quantify the extent to which the variability in neural firing rate was determined by feedback valence and gauge state. We then quantified the time-resolved proportion of cells coding for the task variables (p<0.05 for at least 4 successive bins) and the wPEV in populations of units with short or long TAU (median split by cell type). To assess differences between short and long TAU populations we built null distributions by permuting 1000 times the TAU group allocation of units. We then compared the position of the true data relative to the cumulative distribution of the permutations and set a statistical threshold to a=0.05. The ‘single-unit’ analyses overall revealed similar dynamics for RS long TAU units in MCC and for FS involvement in Feedback than the population coding shown in the main Figure 4 (see Figure 4 supplement 1). It thus suggests that the glmm population analysis conserved some freedom (thanks to Unit as a Random term) to fit the different changes in firing rates. The single-unit analysis however, while revealing (as for the glmm) a higher contribution of long TAU RS MCC units at feedback onset, did not show a bias for short TAU RS for feedback encoding during inter-trial. Thus the two approaches provided complementary outcomes that we now summarize in the Results section (see Result section Temporal signatures are linked to cognitive processing):
“We used data from the whole recordings (all periods) and classed both FS and RS units as either short or long TAU units using a median split. We used a time-resolved generalized mixed linear models (glmm; Figure 4a) to reveal notable dissociations between these populations that we complemented using a more classical approach at the single unit level, using Poisson glm and weighted proportion of variance explained (wPEV; Figure 4 – supplement 1).
Early phases of feedback encoding recruited MCC long TAU populations for both FS and RS units (Figure 4a, upper graphs). This discrepancy was confirmed by a difference between early feedback coding in short and long TAU population at the single unit level (Figure 4 —figure supplement 1). Interestingly FS units in the MCC were mostly engaged in the first second after feedback onset, with a strong bias toward encoding negative feedback (Figure 4a, upper left, positive estimates). Effects were more transient and involved short TAU units in the LPFC (Figure 4a).
During the intertrial interval, feedback valence was represented in different directions between short and long TAU RS populations in MCC, coding being positive for short TAU populations (higher activity for incorrect feedback) and becoming negative for long TAU populations (higher activity for correct feedback). Conversely, solely the population of MCC long TAU RS units coded for the gauge during the intertrial period (Figure 4a, lower graphs). Single unit analyses confirmed the higher contribution of the long TAU population to gauge encoding (Figure 4 – supplement 1).”
Altogether, the previous and new analyses point to differentiations between short/long neuronal timescales and the encoding of behavioral information.
4. Similarly, in Figure 4 the authors suggest that the information is coded differently in the short and long tau neurons. However, they do not perform any statistical test to directly compare these two populations. One option would be to perform a permutation test, where the neurons are randomly allocated into the 'High TAU' or 'Low TAU' group. A similar comment applies to the different groups of neurons qualitatively compared in panel Figure 4C.
We thank the reviewer for this suggestion. To assess how short and long TAU population codings for a given area and cell type differ from null, we have now constructed null distributions of coding (zvalues) by permuting TAU group allocation of units. We performed the permutations for each area and cell type 100 times, this allowed us to conserve differences in sampling (e.g. the population of MCC RS with long TAU is larger than the short TAU one).
We then compare the position of the true data relative to the cumulative distribution of the permutations and set a statistical threshold a=0.05. This is shown now in panel (a) in figure 4, significant statistics being represented by rasters below each curve.
The permutation tests revealed results complementary to the more direct glmm modeling and to wPEV analyses. Glmm and wPEV show that MCC FS cells encode feedback early after its onset, and earlier for long TAU FS than short ones. Note however that the glmm results do differ from the null distributions. The effects of gauge encoding by MCC RS cells is observed for long TAU compared to short Tau cells in a period just before lever onset (glmm and wPEV). Comparing glmm effects to null distribution reveal effects outside of this period.
5. The authors make an interesting and well-supported case for why changing the AHP and GABA-B parameters in their model may be one mechanism which is sufficient to explain the differences in temporal signatures they observed between MCC and LPFC experimentally. However, I think in places the conclusions they draw from this are overstated (e.g. "This suggests that GABA-B inhibitory – rather than excitatory – transmission is the causal determinant of longer spiking timescales, at least in the LPFC and MCC."). There are many other biophysical differences between different cortical regions – which are not explored in the authors' modelling – which could also account for the differences in their temporal signatures. These could include differences in extra-regional input, the position of the region in an anatomical hierarchy, proportion of excitatory to inhibitory neurons, neurotransmitter receptor/receptor subunit expression, connectivity architecture etc. I think the authors should tone down the conclusions a little, and address some more of these possibilities in more detail in their discussion.
We thank the reviewer for pointing out additional mechanisms which could play a role in the differences in temporal signatures. Regarding excitation transmission parameters, exhaustive 1D and 2D explorations (i.e. AMPA and NMDA maximal conductance and time constants) did not allow to find a space region accounting for network activity with realistic statistics and the whole pattern of temporal signatures (timescale and latency of the autocorrelogram in the LPFC and MCC). This was the same for all other parameters or combinations of parameters tested, except for the GABA-B / AHP plane, where it was possible to find such regions in a robust manner. This contrast was striking in our explorations. Nevertheless, we have toned down our conclusion (“Local molecular basis of frontal temporal signatures’ subsection) and adopted a more nuanced formulation. We have also introduced all of the other dimensions rightfully proposed (extra-regional input, position in the hierarchy, etc.; “Local molecular basis of frontal temporal signatures” subsection).
6. For the Hidden Markov Model, I think there are a couple of really important limitations that the authors only touch upon very briefly. Firstly, the authors are performing a population-level analysis on neurons which were not simultaneously recorded during the experiment (only mentioned in the methods). This really affects the interpretation of their results, as presumably the number of states and their duration is greatly influenced by the overall pattern of population activity which the authors are not able to capture. At this stage of the study, I am not sure how the authors can address this point.
We applied Hidden Markov Models (HMMs) to the activity of populations of neurons that were not recorded simultaneously (i.e. pseudo-populations) to analyze overall LPFC and MCC dynamics in terms of neural states for three reasons. First, HMMs were impractical when limited to the groups of simultaneously recorded monkey neurons, as that contained too few neurons (3 or 4 in general, 5 at maximum); therefore, working on pseudo-populations was the best way to improve the amount of data employed and to maximize the information extracted in terms of neural states. Second, whereas working on pseudo-populations degrades quantitative estimations of states (by ignoring a part of correlated activity between neurons), this effect in principle affects absolute quantitative state estimations in the LPFC and MCC equally and likely does not affect the conclusions drawn here, as they are based on relative comparisons between order of magnitudes of state estimates in the LPFC and MCC. Third, previous studies have shown that the amount of information that can be time-decoded in associative areas (such as in the LPFC) is globally similar when using pseudo-populations or simultaneously recorded neurons (Anderson et al., Exp Brain Res 2007; Meyers et al., J Neurophysiol. 2008). We acknowledge that these considerations were absent in the original manuscript and have now introduced them in the Materials and methods (“Hidden Markov Model (HMM) analysis” subsection).
In our results, a similar trend was observed in both monkey and model neural activity: the duration of states detected by HMM was larger in the MCC, compared to the LPFC, with a MCC/LPFC ratio of mean state durations of ~30 in model populations and of ~7.5 in monkey pseudo-populations. To better interpret and strengthen our results, with respect to the question raised by reviewer’s #1, we have now additionally applied HMM analysis to neural activity in the model, using pseudo-populations built of groups of neurons drawn from LPFC and MCC simulations of distinct random neural networks, with group statistics (number of recording sessions and numbers of simultaneously recorded excitatory and inhibitory neurons within each recording session) similar to those of monkey LPFC and MCC data (Figure 7 – supplement 1). The rationale here is to examine whether a higher state duration MCC/LPFC ratio is still present, i.e. assess whether the model still predicts and accounts for higher state durations in monkey MCC pseudo-populations, when using model pseudo-population activity.
We found that, although the MCC/LPFC ratio is decreased when switching to pseudopopulation activity, state duration is still higher in the MCC, compared to the LPFC, with a ratio of ~1.5. In order to better understand the origin of this difference, we performed controls, where we randomly shuffled spikes according to the two orthogonal dimensions of spiking data, i.e. neurons and time, separately from one another. Shuffling spikes by randomly reattributing each individual spike of each neuron to a randomly chosen neuron (neuron shuffle control, repeated 10 times) also gave a similar ratio of ~1.92, showing that the neural identity of spikes did not account for the MCC/LPFC state duration ratio. However, randomly shuffling the timing of individual spikes within each neuron’s activity, or the time shuffle control, led to a ratio of ~1. The combination of both results indicated the MCC/LPFC difference originated from the global temporal structure of the collective network activity (rather than per-neuron temporal structure). To ensure the global temporal structure truly corresponded to spike timing rather than other trivial confounding variables, we controlled for variables we thought could affect the global temporal structure other than spike timing, e.g. the number or the firing frequency of LPFC/MCC neurons. Using the same number of neurons of each type for both LPFC and MCC model pseudo-populations areas (those of the LPFC) gave a ratio of ~1.52, while removing randomly chosen spikes from LPFC activity to lower its mean firing frequency to a value similar to that of the MCC gave a ratio of ~2.03. These controls evidenced that the ~1.5 MCC/LPFC ratio emerged from, and represented, a robust trait of model LPFC and MCC activity temporal structure and not a side effect of possible confounding factors (frequency, neural discharge identity, group statistics). Furthermore, the MCC/LPFC difference is structurally stable, showing strong robustness even to control procedures which strongly degrade model data.
Altogether, these results indicate that, while information is degraded when pseudo-population is used, a part of the information of neural activity’s temporal structure is preserved, resulting in a similar – although decreased – relative state duration ratio between MCC and LPFC. A fortiori, this additional analysis strengthens our results, as it indicates that MCC to LPFC state duration ratios should be higher in real monkey neural populations than the 7.5 ratio found here in pseudo-populations activity (given the 30 vs 1.5 ratios obtained with population vs pseudo-population model activity, respectively). In conclusion, these results indicate that using pseudo-populations in experimental data vs populations in the model is not an issue in drawing conclusions on the fact that neural states display relatively larger durations in the MCC, compared to the LPFC, in both monkey and model activity.
We have now added these extended results as supplementary materials (Figure 7 —figure supplement 1a) to better support our conclusions. We also explicitly mention in the Results that HMM are applied to model populations and monkey pseudo-populations (“Metastable states underlie LPFC and MCC temporal signatures” subsection), and point to specific explanations of this question in Materials and methods (“Hidden Markov Model (HMM) analysis” subsection) and Supplementary Materials (Figure 7 —figure supplement 1). We also mention this issue in the Discussion (“Frontal temporal signatures uncover metastable dynamics” subsection).
Secondly, the experimental data is compared to the network model which is not performing any specific task (i.e. without temporal structure). The authors suggest this may be the reason why their predictions for the state durations (Fig7b) are roughly an order of magnitude out. Presumably, the authors could consider designing a network model which could perform the same task (or a simplified version with a similar temporal structure) as the subjects perform. This would be very helpful in helping to relate the experimental data to the model, and may also provide a better understanding of the functional importance of the metastability they have identified in behaviour.
Building a biophysical neural network realizing the monkey behavioral task, even employing engineer-based learning plasticity rules, is out of scope of the present study. However, to address the concern raised, we devised a specific version of the HMM analysis to assess the possible role of the temporal structure of the task in the difference observed in the MCC compared to the LPFC in experimental data. Specifically, we simulated neural spiking activity of each neuron in LPFC and MCC areas as a non-homogeneous Poisson model (NHPM) whose instantaneous firing rate was given by the post-stimulus time histogram (PSTH) of that neuron upon each epoch type across all possible task events (e.g., reward, go signal, etc.). This simulated activity was only dependent on the coding of task variable averaged across trials, and did not contain more specific fine-grained activity temporal structure within each trial. Applying HMM analysis to such transformed population activity showed no difference of state durations between LPFC and MCC. This indicated that task variable coding was not responsible for the state duration difference between LPFC and MCC but relied on the precise temporal structure of spiking, independent of task epochs. This also indicated that accounting for neural states in monkeys with a biophysical model devoid of task structure was not an issue per se. We have now added this result as Supplementary Materials (Figure 7 —figure supplement 1b) to show that task coding is not essential in understanding the difference in state duration between MCC and LPFC in monkeys and can be accounted for by a model devoid of task temporal structure.
7. It is not clear to me how many neurons the authors included in their dataset, as there appear to be inconsistencies throughout the manuscript (Line 73, Figure 1A-b: MCC = 140, LPFC = 159; line 97-98: MCC = 294, LPFC=276; Figure 2: MCC = 266, LPFC = 258; Methods section line 734-735 and Fig2S2: MCC = 298, LPFC = 272). While this is likely a combination of typos and excluding some neurons from certain analyses, this will need to be resolved. It will be important for the authors to check their analyses, and also add a bit more clarity in the text as to which neurons are being included/excluded in each analysis, and justify this.
We thank the reviewer for this remark, as the reported numbers were indeed confusing. This is because neurons might not meet the criteria used to ensure that analyses contained enough data when performing specific tests, like limiting the extraction of Taus and latencies on specific behavioral periods. We have checked every sample size mentioned in the text and tried to better explain why and when those changes occurred.
Reviewer #2:
The paper investigates the temporal signatures of single-neuron activity (the autocorrelation timescale and latency) in two frontal areas, MCC and LPFC. These signatures differ between the two areas and cell classes, and form an anatomical gradient in MCC. Moreover, the intrinsic timescales of single neurons correspond with their coding of behaviorally relevant information on different timescales. The authors develop a detailed biophysical network model which suggests that after-hyperpolarization potassium and inhibitory GABA-B conductances may underpin the potential biophysical mechanism that explains diverse temporal signatures observed in the data. The results appear exciting, as the proposed relationship between the intrinsic timescales, coding of behavioral timescales, and anatomical properties (e.g., the amount of local inhibition) in the two frontal areas is novel. The use of the biophysically detailed model is creative and interesting. However, there are serious methodological concerns undermining the key conclusions of this study, which need to be addressed before the results can be credited.
We thank the reviewer for their positive appreciation of our study. We address below all the methodological concerns.
1. One of the key findings is the correspondence between the intrinsic timescales of single neurons and their coding of information on different behavioral timescales (Figure 4). However, the method for estimating the intrinsic timescales has serious problems which can undermine the finding.
1.1. The authors developed a new method for estimating autocorrelograms from spike data but the details of this method are not specified. It is stated that the method computes the distribution of inter-spike-intervals (ISIs) up to order 100, which was "normalized", but how it was normalized is not described. The correct normalization is crucial, as it converts the counts of spike coincidences (ISI distribution) into autocorrelogram (where the coincidence counts expected by chance are subtracted) and can produce artifacts if not performed correctly.
We thank the reviewer for pointing out the lack of clarity. We now make a more explicit and detailed explanation of the normalization procedures and the methods for spike autocorrelation in general.
We normalized by calculating the probability distribution function. This doesn’t affect the value of TAU exponentially fitted. In any case, there is a confusion between our measure (spike autocorrelation, calculated on ISIs) and auto-correlation coefficient (calculated across trials). The latter is the classical definition used by (Murray et al., 2014) and requires a specific normalization, i.e. subtracting the mean (across trials) and dividing by the standard deviation (across trials). Our method is based on the classic spike autocorrelogram computing the distribution of lags from each spike and the surrounding emitted spikes.
To clarify the approach, the methods section has been rewritten. It now reads: "To do so, we computed the lagged differences between spike times up to the 100th order, i.e. the time differences between any spike and its n successors (up to n = 100) at the unit level. The lagged differences were then sorted in \Δ_t=3.33ms bins from 0 to 1000ms. Given N the count of lagged differences within bins, we compute the probability distribution function as N/(\sum(N)*\Δ_t). …" in the Spike Autocorrelogram Analysis section of Methods.
1.2. The new method, described as superior to the previous method by Murray et al., 2014, appears to have access to more spikes than the Murray's method (Figure 2). Where is this additional data coming from? While Murray's method was applied to the pre-cue period, the time epoch used for the analysis with the new method is not stated clearly. It seems that the new method was applied to the data through the entire trial duration and across all trials, hence more spikes were available. If so, then changes in firing rates related to behavioral events contribute to the autocorrelation, if not appropriately removed. For example, the Murray's method subtracts trail-averaged activity (PSTH) from spike-counts, similar to shuffle-correction methods. If a similar correction was not part of the new method, then changes in firing rates due to coding of task variables will appear in the autocorrelogram and estimated timescales. This is a serious confound for interpretation of the results in Figure 4. For example, if the firing rate of a neuron varies slowly coding for the gauge size across trials, this will appear as a slow timescale if the autocorrelogram was not corrected to remove these rate changes. In this case, the timescale and GLM are just different metrics for the same rate changes, and the correspondence between them is expected. Before results in Figure 4 can be interpreted, details of the method need to be provided to make sure that the method measures intrinsic timescales, and not timescales of rate changes triggered by the task events. This is an important concern also because recent work showed that there is no correlation between task dependent and intrinsic timescales of single neurons, including in cingulate cortex and PFC (Spitmaan et al., PNAS, 2020).
The reviewer’s comments are very important and suggest that we did not explain the method clearly enough. In the new version we have tried to be more explicit.
In short, the method is based on the distribution of spike emitted after each spike, what has been called in the neurophysiological literature spike autocorrelogram. As explained in the method and main text, the spike autocorrelogram allows extraction of 2 temporal signatures: the latency at the peak, and the decay (as measured by an exponential fit). Author response image 3 illustrates these two measures for 1 example unit:
Contrary to the Murray method it is not based on the autocorrelation of binned spike counts measured on a particular time epoch and across trials. As the reviewer understood, the method requires a spike time timeseries sufficiently long to extract a non noisy spike autocorrelogram. It does not rely on a trial structure imposed by a task. One can measure the temporal signature on data recorded at rest. Importantly, and in contrast to studies using the Murray method, our analyses allowed us to include the large majority of units in our database, and even to seperate by unit types (FS and RS). Note that since in many studies the success rates of fitting an exponential were less than 55%, it followed that about half of the recorded data were excluded from analyses. We think that it could lead to serious sampling biases if for example one dissociates unit types as we did in this study.The reviewer is concerned about a potential relationship between task-induced variations of firing rate and the timescale measured on the spike autocorrelogram, and that this might explain the link observed between encoding of task variables and TAU median-splitted populations of units. As explained above (response to editor) we have made several controls to answer this concern.
We first asked whether TAU and firing average rate are related across units. Because there is no absolute method to measure firing rate, we tested measures for several periods of data (whole recording, feedback epoch, etc..). We found little if any correlation.
In addition, we would like to emphasize the following:
– In Stoll et al. we observed using glmm on spike counts that the gauge size was not encoded continuously during the task (i.e. across trials and during trials of the main task), but only at the time of decisions (see figure 5 in Stoll et al. Nat Comm. 2016). Thus gauge-related changes in firing rates were not structured in a way that it would impact the spike autocorrelograms.
– The difference in time scales between long and short timescale populations are of the order of a few hundreds milliseconds. In contrast the trial structure of the task, the gauge changes and the feedback deliveries, evoked by the reviewer is of the order of several tens of seconds. The two are thus incompatible and it is unlikely that task-related fluctuations explain differences in timescales measures with our methods.
Besides, as mentioned in our response to Reviewer 1 comment 6, we devised a specific version of the HMM analysis to assess the possible role of the temporal structure of the task in the difference observed in the MCC compared to the LPFC in experimental data (see Figure 7 —figure supplement 1). Specifically, we simulated neural spiking activity of each neuron in LPFC and MCC areas as a non-homogeneous Poisson model (NHPM) whose instantaneous firing rate was given by the post-stimulus time histogram (PSTH) of that neuron upon each epoch type across all possible task events (e.g., reward, go signal, etc.). This simulated activity was only dependent on the coding of taskvariable across trials but did not (as opposed to in original data) contain the specific fine-grained activity temporal structure within individual trials (i.e. at the level of individual ISIs) that is internally generated by frontal areas. Applying HMM analysis to such transformed population activity showed that the ratio between LPFC and MCC state duration distributions was lost. This indicated that task variable coding was not responsible for the difference of temporal structure between LPFC and MCC spiking activity, but that it relied on the precise temporal structure of spiking, independent of task epochs. So it is unlikely that the slow timescales arise from some slow firing rate variation due to the structure of the task. This makes sense, as the duration of task epochs was statistically similar when recording in the LPFC and MCC (i.e. the monkey performed the same task).
Regarding the Spitmaan et al. study, one must note that this study, like all other preceding investigations of timescales, could not take into account cell type differences which, according to our results, are a key element discriminating between task dependent properties and temporal signatures. This element put aside, our analyses did not directly test a correlation between estimated behavioral timescales and estimated intrinsic timescales. We simply observed the functional dissociation in cell type and functional properties on median split cell populations. To follow the reviewer’s concern we now mention this in the discussion:
“Spiking timescales distributions have been related to persistent activity, choice value and reward history in the LPFC and MCC (Bernacchia et al. Nat. Neurosci, 2011; Cavanagh et al., Nature Comm 2018; Meder et al., Nature Comm 2017; Wasmuht et al., Nat Comm. 2018), but in a recent study no correlation was observed between task dependent and intrinsic timescales at the unit level (Spitmaan et al., 2020). In all those studies, cell types were not considered. Here, we could not estimate taskrelevant timescales to correlate with TAU, but we found that the spiking timescales of MCC RS units increased on average during periods of engagement in cognitive performance”.
2. The balanced network model with a variety of biophysical currents is interesting and it is impressive that the model reproduces the autocorrelation signatures in the data. However, we need to better understand the network mechanism by which the model operates.
2.1. The classical balanced network (without biophysical currents such as after-hyperpolarization potassium) generates asynchronous activity without temporal correlations (Renart et al., Science, 2010). The balanced networks with slow adaptation currents can generate persistent Up and Down states that produce correlations on slow timescales (Jercog et al., eLife, 2017). Since slow after-polarization potassium current was identified as a key ingredient, is the mechanism in the model similar to the one generating Up and Down states, or is it different? Although the biophysical ingredients necessary to match the data were identified, the network mechanism has not been studied. Describing this network mechanism and presenting the model in the context of existing literature is necessary, otherwise the results are difficult to interpret for the reader.
In the present study, the network model displays collective dynamics that peregrinate between distinct metastable states, each of which implying (and captured by the HMM as) a subset of strongly active neurons. At any time, as consecutive states occur, activity is present in the network, thus being permanent, i.e. with no interruption or ‘DOWN’ state. Such dynamics therefore strongly depart from UP and DOWN dynamics in a recurrent neural network (Jercog et al., eLife, 2017), in which all neurons are either all active or all silent. Moreover, while our network model is aimed at reproducing collective dynamics in frontal areas of awake monkeys, that of Jercog et al. 2017 is designed to capture spiking in anesthetized rodent cortical networks. Engel et al., 2016 have shown that primate visual networks can switch between two states that correspond to ON and OFF attentional states of awake monkeys. These states may resemble UP and DOWN states in other preparations/animals/areas. However, they differ from UP and DOWN states in that OFF state activity corresponded to a weak but non-null activity, whereas OFF states in anesthetized animals or in vitro are defined by the absence of any spiking. Therefore, data from Engel et al. rather likely correspond to switching between two metastable states fulfilling specific computations required by a binary attentional demand. Besides, the main parameters setting the timescale and duration of states in the present study are actually the maximal conductance of GABA-B receptors and the heterogeneity of synaptic weights. By contrast, the AHP maximal conductance is essential in setting the latency of the autocorrelogram, the other dimension of the temporal signature (as well as limiting the amount of fast transitions). AHP conductances have been implied in many possible computational functions such as network decorrelation (Renart et al., Science, 2010), the complexity of network dynamics (Cartling, International Journal of Neural Systems, 1996) or the linearization of the input-output function of neurons (Wang, J Neurophysiol, 1998). While possibly arising in different functional contexts or cooperating together in cortical networks, such diversity of AHP-dependent computational functions may emerge from differences in gating dynamics of AHP currents. Here, for instance, we consider fast AHP (20 ms time constant), while Jercog et al., 2017 consider slow AHP conductance (500 ms). We now clarify these elements in the Discussion (“Frontal temporal signatures uncover metastable dynamics” subsection).
2.2. Does the model operate in a physiologically relevant regime where the firing rates, Fano factor etc. are similar to the data? It is hard to judge from Figure 5b and needs to be quantified.
We agree that this comparison was lacking. We now provide (see Table 2) a quantitative evaluation of whether spiking statistics is similar in monkeys and model data in LPFC and MCC areas. We find that in both cases spiking frequency is below 5 Hz, CV above 1, and CV2 and LV (i.e. local (in the temporal domain) measures of ISI variability) and Fano factors around 1. An exact quantitative match between model and empirical was not sought for being (1) extremely complex to achieve, given the large amount of other essential constraints applied to the model (e.g. connective patterns, AMPA/NMDA ratios, E/I balance, temporal timescales and latencies of both excitatory and inhibitory neurons) and (2) not crucial, as the exact value of firing rates (and second order spiking moments) is not a determinant of timescales (see Figure 7 —figure supplement 1) but a function of biophysical parameters (GABA-B, AHP). In this context, the general quantitative consistency between spiking dynamics in model and monkey data indicates that our recurrent networks qualitatively fit with physiological statistics of frontal areas in the awake regime.
2.3. The latency of autocorrelation is an interesting feature in the data. Since the model replicates this feature (which is not intuitive), it is important to know what mechanism in the model generates autocorrelation latency.
In the model, we find that the principal mechanism responsible for the latency of the autocorrelation is the maximal conductance of the AHP current in excitatory neurons. Indeed, in these neurons, AHP currents are spike-evoked and their presence increases refractoriness of neurons, i.e. first-order latency, which is a major contributor of latency. This was detailed in supplementary material but has now also been stated in the Discussion (“Frontal temporal signatures uncover metastable dynamics” subsection).
3. HMM analysis is used to demonstrate metastability in the model and data, but there are some technical concerns that can undermine these conclusions.
3.1. HMM with 4 states was fitted to the data and model. The ability to fit a four-state HMM to the data does not prove the existence of metastable states. HMM assumes a constant firing rate in each "state", and any deviation from this assumption is modeled as state transitions. For example, if some neurons gradually increase/decrease their firing rates over time, then HMM would generate a sequence of states with progressively higher/lower firing rates to capture this ramping activity.
This is a interesting point and we thank reviewer #2 for this remark. Network neural activity in behaving monkeys displays a dynamical complexity that exceeds that of HMMs, with non-stationary, continuous variations of activity operating at many timescales. Thus, analyzing neural data with HMMs therefore implies a form of simplification due to the representation of non-stationary continuous activity dynamics into stationary discrete states. However, some evidence indicates that HMMs constitute robust models (compared to PCA and PSTH analysis, e.g.), when applied to transient (i.e. non stationary) forms of coding (Jones et al., PNAS, 2007). Actually, HMMs have long been applied in frontal areas in contexts where neural activity is not stationary but evolves with ongoing task constraints and cognitive demands (perceptual events, decision processes, actions, e.g.; La Camera et al., Curr Opin Neurobiol 2019). These studies have evidenced that although non-stationary in essence, neural activity switches between well-separated states of activity – where it dwells most of the time –, within which firing rates are approximately stationary. A large corpus of theoretical studies has shown that states may correspond to attractors such as stable activity state within strongly connected Hebbian assemblies (Hopfield, PNAS, 1982; Brunel et Wang, Journal of Computational Neuroscience 2001; LaCamera et al., Curr Opin Neurobiol 2019). It has also largely been demonstrated that in highdimensional chaotic dynamics (as found in awake cortical recurrent networks), these attractors are not stable in all directions (i.e. they form “pseudo-attractors” or “attractor ruins”), which accounts for transitions between (metastable or quasi-stationary) states (see chaotic itinerancy (Tsuda, Behav Brain Sci 2001), dynamics along heteroclinic channels (unstable directions; Rabinovich et al. Science + PloS CB, 2008) or noise transitions between connectivity-related attractors (La Camera et al., Curr Opin Neurobiol 2019 or Miller, Curr Opin Neurobiol, 2016, e.g.)). In this context, HMM are interesting, as they allow highlighting the property of residing in preferential regions (e.g. Hebbian attractor) (as opposed to erratic wandering with uniform probability) in the neural state-space. Also, they are of particular precision in the temporal dimension, compared to traditional methods (Seidemann et al., J Neurosc 1996; Jones et al., PNAS, 2007).
Besides, we understand that systematic drifts in neuron’s firing frequencies at large time scales may bias HMMs towards modeling drifts as state transitions or bias the relative importance (duration) of specific states at different times. This bias might be important in the monkey recorded data HMM analyses (the neural frequency in our model networks being stationary by design). To specifically answer the issue and disambiguate the possibility that the observed states are reflecting drifts in data frequency, we assessed whether average neural frequency drift and average HMM state durations are correlated (Figure 7 —figure supplement 2). To do so, we repeated the HMM analysis on temporal segments of 50s of the 600s monkey data (i.e. the HMM was performed on spike data from 0-50s, 50100s, etc. all the way up to 550-600s). Moreover, this procedure was repeated for different segment durations (i.e. segments of 50, 60, 75s, 100, 120, 150, 200 300, 600s), to assess many different time scales over which neural frequency drifts might strongly bias HMM state durations, if at all. In each case, the average neural frequency drift was measured as the absolute difference in neural frequency between the first and second half of the data segment (e.g. 0-25 vs 25-50s for first 50s segment), averaged across all neurons.
We found that average neural frequency drifts were correlated to average HMM state durations within each cortical area when pooling all segment durations (LPFC: rho ~0.54, p-val ~9.9*10-5; MCC: rho ~0.38, p-val ~8.3*10-3; only a few correlations survived when looking at individual segment durations; LPFC, 60s: rho ~-0.68, p-val ~4.3*10-2; LPFC, 100s: rho ~0.85, p-val ~3.3*10-2; MCC, 120s: rho ~0.95, p-val ~1.3*10-2). However, state duration and frequency drift were never correlated when pooling data from both cortical areas, be it within or across temporal subdivisions. This was because the drifts in frequency were equivalent in the LPFC and MCC, whereas state durations were much longer in the MCC, compared to the LPFC. Thus, this indicates that average neural frequency drift does not account for the observed difference in state durations between LPFC and MCC.
This conclusion was confirmed by additional analysis of the relation between frequency drift and state durations on the 600s segment. Dividing neurons into two halves – “most” or “least” drifting – according to the amplitude of their frequency drifts across time (calculated similarly as before for each neuron on data from 0-600s), we then analyzed the spiking activity of each group via HMM. The frequency drift was ~6.6x higher (in LPFC and MCC) in the “most drifting” neurons, compared to “least drifting” neurons, indicating a strong difference in frequency drift between both groups across areas. However, the HMM state duration only increased by ~1.75x in “most drifting” vs. “least drifting” neurons across areas (~1.7x for LPFC, 1.8x for MCC). In contrast, the ratio of MCC vs. LPFC average state duration across neuron groups was much larger than 1.75x, being instead around ~5.75x (~5.5x for “least drifting” and ~6x for “most drifting” neurons). This confirms that neural frequency drift was not a major cause in the difference between LPFC and MCC average state durations, similarly to the previous analysis. These analyses are now provided in Figure 7 —figure supplement 2 and Table 3 in the article.
In addition, metastability implies exponential distributions of state durations, which was not verified.
We have now verified that distributions of state durations are exponential. To do so, we performed linear regressions on the logarithm of the count probability (binned in 50ms bins) of HMM state durations for each of the 100 HMM analyses, and reported the slope coefficient p-value and R2, for both cortical areas in both monkey and model data. We observed overall that R2 scores were above ~0.6, with significantly non-zero slopes (p<0.05) in all monkey HMM analyses, in all MCCm models and in the immense majority of PFCm models (98/100). Thus, state duration distributions appear exponential, lending credence to the metastable nature of HMM states. We have now added this analysis as Figure 7 —figure supplement 3.
No model selection was performed to determine the necessary number of states. Therefore, the claims of metastable dynamics are not supported by the presented analysis.
We understand this criticism and have now performed model selection to determine the necessary number of states. To do so we evaluated Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) values for different numbers of states in both cortical areas and in both monkey and model data. We observed that AIC and BIC did not substantially change when varying the number of states in the range 2-10. Other methods for determining the number of states were explored (pseudoresidual error analysis, visual analysis and autocorrelation of HMM-generated spike rasters, and analysis of the structure of individual states) and none indicated a strong preference for any given number of states (not shown). More importantly, even when changing the number of states, HMM state duration results did not substantially differ, compared to those with 4 states. While LPFC and MCC state durations slightly decreased with more states (since the same time period is divided into a greater number of states), MCC state durations remained 1-2 orders of magnitude longer than LPFC state durations. Thus, choosing a number of states in the range 2-10 was equivalent according to this analysis. Since (1) a low number of states is more parsimonious in terms of data interpretation (Pohle et al., 2017) in general, (2) the task structure contains a low number of possible states in terms of actions (four), reward on the last trial (incorrect trial, first correct trial, correct trial after previous correct trials) and behavioral states (exploration, exploitation), and (3) increasing the number of HMM states substantially increased computational time of HMMs, we kept the low number of states chosen, in this instance 4 states. We have now added this analysis as Figure 7 —figure supplement 4.
3.2. HMM was fit to a continuous segment of data lasting 600s, and the data was pooled across different recording sessions. However, different sessions have potentially different trial sequences due to the flexibility of the task. How were different trial types matched across the sessions? If trial-types were not matched/aligned in time, then the states inferred by the HMM may trivially reflect a concatenation of different trial types in different sessions. For example, the same time point can correspond to the gauge onset in one session and to the work trial in another session, and vice versa at a different time. If some neurons respond to the gauge and others to the work, then the HMM would need different states to capture firing patterns arising solely from concatenating the neural responses in this way. This confound needs to be addressed before the results can be interpreted.
We previously showed that the temporal structure of the task had no effect on the LPFC/MCC state duration difference (see Reviewer #1, major point #6), and that pseudo-populations preserve global spike-timing structure (rather than spike-timing structure within individual recording sessions). As such, trial types do not affect HMM state duration, and sessions were thus aggregated without any trial-type alignment.
We have now assessed this issue by two complementary analyses.
First, we have showed that the temporal structure of the task had no effect on the LPFC/MCC state duration difference (Figure 7 – supplement 1b). To do so, we have simulated neural spiking activity of each neuron in LPFC and MCC areas as a non-homogeneous Poisson model (NHPM) with instantaneous firing rate given by its post-stimulus time histogram (PSTH) during behavioral epoch type between all possible task events (e.g., reward, go signal, etc.). This simulated activity only depended on task-variable coding and did not contain more specific fine-grained spiking temporal structure. Applying HMM analysis to this surrogate activity showed no difference of state durations between LPFC and MCC. This indicated that task variable coding was not responsible for the state duration difference between LPFC and MCC but relied on the precise temporal structure of spiking, independent of task epochs. This also indicated that accounting for neural states in monkeys with a biophysical model devoid of task structure was not an issue per se.
Besides, we have also now applied HMM analysis to neural activity in the model, using pseudopopulations built of groups of neurons drawn from LPFC and MCC simulations of distinct random neural networks, with group statistics (number of recording sessions and numbers of simultaneously recorded excitatory and inhibitory neurons within each recording session) similar to those of monkey LPFC and MCC data (Figure 7 – supplement 1a) and showed that pseudo-populations gathered from different sessions preserve global spike-timing structure (rather than spike-timing structure within individual recording sessions).
A more detailed presentation of these supplementary analyses can be found in response to reviewer 1 (major point 6).
Altogether, these new analyses indicate that trial types do not affect HMM state duration, and that sessions can thus be aggregated without any trial-type alignment.
https://doi.org/10.7554/eLife.63795.sa2Article and author information
Author details
Funding
Fondation pour la Recherche Médicale (DEQ20160334905)
- Emmanuel Procyk
Agence Nationale de la Recherche (ANR-10-SVSE4-1441)
- Emmanuel Procyk
Agence Nationale de la Recherche (ANR-16-NEUC-0006-01)
- Bruno Delord
Agence Nationale de la Recherche (ANR-11-LABX-0042)
- Emmanuel Procyk
Fondation pour la Recherche Médicale (FDT201904008187)
- Vincent Fontanier
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank C Nay for administrative support and C Wilson and J Naudé for helpful discussions and proofreading the manuscript.
Ethics
All procedures followed the European Community Council Directive (2010) (Ministère de l'Agriculture et de la Forêt, Commission nationale de l'expérimentation animale) and were approved by the localethical committee (Comité d'Ethique Lyonnais pour les Neurosciences Expérimentales, CELYNE, C2EA683 #42).
Senior and Reviewing Editor
- Timothy E Behrens, University of Oxford, United Kingdom
Version history
- Preprint posted: August 21, 2020 (view preprint)
- Received: October 7, 2020
- Accepted: May 27, 2022
- Accepted Manuscript published: May 30, 2022 (version 1)
- Version of Record published: June 15, 2022 (version 2)
Copyright
© 2022, Fontanier, Sarazin et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
-
- 833
- Page views
-
- 260
- Downloads
-
- 1
- Citations
Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Cell Biology
- Neuroscience
The amyloid-beta (Aβ) plaques found in Alzheimer's disease (AD) patients' brains contain collagens and are embedded extracellularly. Several collagens have been proposed to influence Aβ aggregate formation, yet their role in clearance is unknown. To investigate the potential role of collagens in forming and clearance of extracellular aggregates in vivo, we created a transgenic Caenorhabditis elegans strain that expresses and secretes human Aβ1-42. This secreted Aβ forms aggregates in two distinct places within the extracellular matrix. In a screen for extracellular human Aβ aggregation regulators, we identified different collagens to ameliorate or potentiate Aβ aggregation. We show that a disintegrin and metalloprotease ADM-2, an orthologue of ADAM9, reduces the load of extracellular Aβ aggregates. ADM-2 is required and sufficient to remove the extracellular Aβ aggregates. Thus, we provide in-vivo evidence of collagens essential for aggregate formation and metalloprotease participating in extracellular Aβ aggregate removal.
-
- Developmental Biology
- Neuroscience
Development of the nervous system depends on signaling centers – specialized cellular populations that produce secreted molecules to regulate neurogenesis in the neighboring neuroepithelium. In some cases, signaling center cells also differentiate to produce key types of neurons. The formation of a signaling center involves its induction, the maintenance of expression of its secreted molecules, and cell differentiation and migration events. How these distinct processes are coordinated during signaling center development remains unknown. By performing studies in mice, we show that Lmx1a acts as a master regulator to orchestrate the formation and function of the cortical hem (CH), a critical signaling center that controls hippocampus development. Lmx1a co-regulates CH induction, its Wnt signaling, and the differentiation and migration of CH-derived Cajal–Retzius neurons. Combining RNAseq, genetic, and rescue experiments, we identified major downstream genes that mediate distinct Lmx1a-dependent processes. Our work revealed that signaling centers in the mammalian brain employ master regulatory genes and established a framework for analyzing signaling center development.