Many complex systems give rise to events that are clustered in space and time, thereby establishing a correlation structure that is governed by power law statistics. In the cortex, such clusters of activity, called “neuronal avalanches,” were recently found in local field potentials (LFPs) of spontaneous activity in acute cortex slices, slice cultures, the developing cortex of the anesthetized rat, and premotor and motor cortex of awake monkeys. At present, it is unclear whether neuronal avalanches also exist in the spontaneous LFPs and spike activity in vivo in sensory areas of the mature brain. To address this question, we recorded spontaneous LFPs and extracellular spiking activity with multiple 4 × 4 microelectrode arrays (Michigan Probes) in area 17 of adult cats under anesthesia. A cluster of events was defined as a consecutive sequence of time bins Δt (1–32 ms), each containing at least one LFP event or spike anywhere on the array. LFP cluster sizes consistently distributed according to a power law with a slope largely above –1.5. In two thirds of the corresponding experiments, spike clusters also displayed a power law that displayed a slightly steeper slope of −1.8 and was destroyed by subsampling operations. The power law in spike clusters was accompanied with stronger temporal correlations between spiking activities of neurons that spanned longer time periods compared with spike clusters lacking power law statistics. The results suggest that spontaneous activity of the visual cortex under anesthesia has the properties of neuronal avalanches.
Neurons in primary sensory cortices display firing even in the absence of sensory stimulation. Previously, this spontaneous activity was considered to be noise, discharges of single neurons being stochastic and uncorrelated (Shadlen and Newsome 1998). Recently, this view was challenged by several studies using voltage-sensitive dye imaging (Arieli et al. 1995) and intracellular recordings (Bringuier et al. 1999), suggesting that spontaneous neuronal activity is coherent and correlated within a large cortical area. Activity in different brain areas is linked through a cascade of synaptic inputs that propagates in a wave-like fashion from one cortical site to another (see Wu et al. 2008 for a review). Such activity propagation was shown in both anesthetized and awake animals (Ferezou et al. 2007; Petersen et al. 2003; Xu et al. 2007) Notably, these waves exhibit different sizes (Petersen et al. 2003) and organize in diverse spatiotemporal patterns (Tsodyks et al. 1999), which resemble underlying functional maps (Kenet et al. 2003). It was also shown that these spontaneous waves influence the response to sensory inputs (Ferezou et al. 2007; Petersen et al. 2003) and can account for the variability of evoked activity (Arieli et al. 1996), indicating a pivotal role of spontaneous activity in processing sensory information.
Recently, another type of correlated spontaneous neuronal activity was described in cortical tissue. Beggs and Plenz (2003) found that negative deflections in local field potential signals (nLFP) can propagate in neuronal cultures and acute cortical slices in a nonwave (i.e., noncontiguous) fashion. These cascades of nLFPs, dubbed neuronal avalanches, form spatiotemporal clusters of synchronized activity interrupted by periods of silence and, as shown later in vivo (Gireesh and Plenz 2008), can coexist with theta and beta/gamma oscillations. Furthermore, nLFPs can propagate across many millimeters of the cortex and display long-range temporal correlations, thereby establishing complex spatiotemporal patterns (Beggs and Plenz 2004). Neuronal avalanches emerge during the earliest time of the development of superficial layers in cortex (Gireesh and Plenz 2008; Stewart and Plenz 2008), and their emergence requires a balance of excitatory and inhibitory transmission as well as the presence of the neuromodulator dopamine (Beggs and Plenz 2003; Gireesh and Plenz 2008; Stewart and Plenz 2006). Most importantly, the sizes of neuronal avalanches distribute typically according to a power law with slope −1.5 and show long-range temporal correlations (Beggs and Plenz 2003; Plenz and Thiagarajan 2007). The recent discovery of a power law and long-range temporal correlations in the motor cortex of awake and resting monkeys indicates that neuronal avalanches are not restricted to in vitro and anesthetized in vivo preparations (Petermann et al. 2009). These findings bring the occurrence of avalanches in neural tissue close to the theory of self-organized criticality (SOC), which links power law statistics of events sizes with cascading systems (Bak et al. 1988). Recently, in vitro experiments showed that neuronal networks poised at the critical state display a maximal dynamic range of responses, which disappears when the balance of excitation and inhibition is altered (Shew et al. 2009). This provided first experimental evidence for a possible functional role of critical dynamics such as found in SOC in living neuronal networks.
The LFP signals used in the previous avalanche studies lump together activity of many neurons (mainly synaptic) within a large field of integration. Therefore attempts have been made to extend the investigation of neuronal avalanches to spiking events, and previous studies indeed reported a power law organization similar to that obtained from the earlier LFP studies in vitro (Mazzoni et al. 2007; Pasquale et al. 2008). However, another recent attempt to detect avalanches in the spiking activity in cat parietal cortex in vivo during either awake state or slow wave sleep failed when using small electrode arrays (8 electrodes, aligned linearly) and a large interelectrode distance (1 mm between neighboring electrodes; Bédard et al. 2006). In this study, we made highly parallel recordings of spontaneous spiking activity and LFPs in the visual cortex (area 17) of four anesthetized cats by using 16 or 32 electrodes simultaneously (4 × 4 arrays) with small interelectrode distances (200 μm between neighboring electrodes). This allowed us to test for evidence of neuronal avalanches in LFP events and spikes simultaneously. Our results suggest that spontaneous activity of the visual cortex under anesthesia has statistical properties similar to the neuronal avalanches that were initially described by Beggs and Plenz (2003) and thus may impact the processing of sensory information in a way similar to that of waves of synaptic activity.
Four cats were initially anesthetized with ketamine (Ketanest, Parke-Davis, 10 mg/kg, im) and xylazine (Rompun, Bayer, 2 mg/kg), and the anesthesia was maintained with a mixture of 70% N2O-30% O2 and halothane (1.0%). Tracheotomy was applied, and the animal was fixated in a stereotactic frame. The animals were artificially ventilated, and after craniotomy, the skull was connected to a metal rod, and the halothane level was reduced to 0.5–0.6%. After ascertaining stability of anesthesia to prevent vegetative reactions to somatic stimulation, pancuronium bromide (Pancuronium, Organon, 0.15 mg/kg/h) was applied to obtain paralysis. Glucose and electrolytes were provided by a gastric catheter, and the end tidal CO2 and rectal temperature were maintained between 3 and 4% and 37 and 38°C, respectively. The value of 0.5–0.6% halothane was held constant throughout the experiment except for potentially painful procedures (e.g., intramuscular injection of antibiotic). In this case, we increased the level of halothane to 1.2% 10 min before the procedure and returned immediately back to 0.5–0.6%, which was followed by a period of ≥20 min without new recordings. The nictitating membrane was prepared with neosynephrine, the pupils were dilated with atropine, and the eyes were covered by contact lenses with artificial pupils for protection from desiccation. All procedures abided to the German law for the protection of animals and were supervised by a veterinarian.
We recorded spontaneous activity in area 17 with one or two 16-channel silicon-based microelectrode arrays (Neuronexus Technologies) (Fig. 1). Each array consisted of four 3 mm long shanks, with a profile of 100 × 10 μm at its widest point. Each shank had four electrode contacts. The separation between the neighboring contacts was 200 μm in both directions, i.e., along and across the shanks. This symmetric 4 × 4 arrangement allowed for the maximum distance spanned by the centers of the contacts to equal 600 μm along each dimension and 850 μm diagonally. Thus the recording area of one array spanned ∼0.6 × 0.6 mm2. This array spacing is similar to that used by Beggs and Plenz (2003). Each electrode contact covered an area of 1,250 μm2 and had impedance of 0.3–0.5 MΩ at 1,000 Hz. The arrays were inserted into the cortex always in the same hemisphere such that they penetrated the surface approximately perpendicularly. We recorded mostly the activity from superficial layers; however, we did not identify these layers. All neurons recorded by the same array had overlapping receptive fields. In three cats (cats 1, 3, and 4), the analyzed spontaneous activity was taken from the 1 s period before visual stimulation. For each trial, 1 s recordings of spontaneous activity preceded 4–5 s of visual stimulus presentation (sinusoidal gratings) and were followed by 1 s of recording. Intertrial periods ranged between 1 and 2 s, depending on the experiment. In one cat (cat 2), no visual stimulus was presented, and spontaneous activity was analyzed for the whole trial duration (10 s). The total duration of spontaneous activity obtained in each of the recordings ranged between 100 and 720 s. For each cat, we investigated activity in two different datasets (recordings). In three cats (cats 2, 3, and 4), the two datasets were recorded by the same array but at different time points, the interval between recordings varying between 3 min and 67 h. In one cat (cat 1), the two datasets were recorded simultaneously using two different arrays. Overall, we inserted five different arrays and recorded in total eight datasets, two from each cat. During the recordings, the eyes of the cats were open, and in front of the eyes, a blank (black) computer screen was located. The ambient illumination was low as the lights in the room were strongly dimmed.
Extraction of extracellular spikes and nLFP events
Signals were amplified 1,000× and filtered between 500 Hz and 3.5 kHz for extracting multiunit activity (MUA) and between 1 and 100 Hz for extracting LFPs. MUA signals were sampled with a frequency of 32 kHz, which allowed the later application of off-line spike-sorting techniques. LFP signals were sampled with 1 kHz. The signals were transmitted to an A/D converter and recorded by a customized LabView program running on a PC. Action potentials (spikes) were detected with a two-sided threshold discriminator adjusted manually to yield a signal-to-noise ratio of ∼2:1. For each detected action potential, the time of the event (timestamp) was recorded together with the spike waveform over a duration of 1.2 ms.
After applying spike-sorting techniques based on principal component analysis, we extracted between 43 and 97 units per array (67 ± 17), the majority of which were identified as multi-units. We combined both the multi- and single-unit activity in all analyses to reduce the potential subsampling of spike clusters, which can obstruct the presence of a power law (Priesemann et al. 2009). The firing rates were computed for every unit and recording, and the values ranged between 0.01 and 69.8 Hz, and, as expected, showed a skewed distribution with an overall mean of 7 Hz and median of 2.2 Hz, resulting in an SD of 11.2 Hz. The number of spikes per recording varied between 55,037 and 251,810, with a longer recording duration also leading to the detection of more spikes. In the recordings with power law statistics, the average spike count per second was on average considerably higher [590 ± 112 (SD) spikes/s] than in recordings without a consistent power law (389 ± 179 spikes/s). However, in cat 3, array I–a, in which we also found a power law, the spike density was low (131 spikes/s), and in cat 4, array I–b, in which no power law was detected, the spike count was high (586 spikes/s). LFP events were extracted by applying a threshold of 2 SD of the noise calculated for each electrode, whereas the polarity of this threshold was determined from the polarity of the spike-triggered LFP average. In two recordings, the electrodes showed a spike-triggered LFP average with a negative deflection at the time of spike occurrence and hence we used a threshold of −2 SD. For the remaining six recordings, the unit-triggered LFP average at single electrodes was positive, resulting in a threshold of +2 SD. To determine the time of an event, we detected the point at which the LFP signal reached its maximum between the two threshold crossings at which the signal first exceeded and then returned back below the threshold.
Extraction of event clusters
The spatiotemporal organization of spike and LFP events was studied in the context of neuronal avalanches. We first identified spatiotemporal event clusters in the spiking data based on the same principles as those used for the investigation of avalanches in spontaneous LFP events in earlier studies (Beggs and Plenz 2003; Stewart and Plenz 2006). Spikes were first grouped into small time bins of width Δt, which, across most of analyses and datasets, ranged between 0.5 and 32 ms. At Δt = 32 ms, most recordings either showed a flat distribution of spiking event sizes (i.e., the probability of small and large spiking events was equal) or the number of extracted events was insufficient to compute a distribution. The bins containing at least one spike we refer to as active and to those without spikes as blank. A spatio-temporal cluster begins with a blank bin, is followed by a sequence of active bins, and ends with another blank bin. In other words, sequences of active bins were delineated by periods of no activity, as shown in Fig. 1 (zoom-in). The first and last event clusters of each segment (1 and 10 s in duration) were discarded from further analysis, because they were likely to reflect an incomplete cluster. All the analyses were restricted to neuronal spiking activity obtained from single arrays, and recordings were included into analysis only if all 16 channels of the array detected neuronal spiking activity with signal-to-noise ratio ≥2.
By applying this algorithm, we identified ≤132,400 (47,679 ± 39,684) spike clusters per recording (with Δt = 1 ms). This number depended strongly on the bin size (Δt). If Δt was small, the total number of cluster increased by splitting larger cluster into several smaller ones. Similarly, with large values of Δt, small clusters were combined into larger ones, resulting in fewer clusters. Previous work by Beggs and Plenz (2003) showed that a reasonable estimate of the optimal bin width can be obtained from the gross average of interspike interval (ISI)array distribution, avgΔt, which estimates the average time between successive spikes occurring anywhere in the array (see also Stewart and Plenz 2006 for a detailed explanation of the methods). ISIarrays were calculated always with a resolution of 0.1 ms. The ISIarray distribution quickly decayed to negligible values within ∼50 ms (see Fig. 4B), and there was no requirement to impose any additional cut-off as based on maximal extend of correlations (cf. for the acute slice; Stewart and Plenz 2006). Histograms of ISIs calculated for individual units (ISIunit) were computed always with a resolution of 1 ms and plotted in log-log coordinates. Both ISIarray and ISIunit distributions were fitted with both power law and exponential functions (ISIarray range: 2–30 ms; ISIunit range: 2–500 ms) for statistical purposes.
For each cluster, the cluster size was obtained by counting either 1) the number of spikes across all units and active electrodes or 2) the number of electrodes at which at least one spike was detected. For each array, histograms of cluster sizes were obtained using linear binning, normalized to probability density distributions, and plotted in log-log coordinates for further analysis. Fits of exponential and power-law functions were calculated starting from a minimal cluster size of 3 spikes and ending with ≤40 spikes for Δt values between 1 and 8 ms in four recordings and between 6 and 16 ms in one recording, representing the Δt range in which we detected a power law. The resulting R2 values were averaged across all Δts. In the case that the estimated density plot was not continuous because of an insufficient number of data entries, the maximal cluster size was <40, in which case we truncated the rightmost end of the corresponding plot. This fitting range was optimal for the analysis of power law statistics, because cluster sizes with ≤2 spikes were not always distributed on a straight line (e.g., with Δt ≤2 or ≤3 ms), whereas the distributions of cluster sizes >40 spikes were either curved too (e.g., when Δt ≤2 ms) or sometimes undersampled (when Δt reached values of ∼8 ms), reducing the R2 values of the linear and exponential fits. Preliminary analysis showed no difference in the appearance of a power law and its concomitant exponent across the two definitions of cluster sizes, and thus to characterize the cluster size distributions, we restricted our analysis to the number of spikes, i.e., to definition 1. At least two spikes needed to occur within a sequence of active bins to be counted as a cluster, which excludes from the analysis single, isolated spikes. Spike clusters were considered to be neuronal avalanches if their corresponding size distribution obeyed a power law (i.e., the distribution in a log-log plot followed a straight line and power law fitted better than exponential function). The lifetime of a spike cluster was defined as the length of the uninterrupted sequence of active bins, i.e., the count of bins within the cluster multiplied by Δt. Lifetime distributions were plotted in log-log coordinates, and both power law and exponential fits were computed for lifetimes between 1 and 40 active bins and for bin sizes between 1 and 8 ms (6 and 16 ms for 1 dataset; always in 1 ms steps).
As described for spike clusters, spatiotemporal LFP clusters were also obtained by concatenating active time bins bracketed by at least one empty time bin on each side. For LFP clusters, Δt ranged from 1 to 32 ms, and the LFP cluster size was defined as the number of electrodes for which an LFP event was detected. The histograms of LFP cluster sizes were obtained with logarithmic binning, normalization, and studied in log-log coordinates using linear regression analysis.
Normalized auto-correlation histograms (ACHs) and cross-correlation histograms (CCHs) of spontaneous spiking activity were computed first for each unit and each pair of units, respectively. We averaged the histograms across all units (ACH) and all possible pairs of units (CCH). Consequently, we obtained only one ACH and one CCH per dataset. Normalization was achieved by replacing coincidence counts with Pearson's r computed on the trains of 1s (spike) and 0 s (nonspike) binned to 1 ms precision. The width of an ACH or CCH was estimated at half-height, i.e., at the midpoint between its baseline and the maximum. In all cases, the investigated auto-correlation window was ±50 ms and cross-correlation window was ±100 ms.
To test whether subsampling of spiking events can account for the absence of a power law in some of our recordings, we randomly removed spikes from the recordings with power law and recomputed distributions of cluster sizes (also known as spike-train thinning). Spiking events (i.e., actions potentials) were removed randomly and independently from each unit such that, as a result, only 50, 25, or 10% of the original data were sampled. This subsampling procedure was repeated 10 times for each size of the subsample. Power law and exponential functions were fitted for cluster size distributions (see above) and for Δts (2–8 ms for 4 datasets and 6–16 ms for 1 dataset), where we had previously found a power law. We subsequently averaged across all 10 subsampling repetitions, all investigated Δts, and all tested recordings. To avoid any bias toward a better power law fit because of a distribution that was bimodal, we truncated the tales of all bimodal distributions such that only the single-curved part of these distributions were fitted. Consequently, we obtained two values per subsampling size, one for the average power law fit and one for the average exponential fit, which we compared across different levels of subsampling.
An example of the spontaneous activity recorded from 72 units in parallel is shown in Fig. 1. In the same figure we also show the procedure for extracting the spatiotemporal spike clusters, i.e., candidate avalanches, defined as sequences of active bins delineated by blank bins (see the zoom-in). The same procedure was used to extract LFP event clusters.
LFP cluster size distributions
The cluster sizes of LFP events from spontaneous activity in slice cultures, acute slices, and in vivo from young, anesthetized rats and awake monkeys have been found to distribute according to a power law (Beggs and Plenz 2003; Gireesh and Plenz 2008; Petermann et al. 2009; Stewart and Plenz 2006). In this study, a similar power law distribution was found for LFP cluster sizes recorded during spontaneous activity in primary visual cortex of the anesthetized cat. The density of LFP cluster sizes, P(n), when plotted in log-log coordinates, followed a straight line up to the cluster size of n = 16 (Fig. 2A), which was the largest size that is expected with a 4 × 4 array provided that, within a cluster, LFP events are unlikely to recur, i.e., rarely return back to the electrode at which they once already occurred. This distribution is characteristic of the power law function, P(n)∼nα, with a cut-off at the maximal cluster size, and for which the exponent α indicates the slope in the log-log plot. As shown originally for LFP events in vitro (Beggs and Plenz 2003), the power law in vivo, as well as the cut-off values for LFP cluster sizes, remained also robust across different bin sizes, Δt (in 7 of 8 recordings).
The exponent of the power law allows for a distinction between different types of dynamical systems that may give rise to avalanche-like behavior (i.e., each dynamical system displays a unique exponent) and thus restricts potential mechanisms that may be responsible for the generation of the power law (see Plenz and Thiagarajan 2007). Previous estimates of the exponent α for neuronal avalanches both in vitro and in vivo showed a monotonical increase of α from about −2.2 to −1.2 with an increase in Δt (Beggs and Plenz 2003; Gireesh and Plenz 2008; Petermann et al. 2009; Plenz and Thiagarajan 2007;Stewart and Plenz 2006), and when Δt was chosen as the average time delay between successive events on the array, α was found to be close to −1.5. Similarly, in this study, the slope monotonically increased with Δt from −2 to about –1.2 (Fig. 2B; n = 7 recordings). However, the slope remained largely above –1.5 for Δt >2 ms. Therefore the slopes found in this study in vivo were more shallow from those found previously in vitro (Beggs and Plenz 2003), even though the same interelectrode distance of 200 μm was used in both cases.
Computer simulations of cascading neuronal activity suggested that the distributions of lifetime for clusters of short duration (i.e., less than ∼10Δt) should follow a power law too (Eurich et al. 2002), with an exponent close to −2.0 (Teramae and Fukai 2007; Zapperi et al. 1995). Previous in vitro studies on LFPs reported results consistent with these predictions, showing a power law for the initial part of the lifetime distributions (Beggs and Plenz 2003). We obtained different results for the lifetime distributions of LFP clusters in vivo. For all recordings and with Δt = 1 ms, the entire lifetime distribution displayed curvatures, but on the other hand, deviated also from a simple exponential distribution, showing a substantial increase in long lifetimes relative to what would have been expected if LFP events clustered by chance (Fig. 2C). This finding is consistent with earlier reports of lifetime distributions in vitro (Beggs and Plenz 2003), which also differed from an exponential decay. Moreover, in this study, lifetime distributions of LFP events did not scale, i.e., did not collapse, with an increase in Δt (Fig. 2D), a result that would be expected from theory (Eurich et al. 2002) and that has been already shown experimentally in vitro (Beggs and Plenz 2003).
Spike cluster size distributions
Next, we examined whether spike clusters in vivo exhibited a power law. We varied Δt from 0.5 to 32 ms and found characteristic spike cluster distributions, which depended on Δt. An example distribution computed for different bin sizes is depicted in Fig. 3A (cat 1). For small bin sizes (Δt < 1 ms), the cluster sizes were distributed across a steep, curved line, which had a tendency to become shallower and straighter as the bin size increased. Starting with Δt = 1 ms (in some cases 2 ms; Fig. 3A) the distributions showed a power law that stayed ≤7 or 8 ms. At Δt ≥ 9 ms, the distribution became bimodal with curvatures at the initial part and a horizontal tale, i.e., medium-size and large clusters being about equally probable. In the range 1–7 or 8 ms, with the gradual increase in Δt, the slope of the fitted line (i.e., the exponent of the power law) also gradually increased. This can be explained by the increased likelihood to concatenate spike events (see methods). In one dataset, the power law was detected only between Δt = 6 ms and Δt = 16 ms.
Overall, we found power law distributions of avalanche sizes in five of eight datasets independently of whether the recorded segments were short (cats 1 and 3) or long (cat 2), and the results were similar to a previous report in dissociated cultures (Pasquale et al. 2008), whereby the detection of a power law in spiking activity also depended on Δt. For each of these five experiments (recordings), we fitted the resulting size distributions for different Δts (6–16 ms in 1 recording and 1–8 ms in all other recordings) with an exponential and with a power law function, averaged across all tested Δts and compared the goodness of the two fits (i.e., averaged R2 values). A good fit to power law suggests correlated spikes (Bak et al. 1988). A good fit to an exponential function may suggest either independence between neurons analogous to a Poisson process (Shadlen and Newsome 1998) or subsampling of event clusters, whose true size distribution exhibits a power law but is inaccessible because of a small count of recorded action potentials (Petermann et al. 2009; Priesemann et al. 2009). In these recordings, a power law function fitted the data significantly better than an exponential function (t-test, all t > 2.65; all P < 0.001; df ≥ 13; Fig. 3B), for which the distributions of spike cluster sizes are plotted in Fig. 3C using Δt = 3 ms for four recordings and Δt = 9 ms for one of them. One pair of recordings (cat 1) was obtained simultaneously from the same cat from arrays positioned in the same brain area a few millimeters apart. This suggested that the power law was present simultaneously in two spatially separated arrays. Two other recordings using the same array (i.e., the same position in the cortex) were obtained from another cat (cat 2) at two different instances in time (separated by ∼15 min), which suggested temporal stability for at least a short period of time. The overall linear fit to distributions of event sizes in log-log plots was Rpower2 = 0.93 ± 0.01 (SD), which exceeded considerably the quality of the exponential fit (Rexpon2 = 0.85 ± 0.03; see also Fig. 3B). Therefore similarly to the results obtained with LFPs, the spiking events, which sample much smaller portion of neuronal activity, can show power law distributions of spatially distributed events. The cautionary note is, however, that this result is not obtained in every recording.
Exponent of the power law
A representative result for the estimated exponents of the spike cluster distributions is shown in Fig. 4A, where we plotted the changes in the estimated exponent as a function of the bin size, Δt, for four recordings in which we found a power law distribution for Δt = 1–8 ms (the 5th recording with a power law at Δt = 6–16 ms is not shown). In accordance with the findings based on LFPs, the slope grew with the increase in Δt. This result is expected because the concatenation of spikes at large values of Δt necessarily increases the likelihood of large spike clusters. We also estimated the value of Δt that minimizes the decomposition and concatenation of clusters. To this end, we used the methods of Beggs and Plenz (2003) and computed simply the gross average of ISIarrays (see methods; Fig. 4B). For the recordings that exhibited a power law, the resulting value, avgΔt, varied across different arrays but stayed within the range of Δt (close to the lower limit) where we found a power law, with the averages ranging between 1.5 and 5.9 ms (2.84 ± 1.8). Importantly, when we computed the exponents of event size distributions only for these optimal values of Δt, the values of exponents were remarkably similar across all five recordings for which we have previously established reliable power laws. These values were −1.77, −1.78, −1.81, −1.82, and −1.83, with an average of −1.8 ± 0.03 (Rpower2 = 0.98 ± 0.01, Rexpon2 = 0.90 ± 0.01, for the 5 fitted lines). Consequently, when these distributions of event sizes were plotted on the same graph, they largely overlapped (Fig. 4C). These values of the exponents were lower than those found either in this (Fig. 2B) or in previous studies for LFPs (Beggs and Plenz 2003; Gireesh and Plenz 2008) or those found in previous studies for spike clusters (Mazzoni et al. 2007; Pasquale et al. 2008).
Absence of power law statistics in event size distributions
As mentioned, in three recordings (of 8), we did not find sufficient evidence that the sizes of spike events distributed according to a power law (Fig. 3B). One recording without a power law (cat 3) was made almost 3 days (67 h) after another one that has shown a power law. The remaining two datasets were obtained from the same cat (cat 4) about 1 day apart (28 h). This was despite the presence of a power law in the LFP events in all cases (Fig. 2B). This suggests that the power law statistics might be more robustly found at the level of the LFP compared with unit activity, which is more prone to subsampling. Similarly to the datasets with a power law, for small Δt, the distribution of event sizes was steep and curved (although sometimes only to a small degree). Likewise, an increase in Δt made the distributions of both classes of recordings shallower. The main difference between the two groups was detected with the increase in Δt. Only those five recordings classified as containing a power law exhibited consistently a straight line in a log-log plot, whereas the distributions of the remaining three became gradually more and more curved as Δt increased.
These three recordings were fitted with an exponential and a power law function for bin sizes 1–8 ms. Overall, the quality of fit for exponential functions did not differ from that of a power law, when different values of Δt were taken into account (t-test, all t < 1.71, all P > 0.06, df ≥ 13). As expected, for small bin sizes, the power law fitted the distributions of cluster sizes in two recordings better than did the exponential function (Rpower2 = 0.97; Rexpon2 = 0.88 with Δt = 1 ms), whereas the result was reversed for higher bin sizes (Rpower2 = 0.95; Rexpon2 = 0.98 with Δt = 8 ms). As a result, averages of the goodness of power law and exponential fits computed across all Δts were approximately equal for these two recordings. In the third recording, an exponential function described the distribution of event sizes better than a power law across all bin sizes (cat 4 I-b; on average, Rpower2 = 0.92 ± 0.04 vs. Rexpon2 = 0.95 ± 0.07). Note that the qualities of the exponential and power law fits (i.e., averaged R2 values) are higher in these three recordings than in those with consistent power law. This is because the average duration of the recording was much longer in the former case, increasing the total spike count and, by that, the quality of the fit.
Distributions of ISIs
We also studied the presence of a power law in the distributions of ISIs. These distributions were computed either for each unit separately (ISIunit) or for all the units belonging to an entire array (ISIarray). The ISIunit distributions did not exhibit a straight line in a log-log plot in either of our recordings (Fig. 5B). However, for the datasets that exhibited a power law in the distribution of cluster sizes, a power law function fitted also the distribution of ISIunit better than did the exponential function (Rpower2 = 0.93 ± 0.01; Rexpon2 = 0.85 ± 0.03). The opposite was the case for the datasets without a power law in cluster sizes (Rpower2 = 0.91 ± 0.01; Rexpon2 = 0.94 ± 0.01). Therefore, although the distributions of ISIunit were always curved, the curvature was the smallest in the recordings in which the distributions of cluster sizes obeyed a power law.
In contrast, ISIarray distributions were much more consistent with the analysis of event size distributions. In the recordings in which spike cluster sizes distributed according to a power law, the ISIarray distributions exhibited a straight line (Fig. 5C; Rpower2 = 0.96 ± 0.01; Rexpon2 = 0.88 ± 0.03) and vice versa; if a power law was absent in the distribution of cluster sizes, a straight line was also missing in the ISIarray distribution (Fig. 5D; Rpower2 = 0.92 ± 0.01; Rexpon2 = 0.97 ± 0.02). Thus the analysis of ISIarray distributions suggests conclusions similar to those reached after the analysis of the distributions of events sizes.
The distributions of spike cluster lifetimes showed properties similar to those of the cluster size distributions. The result was dependent on the bin size and whether a power law was detected in the size distribution. In all recordings, the distribution of lifetimes was steep and curved at small bin sizes (≤1 ms; see Fig. 6A for Δt = 1 ms) and became gradually shallower with the increase in Δt (Fig. 6B). In those five recordings with a power law distribution of cluster sizes, the lifetime distributions became straighter with increasing bin sizes, but nevertheless, they never exhibited fully a power law (Fig. 6B). The distributions were closest to a power law with bin sizes 6–7 ms (16 ms for 1 recording)—values similar to the maximum Δt for which a power law of cluster sizes could still be detected. For Δts ≥ 7 ms, the lifetime distributions became bimodal. In the three recordings with absence of a power law distribution in cluster sizes, the increase in Δt had the opposite effect on the lifetime distribution. Much like the cluster sizes, these distributions increased also in curvature with larger Δts (data not shown). We also quantified these results: For all recordings in which cluster sizes exhibited a power law, the lifetime distributions for Δts = 1–8 ms (6–16 ms in 1 case) were approximated more accurately by a power law function (R2 = 0.95 ± 0.01) than by an exponential function (R2 = 0.89 ± 0.2). Conversely, in the absence of a power law in cluster sizes (recordings marked with stars in Fig. 6A), the lifetime distributions were described better by an exponential function than by a power law (Rpower2 = 0.92 ± 0.02; Rexpon2 = 0.97).
Furthermore, recordings with a power law also showed longer lifetimes, as can be seen in Fig. 6A. Lifetime distributions with a power law in cluster sizes were less steep and intercepted the abscissa at a later point (23 ± 4 ms) compared with the lifetime distributions of recordings lacking power law statistics (13 ± 3 ms; marked with stars). The only exception was array I-a of cat 3, which showed a power law in the size distribution (only for Δt ≥ 6 ms) but intercepted the abscissa at a value similar to recordings without a power law (13 ms).
In critical systems, events that occur within cascades can be correlated over long distances and often display long-range temporal correlations (Bak et al. 1988). The same was shown in in vitro experiments in which the existence of a power law in neural activity was associated with long-range correlations across neuronal events (Beggs and Plenz 2003, 2004; Mazzoni et al. 2007). To assess the role of correlations in our in vivo experiments, we computed overall ACHs and CCHs by averaging the individual ACHs and CCHs calculated for all possible units and pairs of units belonging to one array (see methods for more details).
In all cases, ACHs and CCHs showed center peaks that indicated correlations (auto- or cross-) higher than expected by chance. The value of the ACH at 0 time delay is by definition always 1.0 and hence was not used for the analysis. Notably, these correlations differed across arrays with and without power law statistics. Correlations were on average stronger (e.g., the center peaks were higher) in arrays with a power law in cluster sizes (r = 0.019 ± 0.008 and 0.007 ± 0.001 for ACHs and CCHs, respectively) as opposed to cases where a power law was absent [r = 0.007 ± 0.002 and 0.005 ± 0.001 for ACHs and CCHs, respectively; for ACH, t(6) = 2.34, P = 0.027; for CCHs, t(6) = 1.85, P = 0.057]. Example ACHs and CCHs for a power law and a non–power law case are shown in Fig. 7.
Moreover, we found a close relationship between the width of the center peak and the presence of a power law. The recordings in which the spike counts showed power law statistics had significantly wider center peaks, both in ACHs and CCHs, than the recordings in which the power law was absent. The width at the half height of the peak was 14 ± 3 ms for ACHs and 18.6 ± 7.4 ms for CCHs in arrays with power law and 8 ± 2 ms for ACHs and 3.7 ± 0.6 ms for CCHs in arrays without power law [for ACH, t(6) = 2.63, P = 0.02; for CCHs, t(4.1) = 4.51, P < 0.01]. Therefore consistent with the longer lifetimes in the recordings with power law in event sizes, the presence of a power law was associated not only with stronger correlations but also with correlations that spanned considerably longer temporal distances.
Subsampling can destroy power law statistics
As shown in previous studies (Petermann et al. 2009; Priesemann et al. 2009), subsampling, i.e., recording an insufficient number of neuronal events, can mask power law statistics of underlying neuronal dynamics. To test whether subsampling could explain the absence of a power law in our data, we randomly removed spikes from the recordings, the effect of which on the distribution of cluster sizes is shown in Fig. 8A for Δt = 3 ms. The straight line of the fully sampled recording gradually turned into a larger and larger curvature as further spikes were removed from the dataset. To quantify these results, we fitted power law and exponential functions for each subsampling class and for each Δt, where a power law was detected in the original data (in 4 datasets Δt = 2–8 ms, in 1 dataset Δt = 6–16 ms). The resulting fits were averaged across all Δts and all recordings with previously established power law distributions for a given level of subsampling, such that we obtained one power law fit and one exponential fit per subsampling size. Subsequently, the two types of fits were compared, and the results are shown in Fig. 8B. The analysis showed that fully sampled data are fitted better with a power law than with an exponential function [t(84) = 8, P < 0.001]. This difference reduces gradually as data are subsampled to a higher degree, the difference being only marginally significant with the sampling of 10% of the original data [t(643) = 1.74, P = 0.08]. These results show that the power law in cluster size distributions breaks down with subsampling, i.e., with random removal of spikes from the original dataset. However, we also found that the presence of a power law was not tightly linked to the average spike density (i.e., spike count per second in a given recording). Although, on average, the spike count per second was higher in the recordings with power law than in those with curved distributions of event sizes, the lowest spike density of all was found in a recording that also showed power law statistics. Thus power law can also be found with a low level of ongoing activity.
Identifying a critical state of a system, as it has been defined by Bak and others (Bak and Paczuski 1995; Bak et al. 1988; Jensen 1998), is a nontrivial task. Finding a power law in a distribution of events is apparently not sufficient proof of criticality, and other measures are necessary. Much importance has been assigned to long-range correlations between events and to scale invariance, i.e., persistence of the power law across different spatial and temporal scales. First, evidence suggesting the existence of critical dynamics in neural tissues, characterized as neuronal avalanches, was reported from the statistics of the spontaneous LFP activity in cultured neuronal networks and in acute cortical slices (Beggs and Plenz 2003, 2004). In particular, a power law distribution of the event sizes was found, irrespective of whether these sizes were defined as the number of LFP events or the sum of the amplitudes of these events. In addition, the distributions of the events' lifetimes partially obeyed a power law. These power law functions exhibited reliably exponents of −1.5 for event sizes and −2.0 for lifetime distributions, the empirically obtained values matching closely the theoretical predictions (Eurich et al. 2002; Zapperi et al. 1995). Also, the authors described long-range temporal correlations and robustness of power law statistics, when both spatial and temporal scaling operations were applied. The same results were recently obtained from ongoing activity in awake monkeys, suggesting the presence of critical dynamics in the nonanesthetized brain (Petermann et al. 2009). A power law organization of nLFP clusters was accompanied by correlations spanning several seconds in time. This organization also withstood scaling operations in nLFP amplitude threshold, which systematically removed ≤95% of nLFPs.
Our study attempted to expand these findings in two directions. First, we studied the properties of the spontaneous activity in the visual cortex in vivo, and second, we analyzed much more sparse signals based on neuronal spiking activity rather than relying on the LFP. Our results provided support for the conclusions offered by Beggs and Plenz (2003, 2004) and Petermann et al. (2009), but also provided additional information because of several differences between the in vitro and our in vivo results and also between LFP-based and spike-based analyses. We found that the spontaneous spiking activity in cat visual cortex can, under anesthesia, exhibit different types of distributions of spiking events. The distributions in some cases follow a power law, but in other cases, lack such statistics. Importantly, when the power law property was detected, it was found only across a certain range of values for the scaling parameter Δt, which is in contrast to our LFP analysis, where the power law was invariant to the chosen Δt, a scale-invariance characteristic for neuronal avalanches (Beggs and Plenz 2003). This result is similar to a previous study in vitro, in which power law statistics also depended on Δt (Pasquale et al. 2008). When in our data the value of Δt was chosen such that it was most optimal for minimizing concatenation and decomposition of events, the exponent of the power law distribution was highly consistent across different experiments. In these recordings, in which power law statistics was found for event sizes, lifetime distributions were the closest to, although did not fully exhibit, the power law statistics. A similar case was with the ISI distributions computed for an entire array (ISIarray), which also exhibited a power law in cases in which the event sizes did the same. Therefore we found evidence supporting the idea that the spontaneous activity in anesthetized brains is driven by processes that manifest power law statistics consistently across different measures and that are possibly based on cascades of neuronal events or so-called neuronal avalanches.
The exponent of avalanche size distribution was, however, somewhat smaller than that reported in previous studies and expected theoretically. We find consistently the value around −1.8 (as opposed to −1.5 reported previously). On the other hand, we found an exponent larger than −1.5 for LFP clusters for most of the investigated values of Δt (Fig. 2B). It is presently not clear how these differences in slope are related to the type of used signals (spikes vs. LFP), why these results differed from the studies in vitro, and whether anesthesia played a role.
In three of eight recordings we did not observe a power law, but instead, a distribution of event sizes that exhibited an exponential-like distribution. The reasons for this variability in the results also need to be identified. One possible explanation of this finding is subsampling: The invariant presence of power law in LFPs suggests that avalanche-like behavior may exist also in all recordings of spiking activity but, because of the relatively small number of spikes recorded simultaneously, this property may not be detected easily or always (Priesemann et al. 2009). Our subsampling analysis points in a similar direction. Subsampling because of a too small number of electrodes (8) may be the reason that another recent study that recorded spiking activity in vivo did not find a power law distribution of event sizes (Bédard et al. 2006). However, there are also other possible explanations; as in this study, the electrode tips spanned much larger distances (1 mm), and the electrodes were aligned linearly.
At present, we cannot distinguish the subsampling explanation of the lack of power law from changes in cortical states. Spontaneous changes in the depth of anesthesia (Herculano-Houzel et al. 1999), possibly the result of a change in the concentration of acetylcholine (Rodriguez et al. 2004), may be a possible factor that underlies these changes in the properties of neuronal activity. Our findings that the shapes of auto- and cross-correlograms differed in recordings with and without power laws possibly suggest a role of state changes. This, however, does not explain why changes in power law would be observed only in spiking and not in LFP activity. A complete explanation will require accounting for this feature of our results. A similar problem is with a third class of explanation that may hypothesize changes in the statistical properties across different cortical layers. Again, if responsible for the lack of power law, one would expect this factor to affect LFP and spiking activity equally. This conclusion is supported by the finding that one of our arrays showed the differences in the presence of power law and correlations at two different time points, without moving the array to another position.
In this study, we investigated one more prediction. It is expected that cross-correlation analysis of avalanche activity shows correlations between events on a broad temporal scale (Beggs and Plenz 2003). Thus if avalanches are found in neuronal spiking activity, cross-correlation should not show a narrow center peak characteristic for neuronal synchrony associated with beta/gamma oscillations. Instead, one should observe a broad peak that reflects correlations that occur simultaneously at different time scales. This is what we found. Auto-correlation showed a similar internal structure of spike trains.
The functional implications of a putative critical state in the sensory areas remain elusive. Theoretical work (Kinouchi and Copelli 2006) reported that neuronal networks set at criticality display optimal sensitivity to stimuli, a finding that was experimentally confirmed using organotypic cultures (Shew et al. 2009). The preparation showed a maximum range of responses to stimulation by electric currents when nLFP clusters were distributed according to a power law. Possibly, the presence of power law in the visual cortex has a similar impact on processing visual stimuli.
Mechanisms other than self-organized criticality can generate a power law distribution of event sizes (Newman 2005; Plenz and Thiagarajan 2007). Touboul and Destexhe (2010) attributed the occurrence of a power law in nLFP clusters to stochastic dynamics, although with orders of magnitude steeper slopes. Our data leave open the possibility that the mechanisms that generate power laws in spiking and LFP activity are different, because we find different properties in the two types of responses. Further studies will be needed to identify the mechanisms responsible for power law distributions, or lack thereof, in different measures of neuronal activity.
In conclusion, this study provides evidence, for the first time, that not only LFP and spiking activity in vitro can exhibit power law statistics of event sizes, but the same can be the case for neuronal spiking activity recorded in vivo. These results suggest the possibility that neuronal avalanches are a common component of spontaneous brain dynamics and thus may also have implications in understanding how sensory inputs are represented and processed.
This work was supported by the Hertie Foundation, the Max-Planck Society, a grant from Deutsche Forschungsgemeinschaft number NI 708/2–1, the Division of the Intramural Research Program of the National Institute of Mental Health, and the European Community's Seventh Framework Program (FP7) under grant agreement number PITN-GA-2009-237955.
No conflicts of interest, financial or otherwise, are declared by the authors.
We thank M. Wibral and V. Priesemann for useful discussions and for providing sand-pile simulation data on which we tested our computer code for data analysis.
- Copyright © 2010 the American Physiological Society