We first examined the temporal evolution of the spontaneous activity across the $N_E$ 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, $\overrightarrow{\sigma }=\left[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_N\right]$, where ${\sigma }_i=+1$ if the i-th neuron fired a spike in that time bin and ${\sigma }_i=-1$ otherwise. The collective activity was determined by the probability distribution $P\left(\overrightarrow{\sigma }\right)$ over all $2^N$ possible binary patterns. The MEM fits $P_{\text{data}}\left(\overrightarrow{\sigma }\right)$ by finding a distribution $P_{\text{MEM}}\left(\overrightarrow{\sigma }\right)$ that maximizes its entropy under the constraint that the activation rates (${\displaystyle ⟨{\sigma }_i⟩}$) and the pairwise correlations (${\displaystyle <{\sigma }_i{\sigma }_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\left(\overrightarrow{\sigma }\right)\propto e^{-E\left(\overrightarrow{\sigma }\right)}$, where $E\left(\overrightarrow{\sigma }\right)$ is the energy of the pattern $\overrightarrow{\sigma }$, given by: $E\left(\overrightarrow{\sigma }\right)=-\sum _{i=1}^N\left[h_i{\sigma }_i+\frac{1}{2}\sum _{j=1}^N{J}_{ij}{\sigma }_i{\sigma }_j\right]$ (Schneidman et al., 2006; Tkačik et al., 2015). The model parameter $h_i$ represents the intrinsic tendency of neuron i towards activation (${\sigma }_i=+1$) or silence (${\sigma }_i=-1$) and the parameter $J_{ij}$ represents the effective interaction between neurons i and j. Once we learned the parameters ${\displaystyle \mathrm{\Omega }=\left\{\mathit{h},\mathit{J}\right\}}$ 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 $N_E$ 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 $N_{pop}$ 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 2¹⁰ = 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\times N_E$ models, each composed of 10 units. Before studying the evolution of the model parameters ${\displaystyle \mathrm{\Omega }\left(t\right)}$ 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 ($D_{JS}$, 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 ${\displaystyle \mathit{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/D}_{JS}$) 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).

Figure 2

Download asset Open asset

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 $\gamma$; see Materials and methods) between the values in epoch $t$ and those in epoch ${\displaystyle t+\mathrm{\Delta }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 ${\displaystyle t+\mathrm{\Delta }t}$. Specifically, we calculated $⟨D_{JS}⟩\left(\mathrm{\Delta }t\right)$, given by the average Jensen-Shannon divergence between the distribution of data binary patterns in epoch $t$, that is $P_{\text{data,}t}$, and the distribution of binary patterns predicted by the MEM constructed using the data in epoch ${\displaystyle t+\mathrm{\Delta }t}$, that is ${\displaystyle P_{\text{MEM,}t+\mathrm{\Delta }t}}$ (see Materials and methods). We found that $⟨D_{JS}⟩\left(\mathrm{\Delta }t\right)$ increased as a function of $∆t$, indicating that the collective activity changed during the recording session, so that the model parameters ${\displaystyle \mathrm{\Omega }\left(t\right)}$ at epoch $t$ did not predict the collective activity at epoch ${\displaystyle t+\mathrm{\Delta }t}$ (Figure 3B). Indeed, the model parameters substantially changed over time, with a rapidly decaying similarity (Figure 3C, orange and purple traces).

Figure 3

Download asset Open asset

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 ${\displaystyle P_{\text{MEM}}\left(\overrightarrow{\sigma }|\mathrm{\Omega }\right)}$ changes with respect to changes in the parameters ${\displaystyle \mathrm{\Omega }}$. 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 $\gamma =$ 0.882 ± 0.002 for the FIM, $\gamma =$ 0.732 ± 0.003 for the firing rates, $\gamma =$ 0.551 ± 0.004 for the biases, $\gamma =$ 0.364 ± 0.004 for the correlations, $\gamma =$ 0.234 ± 0.003 for the couplings (F_4,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.

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

Download asset Open asset

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 ${\displaystyle \mathrm{\Omega }\left(t^{\prime }\right)}$, calculated at time $t\text{'}$, into the eigenvectors of the FIM at time $t$, denoted ${\mathit{\nu }}_{t,1},{\mathit{\nu }}_{t,2},\dots ,{\mathit{\nu }}_{t,k},\dots$, where $k$ is the rank of the eigenvector (Figure 4C). For each dimension, or eigenvector, we obtained a distribution of projections of parameters ${\displaystyle \mathrm{\Omega }\left(t^{\prime }\right)}$ (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.

To test this hypothesis, we first measured the cortical state in each epoch $t$ using silence density, $CS\left(t\right)$, 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\left(t\right)$. 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 ${\displaystyle D_{JS}\left(t,t^{\prime }\right)=D_{JS}\left(P_{\text{data,}t};P_{\text{data,}{t}^{\prime }}\right)}$, significantly co-varied with the changes in cortical state, given by ${\displaystyle d=\left|CS\left(t\right)-CS\left(t^{\prime }\right)\right|}$ (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

Download asset Open asset

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, $\left|r_{cs}\right|$, between the cortical state $CS\left(t\right)$ and the activity observables. We found that $\left|r_{cs}\right|$ 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 $s_{\text{ne}}$ 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 (${\displaystyle \mathit{h}}$) and couplings (${\displaystyle \mathit{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 ($\left|r_{cs}\right|$) 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.

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\left(t\right)$. 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

Download asset Open asset

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.

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 $h_i$, and pairs or links were associated to parameters $J_{ij}$). 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 $J_{ij}$ (AUC = 0.541–0.766, p < 0.001, for all rats; Figure 7C,G).

Figure 7 with 1 supplement see all

Download asset Open asset

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).

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).

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.

Long continuous recordings (mean, ∼2 h) were divided into $N_E$ 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.

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 (${\sigma }_i=+1$) or did not (${\sigma }_i=-1$) generate one or more spikes. The state of the neural ensemble is described by a binary pattern $\overrightarrow{\sigma }=\left[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_N\right]$, and thus the collective activity is described by the probability distribution $P\left(\overrightarrow{\sigma }\right)$ over all 2^N possible binary patterns. We estimated $P\left(\overrightarrow{\sigma }\right)$ using a Maximum entropy model (MEM). The MEM finds $P\left(\overrightarrow{\sigma }\right)$ 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 (${\displaystyle <{\sigma }_i>}$) and the pairwise correlations (${\displaystyle <{\sigma }_i{\sigma }_j>}$) are preserved. The maximum entropy distribution $P\left(\overrightarrow{\sigma }\right)$ that is consistent with these expectation values is given by the Boltzmann distribution (Schneidman et al., 2006; Tkačik et al., 2015):

where $E\left(\overrightarrow{\sigma }\right)$ is the energy of the pattern $\overrightarrow{\sigma }$, given by:

and ${\displaystyle Z=\sum _{\{\overrightarrow{\sigma }\}}{e}^{-E(\overrightarrow{\sigma })}}$ is the partition function.

The model parameter $h_i$, called intrinsic bias, represents the intrinsic tendency of neuron i towards activation (${\sigma }_i=+1$) or silence (${\sigma }_i=-1$) and the parameter $J_{ij}$ represents the effective interaction between neurons i and j. The estimation of the model parameters ${\displaystyle \mathrm{\Omega }=\left\{\mathit{h},\mathit{J}\right\}}$ 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 $N_{\text{pop}}$ 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\times N_E$ models of 10 units. We were interested on the evolution of the model parameters over time, that is ${\displaystyle \mathrm{\Omega }\left(t\right)}$. 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 ${\displaystyle \mathrm{\Omega }=\left[h_1,h_2,\dots ,h_N,J_{12},J_{13},\dots \right]}$.

The MEM parameters ${\displaystyle \mathrm{\Omega }=\left\{\mathit{h},\mathit{J}\right\}}$ were iteratively adjusted to minimize the absolute difference between the empirical activation rates ($⟨{\sigma }_i⟩$) and correlations ($⟨{\sigma }_i{\sigma }_j⟩$) and those (${⟨{\sigma }_i⟩}_{\text{model}}$, ${⟨{\sigma }_i{\sigma }_j⟩}_{\text{model}}$) predicted by the model through Metropolis Monte Carlo simulations (100,000 samples). Specifically, each iteration is given by: $h_i^{new}=h_i^{old}-\alpha \left({⟨{\sigma }_i⟩}_{\text{model}}-⟨{\sigma }_i⟩\right)$, and $J_{ij}^{new}=J_{ij}^{old}-\alpha \left({⟨{\sigma }_i{\sigma }_j⟩}_{\text{model}}-⟨{\sigma }_i{\sigma }_j⟩\right)$, where α is the learning rate ($\alpha =$ 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).

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

Where $P_{\text{MEM}}$ was given by the Boltzmann distribution of the model, $P_{\text{data}}$ was estimated from the $N$-dimensional binary patterns observed in the data, and:

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 (${\displaystyle <{\sigma }_i>}$) are preserved (i.e., only $\mathit{h}$ is optimized; first-order model). In this case, the pattern energy is given by: $E\left(\overrightarrow{\sigma }\right)=-\sum _{i=1}^N{h}_i{\sigma }_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 $P_1$ of the first-order model (i.e., independent-MEM) and the distribution of the actual data, $D_1=D_{KL}\left(P_1;P_{\text{data}}\right)$, with the Kullback-Leibler divergence between the distribution $P_2$ of the second-order model and the distribution of the actual data, $D_2=D_{KL}\left(P_2;P_{\text{data}}\right)$. Specifically, the Kullback-Leibler ratio is defined as:

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

Because in the MEM, all the information about the collective activity is contained in the probability distribution of the binary patterns, $P\left(\overrightarrow{\sigma }\right)$, one can define the model parameter space as ${\displaystyle P\left(\overrightarrow{\sigma }|\mathrm{\Omega }\right)}$. 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 ${\displaystyle \mathrm{\Omega }}$ and ${\displaystyle \mathrm{\Omega }+\delta \mathrm{\Omega }}$, 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:

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

The FIM represents the curvature of the log-likelihood of the model, ${\displaystyle \log P\left(\overrightarrow{\sigma }|\mathrm{\Omega }\right)}$, 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:

with 1 $\le k,l\le$ 55 and ${\displaystyle \overrightarrow{x}=\left[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_N,{\sigma }_1{\sigma }_2,{\sigma }_1{\sigma }_3,\dots \right]}$.

The FIM was calculated for every neuronal ensemble at every 100-s epoch and it was decomposed into eigenvectors, noted ${\mathit{\nu }}_{t,1},{\mathit{\nu }}_{t,2},\dots ,{\mathit{\nu }}_{t,k},\dots$, 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 $s_{\text{ne,}i}=\frac{1}{N_E}{\sum }_t\left|{\nu }_{t,1}\left(i\right)\right|$, with $1\le i\le N$.

We next constructed a sensitivity measure for the entire population of $N_{\text{pop}}$ neurons. For this, we defined the set of all single neuron indices and all pairs of neurons $\mathit{I}=\left\{1,\dots ,N_{\text{pop}},\left(\mathrm{1,2}\right),\left(\mathrm{1,3}\right),\dots \right\}$. This set has $L=N_{\text{pop}}+\frac{N_{\text{pop}}(N_{\text{pop}}-1)}{2}$ elements. For each element $j$ of $\mathit{I}$, we defined the sensitivity $s_j$ as the average of $s_{\text{en,}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, $s_{\text{ne}}$ denotes the sensitivity within an ensemble of $N=10$ neurons and has 55 elements, and $s$ denotes the sensitivity within the entire population of $N_{\text{pop}}$ neurons and has $L$ elements. This allows comparison of $s$ with statistics derived from the population of $N_{\text{pop}}$ 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 ($a_{t,1},\dots ,a_{t,55}$):

As previous, from $s_{\text{ne}}^w$ one can construct a weighted sensitivity $s^w$ 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.

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 ${\displaystyle t+\mathrm{\Delta }t}$. For example, let $\overrightarrow{r}\left(t\right)$ the average firing rates of the neurons during the epoch $t$, the similarity measure is given by:

Where $N_E$ is the number of epochs and $\rho$ is the Pearson correlation coefficient. In the case of FIM, the matrix was vectorized to calculate $\rho$.

To evaluated how well the data in a given epoch $t$ could be explained by the MEM constructed using the data at time ${\displaystyle t+\mathrm{\Delta }t}$. Specifically, we defined the similarity measure $⟨D_{JS}⟩\left(\mathrm{\Delta }t\right)$, given by the average Jensen-Shannon divergence between the distribution of data binary patterns in epoch $t$, that is $P_{\text{data,}t}$, and the distribution of binary patterns predicted by the MEM constructed using the data in epoch ${\displaystyle t+\mathrm{\Delta }t}$, that is ${\displaystyle P_{\text{MEM,}t+\mathrm{\Delta }t}}$. This measure is given as:

In other words, $1/⟨D_{JS}⟩\left(\mathrm{\Delta }t\right)$ quantifies how well, on average, the model with parameters ${\displaystyle \mathrm{\Omega }\left(t+\mathrm{\Delta }t\right)}$ represents the data from epoch $t$.

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

where $r_{spon}$ is the pre-stimulus average spike count, calculated in the 0.5-s pre-stimulus interval, and $r_{stim}$ 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.

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 $N_E$ 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:

where $p\left(kl\right)$ 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$.

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 (${\sigma }_i$) and the number of coactive neurons (i.e., with ${\sigma }_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.

We used the receiver operating characteristic curve (ROC) to evaluate the separation between the distributions of observables from sloppy and stiff classes. Let $X_{\text{sloppy}}$ and $X_{\text{stiff}}$ 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\left(c\right)$, is build by plotting the probability of ${\displaystyle P\left(X_{\text{sloppy}}>c\right)}$ against the probability of ${\displaystyle P\left(X_{\text{stiff}}>c\right)}$, for each all $c$. The area under the ROC curve (AUC) is a measure of separation between ${P(X}_{\text{sloppy}})$ and ${P(X}_{\text{stiff}})$, and it is given by:

AUC ranges between 0 and 1, with AUC = 0 if ${P(X}_{\text{sloppy}})$ and ${P(X}_{\text{stiff}})$ are completely separated and ${\displaystyle X_{\text{sloppy}}>X_{\text{stiff}}}$, AUC = 1 if ${P(X}_{\text{sloppy}})$ and ${P(X}_{\text{stiff}})$ are completely separated and ${\displaystyle X_{\text{stiff}}>X_{\text{sloppy}}}$, and AUC = 0.5 if ${P(X}_{\text{sloppy}})$ and ${P(X}_{\text{stiff}})$ 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.

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 ${\displaystyle \mathrm{\Omega }\left(t\right)}$, through Monte Carlo simulations of the model to obtain 5000 binary patterns. Third, we repeated the Monte Carlo simulations $N_E$ times. Finally, for each of the $N_E$ pieces of surrogate data, we estimated new MEM parameters, ${\displaystyle {\mathrm{\Omega }}^{\prime }}$, 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.

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 $⟨{\sigma }_i⟩$ and $⟨{\sigma }_i{\sigma }_j⟩$ are given by:

The estimated firing rates and the covariances are given as $r_i=\left(⟨{\sigma }_i⟩+1\right)/\left(2dt\right)$ and $C_{ij}=⟨{\sigma }_i{\sigma }_j⟩-⟨{\sigma }_i⟩⟨{\sigma }_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

Download asset Open asset

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

where ${\displaystyle \mathrm{\Delta }{r}_{t,i}={\left(r_i^{\text{data,t}}-r_i^{\text{model,t}}\right)}^2}$ and ${\displaystyle \mathrm{\Delta }{C}_{t,i,j}={\left(C_{ij}^{\text{data,t}}-C_{ij}^{\text{model,t}}\right)}^2}$, and ${\displaystyle {\sigma }_{\mathrm{\Delta }r}}$ and ${\displaystyle {\sigma }_{\mathrm{\Delta }r}}$ denote the standard deviations of ${\displaystyle \mathrm{\Delta }{r}_{t,i}}$ and ${\displaystyle \mathrm{\Delta }{C}_{t,i,j}}$, respectively —this normalization is needed because ${\displaystyle \mathrm{\Delta }{r}_{t,i}}$ and ${\displaystyle \mathrm{\Delta }{C}_{t,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 $N_{\text{pop}}$ 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 $\left|r_{cs}\right|$ were higher than correlations between estimation errors and $\left|r_{cs}\right|$ (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 $\left|r_{cs}\right|$ and those between estimation errors and $\left|r_{cs}\right|$, 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.

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, ${\displaystyle \overrightarrow{r_e}=\left[r_{e,1},r_{e,2},\dots ,r_{e,N_{pop}}\right]}$, 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 ${\lambda }_{e,i}\left(t\right)$ of the point process, that is ${\lambda }_{e,i}\left(t\right)$ was a step function with ${\lambda }_{e,i}\left(t\right)=r_{e,i}$, for times $t$ within active periods, and ${\lambda }_{e,i}\left(t\right)=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

Download asset Open asset

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, ${\displaystyle {⟨C⟩}_e}$, from the original and the surrogate data were highly correlated (Appendix 1—figure 2B,E), although ${\displaystyle {⟨C⟩}_e}$ 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, ${\displaystyle C_{ij}}$, 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 ${\displaystyle C_{ij}}$ nor BC were predicted by the NHPP: the average fraction of explained variance was equal to R² = 0.148 ± 0.07 and R² = 0.184 ± 0.07 for ${\displaystyle C_{ij}}$ 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.

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

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.

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).

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)

© 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.