Statistical analyses support power law distributions found in neuronal avalanches

Abstract

The size distribution of neuronal avalanches in cortical networks has been reported to follow a power law distribution with exponent close to -1.5, which is a reflection of long-range spatial correlations in spontaneous neuronal activity. However, identifying power law scaling in empirical data can be difficult and sometimes controversial. In the present study, we tested the power law hypothesis for neuronal avalanches by using more stringent statistical analyses. In particular, we performed the following steps: (i) analysis of finite-size scaling to identify scale-free dynamics in neuronal avalanches, (ii) model parameter estimation to determine the specific exponent of the power law, and (iii) comparison of the power law to alternative model distributions. Consistent with critical state dynamics, avalanche size distributions exhibited robust scaling behavior in which the maximum avalanche size was limited only by the spatial extent of sampling ("finite size" effect). This scale-free dynamics suggests the power law as a model for the distribution of avalanche sizes. Using both the Kolmogorov-Smirnov statistic and a maximum likelihood approach, we found the slope to be close to -1.5, which is in line with previous reports. Finally, the power law model for neuronal avalanches was compared to the exponential and to various heavy-tail distributions based on the Kolmogorov-Smirnov distance and by using a log-likelihood ratio test. Both the power law distribution without and with exponential cut-off provided significantly better fits to the cluster size distributions in neuronal avalanches than the exponential, the lognormal and the gamma distribution. In summary, our findings strongly support the power law scaling in neuronal avalanches, providing further evidence for critical state dynamics in superficial layers of cortex.

Figure 1. Avalanche size distributions analyzed in the present study.

A. Average in vitro cluster size distribution in organotypic cortex slice cultures (60 electrodes, 7 cultures, n = 53,443 avalanches on average) in double-logarithmic (upper panel) and linear scale (lower panel). The data set was taken from . Inset: view of a culture on a 88 electrode array (scale bar, 1 mm). B. Average in vivo cluster size distribution from rat somatosensory cortex under urethane anesthesia (27–31 electrodes, 7 recordings, n = 22,321 avalanches on average). Data was taken from . Inset: view of the insertion sites for an 84 array (triangles) in cortical layer 2/3 (vertical scale bar, 1 mm). C. 43-min recording for monkey X (low-density microelectrode array with 32 electrodes in the left primary motor cortex, n = 45,574 avalanches). Data was taken from . D. 30-min recording for the second monkey (monkey Y, high-density microelectrode array with 91 electrodes in the left premotor cortex, n = 24,877 avalanches). Insets in C and D show the location of the multielectrode arrays (scale bar, 10 mm). The size of the arrays (dark squares) is not shown in the actual scale. The number of electrodes in the individual arrays is indicated by arrows in the log-log plots (A–D).

Figure 2. Collapse of rescaled cluster size distributions in neuronal avalanches.

A. Depiction of the rescaling approach for synthetic PMFs for maximum sizes N = 8, 16, 32, 64 (left). The system size, N, corresponds to the number of electrodes included in the analysis. Cluster sizes s were normalized by the system size N (s s/N) and the renormalized probability was obtained according to P(s) P(s)/A(N), resulting in a collapse of the cluster size distributions (right). Here, the definition of A(N) with upper bound N was used (Eq. 16). The vertical arrow indicates the system size (scaled to unity). B. Collapse of rescaled cluster size distributions for average in vitro distributions (n = 7), average in vivo distributions under anesthesia (rat, n = 7), and the two awake monkeys with low- and high-density array, respectively (from left to right). Note that the maximum cluster size for all data sets increases with N with the distribution showing a clear cut-off beyond the system size (s/N = 1). The exponent for the empirical distributions was fitted individually for each system size N (see Materials and Methods ).

Figure 3. Estimation of slope parameter for the in vitro data (n = 7 cultures), the in vivo data under anesthesia (n = 7), and the in vivo recordings in awake monkeys (n = 2).

A. Shown are the average slope parameters and the standard deviations (error bars). The three different estimation methods are: LS least-square estimation with logarithmic binning, KS Kolmogorov-Smirnov statistic, and ML maximum likelihood estimation (see Materials and Methods ). Estimated values of were not statistically different: in vitro, F(2,18) = 0.19, p = 0.827; in vivo (anesthesia), F(2,18) = 0.124, p = 0.884; in vivo (awake), F(2,3) = 0.21, p = 0.821 (one-way ANOVA). Values of were estimated for the entire range of cluster sizes, i.e., from avalanches that included only one electrode to clusters that spanned the entire multielectrode array. B. Autocorrelation of the avalanche sizes for monkey X and Y as a function of the avalanche lag. The autocorrelation showed a fast decay within 10 avalanches (arrow). The shaded areas indicate the autocorrelation (3 SD) for randomly permuted cluster sizes for both monkeys. C. Average values obtained by ML estimation as a function of avalanche lag for monkey X and Y (red and green line, respectively). Error bars denote the standard deviation across the decorrelated sub-sets. The gray lines show meanSD of ML parameter estimates for sample-size matched data from the original sequence of cluster sizes.

Figure 4. Model comparison using the LLR test.

A. Model fits obtained by ML estimation for the power law (red) and the exponential model (green) for the cluster size distribution in monkey X. B. LLR values for increasing threshold, z. Error bars denote the SD across the decorrelated sub-sets of the data. All LLRs were positive and statistically different from zero (p<0.0001). The avalanche lags ranged from 10 (z = −1.5) to 2 (z = −5). C and D. The same for the comparison between the exponentially truncated power law (red) and the lognormal distribution as the alternative model (green). The insets show detailed views of the distributions, corresponding to the respective gray rectangles. The square symbols in (D) indicate the LLR values that were statistically different from zero (p<0.01).

Figure 5. Model comparison based on the KS statistic.

A. PMFs of the two-parameter models for the avalanche size distribution in monkey X. The inset shows a detailed view that corresponds to the gray rectangle. B. Corresponding CDF fits for the same size distribution (i.e., monkey X). The insets show detailed views of the distributions, corresponding to the respective gray rectangles. C. Average KS distance of the model distributions for all data sets (n = 16, which includes 7 data sets recorded in vitro, 7 in vivo under anesthesia, and 2 in vivo awake). Error bars denote the standard deviation. The single-parameter power law and the power law with exponential cut-off yielded significantly better fits to the data than the gamma, the lognormal, or the exponential distribution (Kruskal-Wallis test and Tukey-Kramer multiple comparison, p<0.0001).