Abstract
The human brain is fragile in the face of oxygen deprivation. Even a brief interruption of metabolic supply at birth challenges an otherwise healthy neonatal cortex, leading to a cascade of homeostatic responses. During recovery from hypoxia, cortical activity exhibits a period of highly irregular electrical fluctuations known as burst suppression. Here we show that these bursts have fractal properties, with power-law scaling of burst sizes across a remarkable 5 orders of magnitude and a scale-free relationship between burst sizes and durations. Although burst waveforms vary greatly, their average shape converges to a simple form that is asymmetric at long time scales. Using a simple computational model, we argue that this asymmetry reflects activity-dependent changes in the excitatory–inhibitory balance of cortical neurons. Bursts become more symmetric following the resumption of normal activity, with a corresponding reorganization of burst scaling relationships. These findings place burst suppression in the broad class of scale-free physical processes termed crackling noise and suggest that the resumption of healthy activity reflects a fundamental reorganization in the relationship between neuronal activity and its underlying metabolic constraints.
Introduction
The vulnerability of human cerebral cortex to hypoxia is of paramount importance. Any suspension of the brain's continuous metabolic inflow leaves it susceptible to damage unless supply promptly resumes. Such interruptions occur relatively frequently during complications at birth. Following an initial quiescent period, scalp electroencephalography (EEG) acquired in this critical time window exhibits an abnormal pattern termed “burst suppression,” characterized by the erratic appearance of high-amplitude irregular discharges across most cortical regions punctuating an otherwise low-amplitude background (Niedermeyer et al., 1999). Individual bursts vary greatly in magnitude and shape, ranging from very brief fluctuations barely surpassing amplifier and physiological noise to high-amplitude waveforms lasting several seconds. This pattern of pathophysiological processes, arising from compromise of blood-oxygen supply to the neonatal brain, is termed hypoxic-ischemic encephalopathy (HIE) and is the most frequent severe neurological morbidity in newborns. Outcomes range from complete recovery to permanent neurodevelopmental disability, or even death (Grigg-Damberger et al., 1989; Volpe, 2008). Complete recovery occurs only in cases where burst suppression rapidly resolves and normative, continuous EEG activity resumes within a few hours (Hellström-Westas et al., 1995; ter Horst et al., 2004). Despite its prevalence, the underlying mechanisms of burst suppression in HIE remain poorly understood.
Erratic, burst-like discharges occur in a wide range of physical systems. Examples include the noise emitted by slowly crushed plastic (Kramer and Lobkovsky, 1996), Barkhausen noise emitted by ferromagnets (Spasojević et al., 1996), and violent clusters of earthquakes driven by tectonic plate movements (Sethna et al., 2001). Despite their appearance in such diverse settings, the fluctuations of this “crackling noise” exhibit universal characteristics (Sethna et al., 2001) such as power-law distributions of fluctuation size, implying that they are “scale-free.” More deeply, recent studies have shown that average shapes of crackling noise fluctuations converge toward a simple, universal scaling function at all time scales, despite tremendous diversity of fluctuation morphologies (Baldassarri et al., 2003; Papanikolaou et al., 2011). Analyzing these features yields new insights into underlying mechanisms.
Crackling noise occurs when interactions in a system are in a critical balance between amplification and dissipation (Colaiori et al., 2004). Given the critical balance of excitatory and inhibitory activities in the cortex (Shu et al., 2003), crackling noise might be expected to arise in the brain. Scale-free fluctuations in the spatial extent of neuronal activity, also known as neuronal avalanches, have indeed been reported in cultured tissue (Beggs and Plenz, 2003; Friedman et al., 2012) and in vivo primate recordings (Petermann et al., 2009; Shew et al., 2009). However, despite recent progress (Sporns et al., 2004; Stam and de Bruin, 2004; Kitzbichler et al., 2009; He et al., 2010; Haimovici et al., 2013; Shriki et al., 2013), the existence of such avalanches in humans remains controversial (Dehghani et al., 2012). Here, we show that the human neonatal cortex exhibits a striking example of crackling noise during recovery from an asphyxic insult at birth. Furthermore, the average shape and scaling relationships of the crackling noise fluctuations exhibit a robust reorganization upon the resumption of continuous electrical activity.
Materials and Methods
Data acquisition and preprocessing
We analyzed long-term scalp EEG recordings from all babies admitted to the tertiary level neonatal intensive care unit in Helsinki University Central Hospital due to perinatal asphyxia during a 30 month period from November 2009 to April 2012. Inclusion criteria for these archived data were high-quality EEG data with initial burst suppression within 18 h of birth and sufficient clinical details regarding delivery, drug treatments, other complications immediately following birth, and Apgar scores 1, 5, and 10 min after birth. Burst suppression was identified from contemporaneous clinical records and verified visually by a clinical neurophysiologist (S.V.). Long epochs (40–330 min; Table 1) of relatively artifact-free data were used in this report. Exclusion criteria included the presence of epileptic seizures or excessive high-amplitude, non-neuronal artifacts such as those due to respiration or physical movement. EEG data from a total of 13 neonates (9 male) met these criteria, having presented with either asphyxia at birth (n = 12) or a sudden cardiorespiratory collapse with severe asphyxia immediately following birth (n = 1). In the intensive care unit, clinical management included continuous EEG monitoring and whole-body hypothermia. EEG was recorded using NicOne or Olympic EEG monitors (Natus) at sampling rates 250 and 500 Hz or 256 Hz, respectively, and stored offline for later analysis. We analyzed data from a single biparietal derivation (P3–P4), which is the most common recording configuration in use in the neonatal intensive care setting. Although analysis of the spatial properties of posthypoxic burst suppression would likely be informative, we focus here on the temporal signatures. This approach mirrors treatments of Barkhausen noise in ferromagnets, where a single macroscopic channel is used to record a global measure of the inaccessible microscopic dynamics (Sethna et al., 2001; Zapperi et al., 2005; Papanikolaou et al., 2011). We exported data from the amplifiers into either European Data Format (for NicOne) or ASCII format (Olympic monitor). Preprocessing of data (in MATLAB) was conducted to standardize the sampling rate of all data to 250 Hz, i.e., 500 Hz recordings were downsampled (n = 6) by discarding every second point, and the 256 Hz recording was resampled (n = 1) by using an anti-aliasing (low-pass) filter. The amplitude envelopes of these data were obtained using the magnitude of the Hilbert transform, which was squared to obtain the instantaneous power. Finally, we smoothed the power fluctuations using a 10th-order Savitzky-Golay filter (19 Hz cutoff) to reduce noise on short time scales and assist the convergence of average burst shapes (Papanikolaou et al., 2011). To confirm that the results presented in this report were robust to the choice of filter settings, we repeated all analyses with multiple choices of filter settings for both the raw data and power fluctuations. We also studied the robustness of burst distribution statistics across different choices of these filter settings.
Clinical data summary
Use of these archived, de-identified EEG recordings was approved by the Ethics Committee of the Hospital for Children and Adolescents, Helsinki University Central Hospital.
Statistical characterization
After preprocessing the data, we analyzed each neonatal EEG time series as follows: (1) threshold estimation; (2) burst extraction; (3) calculation of burst areas and durations; (4) fitting burst distributions; (5) calculation of average shapes at different time scales; and (6) quantification of variations in shapes by calculating asymmetry and flattening. Each of these steps is described in detail, below.
Threshold estimation.
We estimated a unique threshold for each time series as that which maximized the number of identified bursts. This heuristic is straightforward and avoids arbitrary or subjective choices. Very low thresholds—such as below the noise level—yield very few bursts. As the threshold rises above the noise floor, suprathreshold bursts emerge, although for small thresholds, many of these will be artificially merged. Very high thresholds fail to identify small bursts and, in the extreme, fail to identify any bursts at all. Hence the number of bursts is a simple unimodal function of threshold with a single maximum (Fig. 1A). We hence used this maximum—unique to each dataset—as the choice of threshold. To reduce computation time, we restricted the set of thresholds tested to 50 quantiles of the data; results are insensitive to the precise number as the relationship between threshold and number of bursts is relatively smooth in all our datasets. Moreover, we verified that all results are insensitive to the precise choice of threshold. Quantitative results, such as the scaling exponent and burst shape, are stable with respect to reasonable changes in threshold. Qualitatively, the existence of power-law scaling and near-universal burst shapes is strongly robust to threshold choice.
Data analysis. A, Data-driven heuristic method to estimate thresholds and identify bursts. The threshold chosen is that which maximizes the number of bursts taken over a set of 50 quantiles of the instantaneous power. B, Example of a single burst spanning 95 data points (380 ms) above threshold (red line). C, Schematic burst with its duration (time above threshold) and size (area under the curve above threshold, shaded). D, Exemplar upper CDF of burst size with fits to Pareto (red), log-normal (blue), Weibull (magenta), exponentially truncated power-law (green), and exponential (light gray) distributions. The exponential distribution's poor fit here is typical of all our data. BA, Burst area. E, Example of a single burst spanning 50 data points (200 ms) above threshold (red). To estimate precise burst start and end points, threshold crossings (red circle) were estimated by linear interpolation between the sub- and superthreshold points, as illustrated in the zoomed inset. F, Two bursts of different original durations (from B and E) interpolated onto the same uniform time axis and rescaled to have unit duration and area.
The present approach for threshold detection yields a single data-driven value and hence allows independent replication, thereby overcoming the established ambiguities associated with visual burst detection (Palmu et al., 2010).
Burst extraction.
A burst is defined as all successive suprathreshold data points, beginning when the signal crosses the threshold from below and ceasing when the signal next crosses the threshold from above. Figure 1B shows an example thresholded burst. This method also yields interburst intervals corresponding to data points between threshold crossings. Bursts of duration <40 ms were not used for further analysis because these bursts are only coarsely sampled by at most 10 samples at 250 Hz, and thus likely susceptible to noise.
Calculation of burst areas and durations.
After thresholding and identifying discrete bursts, we then examined burst size (suprathreshold area under the curve) and burst duration (time between successive threshold crossings), as illustrated in Figure 1C.
Fitting theoretical cumulative distribution functions to burst distributions.
For each set of fluctuation statistics, we calculated the upper cumulative distribution function [i.e., P(X > x)] and fitted candidate distributions to these using the method of maximum likelihood. We tested power law (the Pareto distribution), power law with exponential cutoff, log-normal, stretched exponential (the Weibull distribution), and exponential distributions (Fig. 1D). Following established techniques, each theoretical distribution was fitted to the right-hand side of the empirical cumulative distribution function (CDF), above a lower bound that minimizes the Kolmogorov–Smirnov goodness-of-fit statistic between the power-law model and the data (Clauset et al., 2009). This method of determining the range of the fit from the data strikes a balance between fitting too wide a range (i.e., outside the power-law regime) and too narrow a range (i.e., unnecessarily throwing away data). We used the same fitting range for all five candidate distributions. We estimated a p value for the fitted power law by comparing to 2500 synthetic datasets drawn from a true power law. This estimated the likely deviation between the data and the fitted power law due to finite sampling of a true power law. The p value was taken to be the fraction of synthetic datasets that deviated from the power law by at least as much as the data. A value of p > 0.1 indicates plausibility of the power-law hypothesis (Clauset et al., 2009). We compared the fitted power law with alternative distributions using likelihood ratio tests. Significant deviation of the likelihood ratio from zero was tested using Vuong's methods (Vuong, 1989). For the nested hypothesis of power law versus power law with cutoff (the latter family includes the former), the null hypothesis was that the power law is the best-fitting distribution. For all other tests, the null hypothesis was that both distributions are equally far from the true distribution. In all infants, the exponential distribution was clearly an implausible fit—the CDF in Figure 1D is typical of all datasets—as verified by log-likelihood ratios at least an order of magnitude higher than those of all other fits.
For additional verification of our power-law fits, we also fitted strictly truncated power laws to the burst distributions (Deluca and Corral, 2013). We then used these fits to identify self-consistent ranges of power-law scaling in areas and durations. The strictly truncated power law is essentially the Pareto distribution restricted to a finite range of data bounded by lower and upper cutoffs (cf. the Clauset et al., 2009 method that fits the entire tail above a lower cutoff). These cutoffs were estimated from the data by examining a 100-by-100 grid of logarithmically spaced ranges. For each candidate range, we fitted the strictly truncated power law using maximum likelihood estimates, and generated 50 synthetic datasets. Following the method of Deluca and Corral (2013), we chose the fitted range as the widest (maximizing the ratio of upper and lower cutoffs) that yielded a plausible fit, defined as requiring at least 20% of the synthetic datasets to deviate from the fit by at least as much as the data.
We also analyzed the scaling relationship between burst durations and burst areas. We binned the bursts by duration into 50 logarithmically spaced bins and calculated within each the median duration and median area. Scaling relationships were insensitive to the specific choices of the binning. Plotting these data revealed a straight line in log-log coordinates, indicating a power-law relationship between the two quantities (not to be confused with a power-law distribution of a random variable). We estimated the scaling exponent (slope of the linear relationship in log-log coordinates) using a linear least-squares fit.
Calculation of average shape at different time scales.
To gain deeper insight into the nature of neuronal activity occurring in burst suppression, we next studied the average shape of bursts across a hierarchy of temporal scales. Burst shape analysis is a sharper probe than analysis of size distribution alone, because systems with the same power laws can exhibit different average shapes (Sethna et al., 2001). Scale-dependent deviations in average shapes, such as flattening and asymmetrical skewing (Spasojević et al., 1996; Baldassarri et al., 2003; Colaiori et al., 2004; Zapperi et al., 2005), have been observed in both empirical data and theoretical models, yielding new insights into the underlying mechanisms. To improve the signal-to-noise ratio for the average shape analysis, we pooled bursts over narrow ranges of durations and averaged within these bins (Papanikolaou et al., 2011), rather than averaging only over bursts of exactly the same duration. Hence, the durations were first partitioned onto a finite grid and the average shapes estimated over all bursts within each bin. Here we used a logarithmically spaced partition spanning the observed fluctuation ranges with bin edges at durations of 160 ms, 320 ms, 640 ms, 1.28 s, 2.56 s, and 5.12 s. Durations in the figure captions denote the left endpoints of bins.
To average bursts of different durations within each bin, we first interpolated all bursts onto the lowest duration within the bin. Due to the finite temporal resolution of the data, each threshold crossing occurs between two samples rather than on a single point. We estimated the precise endpoints of each burst by linearly interpolating between the sub- and suprathreshold points on either side of the threshold crossing. Thus, we set each burst to begin and end at the exact threshold, estimated to first order (Fig. 1E). Next we rescaled the bursts to have unit area and unit duration (Colaiori et al., 2004; Fig. 1F). Finally, we averaged together all bursts in each bin to obtain the average shape for that bin.
Quantification of variations in shapes by calculating asymmetry (skewness) and flattening (kurtosis).
To quantify asymmetry and flattening of the average shapes, we estimated measures of skewness (Zapperi et al., 2005) and kurtosis, respectively, as a function of duration. These measures are not to be confused with the skewness and kurtosis (i.e., third and fourth standardized moments) of the time series itself. Skewness Σ is given by
and kurtosis K is given by
where 〈y(t, T)〉 is the average burst shape of duration T,
To estimate skewness and kurtosis trends, we used average shapes calculated on a finer partition of durations than was used for displaying the average shapes (where we used fewer durations for clarity). The linear trends as a function of duration were estimated by the ordinary least-squares method.
Results
Scale-free distribution of burst sizes
We analyzed scalp-EEG data acquired for clinical purposes from 13 term neonates following asphyxia at birth. These data showed the classic temporal sequence of cortical activity during recovery from HIE, namely a period of complete inactivity during which the EEG contains only non-neuronal noise (Fig. 2A), followed by burst suppression (Fig. 2B,C) before the resumption of continuous activity (Fig. 2D). The burst suppression period exhibits highly variable electrical bursts erratically punctuating a low-amplitude, “suppressed” EEG trace (Fig. 2B,C). Whereas the periods of suppression contain low-voltage amplifier noise and non-neuronal artifacts, the bursts reflect intermittent, spontaneous activity in large pools of neurons (Niedermeyer et al., 1999).
From burst suppression to continuous cortical activity in the asphyxic newborn. A, Cortical inactivity immediately after birth asphyxia. B, Sparse occurrences of burst suppression activity. C, Dense burst suppression activity. D, Transition toward normal continuous cortical activity.
Detailed visual inspection of burst power (amplitude squared) revealed irregular size and structure extending across several orders of magnitude (Fig. 3A–C). Therefore, we first analyzed the distribution of burst sizes. Individual bursts were identified by thresholding the data using a simple data-driven heuristic to identify the appropriate threshold (Fig. 1). Cumulative distributions of burst sizes were then derived from the resulting suprathreshold bursts. These distributions were long-tailed, with broad power-law scaling regimes spanning between 2 and 5 orders of magnitude before exponential truncation at their far right tails (Fig. 3D–F). Truncation is common to all natural systems where energy and size are bounded. The complete collection of distributions from all 13 neonates is provided in Figure 4. Approximately one third of our data also exhibited a “knee” before the truncation (denoted by asterisks), which suggests the presence of fluctuations that are larger and more numerous than the background trend. Burst duration CDFs also follow truncated power laws, albeit over a more restricted 1–2 orders of magnitude (Fig. 5).
Scale-free distribution of burst size. A, A 25 min sample recording of burst suppression in neonatal scalp EEG data. B, C, Successive magnifications of instantaneous power show erratic bursts punctuating a low-power noise floor revealing irregular morphology of bursts across several orders of magnitude. D–F, Exemplar upper cumulative distributions of burst area (solid circles) showing lengthy power-law regimes. Solid lines show best-fitting candidate long-tailed distributions: power law (red), Weibull (blue), log-normal (magenta), and exponentially truncated power law (green). Insets show scaling relationship between median burst size and burst duration in double-logarithmic coordinates. G, Scaling relationship between burst area and burst duration for burst suppression (BS, blue) and after the resumption of continuous cortical activity (CTS, red). Crosses denote data points, circles denote means, and lines denote means ± SEMs.
Burst area distributions for all burst suppression recordings. A–M, Upper CDFs of burst area (BA, black) for all infants 1–13, respectively, with fits to Pareto (red), log-normal (blue), Weibull (magenta), and exponentially truncated power-law (green) distributions. Asterisks denote “knees” showing an excess of bursts above the background trend. Insets, Scaling relationships between burst durations and burst areas (black; points are medians after binning both axes), with least-squares linear fit in double-logarithmic coordinates (red). BD, Burst duration.
Burst duration distributions for all burst suppression recordings. A–M, Upper CDFs of burst duration (BD; black) for all infants 1–13, respectively, with fits to Pareto (red), log-normal (blue), Weibull (magenta), and exponentially truncated power-law (green) distributions.
To formally compare candidate probability distributions for these data, we used a hierarchical approach that balances model likelihood with model complexity (Clauset et al., 2009). We compared the widely observed exponential distribution against the standard power-law distribution (also known as the Pareto distribution) and three other heavy-tailed distributions: the power law with exponential cutoff, log-normal, and stretched exponential (or Weibull) distributions (Figs. 3D–F, 4). The likelihoods of all four heavy-tailed distributions easily exceeded that of the exponential distribution. Among the heavy-tailed candidates, model likelihood was highest for the exponentially truncated power-law distribution in 12 of 13 recordings. In the remaining case, the likelihood of the stretched exponential was marginally higher than the truncated power law, but model comparison showed this difference to not be significant. However, when penalized for extra complexity, the simple power-law distribution was more likely than the truncated power-law in three of the datasets (Table 2). In many instances, the match between the empirical and theoretical cumulative distribution functions was extraordinary. In several cases, the fitted distributions showed a close fit to the data >5 orders of magnitude. For the durations, the exponentially truncated power law was also the best-fitting distribution, with the simple power law favored in three datasets (Table 3). Power-law scaling exponents for the fitted area distributions ranged from −2.32 to −1.18 (mean, −1.59; SE, 0.11; 95% CI, −1.38 to −1.81; Table 4), consistent with classic scale-free statistics in critically balanced complex systems (Colaiori et al., 2004). Exponents for the fitted duration distributions ranged from −3.32 to −1.35 (mean, −1.97; SE, 0.16; 95% CI: −2.28 to −1.66; Table 4).
Burst size CDF model-likelihood tests
Burst duration CDF model-likelihood tests
Burst statistics
Another form of power-law distribution is the strictly truncated power law (Deluca and Corral, 2013), which is restricted to a range of data between both lower and upper cutoffs—a narrower range than the tail distributions fitted above. We present these fits overlaid on probability density functions (PDFs) of the burst areas (Fig. 6) and durations (Fig. 7). These PDFs also exhibit striking power-law scaling regimes, complementary to the CDF findings. Burst area fits span up to 7 orders of magnitude, and burst duration fits span up to 3 orders of magnitude. These fits provide additional support for the existence of power-law distributions in our data.
Burst area probability densities for all burst suppression recordings. A–M, Logarithmically binned PDFs of burst area (BA; black) for all infants 1–13, respectively, with fits to strictly truncated power-law distributions (orange). Fits are maximum likelihood estimates and so do not depend on the PDF binning.
Burst duration probability densities for all burst suppression recordings. A–M, Logarithmically binned PDFs of burst duration (BD; black) for all infants 1–13, respectively, with fits to strictly truncated power-law distributions (orange). Fits are maximum likelihood estimates and so do not depend on the PDF binning.
The relationship between median burst size and burst duration also exhibited a striking power law in all neonates, showing a linear regime in double-logarithmic coordinates that in most cases spanned >4 orders of magnitude of burst size (Figs. 3D–F, Fig. 4 insets). Notably, the power-law scaling relationship between size and duration extended to very large scales without truncation in all neonates.
Lengthy EEG recordings acquired immediately after resolution of burst suppression were also available in 10 of these neonates (see Materials and Methods). We extracted bursts using the same thresholding method as for the burst-suppression recordings. We compared the indices of scale-free activity in burst suppression to those derived from the EEG recordings acquired immediately after resumption of continuous activity in the same neonates. Power fluctuations continued to meet statistical criteria for truncated power-law distributions. However, the slope of the scaling relationship between burst size and burst duration showed a significant increase following the transition to continuous EEG (Fig. 3G, p = 0.017, df = 9, t = 2.9, paired t test). Hence, the scaling relationship between burst duration and burst area changes upon cessation of burst suppression and the resumption of continuous cortical activity.
To study the robustness of our results to variations in the methodological details, we tested a range of thresholds and filter settings around the chosen values. For twofold changes in the threshold, the overall shape of the distributions remains largely unchanged, only slightly affecting the smallest bursts (Fig. 8A). The main finding of lengthy scaling regimes is insensitive to the particular choice of threshold. For changes in threshold of ±10%, we extracted the exponentially truncated power-law exponent for burst areas. For 11 of 13 of the burst suppression datasets in Figure 4, the 20% range in threshold yields a <1% change in exponent, showing that our results are robust to choice of threshold. Of the remaining 2 of 13 sets (Fig. 4B, E), the exponent in B changes by 7% and that in E by 20%—in both cases because the threshold lies near a jump in the lower cutoff as determined by the Clauset et al. (2009) algorithm. Note that these two cases both exhibit large bursts above the background trend. Similarly, varying the filter cutoff frequency has little effect on the distributions (Fig. 8B). As is intuitively expected, stronger filtering (lower cutoff frequencies) primarily eliminates the smallest bursts but leaves the bulk of the distribution unchanged.
Robustness of burst distributions to variations in thresholds and filter settings. A, Burst area CDFs for half (blue) and double (red) the chosen threshold (black). B, Burst area CDFs for smoothing filters with cutoffs at 12 Hz (red), 19 Hz (black), and 37 Hz (blue).
Scaling functions of average burst shape
To gain deeper insight into the nature of neuronal activity occurring in burst suppression, we next studied the average shape of bursts across a hierarchy of temporal scales. Scale-specific shapes and their interscale relationships give insight into the mechanisms underlying burst generation (Zapperi et al., 2005; Papanikolaou et al., 2011). Perfect scale invariance occurs when all shapes collapse onto a simple inverted parabola and suggests a critical undamped stochastic process equivalent to a Brownian walk (Baldassarri et al., 2003; Colaiori et al., 2004). Flattening at longer time scales arises in the presence of a simple, constant damping effect. Shape asymmetry, which is often observed in complex physical systems (Sethna et al., 2001; Zapperi et al., 2005; Papanikolaou et al., 2011), has been attributed to more complex, state-dependent effects (see Phenomenological models of burst generation, below).
Average shapes in all of our burst suppression data showed clear leftward asymmetry, increasing in extent toward longer time scales. An exemplar (Fig. 9A) shows a symmetrical, inverted parabolic shape at the shortest time scale (0.16–0.32 s). Thereafter, at successively longer time scales, the average shapes skew progressively to the left. This trend was apparent in formal estimates of shape skewness across successive time scales (Fig. 9B). In contrast, shape kurtosis was approximately scale-invariant (Fig. 9C). Average shapes are also increasingly noisy at longer time scales. This simply reflects the smaller number of bursts in the longer time bins.
Average shape of bursts and their skewness and kurtosis. Top row: Burst suppression (BS). A, Exemplar (single subject) average burst shape 〈y(t,T)〉 during burst suppression across a hierarchy of time scales. Colors correspond to successive durations T (red = 160 ms, yellow = 320 ms, green = 640 ms, cyan = 1.28 s, blue = 2.56 s). B, Trend in skewness (Σ). In this case, the trend is positive, corresponding to the increasingly leftward asymmetry of the average shapes. C, Kurtosis (K) remains constant. Lines (in red) are linear least-squares fits. D–F, Grand average shape for burst suppression activity and corresponding skewness and kurtosis. Bottom row: Continuous EEG (CTS). G–I, Exemplar (single subject) average burst shape following resumption of continuous neocortical activity. Positive skewness is significantly muted at long time scales. J–L, Grand average shape for continuous normal activity data with corresponding skewness and kurtosis. M, Slope of leftward skew [Σ1 (S−1)] as a function of burst duration for burst suppression (blue, left) and after the resumption of continuous cortical activity (red, right). Crosses denote data points, circles denote means, and lines denote means ± SEMs.
We then pooled data from all neonates in the sample and studied the grand average of these scale-specific shapes over all burst suppression recordings (Fig. 9D–F). These grand average shapes exhibited a strong and statistically significant leftward skew at longer time scales (p = 0.012, df = 9, t = 5.3, one sample t test), consistent with the effect seen in the individual recordings. Due to the pooling of data—and hence the increase in the number of bursts—the longer average shapes were smoother in appearance in these grand averages than in shapes derived from individual subjects.
Visual inspection of the average fluctuation shapes following the resumption of continuous activity (Fig. 9G–I) suggested a marked reduction in leftward asymmetry. Indeed, the grand average shapes (Fig. 9J–L) showed only a slight leftward asymmetry at the longer time scales. We quantified this trend by calculating the slope of the linear least-squares fit to skewness versus duration. A paired (repeated-measures) comparison of this burst shape asymmetry effect in the 10 neonates for whom we had both burst suppression and continuous EEG revealed this difference in skewness to be significant (Fig. 9M, p = 0.007, t = 3.47, df = 9, paired t test).
In summary, the transition from burst suppression to continuous EEG was accompanied by changes both in the scaling relationship between area and duration, and in the average shape of electrical fluctuations. Leftward skewing at long time scales was markedly muted following the resumption of continuous EEG. Trends in burst shape are summarized in Table 4.
Interrelating scaling relationships
In critical systems, the various scaling exponents are often interrelated (Sethna et al., 2001; Friedman et al., 2012). We sought to identify a scaling relationship between the exponents for area distributions (a), duration distributions (d), and the duration versus area exponent (c). Theory predicts that if the distribution of burst areas scales as ∼s−a, the distribution of durations scales as ∼T−d, and the average size scales with duration as 〈s(T)〉∼T1/c, then the exponents obey the relationship c = (a − 1)/(d − 1). To test whether this relationship holds in our data, we estimated the three exponents over self-consistent ranges of the data identified using the Deluca and Corral (2013) method. Specifically, we used the fitted range of areas to identify the corresponding range of durations, then fitted the duration distribution and duration versus area relationship to obtain a set of three exponents for each recording. Comparing the predicted and measured exponents, we found good agreement within the uncertainties of the estimates (Fig. 10A).
Scaling interrelations. A, Scaling exponents for area (a), duration (d), duration versus area (c), and the predicted relationship (a − 1)/(d − 1). Uncertainties are 95% confidence intervals derived for a and d from the SD of the maximum likelihood estimator (Deluca and Corral, 2013), for c from the least-squares fit, and for (a − 1)/(d − 1) by propagating the uncertainties from a and d. Points are plotted with small horizontal offsets for clarity. B, C, Collapse of average shapes after rescaling raw averages by a function of duration and the duration versus area exponent, T1−1/c, for exemplar burst suppression (BS; B) and continuous (CTS; C) subjects. Colors correspond to those in Figure 9, and we estimated c by fitting over all bursts used in the averages (i.e., 0.16–2.56 s).
Beyond exponent interrelations, another important signature of criticality is the existence of a scaling function that relates the duration versus area exponent to the average burst shape (Sethna et al., 2001; Friedman et al., 2012). Specifically, it is expected that 〈y(t,T)〉 ∼ T(1/c − 1)F(t/T), where F(t/T) is the scaling function describing the average shape independent of scale. Indeed, we found such a scaling relation in our data (Fig. 10B,C). There is some heterogeneity in our dataset: not every recording collapses with this method over the same range of durations. Full exploration of the heterogeneity of scaling regimes will be the subject of future work.
Phenomenological models of burst generation
The processes underlying bursting phenomena have been successfully elucidated in numerous physical systems through the use of stochastic models (Colaiori, 2008). We analyzed several candidate stochastic models to elucidate the basic mechanisms underlying the average burst shapes in our data. In particular, we sought insights into the appearance of leftward skew at long time scales, but did not aim to develop a detailed model capable of reproducing all power-law exponents or the full EEG dynamics. In these simple phenomenological models, the activity of the cortex is summarized by a single scalar term, x, which is subject to random fluctuations. The simplest continuous stochastic process models fluctuations in x that are subject to uncorrelated random perturbations, given by
where ξ is an uncorrelated zero-mean Gaussian process (the Wiener process). Equation 3 models a perfectly balanced process, where the probability of incremental growth of a burst equals the probability of incremental decay. Whereas Equation 3 has traditionally been used to study the motion of particles subject to random perturbations—that is, Brownian motion—in the present setting, the value of x corresponds to the voltage measured at the scalp, and the fluctuations represent the sum total of excitatory and inhibitory influences at any given moment. Such a process is consistent with a balance of excitatory and inhibitory processes in underlying cortex (Shu et al., 2003).
We generated sample solutions of this stochastic model by numerical integration using the Heun algorithm (Rümelin, 1982). Analysis of the displacement squared of these sample walks recapitulated the well known result (Colaiori, 2008) that fluctuations generated by a Brownian walk exhibit a power-law distribution in size (Fig. 11A, inset). These distributions are truncated at the far right tail in our simulations by an upper bound imposed by the finite sample length. We also observed that the average shapes of these fluctuations converged at all scales toward a single scaling function, namely an inverted parabola (Fig. 11A), without skewing (Fig. 11B) or flattening (Fig. 11C).
Model average shapes of bursts and their skewness and kurtosis. A–C, Symmetrical average burst shape given by model of Equation 5. Inset shows power-law regime in burst size over several orders of magnitude. D–F, Symmetrical but successively flattened shapes given by model of Equation 6. Inset shows the effect of sample damping values λ (black = 0.05, red = 0.075, blue = 0.1) on the upper cumulative distributions of burst area. G–I, Leftward asymmetric shapes without flattening given by model of Equations 4 and 6. Colors correspond to successive durations T (red = 160 ms, yellow = 320 ms, green = 640 ms, cyan = 1.28 s, blue = 2.56 s). Simulations used a time step of 0.04 ms, with ξ having SD 0.1 (in the same arbitrary units as x), and initial conditions x(0) = 0 and λ(0) = 0.
Although this simple random walk captures the power-law scaling seen in our data, none of the empirical burst suppression recordings showed scale-invariant burst shapes. We hence considered more complex phenomenological models. The Brownian walk generalizes to an Ornstein–Uhlenbeck (mean-reverting) process through the introduction of a linear restoring force,
where λ is a constant. The addition of a small, negative restoring force corresponds to a weakly damped, near-critical system, and is a classic model for the motion of a particle within a harmonic potential of concavity λ. In the present setting, this term embodies a weak constant bias toward greater inhibition than excitation within neuronal populations, albeit not of sufficient strength to prevent the triggering of burst activity.
This damping term leads to an exponential truncation of the right-hand tail of the fluctuation size distribution (Fig. 11D, inset), together with a flattening of average fluctuation shape at successively longer time scales (Fig. 11F) consistent with prior analytical results (Colaiori, 2008). This flattening is indicative of a subcritical process. However, bursts at all scales remain symmetric (Fig. 11E). The burst-shape asymmetry that is striking in much of our data suggests more complex, state-dependent effects. The posthypoxic-ischemic neonatal cortex is in a state of metabolic flux and nutrient depletion. Thus, it is reasonable to assume that neuronal metabolites, such as extracellular calcium (Volpe, 2008; Amzica, 2009; Ching et al., 2012), deplete rapidly during the course of a burst and replenish during the subsequent quiescent phase of burst suppression, thereby creating a history-dependent effect. Resource depletion is also implicated in mechanisms of self-limiting disinhibited activity (Staley et al., 1998; Tabak et al., 2006). If history-dependent effects increase the inhibitory bias, this could be captured by modifying the damping term λ to vary according to recent burst activity. A model of domain wall dynamics in Barkhausen noise (Zapperi et al., 2005; Papanikolaou et al., 2011) has been shown to yield left-skewed bursts at short time scales by modifying the damping term according to the recent velocity of the system variable. Here, noting that there is no obvious physical meaning of a “velocity” effect in neonatal cortex, we instead modified the damping term to depend upon the recent energy expended by the system,
This history-dependent damping term has an exponentially decaying memory with time constant α1 and thus relaxes back to zero in the absence of bursting activity. Equation 5 can be recast as a first-order differential equation,
to be coupled with Equation 4. The time scale constant α2 determines how quickly burst energy, x2, drives up the damping, whereas α1 captures the relaxation of the damping back to zero during suppression. Fluctuations arising from this process show leftward asymmetry, increasing at longer time scales (Fig. 11G). This is consistent with a fast-out/slow-return effect of activity-dependent damping. By setting α1 = 0.1 and α2 = 30, we quantitatively recovered the trend toward positive skewness at longer time scales observed in our burst suppression data (Fig. 11H).
Discussion
Scale-free processes have been reported in a variety of neuronal recordings (Beggs and Plenz, 2003; Stam and de Bruin, 2004; Kitzbichler et al., 2009; He et al., 2010; Friedman et al., 2012; Fransson et al., 2013; Haimovici et al., 2013), although their presence is not without controversy (Dehghani et al., 2012). Distributions in prior recordings rarely scale beyond 2 orders of magnitude and often have scaling exponents well below −2, out of the range that is classically considered critical. We find that electrical fluctuations in neonatal cortex during burst suppression show a striking match to a power-law distribution, with scaling behavior consistent with criticality and spanning up to 5 orders of magnitude. Moreover, our findings in clinical human data are arguably the most striking examples of scale-free neuronal dynamics reported to date. Indeed, such a broad and striking match to a power-law model is rare in nature (Stumpf and Porter, 2012), even in physical experiments using carefully designed and controlled conditions, such as the study of Barkhausen noise in ferromagnetic samples. This regularity is perhaps even more remarkable in light of the fact that our data were acquired as part of routine clinical care, in a noisy setting and far from ideal laboratory conditions.
Whereas the scale-free distribution of burst size is a classic hallmark of crackling noise, recent research has moved toward the characterization of average burst shape (Zapperi et al., 2005; Papanikolaou et al., 2011; Friedman et al., 2012). We note that the average shapes previously reported for both Barkhausen noise (Zapperi et al., 2005; Papanikolaou et al., 2011) and neuronal tissue cultures (Friedman et al., 2012) span approximately half an order of magnitude. We were able to reliably extract average burst shapes across 1.5 orders of magnitude. The leftward skew of the average shape of bursts in our clinical data is consistent with observations in some physical systems, such as Barkhausen noise (Spasojević et al., 1996; Zapperi et al., 2005; Papanikolaou et al., 2011). However, the broader scaling range of our data gives greater confidence in the progressive increase in skewness at long time scales. The marked reduction in this leftward skew following the resumption of continuous activity suggests that this feature of burst shapes is sensitive to recovery from HIE. To explain the origin of skewing, we developed a novel stochastic model incorporating state-dependent nutrient depletion so that recent energy expenditure drives an increase in damping (Eqs. 5, 6). This hence introduces an activity-driven bias toward greater inhibition than excitation that is activated during the evolution of each burst. This model yielded a distribution of average pulse shapes that skewed increasingly leftward at longer time scales, consistent with our data. Our model therefore provides greater insight into the generation of bursts in hypoxic data and, in particular, the interaction between neuronal activity and recovering metabolic resources.
The signatures of criticality in our data suggest that the transition from complete quiescence immediately following the HIE insult, through burst suppression and finally to continuous EEG activity, may reflect a non-equilibrium phase transition in neonatal cortex. We therefore position neuronal activity in neonatal HIE within a broad class of irregular physical processes, known as avalanches (Bak et al., 1987), crackling noise (Sethna et al., 2001), or Barkhausen noise (Spasojević et al., 1996), depending on where they are observed. The current findings suggest a novel means of monitoring recovery from neonatal HIE in the clinical setting, while also yielding key quantitative predictions that could be tested with combined neurophysiological and metabolic recordings in a nonclinical, experimental setting. Such experiments, coupled with detailed physiologically based modeling, will be crucial for explaining the differences in dynamics between burst suppression and the recovery phase. Ideally, such modeling would also uncover a physiological control parameter whose trajectory explains the transitions into and out of burst suppression. The precise physiological mechanisms remain to be determined that explain how the dynamics reorganize from burst suppression, exhibiting lengthy power laws in burst area with duration-dependent average burst shapes, to the recovery period where the burst shapes are symmetrical and scale-invariant. Further analysis of the continuous phase will be important in this regard. Scale-free dynamics in the recovery period are also consistent with the existence of long-range correlations as identified using detrended fluctuation analysis in adult EEG (Linkenkaer-Hansen et al., 2001; Palva et al., 2013). Such analyses will be the subject of future work.
Our finding of scale-free dynamics in burst suppression EEG provides dynamical constraints on computational models of this phenomenon. In particular, models should be poised at the cusp of a phase transition where scale-free dynamics emerge, as found in studies of neural networks (Buice and Cowan, 2007; Petermann et al., 2009; Shew et al., 2009) and branching processes (Beggs and Plenz, 2003; Friedman et al., 2012). In contrast, analysis of burst suppression in propofol anesthesia suggests that it exhibits characteristic oscillatory time scales (Ching et al., 2012), including alpha activity and tightly bounded burst amplitudes, for which we do not find evidence in our HIE burst suppression data. Our findings thus exclude models that generate sequences of alternating bursts and suppression by switching between high-amplitude and low-amplitude oscillatory dynamics, as occurs in models positioned close to critical points such as Hopf bifurcations. Informing detailed neuronal models of the cortex at the microscopic scale with insights from bifurcations in biophysical and phenomenological models at the macroscopic scale has been previously achieved for seizures (Lopes da Silva et al., 2003; Breakspear et al., 2006) and the human alpha rhythm (Freyer et al., 2009, 2011, 2012; Roberts and Robinson, 2012). Extending this approach to model burst suppression in HIE is to be the subject of future work.
Burst suppression arising in the unconscious state following a profound metabolic disturbance in cortex can be contrasted to epilepsy, a constellation of diverse neurological disorders in which paradoxical electrical activity arises during sleep or waking, often throughout life, and which may or may not disturb consciousness. Epileptic seizures typically show strong rhythmic content reflecting pathological synchronization, rather than the erratic alternating periods of bursting and quiescence seen in burst suppression. Seizures have been subjected to numerous computational analyses, which typically suggest they reflect a nonlinear bifurcation (Robinson et al., 2002; Lopes da Silva et al., 2003; Breakspear et al., 2006), offering novel opportunities for seizure control (Kramer et al., 2006; Kim et al., 2009). In contrast, there has been little computational treatment of burst suppression following hypoxia, leaving clinical innovation to be guided principally by trial and error. Similar to studies of epilepsy, it is likely that physiologically based modeling of posthypoxic burst suppression will improve our basic understanding, as has begun for anesthetic-induced burst suppression (Ching et al., 2012), and potentially inform clinical interventions.
Our analysis of single-channel EEG—reflecting the mass action of large neural populations—mirrors the analysis of crackling noise in magnets, where typically a single pick-up coil is used to record the mass action of microscopic ferromagnetic domains (Sethna et al., 2001). Although the analysis of the spatial properties of posthypoxic neonatal burst suppression would likely be informative, we focused on temporal properties because inspection of other electrode locations suggested that bursts typically appeared in all scalp locations. Characterization of the spatiotemporal dynamics of posthypoxia burst suppression would require a far denser array of scalp electrodes than is usually used in neonatal intensive care settings (Odabaee et al., 2013). Recent analysis of invasive recordings of burst suppression induced by propofol anesthesia indeed uncovered intriguing spatiotemporal dynamics (Lewis et al., 2013). However, burst suppression following propofol exhibits characteristic scales during both the active and quiescent phases (Ching et al., 2012) and, thus, is categorically different in nature from that following perinatal hypoxia, suggesting different underlying mechanisms.
Our neonatal EEG data provide a dynamic picture of how an otherwise healthy cortex responds to an abrupt interruption of its metabolic supply. It thus provides a unique window into the coupling between neuronal activity and its metabolic underpinning, a phenomenon of very broad importance (Raichle, 2006). The crackling noise that we presently report likely reflects a putative homeostatic response to metabolic depletion, not unlike the erratic, bursty activity observed in other systems faced with resource scarcity as diverse as animal foraging (Viswanathan et al., 1999), gaze allocation in human vision (Brockmann and Geisel, 2000), and financial market volatility (Plerou et al., 2002). Burst suppression also arises in response to a broad variety of other neurophysiological challenges (Ching et al., 2012). Hence, these clinical neurophysiological data may greatly inform our understanding of scale-free dynamics and energy constraints in the human brain which, in turn, holds potential for advancing clinical diagnosis and management (Iyer et al., 2014).
Footnotes
S.F. was supported by a Career Development Award from the National Health and Medical Research Council of Australia. S.V. was funded by European Community's Seventh Framework Programme European Community FP7-PEOPLE-2009-IOF, Grant Agreement 254235.
The authors declare no competing financial interests.
- Correspondence should be addressed to Prof. Michael Breakspear, QIMR Berghofer Medical Research Institute, 300 Herston Road, Herston, QLD 4006, Australia. Michael.Breakspear{at}qimrberghofer.edu.au
