1. Neuroscience
Download icon

Cortical state transitions and stimulus response evolve along stiff and sloppy parameter dimensions, respectively

  1. Adrian Ponce-Alvarez  Is a corresponding author
  2. Gabriela Mochol
  3. Ainhoa Hermoso-Mendizabal
  4. Jaime de la Rocha
  5. Gustavo Deco
  1. Center for Brain and Cognition, Computational Neuroscience Group, Department of Information and Communication Technologies, Universitat Pompeu Fabra, Spain
  2. Institut d’Investigacions Biomèdiques August Pi i Sunyer (IDIBAPS), Spain
  3. Institució Catalana de la Recerca i Estudis Avançats (ICREA), Spain
  4. Department of Neuropsychology, Max Planck Institute for Human Cognitive and Brain Sciences, Germany
  5. School of Psychological Sciences, Monash University, Australia
Research Article
  • Cited 0
  • Views 713
  • Annotations
Cite this article as: eLife 2020;9:e53268 doi: 10.7554/eLife.53268

Abstract

Previous research showed that spontaneous neuronal activity presents sloppiness: the collective behavior is strongly determined by a small number of parameter combinations, defined as ‘stiff’ dimensions, while it is insensitive to many others (‘sloppy’ dimensions). Here, we analyzed neural population activity from the auditory cortex of anesthetized rats while the brain spontaneously transited through different synchronized and desynchronized states and intermittently received sensory inputs. We showed that cortical state transitions were determined by changes in stiff parameters associated with the activity of a core of neurons with low responses to stimuli and high centrality within the observed network. In contrast, stimulus-evoked responses evolved along sloppy dimensions associated with the activity of neurons with low centrality and displaying large ongoing and stimulus-evoked fluctuations without affecting the integrity of the network. Our results shed light on the interplay among stability, flexibility, and responsiveness of neuronal collective dynamics during intrinsic and induced activity.

Introduction

How biological systems achieve a tradeoff between stability and flexibility is a central question in biology. A candidate explanation for the coexistence of these two features is sloppiness (Machta et al., 2013; Transtrum et al., 2015). In general, sloppiness is a property of complex models exhibiting large parameter uncertainty when fit to data, meaning that different combinations of parameters lead to a similar system behavior, while changes in some few critical parameters, called stiff parameters, significantly modifies it. In this way, biological systems can be either robust to large fluctuations of input/environmental signals which effects are embedded in a high-dimensional subspace of insensitive parameters, or, on the contrary, by tuning some few parameters, configured to be highly sensitive and selective to relevant signals.

Recently, it has been shown that the spontaneous activity of neural circuits presents sloppiness both in vitro and in vivo (Panas et al., 2015), suggesting that collective activity is stabilized by a subset of highly active and stable neurons, while the activity and co-activity of the remaining neurons present larger spontaneous fluctuations without strongly affecting the collective statistics. However, this view is challenged by extensive research showing that the spontaneous cortical activity transits through different synchronized and desynchronized cortical states (Marguet and Harris, 2011; Harris and Thiele, 2011; Luczak et al., 2013; Pachitariu et al., 2015) that represent statistically different collective behaviors (Hahn et al., 2017) with different information processing capabilities (Pachitariu et al., 2015; Engel et al., 2016; Beaman et al., 2017). Moreover, how sensory inputs affect sloppiness is unknown and it is a relevant question to understand how sensory stimuli change the network state in a way that responsiveness and stability are ensured. In the present study, we examined how changes in neural network parameters correlate with spontaneous transitions among cortical states and stimulus-evoked responses.

To answer these questions, we recorded the neuronal spiking activity in the primary auditory cortex (A1) of six anesthetized rats. We analyzed the joint activity of groups of neurons while the cortex spontaneously transited through different synchronized and desynchronized cortical states and intermittently received external acoustic stimuli. We used a statistical model to describe the joint spiking activity with a small number of parameters. We found that the estimated parameters of neuronal ensemble activity presented sloppiness and that sensory inputs and cortical state transitions evolved in different pathways in parameter space. Specifically, we found that cortical state transitions evolve along stiff dimensions, whereas sensory-evoked activity evolves along sloppy dimensions. Finally, we showed that stiff parameters are related to the activity and co-activity of neurons with high centrality within the functional network of the recorded neurons.

Results

We recorded spontaneous and stimulus-evoked population activity from the primary auditory cortex (A1) of urethane-anesthetized rats (n = 6) using multisite silicon microelectrodes (see Materials and methods). The data was composed of activity from Npop well-isolated single units (Npop= 44-147 neurons) and some spike-trains from multi-unit activity (3-103 spike-trains). Unless otherwise specified, the analyses present here focused on single-unit activity only. We analyzed the data during spontaneous activity and in response to acoustic 'clicks' (5-ms square pulses; inter-stimulus interval, 2.5 or 3.5 s). To track the evolution of the neuronal activity, we divided each recording session into NE adjacent epochs of 100 s, each one containing 12–29 stimulus presentations. Within each 100-s epoch the data was separated into spontaneous activity, that is the activity during 1.5-s intervals preceding each stimulus (i.e., 18–43.5 s of spontaneous activity in total for each epoch), and stimulus-evoked activity, that is the activity right after the stimulus onset (Figure 1A–B).

Experiment and analysis designs.

(A) Each recording session was divided into NE adjacent epochs of 100 s. (B) Each epoch contained a series of stimulus presentations. Stimuli consisted on acoustic clicks. For each 100-s epoch we collected the spontaneous activity, that is the activity during 1.5-s intervals preceding each stimulus (red intervals), to build concatenated binary data. (C) Binary data was obtained by discretizing time in bins of dt = 10 ms. Within each time bin, the ensemble activity of N neurons was described by a binary vector, σ=σ1, σ2, , σN, where σi=+1 if the i-th neuron generated a spike (black) and σi=-1 otherwise (white). (D) Maximum entropy models were used to describe the binary patterns of subsets of 10 neurons, during each 100-s epoch. The model parameters Ω={h,J} represents the intrinsic tendency of neuron i towards activation (σi=+1) or silence (σi=-1), noted hi, and the effective interaction between neurons i and j, noted Jij.

Description of spontaneous activity patterns using maximum entropy models

We first examined the temporal evolution of the spontaneous activity across the NE epochs. Because we were interested in the evolution of the statistics of ensemble activity, we described the collective activity of groups of N single-units using a maximum entropy model (MEM) (Schneidman et al., 2006; Shlens et al., 2009; Tkačik et al., 2015) in each epoch (see Materials and methods and Figure 1C–D). These models allowed us to describe the patterned activity with a small number of parameters. To fit the model, time was discretized in bins of dt = 10 ms. Within each time bin, the ensemble activity of N neurons was described by a binary vector, σ=σ1, σ2, , σN, where σi=+1 if the i-th neuron fired a spike in that time bin and σi=-1 otherwise. The collective activity was determined by the probability distribution Pσ over all 2N possible binary patterns. The MEM fits Pdataσ by finding a distribution PMEMσ that maximizes its entropy under the constraint that the activation rates (σi) and the pairwise correlations (<σiσj>) found in the data are preserved in the model. It is known that the maximum entropy distribution that is consistent with these constraints is the Boltzmann distribution, Pσe-E σ , where Eσ is the energy of the pattern σ, given by: Eσ=-i=1Nhiσi+12j=1NJijσiσj  (Schneidman et al., 2006; Tkačik et al., 2015). The model parameter hi represents the intrinsic tendency of neuron i towards activation (σi=+1) or silence (σi=-1) and the parameter Jij represents the effective interaction between neurons i and j. Once we learned the parameters Ω={h,J} using a gradient descent algorithm (see Materials and methods), the expected probability of any pattern is known. For each recording session and for each of the NE epochs, we fitted the model using the spontaneous binarized activity from an ensemble of N = 10 randomly selected single neurons from the entire population of Npop single neurons. We chose N = 10 because 100-s epochs provided around 5000 observed spontaneous patterns, which is a reasonable amount to get an estimate of the distribution of the 210 = 1024 possible patterns. To accurately estimate models of larger N, the epochs ought to be much larger preventing possibility to investigate the temporal evolution of the model along the experiment. We finally repeated the process of randomly choosing N = 10 single units Q times for each experiment (for datasets 3 and 5: Q=10 ensembles, otherwise: Q=20). In summary, for each recording session, we built Q×NE models, each composed of 10 units. Before studying the evolution of the model parameters Ω(t) across epochs (t= 0, 1, 2..), we first evaluated how well the MEM fitted the data.

For each epoch, we used the Jensen-Shannon divergence (DJS, see Materials and methods) to measure the similarity between the probability distribution of the empirical and model binary patterns (Figure 2A–C). We compared this similarity to the distribution of binary patterns predicted from independent-MEMs, for which only the activation rates were preserved (i.e., only h was optimized). We found that the empirical distribution was well approximated by MEMs and that, for all recording sessions, the goodness-of-fit (i.e., 1/DJS) was orders of magnitude higher for MEMs than for independent-MEMs (Figure 2C), leading to excellent model performances (i.e., Kullback-Leibler ratio equal to 0.95 ± 0.03 on average, see Materials and methods).

Fitting maximum entropy models (MEMs) to spontaneous activity patterns.

(A) Comparison between the probability distribution of empirical binary patterns and the probability distribution predicted by MEMs (black dots) and independent models (gray dots), for all epochs and all neuronal ensembles. Every point represents a binary pattern that has appeared in the data at least once. Red line represents the identity line. (B) Jensen-Shannon divergence (DJS) between spiking data and MEMs, and between spiking data and independent models, across time, averaged across neuronal ensembles. Error bars are smaller than the widths of the traces. Data in (A) and (B) correspond to one example rat (#1). (C) Goodness-of-fit (1/DJS) for MEMs and for independent models, averaged over all models (i.e., all ensembles and all epochs), for each rat. Error bars indicate SEM.

Temporal evolution of activity observables, model parameters, and their sensitivity

We next analyzed the temporal evolution of the different spiking data statistics and the model parameters. We first measured the temporal variation of the activity observables (i.e., firing rates and pairwise correlations) by calculating the average Pearson correlation (or similarity γ; see Materials and methods) between the values in epoch t and those in epoch t+Δt (Figure 3A). This similarity rapidly decayed with t, indicating that the observables substantially changed over time. We next examined how much these variations influenced the evolution of the collective activity characterized by the distribution of binary patterns. For this, we evaluated how well the data in a given epoch t could be explained by the MEM constructed using the data at time t+Δt. Specifically, we calculated DJSΔt, given by the average Jensen-Shannon divergence between the distribution of data binary patterns in epoch t, that is Pdata, t, and the distribution of binary patterns predicted by the MEM constructed using the data in epoch t+Δt, that is PMEM,t+Δt (see Materials and methods). We found that DJSΔt increased as a function of t, indicating that the collective activity changed during the recording session, so that the model parameters Ω(t) at epoch t did not predict the collective activity at epoch t+Δt (Figure 3B). Indeed, the model parameters substantially changed over time, with a rapidly decaying similarity (Figure 3C, orange and purple traces).

The sensitivity of model parameters was more stable than activity observables.

(A) Similarity (i.e., Pearson correlation coefficient) of mean firing rates (red) and pairwise correlations (purple) as a function of elapsed time Δt. (B) Jensen-Shannon divergence (DJS) between the distribution of empirical spiking patterns in epoch t and the distribution of binary patterns of the pairwise MEM in epoch t + Δt, averaged over all t. (C) Left: Similarity of Fisher information matrix (FIM) elements, biases (hi), and couplings (Jij) as a function of elapsed time Δt. Right: FIMs at time t = 0 h and t = 1 h. Data in (A), (B), and (C) correspond to one example rat (#1); traces show averages over neuronal ensembles and shaded areas correspond to SEM. (D) Similarity of FIM elements, rates, biases, correlations, and couplings after 1/2 hour (i.e., Δt = 30 min), averaged across all neuronal ensembles, for each rat. Error bars indicate SEM.

As shown in Panas et al. (2015), changes in model parameters can differently contribute to collective activity, since the model can be sensitive to changes in some few combinations of parameters. Following this, we next evaluated the sensitivity of model parameters by calculating the Fisher information matrix (FIM, see Materials and methods) for each neuronal ensemble and each epoch. The FIM measures how much the model log-likelihood PMEM(σ|Ω) changes with respect to changes in the parameters Ω. We first notice that the FIM had the highest stability across time, compared to the data firing rates and correlations and the model parameters (Figure 3C, blue trace). Indeed, the similarity after 1/2 hour was γ= 0.882 ± 0.002 for the FIM, γ= 0.732 ± 0.003 for the firing rates, γ= 0.551 ± 0.004 for the biases, γ= 0.364 ± 0.004 for the correlations, γ= 0.234 ± 0.003 for the couplings (F4,495 = 305.73, p<0.001, one-way ANOVA followed by Tukey's post hoc analysis) (Figure 3D). Altogether these results show that the sensitivity of the model parameters remained relatively stable despite substantial changes in firing rates, correlations, collective activity and the model parameters themselves.

Spontaneous neuronal activity presents sloppiness

Having shown that the sensitivity of model parameters was relatively stable during the recording sessions, we next studied the structure of the FIMs. First, we noted that most elements of the FIM had near-zero values (Figure 4A) indicating that most of the parameters had a small effect on the model log-likelihood. In contrast, a small fraction of elements had values strongly different from zero as revealed by the heavy tail of the distribution of FIM values (Figure 4A). To identify the parameter combinations that had the strongest effect on model behavior, we decomposed the FIM into eigenvectors and classified them according to their eigenvalue (Figure 4B). We observed that, except for some few eigenvalues, most of the FIM eigenvalues were small, corresponding to combinations of parameters that had little effect on model behavior. These unimportant parameter combinations defined the sloppy dimensions of the model. The few eigenvectors with large eigenvalues defined the stiff parameter dimensions along which the model behavior was strongly affected.

Figure 4 with 1 supplement see all
Model parameters presented sloppiness and predominantly evolved along sloppy dimensions.

(A) Distribution of FIM elements, for all epochs, all neuronal ensembles, and all rats. (B) Eigenvalues of the FIM, average across epochs and neuronal ensembles, for an example rat (# 1). Shaded areas represent SEM. Stiff and sloppy dimensions correspond to FIM eigenvectors of lowest and highest ranks, respectively. (C) Projection of Ω(t') into the first three eigenvectors of the FIM from a given epoch t, for all t't. Data from one neuronal ensemble from rat 1. (D) Projection of Ω(t') into the first and the 20th eigenvectors of the FIM from epoch t, noted νt,1 and νt,20, respectively, for all t't. Top inset: distribution of projections into νt,1. Right inset: distribution of projections into νt,20. Note higher variance of projections into νt,20 than into νt,1. Data from one neuronal ensemble from rat 1. (E) Average variance of projections of Ω(t') into the different eigenvectors of the FIM at epoch t (for all t't), for the different rats. Traces represent average over neuronal ensembles and shaded areas represent SEM.

In the following we showed that the temporal evolution of the model parameters occurred predominantly along the sloppy dimensions. For this, we projected the parameters Ω(t), calculated at time t', into the eigenvectors of the FIM at time t, denoted νt,1, νt,2, ,νt,k, , where k is the rank of the eigenvector (Figure 4C). For each dimension, or eigenvector, we obtained a distribution of projections of parameters Ω(t) (Figure 4D). To quantify how much the parameters varied along each eigenvector, we calculated the average variance of each projection as a function of the rank of the eigenvector. We found that the projection variance increased as a function of the eigenvector’s rank for all datasets (Figure 4E). This indicates that the model parameters predominantly evolved along sloppy dimensions (i.e., FIM eigenvectors of highest rank k), while they remained relatively stable along stiff dimensions (i.e., FIM eigenvectors of lowest rank k). Using stationary surrogate data, we controlled that these parameter fluctuations were not fully explained by estimation errors and, furthermore, that parameter fluctuations along sloppy dimensions were those that deviated the most from the stationary case (see Figure 4—figure supplement 1). Nevertheless, we noted that the projection variance into the stiff dimensions, albeit small, was not zero. This means that the model also evolved along parameter dimensions that had a strong impact on the collective activity. We hypothesized that changes in collective behavior, associated to changes in stiff parameters, were related to changes in cortical state.

Cortical state transitions evolve along stiff dimensions

To test this hypothesis, we first measured the cortical state in each epoch t using silence density, CS(t), defined as the fraction of 20-ms time bins with zero population activity, that is no spikes from any neuron (see Materials and methods) (Luczak et al., 2013; Pachitariu et al., 2015; Mochol et al., 2015). To obtain the most accurate estimate of silence density, we used all the spikes from the merge of all the single-units and multi-units in the calculation of CS(t). During the course of the experiment, we observed large fluctuations in silence density, with low and high values associated to desynchronized and synchronized cortical states, respectively (Figure 5A). We found that differences in collective dynamics in different epochs, quantified by DJS(t,t)=DJS(Pdata,t;Pdata,t), significantly co-varied with the changes in cortical state, given by d=|CS(t)CS(t)| (averaged correlation coefficient 0.56 ± 0.08, p < 0.001) (Figure 5B–C). Thus, changes in collective behavior correlated with changes in cortical state.

Figure 5 with 1 supplement see all
Spontaneous neuronal activity evolved along stiff dimensions during cortical state transitions.

(A) Silence density was used to characterize the cortical state. Green inset: low values of the silence density indicate desynchronized cortical states (each row represents the spike train of a single-unit). Red inset: high values of the silence density indicate synchronized cortical states. Data from rat 1. (B) Difference in collective pattern statistics, that is DJS, between different epochs, t and t', as a function of the corresponding difference in silence density, noted d. Each gray dot corresponds to a pair of epochs (t,t'). The solid line indicates the average relation between DJS and d; error bars indicate SD. Data from all neuronal ensembles from rat 1. (C) Correlation coefficient between DJS and d, for all rats. *: p < 0.001. (D) Top: Distribution of the absolute value of the correlations between the cortical state and the activity observables, noted rcs. Bottom: Distribution of parameter sensitivity values for biases (h parameters) and couplings (J parameters), for all models from all rats. (E) rcs vs. sensitivity s of all activity variables (i.e., firing rates and pairwise correlations). Correlation: rc = 0.36; p <0.001. Data from rat 4 (Npop = 72). (F) Correlation between rcs and the sensitivity for each dataset. *: p < 0.001. Error bars indicate correlation 95% confidence interval. (G) rcs for sloppy and stiff variables. *: p < 0.001, paired t-test. Error bars indicate SEM.

We next asked which activity observables, that is the firing rate of each neuron and all pairwise correlations, related more to cortical state transitions. For this, we calculated the absolute correlation, rcs, between the cortical state CS(t) and the activity observables. We found that rcs was broadly distributed between 0 and 0.94, thus some observables correlated more with the cortical state (Figure 5D, top panel). Next, to relate the sensitivity of model parameters (their stiffness) to the activity observables, we measured the sensitivity of a given parameter by its average contribution to the first eigenvector of the FIM and we associated it to the corresponding observable (Panas et al., 2015). We defined the sensitivity sne at the neuronal ensemble level and the sensitivity s at the population level (see Materials and methods). Note that the ranges of the sensitivity of biases (h) and couplings (J) were similar (Figure 5D, bottom panel), and that sensitivities calculated in the first and the second halves of the recording session were highly correlated (correlation coefficient > 0.82, for all rats; average: 0.89 ± 0.03). We found a significant positive correlation between the associated sensitivity (s) and the correlation with the cortical state (rcs) in 5/6 datasets (Figure 5E–F). Thus, the observables that correlated more with the cortical state were those with the highest associated sensitivity. This result led us to separate the activity observables into two classes, called 'sloppy' and 'stiff', based on whether the associated sensitivity (s) was lower or higher than the median of s. We found that stiff variables were significantly more correlated with the cortical state than the sloppy variables (p<0.01 for all datasets, paired t-test; Figure 5G). This relationship was preserved when using an alternative, more general definition of sensitivity that considered the contribution to all eigenvectors of the FIM, instead of the contribution to the first eigenvector only (see Figure 5—figure supplement 1). Altogether, these results indicate that neuronal activity and co-activity preferentially evolved along sensitive (stiff) parameter dimensions during cortical state transitions.

Sensory-evoked activity evolves along sloppy dimensions

The above results indicate that, although intrinsic spontaneous dynamics predominantly evolved along sloppy dimensions (Figure 4F), cortical state transitions were governed by changes in stiff parameters (Figure 5G). We next investigated which parameter dimensions were explored when the neural network was driven by external sensory inputs, that is during stimulus-evoked activity (Figure 6A). We observed that evoked responses (which could be increased or decreased with respect to pre-stimulus baseline firing rate) were larger for sloppy neurons than for stiff neurons (Figure 6B–C). To quantify the responsiveness of each neuron, we calculated the modulation index (MI, see Materials and methods) of each neuron in response to acoustic stimuli. We next calculated the relation between MI, calculated during evoked activity, and the sensitivity s associated to firing rates, calculated during the spontaneous activity as above. We found that the more responsive neurons were those with the lowest associated sensitivity (Figure 6D–E). This indicates that stimulus-evoked neuronal activity evolved mostly along sloppy dimensions. This result was replicated when using a more general definition of sensitivity that considered the contribution to all eigenvectors of the FIM (see Figure 6—figure supplement 1A). Finally, we evaluated the difference, noted ΔMI, between the MI of sloppy and stiff neurons as a function of cortical state CS(t). Specifically, first, the MI values in each epoch were averaged according to different ranges of the silence density. Second, the MI values of sloppy and stiff neurons were compared within each range. We found that ΔMI was maximal during desynchronized activity, and minimal during synchronized activity (Figure 6F, see also Figure 6—figure supplement 1B). Thus, the cortical activity during stimulus response evolved predominantly along sloppy dimensions for the desynchronized cortical state, while, in the synchronized state, the dominance of sloppy fluctuations was reduced, and stiff fluctuations became comparable.

Figure 6 with 1 supplement see all
Stimulus-evoked neuronal activity evolved along sloppy dimensions.

(A) Population responses to acoustic clicks. (B–C) Median-split of sensitivity s was used to separate stiff neurons and sloppy neurons. The mean responses for stiff and sloppy neurons are shown in the case of excited responses (B) and suppressed responses (C). The responses were normalized by the average pre-stimulus activity r0. Shaded areas correspond to SEM. Data in (A), (B), and (C) correspond to one example rat (#1). (D) Modulation index (MI) as a function of associated sensitivity of firing rates (si, with 1 iNpop), for each dataset. Each dot corresponds to a single neuron of the recorded population. The correlation between MI and sensitivity was negative for all datasets (rc: correlation coefficient; p: p-value). Solid lines indicate exponential fits. (E) Correlation between MI(t), calculated in epoch t, and νt,1, for all neuronal ensembles of each of the rats. On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. Asterisks indicate significantly negative medians (p < 0.001, two-sided signed rank test). (F) Difference of the MI of sloppy neurons minus the MI of stiff neurons as a function of cortical state (i.e., silence density), averaged for all rats (black trace; error bars indicate SEM). The gray bars indicate the distribution of silence density values.

Finally, correlations between cortical state fluctuations and sensitivity and between MI and sensitivity could reflect a dependency between sensitivity and model estimation errors. To test this, we evaluated the mean error on the model estimation of the observables and test its interaction with sensitivity, cortical state fluctuations, and MI (see Appendix 1 and Appendix 1—figure 1). We found that model estimation errors correlated with sensitivity, but they could not fully explain neither the positive correlation between sensitivity and cortical state nor the negative correlation between sensitivity and MI.

Stiff parameters were associated to central neurons within the neuronal network

In this section, we further investigate the properties of neurons and pairs of neurons with respect to their associated parameter sensitivity. As above, we separated the neurons and pairs of neurons into two classes, called 'sloppy units/pairs' and 'stiff units/pairs', based on whether the associated sensitivity (s) was lower or higher than the median s (units were associated to parameters hi, and pairs or links were associated to parameters Jij). With this dichotomization, we found that stiff units were significantly more active than sloppy units (Figure 7A). We quantified this by performing receiver operating characteristic (ROC) analysis and used the area under the ROC curve (AUC) as a measure of how well the firing rates distributions of the two classes were separated (AUC = 0.961–0.998, p < 0.001, for all rats; Figure 7G). Stiff neurons were also significantly more correlated among them than sloppy neurons (AUC = 0.615–0.932, p < 0.001, for all rats; Figure 7B,G). The distributions of correlations remained well separated when calculated for the links, that is pairs of neurons with associated parameters Jij (AUC = 0.541–0.766, p < 0.001, for all rats; Figure 7C,G).

Figure 7 with 1 supplement see all
Central and highly active neurons were associated to stiff parameters.

(A) Distribution of firing rates of sloppy and stiff neurons. (B) Distribution of correlations among sloppy neurons and among stiff neurons. (C) The distribution of correlations was also calculated for the links (i.e, pairs of neurons with associated parameters Jij). Note that, in principle, links can be related to pairs composed of one sloppy and one stiff neuron. (D) Distribution of betweenness centrality of sloppy and stiff neurons. (E) Connectivity graph: each node represents a neuron and links represent significant correlations between pairs of neurons. The graph was plotted using force-directed layout, that is using attractive forces between strongly connected nodes and repulsive forces between weakly connected nodes. Left: the nodes were colored as a function of betweenness centrality. Right: the nodes were colored as a function of associated sensitivity s. Note the high overlap between both color labeling methods, indicating that sensitivity was highly predictive of the centrality of the nodes. (F) Distribution of neuron-to-population couplings of sloppy and stiff units. Panels A–F show data from rat 1. (G) Area under the receiver operating curve (AUC) quantifying the separation of distributions of sloppy and stiff classes. All AUC values were significantly higher than 0.5 (p < 0.001).

To further investigate the structure of correlations, we evaluated the centrality of stiff and sloppy neurons within the observed network of neurons. For this we used the betweenness centrality (BC), a measure of node centrality in a graph or network, which in our case was given by the functional connectivity matrix among the recorded neurons (see Materials and methods). The BC measures the extent to which a node in the graph tends to lay on the shortest path between other nodes. Thus, a node with higher BC has more influence over the network, because more information passes through that node. We found that stiff neurons had significantly more centrality in the functional connectivity graph than sloppy neurons (AUC = 0.740–0.831, p<0.001, for all rats; Figure 7D,G). This indicates that stiff neurons were part of the core of the graph, while sloppy neurons were part of the graph periphery, as clearly shown using graph visualization (Figure 7E; Fruchterman and Reingold, 1991). BC values were correlated with firing rates (correlation coefficient: 0.59 ± 0.11), which could suggest that differences in BC between stiff and sloppy neurons were simply a consequence of differences in firing rates. However, using surrogate data that preserved the observed firing rates and produced correlations through global modulations, we found that neither the structure of correlations nor the BC values could be trivially predicted by globally modulated firing rates but they were rather suggestive of functional interactions (see Appendix 1 and Appendix 1—figure 2). Thus, in addition to different firing rates, different correlations and BC values were supplementary features of stiff and sloppy neurons.

Moreover, previous work has shown that cortical neurons differ in their coupling to the population activity, with neurons that activate most often when many others are active and neurons that tend to activate more frequently when others are silent (Okun et al., 2015). Thus, along with centrality, we calculated the neuron-to-population coupling, given by the Pearson correlation between the activity of each neuron i and the number of coactive neurons (excluding neuron i; see Materials and methods). We found that stiff neurons were significantly more coupled to the population activity than sloppy neurons (AUC = 0.603–0.939, p < 0.001, for all rats; Figure 7F,G). In summary, stiff units were more active, more central, more coupled among them, and more coupled to the population activity than sloppy units. The same results were found when using a more general definition of sensitivity that considered the contribution to all eigenvectors of the FIM (Figure 7—figure supplement 1).

Discussion

We here studied the changes in activity caused by intrinsic (i.e. cortical state) and extrinsic (i.e., stimulus-evoked) sources in A1 neuronal ensembles in an estimated parameter space. The parameter space was obtained using the maximum entropy principle, providing a handful number of parameters describing the probability of all possible binary activity patterns. These parameters differed in their impact on collective activity that was sensitive to a few combinations of parameters, called stiff dimensions, but insensitive to many others called sloppy dimensions. Our results suggest that spontaneous cortical state transitions and stimulus-driven activity evolved along different parameter dimensions. Indeed, in one hand, while most of the fluctuations during spontaneous activity evolved along sloppy dimensions, some residual ongoing fluctuations evolved along stiff dimensions, and these fluctuations were correlated with synchronized/desynchronized cortical state transitions. On the other hand, stimulus-induced activity was larger in sloppy dimensions than in stiff dimensions, an effect that was most prominent during the desynchronized cortical state. Note that the observation that both spontaneous and stimulus-driven activities predominantly evolve along sloppy dimensions results from the strong similarity of spontaneous and evoked activity, reported in several previous studies (Arieli et al., 1996; Kenet et al., 2003; MacLean et al., 2005; Luczak et al., 2009). Finally, by classifying the neurons as stiff versus sloppy neurons (i.e., those contributing more or less to the principal stiff dimension) we found that the firing rates and the functional connectivity topology significantly differed between the two classes of neurons. It should be noted, however, that, since sensitivity is a continuous variable, the two classes of neurons that we defined here do not represent two disjoint groups but rather represent two parts of a continuum.

The observation that stimulus-induced activity evolved along sloppy dimensions can have important functional implications. It suggests that a stimulus can modulate the activity of a subset of sloppy neurons without entirely affecting the collective activity. This could be an efficient functional architecture to encode sensory information without perturbing other ongoing or memory-stored processes. Consistent with this view and with previous studies (Margolis et al., 2012; Mizuseki and Buzsáki, 2013; Panas et al., 2015), our results suggest that the integrity of the network is ensured by a core of highly active stiff neurons, which have strong functional connections among them (either through anatomical connections or common inputs), while topologically peripheral sloppy neurons (within the functional connectivity graph) can be largely modulated by external inputs. A similar sub-network of highly active, interconnected neurons has been recently identified in the mice neocortex (Yassin et al., 2010). Importantly, sensory input was not required to drive these cells. Previous studies of complex systems have derived general principles of core/periphery network structures: the network periphery is more variable, evolvable, and plastic than the network core, while the network core facilitates system robustness (Kitano, 2004; Csermely et al., 2013). Thus, we hypothesize that sloppy neurons could also be more affected by synaptic plasticity, allowing for network reconfiguration without loss of stability. Consistent with this, previous work on whole-brain fMRI has observed core stability and peripheral flexibility over the course of learning (Bassett et al., 2013), and recent analyses of functional networks from calcium imaging data recorded in mouse primary auditory cortex revealed a stable core and a variable periphery over time (Betzel et al., 2019). Furthermore, we observed that stimulus responses evolved more pronouncedly along sloppy dimensions in the desynchronized state, while in the synchronized state fluctuations along sloppy and stiff dimensions were comparable (Figure 6F). This supports the view that responses along sloppy dimensions provide information processing benefits, since previous studies have shown that auditory stimuli in rodents (Marguet and Harris, 2011; Pachitariu et al., 2015) and visual stimuli in both rats (Goard and Dan, 2009) and monkeys (Beaman et al., 2017) are better represented in the desynchronized state as compared to the synchronized state.

The properties of spontaneous and induced cortical dynamics observed in the present anesthetized condition are likely to be relevant also during wakefulness. Indeed, several studies reported the existence of synchronized cortical states during wakefulness (for review see Zagha and McCormick, 2014), and global fluctuation resembling transitions between up and down periods during alert or quiescent wakefulness (Petersen et al., 2003; Luczak et al., 2007; Poulet and Petersen, 2008; Zagha et al., 2013; Tan et al., 2014; Engel et al., 2016) or even during task engagement (Sachidhanandam et al., 2013). Moreover, sloppiness has been observed in asynchronous spontaneous activity under light anesthesia (Panas et al., 2015), we thus expect to observe a similar stiff-sloppy architecture in the awake state. However, we believe that the comparison of Fisher information matrices during wakefulness and during different levels of anesthesia could provide valuable information about the principles governing vigilance.

We found that stiff neurons were more linked to the observed neuronal population activity than sloppy neurons. Stiff neurons had higher centrality in the functional connectivity graph and higher coupling to the population activity than sloppy neurons. Previous research showed that neurons differ in their coupling to the population activity, with neurons that activate most often when many others are active, called ‘choristers’, and neurons that tend to activate more frequently when others are silent, called ‘soloists’ (Okun et al., 2015). Our results suggest that stiff and sloppy neurons are chorister and soloist neurons, respectively. In other words, changes in the activity of stiff/chorister neurons lead to changes in collective behavior (i.e., cortical states), while the activity of sloppy/soloist neurons can spontaneously fluctuate or respond to stimuli without strongly affecting the collective behavior. Thus, we believe that the roles of stiff/chorister neurons and sloppy/soloist neurons are important to understand tradeoffs between responsiveness and stability of the network. Furthermore, we here studied the evolution of neuronal activity on the time scale of hours and found that fluctuations on stiff parameter dimensions were the weakest and were related to cortical state transitions, which time scale is in the order of tens of minutes (Hahn et al., 2017; Mochol et al., 2015). Previous studies have reported prominent changes on neuronal activity and tuning properties over days, but with stable decoding performances of population activity (Chestek et al., 2007; Ziv et al., 2013; Panas et al., 2015). However, we hypothesize that learning or adaptation to changing environments could lead to large changes in collective activity. In that case, particular attention could be paid to the influence of high-order areas on the activity of subsets of stiff and sloppy neurons from sensory areas, as top-down regulation might be a mechanism to control the stabilizing network core.

The existence of cortical neurons with different sensitivities (from sloppy to stiff neurons) provides new valuable architectural constrains for models of the brain state and its transitions. Several past studies have modeled the synchronized brain dynamics as transitions between two attractors. Depending on the model specificity those transitions could be noise driven (Mejias et al., 2010; Mochol et al., 2015; Jercog et al., 2017) or caused by some fatigue mechanism (Compte et al., 2003; Hill and Tononi, 2005; Mattia and Sanchez-Vives, 2012). To make the system works in a desynchronized regime it was enough to increase the background input to the network (Bazhenov et al., 2002; Hill and Tononi, 2005; Curto et al., 2009; Destexhe, 2009; Mochol et al., 2015). Given our present results, the models could be extended to include a network core/periphery architecture, a non-homogeneous background input preferentially targeting the network core, and different stimulus spatial distributions. Such a model would provide insights on the interplay between cortical state transitions and sensory representation. Moreover, our findings question the view that the mechanisms by which background and stimulus inputs impact the dynamics are similar, as assumed in the simple bi-stable rate model (Mochol et al., 2015).

Finally, we here described the patterned activity of small (N = 10) neuronal ensembles using MEMs. It is known that MEMs of small sizes can present departures from the observed distribution of summed activities and higher-order correlations (Tkačik et al., 2014). Recent advancements on learning algorithms allow to construct MEMs of ∼100 neurons. However, these models cannot be used in a time-resolved manner, as we did here, due to limited data in each epoch. Small model sizes are thus the cost to pay to study the evolution of collective activity over time in a meaningful time scale (i.e., the one of cortical state transitions).

Materials and methods

Key resources table
Reagent type (species)
or resource
DesignationSource or referenceIdentifiersAdditional
information
Biological sample (Sprague–Dawley rat)Sprague–Dawley rathttps://doi.org/10.1073/pnas.1410509112six rats, 250–400 g
Software, algorithmMatlabMathWorksRRID:SCR_001622All analyses
Software, algorithmKlustakwikhttp://klustakwik.sourceforge.net/RRID:SCR_014480Spike sorting (detection and initial clustering)
Software, algorithmEToShttp://etos.sourceforge.net/Spike sorting (detection and initial clustering)
Software, algorithmKlustershttp://neurosuite.sourceforge.net/Spike sorting (clustering)

Ethics statement

Request a detailed protocol

All experiments were carried out in accordance with protocols approved by the Animal Ethics Committee of the University of Barcelona (Comité d’Experimentació Animal, Universitat de Barcelona, Ref 116/13).

Experimental techniques

Request a detailed protocol

We analyzed the neuronal activity recorded in the primary auditory cortex (A1) of 6 anesthetized rats (Sprague–Dawley; 250–400 g). The experimental procedures and spikes sorting procedures have been previously described in Mochol et al. (2015). Briefly, rats were anesthetized with urethane (1.5 g/kg body weight) and silicon microelectrodes (Neuronexus) with 32 or 64 channels were inserted in deep layers (depth, 600–1,200 μm) of the primary auditory cortex. The spiking activity from single units and multi-units (i.e., neurons that were not well isolated) was simultaneously recorded during spontaneous activity and in response to acoustic ‘clicks’ (5 ms square pulses; interstimulus interval, 2.5 or 3.5 s; see Table 1). In some datasets, double clicks (5 ms square pulses; 50- or 100 ms inter-click interval) were also presented, but, in the present study, we analyzed only the responses to single click. The spiking data is publicly available here: https://github.com/adrianponce/Spont_stim_spiking_A1.

Table 1
Number of neurons (SU: single-units, MU: multi-unit), number of 100 s epochs, number of stimulus presentations in 100 s epochs, and number of neuronal ensembles, for each dataset.
No. of neuronsNo. of 100 s epochsStimulus presentations in 100 s epochsNo. of neuronal ensembles (Q)
Rat 1SU: 81; MU: 316312–1420
Rat 2SU: 147; MU: 137412–1420
Rat 3SU: 44; MU: 307012–2010
Rat 4SU: 72; MU: 1035910–2020
Rat 5SU: 58; MU: 392928–2910
Rat 6SU: 112; MU: 832817–2920

Cortical state

Request a detailed protocol

Long continuous recordings (mean, ∼2 h) were divided into NE 100-s epochs, and cortical state was estimated in each epoch based on spontaneous pooled population activity, that is the merge of single and multiunit spike trains during the 1.5-s intervals preceding each stimulus presentation. Cortical state was quantified using silence density defined as the fraction of 20-ms time bins with no population activity. Silent and active periods were obtained from the merge of consecutive empty and nonempty bins, respectively.

Maximum entropy models

Request a detailed protocol

The spontaneous spiking activity of ensembles of N single neurons was studied using statistical modeling based on maximum entropy principle. The ensemble activity was binarized in non-overlapping time bins of dt = 10 ms, during which neuron i either did (σi=+1) or did not (σi=-1) generate one or more spikes. The state of the neural ensemble is described by a binary pattern σ=σ1, σ2, , σN, and thus the collective activity is described by the probability distribution Pσ over all 2N possible binary patterns. We estimated Pσ using a Maximum entropy model (MEM). The MEM finds Pσ by maximizing its entropy under the constraint that some empirical statistics are preserved. A pairwise-MEM provides a solution under the constraint that the activation rates (<σi>) and the pairwise correlations (<σiσj>) are preserved. The maximum entropy distribution Pσ that is consistent with these expectation values is given by the Boltzmann distribution (Schneidman et al., 2006; Tkačik et al., 2015):

(1) P(σ)=eE( σ ){σ}eE( σ ),

where Eσ is the energy of the pattern σ, given by:

(2) E(σ)=i=1Nhiσi12i=1Nj=1NJijσiσj,

and Z={σ}eE( σ ) is the partition function.

The model parameter hi, called intrinsic bias, represents the intrinsic tendency of neuron i towards activation (σi=+1) or silence (σi=-1) and the parameter Jij represents the effective interaction between neurons i and j. The estimation of the model parameters Ω={h,J} was achieved through a gradient descent algorithm (see below). For each recording session, we constructed models for Q ensembles of N=10 randomly selected single neurons from the entire population of Npop single neurons and learned the model parameters using the spontaneous binarized activity within each 100-s epoch. Thus, for each recording session, we built Q×NE models of 10 units. We were interested on the evolution of the model parameters over time, that is Ω(t). Note that, for a given model, the number of free parameters is the sum of intrinsic biases and effective couplings, N+N(N-1)/2= 55, that is Ω=[h1, h2, , hN,J12, J13,].

Estimation of MEM parameters

Request a detailed protocol

The MEM parameters Ω={h,J} were iteratively adjusted to minimize the absolute difference between the empirical activation rates (σi) and correlations (σiσj) and those (σimodel, σiσjmodel) predicted by the model through Metropolis Monte Carlo simulations (100,000 samples). Specifically, each iteration is given by: hinew=hiold-ασimodel-σi, and Jijnew=Jijold-ασiσjmodel-σiσj, where α is the learning rate (α= 0.1). In our study we stopped the re-estimations once the differences between the empirical and model values are less than a tolerance threshold (0.005) or if this tolerance was not reached within a maximum number of iterations (100).

MEM goodness-of-fit

Request a detailed protocol

The goodness-of-fit of the MEMs was evaluated using the Jensen–Shannon divergence (DJS) between the probability distribution of the empirical and model binary patterns (Marre et al., 2009). DJS is a symmetric version of the Kullback-Leibler divergence (DKL) and is given as:

(3) DJS(Pdata;PMEM)=12DKL[Pdata;(Pdata+PMEM)2]+12DKL[PMEM;(Pdata+PMEM)2],

Where PMEM was given by the Boltzmann distribution of the model, Pdata was estimated from the N-dimensional binary patterns observed in the data, and:

(4) DKL(P1;P2)={x}P1(x)logP1(x)P2(x).

The fitting of MEM (second-order model) was compared to the fit obtained using independent-MEM, that is in which only for which only the activation rates (<σi>) are preserved (i.e., only h is optimized; first-order model). In this case, the pattern energy is given by: Eσ=-i=1Nhiσi.

Furthermore, the performance of the model can be evaluated using the Kullback-Leibler ratio, R (Shlens et al., 2009). This ratio is given by comparing the Kullback-Leibler divergence between the distribution P1 of the first-order model (i.e., independent-MEM) and the distribution of the actual data, D1=DKLP1;Pdata, with the Kullback-Leibler divergence between the distribution P2 of the second-order model and the distribution of the actual data, D2=DKLP2;Pdata. Specifically, the Kullback-Leibler ratio is defined as:

(5) R=D1D2D1.

This ratio can range between 0 and 1, with one giving the highest performance.

Fisher information matrix

Request a detailed protocol

Because in the MEM, all the information about the collective activity is contained in the probability distribution of the binary patterns, Pσ, one can define the model parameter space as P(σ|Ω). We were interested in knowing which parameters, or combination of parameters, have a strong effect on the collective activity. To measure how distinguishable two models, with parameters Ω and Ω+δΩ, are based on their predictions, we used the Fisher information matrix (FIM). Indeed, the Kullback-Leibler divergence between the two models can be written as:

(6) DKL(Ω;Ω+δΩ)=FIMklδΩkδΩl+𝒪(δΩ3),

where 1 k,l 55, and the matrix FIM is given by:

(7) FIMkl={σ}P(σ|Ω)logP(σ|Ω)ΩklogP(σ|Ω)Ωl.

The FIM represents the curvature of the log-likelihood of the model, logP(σ|Ω), with respect to the model parameters. It quantifies the sensitivity of the model to changes in parameters. By calculating the eigenvalues of the FIM, we can determine which combinations of parameters affect the most the model's behavior.

In the case of MEM, the FIM can be easily obtained by using Equations 1, 2, and 7. As a result, the FIM is given by the covariance matrix of observables associated to the parameters which can be calculated from the model through Metropolis Monte Carlo simulations (500,000 steps), that is:

(8) FIMkl=xkxlxkxl,

with 1 k,l 55 and x=[σ1, σ2, , σN, σ1σ2, σ1σ3,].

Sensitivity measures

Request a detailed protocol

The FIM was calculated for every neuronal ensemble at every 100-s epoch and it was decomposed into eigenvectors, noted νt,1, νt,2, ,νt,k, , where k is the rank of the eigenvector and t denotes the epoch. Following Panas et al. (2015), within each neuronal ensemble, we measured the sensitivity of a given parameter by its averaged contribution to the first eigenvector of the FIM, that is the sensitivity of i-th parameter is given by sne,i=1NEtνt,1(i), with 1iN.

We next constructed a sensitivity measure for the entire population of Npop neurons. For this, we defined the set of all single neuron indices and all pairs of neurons I=1,,Npop,1,2,1,3, . This set has L=Npop+Npop(Npop-1)2 elements. For each element j of I, we defined the sensitivity sj as the average of sen,i over the neuronal ensembles that contained the j-th single neuron or the pair of neurons (i.e., those neuronal ensembles for which i maps to j). In other words, sne denotes the sensitivity within an ensemble of N=10 neurons and has 55 elements, and s denotes the sensitivity within the entire population of Npop neurons and has L elements. This allows comparison of s with statistics derived from the population of Npop neurons.

Parameters that contributed less to the first eigenvector could in principle contribute to the other stiff dimensions (those with lower rank k, e.g., k= 2). For this reason, we also considered an alternative definition of sensitivity that considers the weighted contribution to all eigenvectors of the FIM. For each neuronal ensemble and each 100-s epoch t, we defined the weighted sensitivity of the parameter i as the temporal average of its contribution to the eigenvectors of the FIM, weighted by the associated eigenvalues (at,1, , at,55):

(9) sne,iw=1NEt=1NEk=155at,k|νt,k(i)|at,1+at,2+ +at,55.

As previous, from snew one can construct a weighted sensitivity sw at the population level.

Finally, we separated the activity observables into two classes, called “sloppy” and “stiff”, based on whether the associated sensitivity s was lower or higher than the median sensitivity.

Similarity measures

Request a detailed protocol

Temporal variations of model parameters and data statistics were quantified using the average correlation between the parameters/statistics at time t and the parameters/statistics at time t+Δt. For example, let r(t) the average firing rates of the neurons during the epoch t, the similarity measure is given by:

(10) γ(Δt)=1NEΔtt=1NEΔtρ[r(t),r(t+Δt)],

Where NE is the number of epochs and ρ is the Pearson correlation coefficient. In the case of FIM, the matrix was vectorized to calculate ρ.

To evaluated how well the data in a given epoch t could be explained by the MEM constructed using the data at time t+Δt. Specifically, we defined the similarity measure DJSΔt, given by the average Jensen-Shannon divergence between the distribution of data binary patterns in epoch t, that is Pdata, t, and the distribution of binary patterns predicted by the MEM constructed using the data in epoch t+Δt, that is PMEM,t+Δt. This measure is given as:

(11) DJS(Δt)= 1NEΔtt=1NEΔtDJS(Pdata,t;PMEM,t+Δt).

In other words, 1/DJSΔt quantifies how well, on average, the model with parameters Ω(t+Δt) represents the data from epoch t.

Modulation index

Request a detailed protocol

We quantified the responsiveness of the neurons to sensory stimuli through the modulation index (MI) defined as:

(12) MI= |rstimrspon|rstim+rspon,

where rspon is the pre-stimulus average spike count, calculated in the 0.5-s pre-stimulus interval, and rstim is the average spike count calculated from stimulus onset to 0.5 s after stimulus onset. With this definition, strongly increased or suppressed stimulus responses, with respect to pre-stimulus activity, lead to high MI values.

Betweenness centrality

Request a detailed protocol

For each recording session, we analyzed the network defined by the Pearson correlation matrix of the activities of all single units. The centrality of a neuron, or node, within the network was quantified using the betweenness centrality (BC) measure. BC is given by the number of shortest paths that pass through a given node. The correlation matrix was compute for all 100-s epochs and, for each matrix element, we tested whether the mean of the NE correlation values differs from 0 (t test followed by Bonferroni correction), resulting in a binary graph G with entries equal to 1 if correlation were significantly different from zero (corrected p-value < 0.05) and 0 otherwise. The BC for each node of the graph was given by:

(13) BC(i)=kilp(kl;i)p(kl)

where p(kl) is the total number of shortest paths from node k to node l and p(kl;i) is the number of those paths that pass through i.

Neuron-to-population coupling

Request a detailed protocol

To quantify the coupling of each neuron to the activity of the neuronal population, we calculated, for each epoch, the Pearson correlation between the activity of each neuron (σi) and the number of coactive neurons (i.e., with σi=+1) at each time bin (dt = 10 ms) from the neuronal population of single units (without including the neuron i). The neuron-to-population coupling was given by the average of the correlation coefficient across epochs.

ROC analysis

Request a detailed protocol

We used the receiver operating characteristic curve (ROC) to evaluate the separation between the distributions of observables from sloppy and stiff classes. Let Xsloppy and Xstiff be the sloppy variables, that is those variables with associated sensitivity (s) lower than the median s, and the stiff variables, that is those variables with associated sensitivity (s) higher than the median s, respectively. The ROC curve, f(c), is build by plotting the probability of P(Xsloppy>c) against the probability of P(Xstiff>c), for each all c. The area under the ROC curve (AUC) is a measure of separation between P(Xsloppy) and P(Xstiff), and it is given by:

(14) AUC= f(c)dc.

AUC ranges between 0 and 1, with AUC = 0 if P(Xsloppy) and P(Xstiff) are completely separated and Xsloppy>Xstiff, AUC = 1 if P(Xsloppy) and P(Xstiff) are completely separated and Xstiff>Xsloppy, and AUC = 0.5 if P(Xsloppy) and P(Xstiff) are undistinguishable. We used a permutation test (1000 re-samples), in which observables and classes were randomly associated, to assess AUC values that were significantly different from 0.5.

Stationary surrogates

Request a detailed protocol

To construct the stationary surrogates we first randomly selected a reference epoch t. Second, we generated binary data using the MEM estimated from the spiking data at this reference epoch, that isusing parameters Ω(t), through Monte Carlo simulations of the model to obtain 5000 binary patterns. Third, we repeated the Monte Carlo simulations NE times. Finally, for each of the NE pieces of surrogate data, we estimated new MEM parameters, Ω, using gradient descend, and we calculated the corresponding Fisher Information Matrix (FIM) using 500,000 Monte Carlo steps as described above. By construction, the obtained surrogate data were stationary and had the same length of the original spiking data. Thus, parameter fluctuations in the surrogate data were only due to model estimation errors.

Appendix 1

Estimation errors and sensitivity

Correlations between cortical state fluctuations and parameter sensitivity and between MI and sensitivity could reflect a dependency between sensitivity and model estimation errors. To test this, we here evaluated the error on the model estimation of the observables (firing rates and correlations) as a function of sensitivity. For each epoch t, the firing rates and covariances estimated by the MEM can be calculated from the model parameters by means of the Boltzmann distribution. The estimated moments σi and σiσj are given by:

(15) σi={σ}P(σ)σi,
(16) σiσj={σ}P(σ)σiσj.

The estimated firing rates and the covariances are given as ri=σi+1/(2dt) and Cij=σiσj-σiσj (Appendix 1—figure 1A-B). We found that the firing rates and covariances estimated from the data and those estimated by the MEM highly correlated (average correlation coefficient for rates: 0.999 ± 0.001; for covariances: 0.956 ± 0.018, Appendix 1—figure 1C), indicating acceptable fits of the models. Indeed, the model performances, evaluated through the Kullback-Leibler ratio R, were also high (average equal to 0.95 ± 0.03, Appendix 1—figure 1D).

Appendix 1—figure 1
Estimation errors and sensitivity.

(A-B) Comparison between the firing rates and covariances measured in the data and those estimated by the MEM. Data from rat 4. Each point represents the activity of a neuron (in panel A) or co-activity of a pair of neurons (in panel B) from a given neuronal ensemble and at a given epoch. The data of all neuronal ensembles and all epochs are presented. Insets show the resulting differences between statistics estimated from the data and from the model. (C) Correlation between the firing rates and covariances measured in the data and those estimated by the MEM, for each dataset. (D) Average Kullback-Leibler ratios for all learned models of each dataset. Error bars indicate SEM. (E) Correlation between sensitivity s and estimation errors Err. (F) Full squares: correlation between sensitivity and rcs (same values as in Figure 5F). Open circles: correlation between estimation errors and rcs. For each dataset the correlation coefficients were compared, given the observed correlations in panel E using Meng's z-test for dependent correlations, *: p < 0.01. (G) Full squares: correlation between sensitivity of firing rates and MI (same values as in Figure 6D). Open circles: correlation between firing rate estimation errors and MI. *: p < 0.01, Meng's z-test for dependent correlations. In panels (C) and (EG), error bars indicate correlation 95% confidence interval.

We next defined the estimation errors as the averaged squared difference between the values estimated from the data and those estimated by the model:

(17) EE[ri]=1NEt=1NEΔrt,iσΔr,
(18) EE[Cij]=1NEt=1NEΔCt,i,jσΔC,

where Δrt,i=(ridata,trimodel,t)2 and ΔCt,i,j=(Cijdata,tCijmodel,t)2, and σΔr and σΔr denote the standard deviations of Δrt,i and ΔCt,i,j, respectively —this normalization is needed because Δrt,i and ΔCt,i,j had different magnitudes. Using the same procedure as for sensitivity s, we obtained, for each dataset, a variable Err that represented the mean estimation error of the Npop neurons and all pairwise combinations.

Since estimation errors correlated with sensitivity (average: 0.192 ± 0.055; Appendix 1—figure 1E), we asked whether the correlations between sensitivity and cortical state observed in Figure 5F could be fully predicted by correlations between errors and cortical state. We found that correlations between sensitivity s and rcs were higher than correlations between estimation errors and rcs (0.252 ± 0.056 vs. 0.168 ± 0.017 on average; Appendix 1—figure 1F). We used Meng's z-test for dependent correlations to compare the correlations between sensitivity s and rcs and those between estimation errors and rcs, given the correlation between s and errors, and found significant differences in 4/6 datasets (p < 0.05). This was the case for datasets with the highest model performances (datasets 1–4).

Finally, we tested whether estimation errors for firing rates correlated with the modulation index (MI) of the neurons and compared the results with those obtained in Figure 6D. We found that correlations between estimation errors and MI were substantially weaker than correlations between estimation errors and MI (−0.506 ± 0.043 vs. −0.086 ± 0.017 on average; p<0.05 for all datasets, Meng's z-test for dependent correlations; Appendix 1—figure 1G). We concluded that, although estimation errors correlated with sensitivity, they did not fully explain neither the correlation between sensitivity and cortical state (for the datasets with highest model performance) nor the correlation between sensitivity and MI.

Non-homogeneous Poisson process surrogate data

In Figure 7A–D we found that stiff neurons were more active, more coupled among them, and more central than sloppy neurons. We here investigated the possibility that correlations and BC values could be trivially predicted by firing rates. We tested whether the structure of correlation could arise from spontaneous global fluctuations increasing or decreasing the firing of the neurons, instead of functional connectivity among them. For this, we constructed surrogate data that preserved the mean firing rate of the neurons, and the number and duration of silent (i.e., no spikes, empty 20 ms time bins) and active (i.e., non-empty 20 ms time bins) periods in each epoch, but without other structured correlations among the units.

To do this, we generated non-homogeneous Poisson process (NHPP) surrogate data as follows. We first calculated the average spontaneous firing rates of the single neurons during active periods in the original data, re=[re,1,re,2,,re,Npop], given by the number of spikes, in a given epoch e, divided by the total duration of active periods within this epoch. During silent periods the firing rate of all single units was zero. Next, for each unit i and within each epoch e, we generated non-homogeneous Poisson spike trains using the estimated firing rate as the intensity λe,i(t) of the point process, that is λe,i(t) was a step function with λe,i(t)=re,i, for times t within active periods, and λe,i(t)=0, otherwise. The resulting surrogate data preserved the firing rate of each neuron and silence periods of population activity (Appendix 1—figure 2A,E).

Appendix 1—figure 2
Correlations and centrality in non-homogeneous Poisson process surrogate data.

(A) Relation between firing rates in the original data and firing rates in the non-homogeneous Poisson process (NHPP) surrogates. R2 indicates the fraction of explained variance. The black line indicates the identity line. (B) Relation between the mean pairwise correlation in the original data and that in NHPP surrogates. Each point represents the mean correlation (averaged over neurons) in a given epoch. (C) Relation between pairwise correlations, averaged over the NE epochs, in the original data and those in NHPP surrogates. Each point represents a pair of neurons. (D) Relation between BC values in the original data and those in NHPP surrogates (see Materials and methods for calculation of BC). Data in (A), (B), (C), and (D) correspond to one example rat. (E) Explained variance for firing rates, correlations, and BC values, averaged over datasets. Error bars indicate SEM. (F) Correlation between firing rates and BC values, averaged over datasets. Error bars indicate SEM.

Despite being independent during active periods, the units in the NHPP data were correlated due to common global fluctuations across epochs and silent and active periods. Indeed, the mean correlation values in each epoch, Ce, from the original and the surrogate data were highly correlated (Appendix 1—figure 2B,E), although Ce was systematically lower in the surrogate data than in the original data. We tested how well this simple model predicted the structure of correlations, by calculating the pairwise correlations, Cij, and the BC values in the surrogate data, and comparing them to the values obtained in the original data (Appendix 1—figure 2C-D). We found that neither Cij nor BC were predicted by the NHPP: the average fraction of explained variance was equal to R2 = 0.148 ± 0.07 and R2 = 0.184 ± 0.07 for Cij and BC, respectively (Appendix 1—figure 2E). Moreover, the correlation between firing rates and BC values was substantially reduced in the NHPP, with correlation coefficients equal to 0.59 ± 0.11 and 0.19 ± 0.15 for the original and the NHPP, respectively (p < 0.01, t-test). We concluded that correlations and centrality measures were not fully predicted by globally modulated firing rates but were rather suggestive of functional interactions.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47

Decision letter

  1. Michael J Frank
    Senior Editor; Brown University, United States
  2. Brice Bathellier
    Reviewing Editor; CNRS, France
  3. Brice Bathellier
    Reviewer; CNRS, France
  4. Matthias H Hennig
    Reviewer
  5. Mark D Humphries
    Reviewer; University of Nottingham, United Kingdom

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This study shows that auditory cortex dynamics can be decomposed into a small set of critical dimensions and a set of dimensions that rather account for small variations of the dynamics. In addition, remarkably, critical dimensions are not shared between spontaneous and evoked activity, suggesting that the two types of activity come from different processes. This is an exciting result that should motivate important refinements in current models of auditory cortex dynamics.

Decision letter after peer review:

Thank you for submitting your article "Cortical state transitions and stimulus response evolve along stiff and sloppy parameter dimensions, respectively" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Brice Bathellier as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Michael Frank as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Matthias H. Hennig (Reviewer #2); Mark D Humphries (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

This paper uses model-based analysis of extracellular recordings from the rat auditory cortex to provide evidence for the hypothesis that the stability of the population in cortical networks depends on a small number of neurons with highly sensitive parameters, while parameters associated with the majority of neurons are insensitive to changes. Moreover, the authors show that transitions between desynchronised and synchronised states are associated with changes in the sensitive parameters, while during sensory stimulation this is state-dependent. Finally, it is shown that neurons with sensitive parameters have hub-like properties in the network.

Overall these are very interesting and in fact surprising results. The paper is well written and presented, but a few complementary analyses/clarifications are needed to make the important but complex message of this study more accessible and convincing.

Essential revisions:

1) The central conclusions in this paper are based on analysis of the Fisher information matrix (FIM) for pairwise maximum entropy models fit to binarised multi-neuron patterns. One possible caveat here is that the models were fit on relatively small data sets, which could affect the analysis downstream as the insensitive directions in the FIM could just reflect poor parameter estimates. The around 5000 patterns used here to fit the model are very few, and estimates of correlations used to constrain the model may be poor, given cortical activity is quite sparse. Possible biases should be assessed more thoroughly. In particular it is important to check for the impact of estimation errors given the short epochs. One key measure to check is the sensitivity measure, as this defines both core results (spontaneous = stiff; evoked = sloppy). The authors should provide an estimate of the variation in "sensitivity" as a function of model estimation error.

Most importantly:

– Figure 5E-G, check for how correlations between sensitivity and "state" change with errors in model estimation.

– Figure 6D-E, check for how correlations between sensitivity and the modulation index change with errors in model estimation.

2) It is hard to follow what is 'stiff' and what is 'sloppy', as in most of the cases the frontier between the two categories is obviously not sharp. It is clear that there are two ends of the parameter distribution (one end is more stiff and the other more sloppy), but in between, there are many intermediates. The authors must be more explicit about this and make more precise wherever possible what they mean when they define these two categories.

– In Figure 4B they should display the boundary between stiff and sloppy with a shading or a vertical line. Also, there is a kink in the eigenvalue curve. Is it related to model/population size? Are the eigenvalues beyond this kink relevant?

– It should be explained why stiffness of individual neurons is chosen to be their contribution to the first eigenvalue only. Are further eigenvalues not relevant?

– Together it should be clarified in a dedicated paragraph that stiff and sloppy parameters/neurons are two parts of a wide distribution.

https://doi.org/10.7554/eLife.53268.sa1

Author response

Essential revisions:

1) The central conclusions in this paper are based on analysis of the Fisher information matrix (FIM) for pairwise maximum entropy models fit to binarised multi-neuron patterns. One possible caveat here is that the models were fit on relatively small data sets, which could affect the analysis downstream as the insensitive directions in the FIM could just reflect poor parameter estimates. The around 5000 patterns used here to fit the model are very few, and estimates of correlations used to constrain the model may be poor, given cortical activity is quite sparse. Possible biases should be assessed more thoroughly. In particular it is important to check for the impact of estimation errors given the short epochs. One key measure to check is the sensitivity measure, as this defines both core results (spontaneous = stiff; evoked = sloppy). The authors should provide an estimate of the variation in "sensitivity" as a function of model estimation error.

As suggested by the reviewers, in the revised manuscript, we now evaluated the correlation between sensitivity and the error on the model estimation of the firing rates and covariances. This analysis is mentioned in the main article in subsection “Sensory-evoked activity evolves along sloppy dimensions” and described in detail in the Appendix 1 and in Appendix—figure 1.

The estimated firing rates and covariances can be calculated from the model parameters by means of the model's Boltzmann distribution. The firing rates and covariances estimated from the data and those estimated by the MEM highly correlated (average correlation: 0.999 ± 0.001 for rates; 0.956 ± 0.018, for covariances), indicating acceptable fits of the models (Appendix—figure 1A-C). This goes in line with the reported high Kullback-Leibler ratios (0.95 ± 0.03 on average) that quantify model performances (Appendix—figure 1D).

We next defined the estimation errors as the averaged normalized squared difference between the values estimated from the data and those estimated by the model. We found a correlation between the estimation and sensitivity of 0.192 ± 0.055 (Appendix—figure 1E). We next asked whether the correlations between sensitivity and cortical state observed in Figure 5F could be fully predicted by correlations between errors and cortical state. We found that correlations between sensitivity and were higher than correlations between estimation errors and (0.252 ± 0.056 vs. 0.168 ± 0.017, Appendix—figure 1F). We used Meng's z-test for dependent correlations to compare the correlations between sensitivity and those between estimation errors and, given the correlation between and errors, and found significant differences in 4/6 datasets (p < 0.05). This was the case for datasets with the highest model performances (datasets 1–4). Finally, we tested whether estimation errors for firing rates correlated with the modulation index (MI) of the neurons and compared the results with those obtained in Figure 6D. We found that correlations between estimation errors and MI were substantially weaker than correlations between estimation errors and MI (-0.506 ± 0.043 vs. -0.086 ± 0.017 on average; p < 0.05 for all datasets, Meng's z-test; Appendix—figure 1G).

We concluded that, although estimation errors correlated with sensitivity, they did not fully explain neither the correlation between sensitivity and cortical state (for the datasets with highest model performance) nor the correlation between sensitivity and MI.

Most importantly:

– Figure 5E-G, check for how correlations between sensitivity and "state" change with errors in model estimation.

As mentioned above, estimation errors did not fully explain the correlation between sensitivity and cortical state.

– Figure 6D-E, check for how correlations between sensitivity and the modulation index change with errors in model estimation.

As mentioned above, estimation errors did not fully explain neither the correlation between sensitivity and MI.

2) It is hard to follow what is 'stiff' and what is 'sloppy', as in most of the cases the frontier between the two categories is obviously not sharp. It is clear that there are two ends of the parameter distribution (one end is more stiff and the other more sloppy), but in between, there are many intermediates. The authors must be more explicit about this and make more precise wherever possible what they mean when they define these two categories.

– In Figure 4B they should display the boundary between stiff and sloppy with a shading or a vertical line.

Stiff and sloppy dimensions are those associated with low and high rank eingenvectors of the FIM, respectively. In sloppy systems there is no clear cut between stiff and sloppy dimensions, so it is not possible to draw a line separating stiff and sloppy dimensions, thus we added a shading in Figure 4B. We concentrated into the first eigenvector to estimate the sensitivity of parameters, but, as shown below, we now also provide a more general definition of sensitivity that considers the weighted contribution to all eigenvectors (see below).

Also, there is a kink in the eigenvalue curve. Is it related to model/population size? Are the eigenvalues beyond this kink relevant?

In our study, the stiff dimension was defined as being the first eigenvector of the FIM (see next point to a discussion about an alternative definition). The kink in the eigenvalue curve observed in Figure 4B was observed in two datasets, for rat 1 (which is presented in the figure) and rat 2. For the other datasets, this kink was not clearly observed as shown Author response image 1. Thus, this observation is not systematic in our sample. The relation between the kink observed in rats 1 and 2 and the size of the population is not straightforward. One could think that the first eigenvectors relate to the bias parameters of the neurons, but, as shown in Figure 5D (bottom), biases and couplings contributed equally to the first eigenvector. We believe that a discussion about the kink observed in those datasets and its relation to model size is out of the scope, since our results are not based on this part of the spectrum but mainly focus on the first eigenvector. Concerning the relevance of eigenvectors beyond the first one, see next point.

Author response image 1

– It should be explained why stiffness of individual neurons is chosen to be their contribution to the first eigenvalue only. Are further eigenvalues not relevant?

Following the work of Panas et al., 2015, we have measured the sensitivity of a given model parameter by its contribution to the first eigenvector of the FIM averaged across time. However, it could happen that the parameters that contribute less to the first eigenvector do contribute to the other stiff dimensions (those with lower rank, for example 2). To take into account this, in the new manuscript, we added an alternative definition of sensitivity that considers the weighted contribution to all eigenvectors of the FIM. Specifically, for each 100-s epoch and each neuronal ensemble, we calculated the eigenvectors of the FIM (vt,1,…,vt,55) and their associated eigenvalues (at,1,…,at,55). We defined the weighted sensitivity of the parameter as the temporal average of its weighted contribution to all eigenvectors of the FIM, with weights equal to the associated eigenvalues:

sne,iw=1NE∑t=1NE∑k=155at,k|Vt,ki|at,1+at,2+at,55

Just as for the original definition, one can construct a weighted sensitivity at the population level.

Using this more general definition, all the results were preserved. This version of the sensitivity is described in Materials and methods (subsection “Sensitivity measures”) and the results are described in Figure 5—figure supplement 1, Figure 6—figure supplement 1, Figure 7—figure supplement 1, and in the main article in subsections “Cortical state transitions evolve along stiff dimensions”, “Sensory-evoked activity evolves along sloppy dimensions” and “Stiff parameters were associated to central neurons within the neuronal network”.

– Together it should be clarified in a dedicated paragraph that stiff and sloppy parameters/neurons are two parts of a wide distribution.

We explicitly mention this point in the discussion of the revised manuscript. See paragraph one and five in the Discussion.

https://doi.org/10.7554/eLife.53268.sa2

Article and author information

Author details

  1. Adrian Ponce-Alvarez

    Center for Brain and Cognition, Computational Neuroscience Group, Department of Information and Communication Technologies, Universitat Pompeu Fabra, Barcelona, Spain
    Contribution
    Conceptualization, Formal analysis, Investigation, Methodology
    For correspondence
    adrian.ponce@upf.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1446-7392
  2. Gabriela Mochol

    Center for Brain and Cognition, Computational Neuroscience Group, Department of Information and Communication Technologies, Universitat Pompeu Fabra, Barcelona, Spain
    Contribution
    Data curation, Formal analysis, Investigation, Methodology
    Competing interests
    No competing interests declared
  3. Ainhoa Hermoso-Mendizabal

    Institut d’Investigacions Biomèdiques August Pi i Sunyer (IDIBAPS), Barcelona, Spain
    Contribution
    Data curation, Formal analysis, Methodology
    Competing interests
    No competing interests declared
  4. Jaime de la Rocha

    Institut d’Investigacions Biomèdiques August Pi i Sunyer (IDIBAPS), Barcelona, Spain
    Contribution
    Data curation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3314-9384
  5. Gustavo Deco

    1. Center for Brain and Cognition, Computational Neuroscience Group, Department of Information and Communication Technologies, Universitat Pompeu Fabra, Barcelona, Spain
    2. Institució Catalana de la Recerca i Estudis Avançats (ICREA), Barcelona, Spain
    3. Department of Neuropsychology, Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany
    4. School of Psychological Sciences, Monash University, Melbourne, Australia
    Contribution
    Funding acquisition, Methodology
    Competing interests
    No competing interests declared

Funding

European Commission (Flag-Era JTC PCI2018-092891)

  • Adrian Ponce-Alvarez
  • Gustavo Deco

Ministerio de Economía y Competitividad (Juan de la Cierva Fellowship IJCI-2014-21937)

  • Gabriela Mochol

European Regional Development Fund (SAF2013-46717-R)

  • Jaime de la Rocha

Ministerio de Economía y Competitividad (BES-2011-049131)

  • Ainhoa Hermoso-Mendizabal

Ministerio de Economía y Competitividad (SAF2010-15730)

  • Jaime de la Rocha

European Regional Development Fund (SAF2010-15730)

  • Jaime de la Rocha

Ministerio de Economía y Competitividad (SAF2013-46717-R)

  • Jaime de la Rocha

Horizon 2020 Framework Programme (785907 HBP SGA1)

  • Gustavo Deco

Horizon 2020 Framework Programme (785907 HBP SGA2)

  • Gustavo Deco

Ministerio de Economía y Competitividad (PSI2016-75688-P)

  • Gustavo Deco

Catalan Research Group Support (2017 SGR 1545)

  • Gustavo Deco

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

APA and GD received funding from the FLAG-ERA JTC (PCI2018-092891). GD acknowledges funding from the European Union’s Horizon 2020 FET Flagship Human Brain Project under Grant Agreement 785907 HBP SGA2, the Spanish Ministry Project PSI2016-75688-P (AEI/FEDER) and the Catalan Research Group Support 2017 SGR 1545. GM was supported by a Juan de la Cierva fellowship (IJCI-2014–21937) from the Spanish Ministry of Economy and Competitiveness. AHM received support from the Spanish Ministry of Economy and Competitiveness (BES-2011–049131). JR received funding from the Spanish Ministry of Economy and Competitiveness together with the European Regional Development Fund Grants SAF2010-15730 and SAF2013-46717-R. We thank L Hollender for sharing her data.

Ethics

Animal experimentation: All experiments were carried out in accordance with protocols approved by the Animal Ethics Committee of the University of Barcelona (Comité d'Experimentació Animal, Universitat de Barcelona, Reference: 116/13).

Senior Editor

  1. Michael J Frank, Brown University, United States

Reviewing Editor

  1. Brice Bathellier, CNRS, France

Reviewers

  1. Brice Bathellier, CNRS, France
  2. Matthias H Hennig
  3. Mark D Humphries, University of Nottingham, United Kingdom

Publication history

  1. Received: November 1, 2019
  2. Accepted: March 16, 2020
  3. Accepted Manuscript published: March 17, 2020 (version 1)
  4. Version of Record published: March 31, 2020 (version 2)

Copyright

© 2020, Ponce-Alvarez 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

  • 713
    Page views
  • 129
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)