Scale-free avalanches in quiet wakefulnessA necessary feature of critical dynamics is scale-free statistics15 in the form of$$p(S)\propto {S}^{-\tau },$$
(1)
where S is the spatiotemporal size of clustered activity (in this case, calcium transients), also known as an avalanche, and τ is the critical exponent. To define calcium-transience avalanches across the field of view, we first binarized the processed Thy1-jRGECO1a calcium fluorescence signal (here on out referred to as “the calcium signal”) using a per-pixel threshold of ϕ standard deviations above the mean (see also methods for further details). We utilize jRGECO as this sensor has a favorable combination of brightness (high SNR) and decay time28, and (as with other red-shifted calcium biosensors) reduced hemodynamics contamination. An example of the mean binarized signal (ϕ = 1) across the duration of a single recording is shown in Fig. 1a. Above threshold, or “active”, pixels are clustered into avalanches.Fig. 1: Avalanche analysis of quiet wakefulness controls (QW).a The average number of activations for a single QW recording is illustrated (ϕ = 1). The midline overlying the sagittal sinus is removed to prevent cross-hemisphere merging. b An example of activity clustering is illustrated. Any active site within both one frame and a neighborhood of radius r of the active reference pixel (Ref.), such as the red pixel but not the yellow one, is clustered to the same avalanche as Ref. c The mean partial correlation function averaged over all awake recordings (N = 12) is shown, with shaded area denoting standard deviation. The zero crossing (r*) defines the neighborhood used to classify avalanches. d Avalanche size distributions for various thresholds ϕ are shown (in units of standard deviations), with the space shuffled surrogate distribution and the dot-dashed line representing τ. Source data are provided as a Source data file.Cortical dynamics have long-range functional connectivity extending beyond a single pixel. As such, clustering schemes based only on the nearest neighboring pixels may break up causally related avalanches. To account for this we define a radius of direct influence, r. Two active pixels within a distance r, and one time frame of one another, are clustered within the same avalanche (Fig. 1b)29. Pre-existing avalanches crossing this boundary are “merged”, and both are (retroactively) relabeled to be the same avalanche. The radius r is estimated from the first zero crossing of the partial correlation function (PCF, see also methods), occurring at r* = 8 ± 2 (Fig. 1c). Once clustered, the normalized avalanche size (S) is a unit-less value defined as the total number of active pixels participating in the avalanche (repeat activation included), divided by the number of cortical pixels per hemisphere (\({{{\mathcal{O}}}}(1{0}^{4})\)).Probability densities were estimated via logarithmic binning of avalanche sizes for all quiet wakefulness (QW) recordings (N = 12). For ϕ = 1, the avalanche rate across recordings was 13 ± 6 avalanches/second, and a total of 7(2) × 102 avalanches/recording. QW avalanche size distributions followed scale-free statistics with an estimated mean exponent, averaged across recordings, τ = 1.2(1), and a typical dynamic range of \(\Delta={S}_{\max }/{S}_{\min }\approx {{{\mathcal{O}}}}(1{0}^{4})\). Figure 1d presents the avalanche size distribution across multiple thresholds obtained by concatenating across all recordings, which is representative of individual recordings (see S1).We also analyzed surrogate distributions generated by random shuffling of pixels in space, and random cyclic permutations of pixel signals in time17,30. The former, shown as the dashed curve in Fig. 1d, preserves the tail (due to finite size effects) but changes the head of the distribution. The latter consistently generated never-ending, system-sized avalanches, with smaller events being almost entirely non-existent (not shown for figure clarity). This analysis suggests both spatial-temporal correlations are important aspects of the avalanche distribution p(S).In addition to avalanche sizes, we analyzed avalanche durations, T, defined as the number of seconds (or frames) that an avalanche took to occur. Their range of values is much smaller than that of avalanche sizes, covering only about two orders of magnitude. Nevertheless, their distributions, p(T), also followed power-law statistics typically over one order of magnitude with an exponent α = 1.7(2) (more on durations in S2). If criticality holds, α and τ uniquely determine the relation between the average avalanche duration for a given size, 〈T 〉(S), and S16. Specifically, 〈T 〉(S) ∝ Sγ with$$\gamma=\frac{\tau -1}{\alpha -1}\equiv \gamma (\alpha,\tau ).$$
(2)
γ can be directly estimated via robust linear fit to 〈T 〉(S) (see S2), and compared to γ(α, τ), which is calculated from the above equation. The independent equivalency of γ and γ(α, τ) can be considered a test for criticality31. We found that in ten out of twelve recordings, the equivalence between γ and γ(α, τ) held within statistical uncertainties of one σ. Of the two recordings where the scaling relation only held within three σ, we found that it was because the size exponent τ was very close to 1, where the scaling relation given by Eq. (2) becomes problematic and systematic uncertainties are large. Averaged over the twelve recordings we found γ = 0.31(4), and γ(α, τ) = 0.34(9). The exponents mentioned above, as well as those obtained in the following sections, are all summarized in Table 1.Table 1 Estimated critical exponents averaged across recordings belonging to the different conditions, with the number of recordings in bracketsAvalanche statistics in the presence of anestheticsLow-dose anesthesiaNext, we analyzed low-dose anesthesia recordings corresponding to isoflurane 1% (N=6), ketamine 10 mg/kg (N=8), and pentobarbital 12.5 mg/kg (N=6). Examples of the cortex-wide average calcium trace in a window of 20 s is shown in Fig. 2a for each anesthetic and QW. Figure 2b shows no significant changes to the PCF in any case. At ϕ = 1, avalanche size distributions follow a power-law distribution (Fig. 2c). The average exponents for isoflurane and ketamine were estimated to be τ = 1.24(4) and τ = 1.08(7), respectively. Pentobarbital recordings had a lower average exponent (τ = 1.0(1)), and in three recordings τ < 1. For all three anesthetics, the fitting domain was reduced to \(\Delta \approx {{{\mathcal{O}}}}(1{0}^{3})\).Fig. 2: Avalanche analysis of low-dose anesthetic recordings.a Examples of (z-scored) traces of the total calcium activity. b The mean estimated PCF for isoflurane 1% (N = 6), ketamine 10 mg/kg (N = 8), and pentobarbital 12.5 mg/kg (N = 6) with quiet wakefulness (QW) for reference. The vertical line shows the first zero crossing for QW. The shaded area for QW denotes uncertainty standard deviation (omitted from others for clarity). c Avalanche size statistics along with QW for reference. The dot-dashed reference line is the QW reference. d KS distance between individual recordings, against the concatenated QW distribution. Axis limits expanded to allow for direct comparison with surgical plane anesthesia in Fig. 3. Source data are provided as a Source data file.Per recording, avalanche durations had an average exponent of α = 1.5(2) for isoflurane, α = 1.9(3) for ketamine, and α = 1.7(2) for pentobarbital. The scaling relation given by Eq. (2) held for all testable (i.e., τ > 1) isoflurane and pentobarbital recordings within two σ, while it held in five out of seven ketamine recordings. For the other two recordings, τ was very close to 1 where systematic uncertainties can become very large as mentioned above. The aforementioned exponents are shown in Table 1.We also calculated the Kolmogorov–Smirnov distance (DKS) of avalanches on a per-recording basis against the concatenated QW statistics as a binning-independent measure of (dis)similarity (Fig. 2d) (see methods). Per-recording QW recordings were tested for consistency and had a mean value of 0.07(4). Isoflurane exhibited similar behavior (〈DKS〉 = 0.09(3)), whereas ketamine recordings showed larger deviations (〈DKS〉 = 0.15(5)), which is consistent with what is observed in Fig. 2c. Pentobarbital recordings had a KS distance of 〈DKS〉 = 0.2(1), but also showed higher variability.To test statistical significance, we performed a two-sample Wilcoxon rank-sum test (rank-sum, MATLAB) with Bonferroni–Holm correction between both avalanche size exponents and DKS values of QW and all low-dose cases. Specifically, for the critical exponents τ, we found that all low-dose cases were not statistically different from QW at the 95% significance level. The rank-sum test between the distributions of DKS values suggested that only the low-dose pentobarbital is statistically different (p = 0.002).Surgical plane anesthesia: predominant dynamicsExamples of the total calcium activity under surgical plane anesthesia levels are shown in Fig. 3a. Isoflurane and ketamine each produced two distinct dynamical modes. Isoflurane 2% predominantly produced smallscale fluctuations interrupted by aperiodic cortex-wide bursts (7 recordings out of N = 10 total), whereas ketamine 100 mg/kg was typically (6 out of N = 7) associated with wave-like activity. In these cases, which we call the “predominant” case, the changes in calcium dynamics can be clearly differentiated from QW by way of phase-space representations obtained by Hilbert transform (S3). The other recordings, which could not be differentiated from QW, are discussed in the next section.Fig. 3: Avalanche analysis of predominant dynamics for surgical plane anesthesia.a Examples of (z-scored) traces of the total calcium activity. b The mean estimated PCF for the isoflurane 2% (N = 7), ketamine 100 mg/kg (N = 6) and pentobarbital 80 mg/kg (N = 4) with quiet wakefulness (QW) for reference. The vertical line shows the first zero crossing for QW. The shaded area for QW denotes uncertainty standard deviation (omitted from others for clarity). c Avalanche statistics along with QW for reference. The dot-dashed reference line is the QW reference. d KS distance between individual recordings, against the concatenated QW distribution. Awake-like recordings (indicated by diamonds) are included for comparison. Source data are provided as a Source data file.As Fig. 3a shows, Pentobarbital 80 mg/kg (N = 4) produced more qualitatively regular dynamics, though not to the same extent as ketamine. Figure 3b shows the mean PCF, averaged across recordings. We found that isoflurane heavily reduced the PCF while ketamine and pentobarbital appeared to produce a minor suppression of the PCF, but more statistics would be needed to make a definitive claim. Figure 3c shows the avalanche size distributions. The wave-like dynamics under ketamine produce a characteristic size which manifest as a large peak near S ≈1, indicating excessive cortex-wide avalanches consistent with wave-like dynamics. Conversely, isoflurane produced a reduction in large avalanche sizes. We also tested increasing the neighborhood radius to r = 14 for the isoflurane analysis due to the elongated PCF in 3b, but found no significant differences. Interestingly, the higher dose of pentobarbital (80 mg/kg) produced a similar result to the lower dose (12.5 mg/kg)—a shift in the exponent to τ = 1.0(1), and a duration exponent of α = 1.8(5). A two-sample Wilcoxon rank-sum test with Bonferroni–Holm correction between QW and pentobarbital for τ indicated that the estimated critical exponents were statistically different (p = 0.004). The scaling relation (Eq. (2)) held within two σ in two of the three testable recordings. Finally, Fig. 3d shows the per-recording KS distance relative to QW. Relative to the low-dose anesthetics, surgical plane levels of ketamine and isoflurane both showed substantially increased deviations away from QW (〈DKS〉 = 0.4(1) for isoflurane and 〈DKS〉 = 0.6(2) for ketamine). Conversely, pentobarbital was consistent with sub-anesthetic results in Fig. 2d (〈DKS〉 = 0.2(1)).Surgical plane anesthesia: awake-like dynamicsInterestingly, we also observed dynamics under ketamine (1 out of 7) and isoflurane (3 out of 10) resembling QW (see S3 and Fig. S6a), despite the surgical plane of anesthesia (S4). We will refer to these as awake-like (AL). AL recordings exhibited statistics different from other in-group recordings. For example, the avalanche rate of the AL ketamine recording was 15 Hz (consistent with QW), as opposed to 4 ± 1 Hz for the other ketamine recordings. In both AL ketamine and isoflurane, both the PCF and avalanche statistics resembled QW (Fig. S6b, c). The estimated exponents were also consistent with QW; τ = 1.17(4) for isoflurane, and τ = 1.05(5) for ketamine (uncertainty here is from MLE as there is only one recording), and had comparable dynamical ranges (\(\Delta \approx {{{\mathcal{O}}}}(1{0}^{3})\)). AL recordings also had the lowest KS distance of surgical plane anesthesia recordings (Fig. 3d), emphasizing the similarity with QW further. Moreover, in all awake-like isoflurane recordings, the scaling relation (Eq. (2)) held within two σ, while it was not in the single ketamine recording because τ was very close to unity. Note that no meaningful statistical comparison between QW and AL could be performed due to the low sample size of AL.Spatial analysisHere we study avalanche initiation/nucleation sites across the cortical surface (see spatial analysis in methods). In all cases, avalanches are preferentially initiated in key locations, as seen in Fig. 4, but are sensitive to pharmacological manipulation. Three regions generated the most avalanches in QW—the primary Somatosensory (SSp) area, particularly areas associated with both the mouth and nose, the retrosplenial (RSP) area, and the secondary somatomotor (MOs) area. The Gini coefficient, G, is also calculated as a measure of focalization32:$$G={(2{n}^{2}\bar{x})}^{-1}{\sum}_{ij}^{n}| {x}_{i}-{x}_{j}|,$$
(3)
where xi is the value of the ith pixel, and \(\bar{x}\) is the average across all pixels. G = 0 indicates a perfectly uniform distribution of activity, while G = 1 indicates entirely localized. For QW, we found G = 0.437. We organize our observations on a per-drug basis.Fig. 4: The average avalanche initiation maps for each case.Color bars indicate the number of times a pixel was classified as an initiation site, averaged over all in-group maps, weighted by the number of observed avalanches with the Gini coefficient G indicated. Awake-like recordings were kept separate as per the previous analysis. The somatomotor (MOs), somatosensory (SSp), and retrosplenial (RSP) areas are also labeled as regions of interest for reference in the first panel. The panel for quiet wakefulness has the Allen Institute atlas outlined for reference.PentobarbitalSpatial maps of for both 12.5 mg/kg and 80 mg/kg pentobarbital resembled QW. However, the heavier doses also showed recruitment of additional sub-regions, which is reflected in G initially decreasing (relative to QW) under low doses, but then increasing in high doses. Along with the mouth and nose sub-regions, we also observe increased activation in upper, lower, and truck areas of the SSp, as well as greater coverage of the MOs.IsofluraneIsoflurane 1% was associated with SSp avalanche initiation similar to QW, but reduced initiation in RSP and MOs. At a surgical plane of anesthesia (2% isoflurane), we observed initiations confined to the hind-limb region of SSp, and a reduction of initiations in facial SSp regions. Spatial maps of AL recordings resembled those of isoflurane 1% (see S3A). Interestingly, despite low-dose baring the most resemblance to QW in terms of avalanches, all isoflurane-associated maps are overall more homogeneous, which is reflected by the Gini coefficient.KetamineKetamine 10 mg/kg recordings displayed activations in the SSp. A suppression of activation in the RSC and MOs was also observed, though not as drastic as in the isoflurane case. Ketamine 100 mg/kg recordings exhibited wave-like dynamics, predominately initiating in the MOs, the result of which leads to the highest observed G among all cases. The single AL recording did not show significant activation in the MOs and instead resembled the 10 mg/kg ketamine map (see S3A).
Existence of multiple transitions of the critical state due to anesthetics
