## Abstract

The framework of criticality provides a unifying perspective on neuronal dynamics from *in vitro* cortical cultures to functioning human brains. Recent findings suggest that a healthy cortex displays critical dynamics, giving rise to scale-free spatiotemporal cascades of activity, termed neuronal avalanches. Pharmacological manipulations of the excitation-inhibition balance (EIB) in cortical cultures were previously shown to result in deviations from criticality and from the power law scaling of avalanche size distribution. To examine the sensitivity of neuronal avalanche metrics to altered EIB in humans, we focused on epilepsy, a neurological disorder characterized by hyperexcitable networks. Using magnetoencephalography, we quantitatively assessed deviations from criticality in the brain dynamics of patients with epilepsy during interictal (between-seizures) activity. Compared with healthy control subjects, epilepsy patients tended to exhibit a higher neural gain and larger avalanches, particularly during interictal epileptiform activity. Moreover, deviations from scale-free behavior were exclusively connected to brief intervals at epileptiform discharges, strengthening the association between deviations from criticality and the instantaneous changes in EIB. The avalanches collected during interictal epileptiform activity had not only a stereotypical size range but also involved particular spatial patterns of activations, as expected for periods of epileptic network dominance. Overall, the neuronal avalanche metrics provide a quantitative novel description of interictal brain activity of patients with epilepsy.

**SIGNIFICANCE STATEMENT** Healthy brain dynamics requires a delicate balance between excitatory and inhibitory processes. Several brain disorders, such as epilepsy, are associated with altered excitation-inhibition balance, but assessing this balance using noninvasive tools is still challenging. In this study, we apply the framework of critical brain dynamics to data from epilepsy patients, which were recorded between seizures. We show that metrics of criticality provide a sensitive tool for noninvasive assessment of changes in the balance. Specifically, brain activity of epilepsy patients deviates from healthy critical brain dynamics, particularly during abnormal epileptiform activity. The study offers a novel quantitative perspective on epilepsy and its relation to healthy brain dynamics.

## Introduction

Epilepsy, a chronic neurological disorder that affects ∼1% of the world's population (World Health Organization, 2015), is associated with an altered excitation-inhibition balance (EIB) in the cortex (Eichler and Meier, 2008; Fritschy, 2008; Badawy et al., 2012). The manifestation of hyperexcitable networks and the underlying epileptic processes are mediated by changes in numerous physiological and dynamical factors of both the excitatory and inhibitory circuits (Wendling et al., 2005; Eichler and Meier, 2008; Fritschy, 2008). Accordingly, it has been suggested that cortical excitability in patients with epilepsy is constantly changing (Badawy et al., 2012). Along with the occurrence of recurrent seizures (Banerjee et al., 2009), interictal (between-seizures) periods are likewise characterized by distinct excitability features (McCormick and Contreras, 2001; Badawy et al., 2012). A magnetoencephalography (MEG) or EEG trace often shows interictal epileptiform activity (IEA), namely, abnormal waveforms not associated with seizure symptoms. However, it is still subject to debate whether interictal spikes forerun, drive, or rather protect against the outburst of seizures (Rogawski, 2006; Staley et al., 2011).

In the past decade, numerous studies have presented supporting evidence that a healthy cortex displays near-critical dynamics; namely, it operates on the border between premature termination and runaway explosive growth of neuronal activity, thereby enabling efficient information propagation (Beggs and Timme, 2012; Shew and Plenz, 2013; Massobrio et al., 2015). A critical state is accompanied by scale-free cascades of activity, termed neuronal avalanches (Plenz, 2012). The sizes of these cascades obey a power law distribution with a specific exponent of α = −3/2. These cascades have been successfully described by a critical branching process (Harris, 1989). Accordingly, the gain of neural systems is reflected through the branching parameter, σ. In critical systems, σ = 1; namely, one event in a cascade leads, on average, to a single future event, whereas subcritical dynamics correspond to a lower gain, and supercritical dynamics to a higher gain.

It has previously been shown, both theoretically (Poil et al., 2012) and experimentally, that avalanche analysis is sensitive to EIB manipulations. When the ratio of excitation to inhibition in cortical cultures was pharmacologically disturbed by the introduction of antagonists of fast glutamatergic or GABAergic synaptic transmission, the systems deviated from criticality (Shew et al., 2009). In systems with reduced GABA inhibition, the power law was destroyed, and the avalanche size distribution became bimodal with more large-sized avalanches than expected from a power law distribution, reflecting supercritical dynamics. Consistent with these findings, sleep deprivation, which was shown to increase cortical excitability in humans (Huber et al., 2013) and is known to trigger epileptic seizures (Kotagal and Yardi, 2008), was also demonstrated by neuronal avalanche analysis to correlate with deviations toward supercritical dynamics (Meisel et al., 2013).

Can the neuronal avalanche analysis provide a sensitive tool to assess the constantly changing excitability of the epileptic brain? In this study, we examined neuronal avalanches in interictal resting state activity of patients with epilepsy, vis-à-vis healthy control subjects, with the aim to monitor deviations from the critical state. The neuronal dynamics of epilepsy patients can potentially deviate from near-critical brain dynamics, even in prolonged interictal periods during which the patients do not experience seizures. Another central aspect is examining whether IEA, as opposed to seemingly “normal” activity within the interictal periods, indicates a difference in the operating distance from criticality. This may provide mechanistic insights into the dynamics of EIB in epileptic networks outside of seizures. Particularly, as interictal recordings are routinely used in epilepsy diagnosis and presurgical evaluation to localize epileptogenic zones (Stefan et al., 2003; Wang G. et al., 2011; Lascano et al., 2012), a better understanding of the neuronal dynamics within interictal periods is of prime importance. Finally, because a substantial body of work in basic cognitive research is based on interictal intracranial recordings from patients with epilepsy, an examination of the validity of generalizing from nonepileptiform brain activity of patients with epilepsy to the general population is in itself a significant goal.

## Materials and Methods

##### Participants.

Patients with refractory epilepsy were referred to the Electromagnetic Brain Imaging Unit at Bar-Ilan University for presurgical assessment. Twenty patients of either sex (adults: *n*_{1} = 12, age_{1} = 24.9 ± 6.2 years; children: *n*_{2} = 8, age_{2} = 9.1 ± 3.2 years) were randomly selected from a large centralized database under the inclusion criteria that the patients underwent an MEG recording while awake and that the reviewing neurologist reported abnormal epileptiform activity in the recorded data. The patients were treated by antiepileptic drugs (detailed clinical information is provided at Table 1) and underwent MEG recordings during the evening, without daytime sleep. The study was approved by the Ethics Committee of the referring hospital, Tel-Aviv Sourasky Medical Center, in accordance with the Declaration of Helsinki. Eighteen age-matched healthy control subjects of either sex (adults: *n*_{1} = 12, age_{1} = 25.4 ± 5.6 years; children: *n*_{2} = 6, age_{2} = 10.4 ± 1.5 years) were also recruited to the study and underwent MEG recordings during varying hours during the day or evening. Informed consent for the MEG recordings was obtained from each subject or parent. Two of the patients with epilepsy were younger (patients 9 and 20, Table 1) than the age-matched control subjects and were therefore excluded from between-group comparisons (epilepsy children group: *n*_{2} = 6, age_{2} = 10.4 ± 2.6 years). Another patient (patient 12, Table 1) had insufficient amount of epileptiform activity for purpose of analyses; thus, this patient and her age-matched healthy control subject were excluded from between brain activity period type comparisons.

##### Data acquisition and preprocessing.

MEG recordings were conducted with a whole-head, 248-channel magnetometer array (4-D Neuroimaging, Magnes 3600 WH) in a magnetically shielded room while participants lay supine with their heads positioned in the MEG helmet. Reference coils, located ∼30 cm above the head, were used to remove environmental noise. Accelerometers (Bruel and Kjaer) attached to the gantry were used to remove vibration noise. Five localization coils were attached to the head of each participant before data acquisition to monitor the head position relative to the MEG sensors. Head localization measurements (1 mm precision) were performed several times during long recordings to ensure that the participant remained still. Head shape and the positions of the coils were digitized using a Pollhemus FASTTRAK digitizer. MEG data were digitized with a sample rate of either 678.17 Hz or 1017.25 Hz and bandpass filtered online at 1–200 Hz or 0.1–400 Hz, respectively. The 50 Hz signal from the power outlet was recorded by an additional channel, and the average power-line response to a power cycle was subtracted from every MEG sensor, enabling the noise and its harmonics to be cleaned without the need to apply a notch-filter (Tal and Abeles, 2013).

Data processing and analysis were performed using MATLAB 2011b (The MathWorks) and Fieldtrip open-source toolbox for Advanced MEG Analysis (Oostenveld et al., 2011). MEG data were first cleaned for line frequency, building vibration, and heartbeat artifacts with in-house open-source software (Tal and Abeles, 2013). Additionally, data from up to two malfunctioning MEG sensors were discarded. Data were downsampled to 678.17 Hz if needed and bandpass filtered offline between 1 and 80 Hz. Epochs contaminated by muscle or jump (in the MEG sensors) artifacts were marked and subsequently excluded from all analyses. Independent component analysis (ICA) was performed (Jung et al., 2000) to ensure the removal of all eye movements, blinks, and remaining heartbeat artifacts. ICA components reflecting such artifacts, as determined by visual inspection of the 2D scalp maps and time course of that ICA component, were rejected, and the remaining components were used to reconstruct the data.

Interictal MEG datasets from the patients with epilepsy were reviewed by a neurologist (M.M.), with expertise in MEG and EEG in epilepsy. He identified interictal epileptiform transients within each dataset, which were accordingly defined as IEA periods. The remaining segments from each dataset were thus defined as non-IEA periods. Notably, manual detection by human experts is the gold standard, against which automatic algorithms are benchmarked (Halford, 2009). Nonetheless, for the purposes of this manuscript, it suffices to consider IEA periods as only having higher propensity to consist of epileptiform activity. The segments were marked and tagged and were later analyzed separately. Furthermore, a heuristic counting of interictal epileptiform discharges (IEDs) was obtained by focusing, within the previously identified IEA, on joint extrema suprathreshold events (i.e., peaks) in at least 8 (∼3%) sensors occurring within 150 ms of the first event. Thus, this criterion determined a minimal duration of discharge width combined with a refractory period (i.e., until the next extremum can be categorized as a separate IED). The chosen parameters were found to accord with the waveforms of epileptiform activity in our datasets (Nowak et al., 2009).

##### Signal discretization, cascade size and duration distributions, and power law statistics.

The signal from each sensor was *z*-scored by subtracting its mean and dividing by the SD. The mean for each sensor was calculated over all non-IEAs or healthy preprocessed periods. Both positive and negative excursions beyond the chosen threshold of 3 SDs from the mean for each sensor were identified. A single event was identified per excursion at the most extreme value (maximum for positive excursions and minimum for negative excursions). Importantly, under these parameter and procedure selections, there was no significant difference in event rate between healthy control subjects and epilepsy patients in non-IEA periods (*p* = 0.4; see Fig. 3*A*); thus, these periods were neither favored nor discriminated against in discretization. The time series of events obtained from each sensor was discretized with time bins of duration Δt. The timescale of the analysis, Δ*t* = 2.95 ms, was twice Δt_{min}, which is the inverse of the data sampling rate (Shriki et al., 2013). A cascade was defined as a continuous sequence of time bins in which there was an event on any sensor, ending with a time bin with no events on any sensor. The number of events on all sensors in a cascade was defined as the cascade size. Cascades from IEA and non-IEA periods were collected separately. Notably, if a cascade occurred during a tagging transition between an IEA and a non-IEA period, it was exclusively assigned to the IEA period. Thus, there was no breaking of evolving cascades due to tagging.

According to the theory of critical branching processes, power law behavior is predicted at the critical state (Harris, 1989). The neuronal cortical sources undergo linear mixing at the sensor level as each MEG sensor linearly sums contributions from multiple cortical sources. In a previous study (Shriki et al., 2013), power law distributions in a simulated network were observed at both the source and sensor level, at criticality. Increasing the degree of overlap in sensitivity of nearby sensors caused a more shallow power law exponent, and a slight underestimation of small cascades, while the branching parameter remained close to 1. At small sensor overlap, the power law exponent was similar at both the sensor and source levels and was close to −1.5. The sensor level analysis is therefore informative to assess criticality, and particularly to assess relative deviations, as is done in this study. Yet, inferring the spatial spread and propagation of a cascade at the sensor level is limited by the mixing of sources.

The fit of the avalanche size and duration distributions to a power law were analyzed as described previously (Clauset et al., 2009; Klaus et al., 2011). The candidate distributions were power law and exponential distributions, both characterized by a single parameter (degree of freedom) and log-normal and exponentially truncated power law distributions, both characterized by two parameters. All distributions were limited to the range between a minimum and a maximum size. Power laws were modeled as follows:
Exponential functions were modeled as follows:
Log-normal functions were modeled as follows:
Exponentially truncated power laws were modeled as follows:
where C_{α}, C_{λ}, C_{μ,σ}, and C_{α,λ} are normalization factors.

The parameters *x*_{min} and *x*_{max} were set to include all observed avalanches.

A maximum likelihood estimation was applied directly to the sample of avalanche sizes (and durations). Assuming independence of avalanche sizes (and durations) and a sample of *n* avalanches, the likelihood of a sample of avalanche sizes given the power law and exponential models, and a parameter α or λ, respectively, is the product of individual probabilities of each avalanche size as follows:
While the log-likelihood is given by the following:
The best fit parameters for the power law and exponential distributions (α̂ and λ̂) were calculated by maximizing the log-likelihood as a function of the parameter. To determine whether a power law or an exponential distribution have a higher likelihood to model the data, the log of the likelihood-ratio (LLR) was taken with the best fit parameters as follows:
Thus, a positive LLR indicates that the power law model is more likely, whereas a negative LLR indicates that the exponential model is more likely; for an LLR of zero, neither distribution is more likely. To determine whether the LLR was significantly different from zero, the *p* value of the LLR was calculated as follows:
where: σ^{2} = with and .

The LLR for the comparison between exponentially truncated power law and log-normal distributions can be calculated analogously (Eqs. 7, 8). The conventional choice of an exponentially truncated power law model is motivated by the predicted cutoff of the power law around system size (i.e., number of sensors in the array). Nevertheless, this model is suboptimal, as the exponent is applied over all effective range, therefore affecting the entire fit and the obtained power law exponent. Since log-normal and exponentially truncated power law have an additional degree of freedom compared with the power law and exponential models, the LLR test would have been difficult to interpret if intermixingly compared. Thus, when comparing distributions for each single subject and period type, only models with the same number of degrees of freedoms were tested against each other (Klaus et al., 2011). Additionally, at the group level, we compared all four competing models using Bayesian Model Selection (BMS) (Penny et al., 2004; Penny, 2012). The BMS determined the most likely model for a given dataset taking into account model complexity. In practice, this was obtained by plugging into the BMS open code [part of the SPM12 package (http://www.fil.ion.ucl.ac.uk/spm/)] estimations of log-evidences for each model over subjects and period types. The log-evidence is a measure of model goodness, which can be decomposed into an accuracy (i.e., model fit or likelihood) and complexity (which serves as penalty) terms, where the best model balances the two terms to give the highest log-evidence. The approximations to the log-evidence used in this study are the Akaike Information Criterion (AIC), and the Bayesian Information Criterion (BIC). The AIC and BIC approximate the complexity with the number of parameters and the number of parameters scaled by the log of the number of observations, respectively. These approximations are simpler and potentially less accurate than the negative free-energy bound on the log-evidence (Stephan et al., 2009) but do not rely on additional assumptions regarding the parameters of the fitted models.

Throughout this article, the reported power law exponent, α, is the one fitted to the simple power law model (Eq. 1), and not to the exponentially truncated power law (Eq. 4). This was chosen because of two complimentary reasons: (1) according to the BMS analysis applied at the group level, the power law model is more consistent with our datasets; and (2) although the exponential decay multiplication factor is conventionally inserted to capture the cutoff of the power law behavior at approximately system size (in our case, size of the sensor array), it affects the entire range, thus distorting the power law exponent itself. Nonetheless, the fitted exponential modulation parameter, λ, is a measure of interest and was a subject for group and period type comparisons.

##### Estimation of deviation from the critical neuronal avalanche size distribution, calculation of branching parameter, and avalanche shape collapse analysis.

The nonparametric measure, κ, quantifies the difference between an experimental cumulative density function (CDF) for cascade sizes and the theoretical reference CDF. As for neuronal avalanches, the probability density function (PDF) of cascade size *x* follows a power law with a slope α = −1.5. The corresponding CDF for cascade size, F^{NA}(β), specifies that the fraction of measured cascade sizes *x* < β is a −0.5 power law function, F^{NA}(β) = (1 − )^{−1}(1 − ) *for l* < *x* < *L*. Therefore, κ is expressed as follows:
where β_{k} are *m* = 10 cascade sizes logarithmically spaced between the minimum and maximum observed cascade size. The measure κ was found to be more accurate in measuring deviation from neuronal avalanches than other nonparametric comparisons of CDFs (e.g., Kolmogorov–Smirnov and Kuiper's tests) (Shew et al., 2009).

The branching parameter, σ, representing the gain of the system, was estimated by calculating the ratio of the number of events in the second time bin of a cascade to that in the first time bin. This ratio was averaged over all cascades for each participant and for IEA and non-IEA periods separately, with no exclusion criteria (Shriki et al., 2013) as follows:
where *N*_{av} is the total number of avalanches in the particular dataset and *n*_{events} represents the number of events in a particular bin.

The values of α and σ vary with the time scale of analysis of Δt. Consistent with Shriki et al. (2013), we find that for Δt close to 3–4 ms the obtained values, are very close to the ones predicted for a critical branching process. Here, the chosen timescale is 2.95 ms. Additionally, the obtained values were only weakly sensitive to changes in the threshold. As was shown previously, for different choices of the threshold, there are Δt values that are consistent with the branching process, and those Δt values are all in the same ballpark (Shriki et al., 2013).

Importantly, because IEA and non-IEA periods have different durations and rates (see Fig. 3*A*) across patients, it was necessary to establish a more valid basis for comparison; we therefore used an identical number of avalanches (*N*_{samples} = 3500), randomly selected from each dataset (i.e., IEA periods, non-IEA periods, and full recording of patients with epilepsy, and healthy control subjects) to construct size distributions and assess the branching parameter. One patient (patient 12, Table 1), who did not have a sufficient IEAs to obtain a substantial amount of avalanches, was excluded from all IEA and non-IEA comparisons (*n* = 19; *n*_{1} = 11, *n*_{2} = 8), as well as from all IEA, non-IEA and healthy control comparisons (*n* = 17; *n*_{1} = 11, *n*_{2} = 6). Noticeably, for all datasets, no significant differences were found between α, κ, and σ based on the full recordings and those based on the 3500 sampled avalanches (*p* was in the range of 0.8–0.99 for all parameters and datasets). The statistical tests used throughout this study were *t* test, when comparing epilepsy patients and healthy control subjects, and one-way ANOVA, when comparing IEA and non-IEA periods from patients to healthy subjects. *Post hoc* comparisons were carried out by pairwise *t* test. Multiple comparisons were Bonferroni corrected. Nonetheless, in several cases, the assumption of homogeneity of variances has been violated. Accordingly, in these cases, the statistical tests used were Welch ANOVA with a Games–Howell test.

At criticality, an avalanche temporal profile is expected to show a characteristic shape which scales with the avalanche duration (Friedman et al., 2012): where S(t, T) is the number of events at time t in an avalanche of duration T, is the expected characteristic shape with a rescaled time axis (i.e., relative to the avalanche duration), and b is the critical exponent.

To examine whether the avalanche shapes will collapse to a single characteristic shape, , the avalanches were divided into groups of same duration and the time axis was scaled by that duration, to obtain a common axis between 0 and 1. Following, a common binning at the highest temporal resolution was used, and each avalanche was interpolated accordingly. Avalanches of the same duration were averaged to obtain a set of relatively smoother temporal profiles (i.e., shapes, one for each duration) of which amplitudes are predicted to scale by a scaling function of the form χ(*T*)*T̃ ^{b}*. The optimal χ̂(

*T*) was found by minimizing the squared difference between the scaled shapes and the mean shape relative to the mean shape. The resulting χ̂(

*T*) was plotted as a function of duration, T, in log-log coordinates, and a linear fit was used to extract the exponent of the anticipated power law. In our datasets, we found that, by accumulating all avalanches collected from all epilepsy patients, the maximal duration that has at least 1000 avalanche samples is 14 bins in non-IEA periods and 24 bins in IEA periods. The lower limit for duration was determined by supporting a meaningful temporal profile, and it was set to 5 bins. Thus, the analysis could only be applied to this relatively narrow range of durations.

## Results

Using the framework of neuronal avalanches, we investigated interictal brain activity in patients with refractory epilepsy and compared the findings with brain activity of healthy control subjects. The neuronal avalanche analysis constitutes the detection of large positive and negative signal deflections (e.g., Fig. 1*A*). Accordingly, by applying an amplitude thresholding operation (threshold of ±3 SD), the discretization of the continuous electromagnetic signals into highly synchronized point events was obtained. Subsequently, the discrete events were clustered into spatiotemporal cascades based on temporal proximity within a time bin, Δt, of duration 2.95 ms. Importantly, the findings for control subjects (*n* = 18) supported these choices of threshold and Δt: the branching parameter, σ, was ∼1 (σ = 1.06 ± 0.15), the distribution of cascade sizes followed a power law with an exponent, α, that was close to −3/2 (α = −1.52 ± 0.08), and the measure, κ, which quantifies the differences between an experimental cascade size CDF and the expected theoretical reference, was ∼1 (κ = 0.97 ± 0.02) (please also see Table 2, where the values of metrics are indicated for adults and children separately). Indeed, a maximum likelihood based analysis demonstrated a significantly higher likelihood of a power law compared with an exponential function for all subjects (see Materials and Methods, *p* < 0.01), as well as a significantly higher likelihood of an exponentially truncated power law compared with a log-normal function in 17 subjects (*p* < 0.05; 1 patient demonstrated *p* > 0.2). At the group level, according to BMS, the expectation of the posterior gives highest likelihood for the power law model according to the AIC log-evidence approximation, and for the exponentially truncated power law according to the BIC approximation. By AIC: power law, 0.79; exponential, 0.06; log-normal, 0.06; exponentially truncated power law, 0.09. By BIC: power law, 0.12; exponential, 0.05; log-normal, 0.08; exponentially truncated power law, 0.75. As expected for healthy subjects, all the values indicate close proximity to a critical branching process, supporting the parameters' choice. Nonetheless, sensitivity to variations in the time scale Δt and in the threshold were small and closely followed previously reported results (Shriki et al., 2013).

Figure 1 illustrates several consecutive neuronal avalanches associated with interictal epileptiform discharges in a particular patient. The high positive and negative amplitude deflections that can be seen in Figure 1*A* were converted into suprathreshold discrete events that are presented on a raster plot (Fig. 1*B*). The time axis (horizontal, Fig. 1*B*) was binned according to Δt (gray vertical lines, Fig. 1*B*), and each uninterrupted consecutive group of events represents a neuronal avalanche.

Figure 1*C*, *D* illustrates the spatiotemporal nature of the neuronal avalanches. They present topographic maps across time, each portraying the inward and outward magnetic fields (i.e., red and blue extrapolations across the 2D layout of MEG sensors) that surround the intracellular current flows of a synchronized and aligned group of cortical neurons (forming a dipole). In this particular case, the measured activity is IEA. Noticeably, the propagation of these current flows in the opposite direction or across different orientations of the folded cortex may change the polarity and spatiotemporal distribution of the measured fields (Fig. 1*D*). However, because our interest lies in the identity and magnitude of activated networks, the extreme of each suprathreshold deflection, whether positive or negative, was marked as a discrete event (symbolized as black dots on the topographic maps), and the combination of such discrete events, in-between silent time bins, as an avalanche (the topographic maps of each avalanche time bins are grouped with a purple curly bracket; Fig. 1*C*,*D*).

### Comparing patients with epilepsy and healthy control subjects

We compared the values of the major neuronal avalanche metrics from the full recordings of epilepsy patients with those of healthy control subjects. The dominant tendency in the sampled epilepsy population was toward substantially higher values than the expected theoretical values for a critical process (Harris, 1989) and, accordingly, higher than the values for healthy control subjects, suggesting an altered EIB in the epilepsy patients [for patients with epilepsy (*n* = 18), α = −1.37 ± 0.12, κ = 1.03 ± 0.05, σ = 1.29 ± 0.23; vs healthy control subjects, *p* ≈ 0.0001 for α and κ, and *p* ≈ 0.001 for σ] (Fig. 2*B*, a 3D phase plot; Table 2, mean ± SD of adults and children). The avalanche size distributions of one patient and one control subject are shown in Figure 2*A*. As demonstrated in Figure 2*A* for a single patient, in some of the patients with epilepsy, there were higher probabilities for large-size avalanches than expected from a −1.5 power law (dashed black line), and accordingly the slope of the avalanche size distributions on log-log coordinates was substantially shallower. Nonetheless, the observed behavior of the avalanche size distributions of all patients was noticeably scale-free, with a cutoff (i.e., a reduction in the frequency of avalanches larger than a particular size compared with a power law behavior) that corresponded to the size of the sensor array. The mean estimated cutoff increases with the size of the (sub-)sensor array (epilepsy patients: *R*^{2} = 0.94, *p* < 0.05, healthy control subjects: *R*^{2} = 0.91, *p* < 0.05), whereas the cutoff occurred at significantly larger size avalanches in patients with epilepsy under the same spatial constraints; whole array: *p* = 0.02, Bonferroni corrected. Additionally, the mean exponent of the exponential modulation, λ, in the exponentially truncated power law model was significantly smaller for epilepsy patients, λ = 0.021 ± 0.012, than for healthy control subjects, λ = 0.039 ± 0.013 (*p* < 10^{−4}). A maximum likelihood based analysis demonstrated a significantly higher likelihood of a power law compared with an exponential function for all subjects (see Materials and Methods, *p* < 10^{−7}), as well as a significantly higher likelihood of an exponentially truncated power law compared with a log-normal function in 7 patients (*p* < 0.05; 3 patients demonstrated a significance tendency *p* < 0.1, 8 patients demonstrated *p* > 0.2). At the group level, according to BMS, the expectation of the posterior gives the highest likelihood for the power law model. By AIC: power law, 0.83; exponential, 0.06; log-normal, 0.06; exponentially truncated power law, 0.05. By BIC: power law, 0.68; exponential, 0.05; log-normal, 0.20; exponentially truncated power law, 0.07. Overall, it seems that, although some patients occupied a higher-values region of the (α, κ, σ) phase space, there was also a substantial overlap with the upper end of the healthy controls' distribution, and therefore the two populations are not completely separable on the basis of these parameters.

The temporal organization of neuronal avalanches has additional special characteristics, which are predicted for a critical branching process (Harris, 1989). Accordingly, the avalanche duration, T, follows a power law P(T) ∼ T^{β} with an exponent β = −2 (Plenz, 2012). A maximum likelihood estimation was applied directly to the sample of avalanche durations to estimate the power law exponent. Whereas for healthy control subjects β ≈ −2 [(*n* = 18), β = −2.02 ± 0.10], for patients with epilepsy the exponent is higher than expected for critical systems [(*n* = 18), β = −1.77 ± 0.14; (*p* < 10^{−9})]. Nonetheless, a significantly higher likelihood of avalanche duration distributions to be modeled by an exponential function compared with a power law was found for 13 patients (*p* < 0.001; 2 patients demonstrated a significance tendency *p* < 0.1, whereas 3 patients demonstrated *p* > 0.1), as well as for 17 control subjects (*p* < 10^{−7}). Additionally, a significantly higher likelihood of a log-normal function compared with an exponentially truncated power law was found for 9 patients (*p* < 0.05), whereas only 2 patients demonstrated a higher likelihood of an exponentially truncated power law compared with a log-normal function (*p* < 0.05; 7 patients demonstrated *p* > 0.1). In the control group, 4 subjects demonstrated a higher likelihood of a log-normal function, whereas 4 subjects demonstrated a higher likelihood of an exponentially truncated power law (*p* < 0.05; 10 subjects showed no preference to either function *p* > 0.1). At the group level, according to BMS, for epilepsy patients the expectation of the posterior gives the highest likelihood for the exponential and log-normal models according to the AIC log-evidence approximation, and of the log-normal according to the BIC approximation. By AIC: power law, 0.07; exponential, 0.50; log-normal, 0.37; exponentially truncated power law, 0.06. By BIC: power law, 0.05; exponential, 0.08; log-normal, 0.80; exponentially truncated power law, 0.07. For the control subjects, the expectation of the posterior gives highest likelihood of the exponential model according to the AIC log-evidence approximation, and of the log-normal and exponential according to the BIC approximation. By AIC: power law, 0.05; exponential, 0.79; log-normal, 0.10; exponentially truncated power law, 0.06. By BIC: power law, 0.05; exponential, 0.42; log-normal, 0.44; exponentially truncated power law, 0.09. Overall, in contrast to theory, neither the power law model nor the exponentially truncated power law seems to model the datasets of healthy control subjects, as well as epilepsy patients. However, the duration distributions, as opposed to the size distributions, have a particularly narrow range of effective power law regimen due to the smaller magnitudes involved. Moreover, as reported previously (Arviv et al., 2015), there is a consistent, across all subjects, underrepresentation of small durations than would have been expected from a power law (see Fig. 4*A*). This underrepresentation might be related to sensor overlap (Shriki et al., 2013). Additionally, the theoretical prediction of power law behavior for the duration distribution assumes a fixed time-step for each step in a cascade. The natural temporal jitter in the propagation of neuronal activity may lead to deviations from the expected power law. Therefore, the fitted exponent, β, should be considered as the linear approximation to the distribution in log-log coordinates. In this sense, these fitted exponents are meaningful, and as shown above, they exhibit a significant difference between epilepsy patients and healthy control subjects.

### Comparing IEA and non-IEA periods of interictal activity

In the interictal brain dynamics of patients with epilepsy, some of the neuronal avalanches originate from IEA periods, as in Figure 1, whereas others originate from nonepileptiform brain dynamics, which we termed non-IEA. The cascades collected from IEA versus non-IEA periods, which exhibit dissimilar dynamics and morphology of measured waveforms, may be characterized by a difference in their distances from the critical point. Therefore, we examined whether the characteristics of IEA and non-IEA in the epilepsy patients differed from the near-critical dynamics of the healthy control subjects. In contrast to the significant difference in the suprathreshold event rate between epilepsy patients and healthy control subjects (*p* = 0.002; Fig. 3*A*, inset), the non-IEA periods of patients with epilepsy have a similar event rate to that of healthy control subjects (*p* = 0.66; Fig. 3*A*). Indeed, the analysis shows that the event rate during IEA periods was significantly higher than the rate during non-IEA periods and the rate in healthy control subjects (*p* < 0.0005, for both). Because Levene's test resulted in the rejection of the null hypothesis of equal variances, the reported *p* values are those of a Welch ANOVA with a Games–Howell test. Additionally, the mean duration (and size) of the avalanches during non-IEA in the patients with epilepsy were significantly (and with a significance tendency) greater than those in healthy control subjects (*p* < 0.001 and *p* = 0.07, respectively), whereas the avalanches during the IEA periods were even larger (*p* < 0.001 and *p* = 0.01 for duration and size, respectively, compared with those in the non-IEA periods and to those in healthy control subjects). Again, because Levene's test resulted in the rejection of the null hypothesis of equal variances, the reported *p* values are those of a Welch ANOVA with a Games–Howell test. In the phase plot presented in Figure 3*D*, there was a considerable overlap between patients with epilepsy during non-IEA periods and healthy control subjects (which, indeed, did not include only the lower end region occupied by healthy control subjects). The differences in metrics of non-IEA periods compared with healthy control subjects were significant for α (*p* = 0.05) and with a significance tendency (*p* = 0.06) for κ; differences for σ were not significant (*p* = 0.45). Importantly, the mean values across epilepsy patients for non-IEA periods (*n* = 17; α = −1.47 ± 0.06, κ = 0.99 ± 0.02, σ = 1.10 ± 0.13), are similar to those for healthy control subjects (*n* = 17; α = −1.52 ± 0.08, κ = 0.97 ± 0.02, σ = 1.06 ± 0.16) and are likewise consistent with the predicted values for a critical branching process. However, the IEA periods were characterized by a substantial deviation from the critical point toward supercritical dynamics (*n* = 17; α = −1.22 ± 0.10, κ = 1.09 ± 0.05, σ = 1.65 ± 0.24) (Table 2 summarizes metrics' values for IEA and non-IEA periods, and healthy control subjects). The values of the avalanche metrics during IEA periods had only a minor overlap in this phase space with non-IEA periods and with the avalanche metrics in healthy control subjects (*p* < 10^{−9}). Remarkably, in each patient with epilepsy, the parameters for the IEA periods were larger than those for the non-IEA periods. The avalanche size distributions of any particular patient demonstrated the above-described divergence of the IEA and non-IEA distributions while still exhibiting scale-free behavior with a cutoff at approximately the size of the sensor array (Fig. 3*C*, a single patient). The mean estimated cutoff across patients increased with the size of the sensor array (IEA, *R*^{2} = 0.97, *p* < 0.01; non-IEA, *R*^{2} = 0.94, *p* < 0.05; Fig. 4*A*), while the cutoff occurred at significantly larger size avalanches during IEA periods under the same spatial constraints; whole array: *p* < 10^{−6} (*p* < 10^{−6}), subsamples of array (Fig. 4*A*, insets): *p* < 10^{−7} (*p* < 10^{−5}), *p* < 10^{−5} (*p* < 10^{−4}), *p* < 0.001 (*p* < 0.05), *p* = 0.001 (*p* < 0.05), Bonferroni corrected, compared with non-IEA periods (and healthy control subjects). Additionally, the mean exponent of the exponential modulation, λ, in the exponentially truncated power law model was significantly smaller for IEA periods, λ = 0.013 ± 0.007, than for non-IEA periods, λ = 0.034 ± 0.010 (*p* < 10^{−6}), healthy control subjects, λ = 0.039 ± 0.013 (*p* < 10^{−7}), while no significant difference was found between non-IEA periods and healthy control subjects (*p* = 0.23). A maximum likelihood based analysis demonstrated a significantly higher likelihood of a power law compared with an exponential function for both IEA and non-IEA periods (see Materials and Methods, *p* < 10^{−7}). For IEA periods, there was a significantly higher likelihood of an exponentially truncated power law compared with a log-normal function in 13 patients (*p* < 0.05; 1 patient demonstrated a significance tendency *p* < 0.1, 3 patients demonstrated *p* > 0.2). For non-IEA periods, a significantly higher likelihood of an exponentially truncated power law compared with a log-normal function was found in 11 patients (*p* < 0.05; 1 patients demonstrated a significance tendency *p* < 0.1, 5 patients demonstrated *p* > 0.1). At the group level, according to BMS, for IEA periods the expectation of the posterior gives the highest likelihood for the power law model. By AIC: power law, 0.83; exponential, 0.05; log-normal, 0.06; exponentially truncated power law, 0.06. By BIC: power law, 0.77; exponential, 0.05; log-normal, 0.11; exponentially truncated power law, 0.07. For non-IEA periods, the expectation of the posterior gives the highest likelihood for the power law model. By AIC: power law, 0.80; exponential, 0.06; log-normal, 0.07; exponentially truncated power law, 0.07. By BIC: power law, 0.48; exponential, 0.06; log-normal, 0.24; exponentially truncated power law, 0.22.

Next, a maximum likelihood estimation was applied to avalanche durations of both IEA and non-IEA periods of each patient to determine the power law exponent. The mean power law exponent fitted to IEA periods (*n* = 17) is β = −1.60 ± 0.13, which is significantly shallower than the value for non-IEA periods, β = −1.90 ± 0.08, and the value for healthy control subjects (*n* = 17) β = −2.02 ± 0.10 (*p* < 10^{−9} and *p* < 10^{−11}, respectively). Moreover, non-IEA periods are significantly different from healthy control subjects (*p* < 10^{−6}), although both demonstrate proximity to β = −2, which is consistent with the predicted value for a critical branching process. Nonetheless, for IEA periods, a significantly higher likelihood of an exponential function compared with a power law was found for 14 patients, whereas a higher likelihood of a power law was found in 1 patient (*p* < 0.05; 2 patients demonstrated *p* > 0.4). Additionally, two patients demonstrated a significantly higher likelihood of a log-normal function compared with an exponentially truncated power law (*p* < 0.005, and 1 patient demonstrated a higher likelihood of an exponentially truncated power law (*p* < 0.05; 1 patient demonstrated a significance tendency *p* < 0.1, 13 patients demonstrated *p* > 0.2). For non-IEA periods, a significantly higher likelihood of an exponential function compared with a power law was found for 16 patients (*p* < 5 × 10^{−4}, 1 patient demonstrated *p* > 0.1). Additionally, in 4 patients, a significantly higher likelihood of a log-normal function compared with an exponentially truncated power law was found for (*p* < 0.05), and 1 patient demonstrated a higher likelihood of an exponentially truncated power law compared with a log-normal function (*p* < 0.05; 12 patients demonstrated *p* > 0.1). At the group level, according to BMS, for IEA periods the expectation of the posterior gives the highest likelihood for the exponential model according to the AIC log-evidence approximation, and for the log-normal according to the BIC approximation. By AIC: power law, 0.08; exponential, 0.73; log-normal, 0.12; exponentially truncated power law, 0.07. By BIC: power law, 0.06; exponential, 0.12; log-normal, 0.73; exponentially truncated power law, 0.09; and also for non-IEA periods. By AIC: power law, 0.06; exponential, 0.78; log-normal, 0.09; exponentially truncated power law, 0.07. By BIC: power law, 0.06; exponential, 0.21; log-normal, 0.60; exponentially truncated power law, 0.13. Overall, neither the power law model nor the exponentially truncated power law seems to model the datasets of patients at IEA and non-IEA periods. The model inversion of duration distributions suffers from the shortcoming described in the previous section. However, the fitted exponents, β, exhibit a significant difference between IEA and non-IEA periods.

Accumulating all avalanches from all patients of IEA periods and non-IEA periods separately, enables a “grand” (across all patients) perspective. Intermixing data from different patients into a common pool has the benefit of increasing statistical strength, although it mounts up variability among patients, particularly at IEA periods. Figure 4*A* demonstrates “grand” avalanche duration distributions of IEA and non-IEA periods. The power law exponent of the IEA curve with the highest likelihood is β = −1.56 and of the non-IEA curve is β = −1.92. Clearly, the curve of IEA periods is shallower and demonstrates a tendency of all patients toward avalanches of longer durations during IEA periods (Fig. 3*B*). Noticeably, in the “grand” cascade duration distributions, as in individual patients' distributions, there is an underrepresentation of cascades of short durations, which may result from sensor overlap.

So far, we have focused on power law distributions, which reflect scale-free dynamics. Although the observed exponents are consistent with those predicted for a critical system, in general, power law behavior may also arise from noncritical dynamics. To go beyond power law distributions, we also looked for universal scaling relationships in the data. Specifically, we applied the avalanche shape collapse analysis (Friedman et al., 2012) (see Materials and Methods), which is an additional way to probe for critical dynamics. This analysis relies on having a large amount of avalanches of several durations, and thus benefits from the “grand” perspective. We note that our data-poll supported these type of analyses only in the grand (across all patients) perspective and could span a limited duration range (5 × Δt up to 14 × Δt in non-IEA, and up to 24 × Δt in IEA; see Materials and Methods). An avalanche shape, *S(t, T)*, is the temporal profile of the number of events in each time bin, *t*, along the avalanche duration, *T*. For a critical system, avalanche shapes are predicted to have a universal shape, which scales as a function of avalanche duration, with a scaling exponent, b (Eq. 11). Avalanche shapes before and after collapse are depicted in Figure 4*B* (non-IEA: bottom left, IEA: bottom right). The dependence of the scaling factor on the avalanche duration is shown in Figure 4*B*, displaying a power law with an exponent of 0.33 for both non-IEA (top left panel) and non-IEA (top right panel). Notably, although the same scaling exponent was found for non-IEA and IEA periods, the avalanche shapes themselves have higher amplitudes at IEA versus non-IEA (Fig. 4*B*, bottom). The scaling of avalanche shapes further strengthens the evidence for critical dynamics. Nonetheless, the scaling exponents derived from the shape collapse analysis are predicted to obey the following relationship (Sethna et al., 2001; Friedman et al., 2012): = *b* + 1. In our data, the right side is ∼1.3 for both non-IEA and IEA periods, whereas the left side is ∼1.9 and 2.7, respectively. Therefore, not only that the two sides of the equations differ, the obtained scaling of avalanche shape collapse does not disclose the deviations of IEA periods from the critical point (α = −1.5, β = −2.0) and from the non-IEA periods. Indeed, the above analyses have shortcomings: (1) intermixing data of several subjects, which may differ in temporal organization of neuronal avalanches and distance from critical dynamics, may obscure our results; (2) the model inversion of the duration distributions suffers from an underrepresentation at short durations and from a narrow effective range in which a power law regimen can exist; and (3) the shape collapse analysis was applied to only a relatively small range of durations and the obtained shapes were not smooth, suggesting sensitivity to the relatively small amount of avalanches used in averaging over these temporal profiles. Nonetheless, the scaling exponent reported here is in agreement with previous results (Arviv et al., 2015). Although noisy at longer durations, the avalanche shape does not seem to be parabolic (Friedman et al., 2012). Rather, it has a piecewise behavior with a flat part between initiation and termination. As reported previously, the initiation and the termination of avalanches are correlated and reflect the gain of the system, whereas the middle part may reflect different dynamical processes (Arviv et al., 2015). The new findings presented here suggest that deviations from criticality, which are reflected in the gain of the system, may not necessarily be reflected in the scaling relationships of avalanche shapes. These findings call for further investigation.

### Deviations from critical dynamics are more pronounced during short time intervals around discharges

Overall, it seems that neuronal-avalanche-based metrics distinguish IEA periods from non-IEA periods, both at the interpatient level and, even more so, without exception at the intrapatient level. Indeed, by heuristically counting the number of IEDs as extrema within IEA periods (see Materials and Methods), we found a significant correlation between the number of interictal epileptiform discharges across 2, 5, and 7 min and α, κ, and σ (e.g., the correlations at 7 min segmentation were as follows: α: <R> = 0.83 ± 0.07 [range: 0.71 to 0.93], *p* < 10^{−3} for all subjects; κ: <R> = 0.75 ± 0.08 [range: 0.61 to 0.91], *p* < 10^{−7} for all subjects; σ: <R> = 0.84 ± 0.09 [range: 0.64 to 0.95], *p* < 0.01 for all subjects, Bonferroni corrected). Therefore, specifically focusing on IED transients is important for examining the degree of deviation from criticality at the short time scale of discharges within the IEA. For this purpose, the recordings of 5 subjects with substantial amount of epileptiform activity were examined at varying interval lengths surrounding discharges. These intervals were defined as new IEAs, ranging from ±200 ms (Fig. 5*B*, brown line) to ±10 ms (Fig. 5*B*, green line) around each identified discharge. Accordingly, the new non-IEAs included all cascades outside the corresponding intervals around the IEDs (Fig. 5*B*, dashed line of same color as the color associated with the corresponding new IEA of each interval length). The cascade size distributions from a representative subject are shown in Figure 5*B* [the orange lines denote analysis according to the definition of the expert (M.M.) of IEA]. At increasingly narrower time intervals around the discharge, the deviation from a scale-free behavior became more pronounced, resembling a bimodal distribution. Thus, deviation from scale-free behavior was concentrated around very short IED transients. Avalanches of a specific size range became more frequent, substantially surpassing the expected probability from a power law distribution. Moreover, when shifting the ±10 ms interval relative to the IED at ±25 ms steps, such that the IEDs will not be at the center, only at intervals proximate and after the IEDs (∼25 ms), the deviation from scale-free behavior remained pronounced, reflecting temporal asymmetry with respect to the IED. The correlation between the probability of large-size avalanches and the temporal proximity to IEDs is also presented in Figure 5*C*. Grouping avalanches by their temporal separation from the nearest IED shows that, at close proximity to an IED, there is a substantial peak in the mean size of avalanches, which remains above baseline for tens of milliseconds (Fig. 5*C*, right). Noticeably, the baseline of mean avalanche size over IED periods is by itself larger than over non-IEA periods. This may be affected by IED clustering, although the analysis was limited to relatively short time intervals, of plus or minus the order of magnitude of the time difference between IED events. This analysis also reveals temporal asymmetry, as the largest mean avalanche size is obtained at a few milliseconds shift in time from the IED. This temporal asymmetry is partially due to the definition of the IED as the peak at the first participating sensor. Defining the IED at the last participating sensor, or even more so at the last suprathreshold event, decreased the asymmetry but did not eliminate it. Temporal asymmetry may suggest that the deviation from criticality is more prolonged after the IED than before and that it reflects some asymmetry in the underlying dynamical processes (fast ignition and slow decay). However, because of the large number of IED events required, these analyses could only be applied to very few subjects. Thus, further study is needed to draw general conclusions about asymmetry from these data and there may be individual differences.

### Neuronal avalanches during IEA periods are more stereotypical

The stereotypical nature of epileptiform activity may be revealed throughout the collected avalanches, not only in terms of their size range but also in terms of the spatial identity of the participating sensors, which may potentially reflect the activation of epileptic networks. One way to quantify this is by examining whether avalanches tend to involve the same groups of sensors or different combinations of sensors. Given the *n* most active sensors, characterized by the highest event rate, we calculated the percentage of avalanches that involved these sensors among all avalanches of size *n* or more. We repeated this analysis for different values of *n* and applied it separately to IEA and non-IEA periods. Figure 6*A* shows that, during IEA periods, there is a much higher prevalence of avalanches involving the most active sensors (*p* < 10^{−8}, for 1, *p* < 10^{−7}, for 2, *p* < 10^{−6}, for 3 and 5, *p* < 10^{−5} for 10, *p* < 10^{−4} for 20, *p* = 0.002 for 30 and with a significance tendency *p* = 0.08 for 50 sensors in combination, Bonferroni corrected). This higher propensity of coactivation of leading sensors (within avalanches) suggests a more spatially stereotypical activity during IEA. Moreover, by focusing on two-sensor combinations within an avalanche, and calculating for each pair of sensors the probability to co-participate in an avalanche during IEA periods versus the probability during non-IEA periods, it became apparent that the corresponding probability matrices displayed different patterns (Fig. 6*C*, left panels, an example from a single subject). Of particular interest are pairs of sensors that have high probability to coparticipate in avalanches. Therefore, we converted the probability matrices into binary matrices, assigning 1 to pairs with a probability that is greater than the mean + 1 SD, and zero otherwise. Then, based on these matrices, we calculated the degree for each sensor (i.e., the number of connections to that sensor) (Rubinov and Sporns, 2010). The mean degree per sensor across patients was significantly higher during IEA periods than during non-IEA periods (IEA, 36.1 ± 7.5; non-IEA, 27.7 ± 2.8; higher at each patient; controls, 28.7 ± 3.3). Accordingly, the degree density per patient (i.e., the percentage of existing connections from possible connections) was significantly higher during IEA periods than during non-IEA periods and significantly higher than in healthy control subjects (IEA, 14.6 ± 3.1%; non-IEA, 11.2 ± 1.1%, *p* < 0.001; higher at each patient; controls, 11.7 ± 1.3%, *p* < 0.005; non-IEA periods do not differ significantly from values for healthy control subjects, *p* = 0.6) (Fig. 5*B*, bottom panel) (Rubinov and Sporns, 2010). Because Levene's test resulted in the rejection of the null hypothesis of equal variances, the reported *p* values are those of a Welch ANOVA with a Games–Howell test. Furthermore, the dispersion of the degree distribution across sensors, measured as the SD of the distribution of each patient, was higher for IEA periods than for non-IEA periods across patients (IEA, 29.1 ± 13.5; non-IEA, 12.0 ± 6.3, *p* < 0.001; higher at each patient; controls: 13.5 ± 5.4, *p* < 0.001; non-IEA do not differ significantly from control, *p* = 0.7) (Fig. 6*B*, top). Testing the homogeneity of variances by Levene's test resulted in the rejection of the null hypothesis of equal variances for 14 of 19 patients (*p* < 0.05). Figure 6*C* (middle panels) shows the degree distribution from a single patient (for clarity, each distribution was divided by its mean). Clearly, the degree distribution is more spread out during IEA and more centered near the mean during non-IEA periods. Plotting the degree associated with each sensor as a topographic map (Fig. 6*C*, right panels) demonstrated the realization of relatively extreme values of both high and low degrees and revealed a relation between sensors of high degree and epileptic networks. As apparent from Figure 6*C* (top right), this patient has a right dorsolateral frontal epileptic focus. In contrast, the topographic maps of non-IEA periods revealed a narrower range of degree. This patient, as most control subjects and other patients during non-IEA periods, exhibited a relatively higher degree at medial frontal and medial posterior sites (Fig. 6*C*).

## Discussion

Accumulating evidence supports the hypothesis that healthy brain dynamics maintain proximity to a critical state (for review, see Massobrio et al., 2015). To examine the sensitivity of metrics of criticality to changes in the balance of excitation and inhibition in epilepsy, we evaluated deviations from critical brain dynamics in patients with refractory epilepsy outside of seizures.

### Neuronal avalanches in epilepsy: an across-subject outlook

We found a tendency among the epilepsy patients toward higher values of the branching parameter (i.e., of the neural gain) and toward shallower avalanche size distributions compared with the healthy control subjects. Within each group, there was a large variability in the neuronal avalanche metrics, which may reflect differences in the individual distance from critical dynamics or variations in the physical constraints of the measurements, such as the size and position of the head within the MEG helmet (the latter option is the less likely, because we found no differences in the calculated neuronal avalanche metrics between children and adults; Table 2). Nonetheless, despite this variability, the differences between the two groups were significant. Antiepileptic drugs, with which all patients in this study were treated, tend to normalize the EIB and may reduce the density of interictal epileptiform elements (Duncan, 1987; Khachidze et al., 2010). Nonetheless, the patients in this study have drug-resistant epilepsy; thus, in these patients, antiepileptic drugs do not fully control seizures. The influence of antiepileptic drugs would diminish the difference between groups. The fact that this difference is still present indicates that the antiepileptic drugs only partially influenced the MEG traces. Assigning the obtained deviations to epilepsy was further supported by examining periods of IEA and non-IEA periods separately. The division between the epilepsy patients and the healthy subjects became much more pronounced when we focused on IEA periods, suggesting a shift toward an excitation-dominated state. Nevertheless, we note that the avalanche analysis and the branching process description, by themselves, cannot directly determine the neuronal populations or mechanisms responsible for the observed changes. For example, the effect of increased synaptic excitation may be equivalent to that of decreased inhibition. However, in general, supercritical dynamics is likely to correspond to an excitation-dominated state, whereas subcritical dynamics is likely to reflect an inhibition-dominated state, whatever the nature of the precise underlying mechanisms. Several theoretical physiological models relate the emergence of seizure activity to mechanisms involving changes in the underlying excitatory and/or inhibitory dynamics (da Silva et al., 2003; Breakspear et al., 2006; Jirsa et al., 2014). Furthermore, differences in network connectivity, short-term depression and facilitation, single neuron excitability, and additional factors can also determine whether the network will be subcritical, critical, or supercritical. Studying more biologically detailed models using the analyses proposed here could constrain the possible mechanisms and provide new insights into the changes that take place.

Despite the higher values of the branching parameter, the avalanche size distributions during IEA periods still displayed power law behavior, albeit with higher exponents. This finding suggests that scale-free behavior can still be maintained even when the balance is perturbed. We found substantial deviations from scale-free behavior only at very short time intervals around the interictal epileptiform discharges themselves. Presumably, during seizures, an excessive and sustained hypersynchronous discharge activity may demonstrate a more durable bimodal behavior, similar to the behavior of cortical tissues under pharmacologically induced shifts in the EIB (Shew et al., 2009). If so, this view may support a continuous outlook between the isolated abnormal waveforms and seizures (Fisher et al., 2014).

The neuronal avalanche metrics during non-IEA periods in the epilepsy patients were close to those expected from a critical branching process (α = −1.5, κ = 1, σ = 1). However, they demonstrated a small but significant positive difference from the brain dynamics of the healthy subjects in terms of a longer mean duration of neuronal avalanches and a higher mean power law exponent. In contrast, there were no significant differences in the event rate or in the branching parameter. Notably, several studies demonstrated that antiepileptic drugs influence the spectrum of background activity (Neufeld et al., 1999; Clemens et al., 2007; Cho et al., 2012), yet such spectral differences may specifically relate to good (and not to poor) drug responsiveness (Kim et al., 2015). Accordingly, the differences in neuronal avalanche metrics between control subjects and non-IEA periods in patients were minimal.

The small differences between non-IEA periods and the brain dynamics of the healthy subjects may not necessarily reflect a divergence from critical dynamics in patients with epilepsy when outside of IEA. A possible source of the observed differences may be a suboptimal division of our datasets between IEA and non-IEA periods. An additional factor that affects the neuronal avalanche analysis is prolonged wakefulness (Meisel et al., 2013). The epilepsy patients in our study had been instructed not to sleep during the day so as to increase the yield of IEDs during the clinic visit in the evening (Pillai and Sperling, 2006). Thus, although these patients were not sleep deprived, there may have been up to a few hour dissimilarity in sustained wakefulness before MEG recordings between patients and control subjects, which may account for the small difference between the neuronal avalanche distributions. Sustained wakefulness was shown to affect neuronal avalanche metrics, demonstrating, for example, an increase in the branching parameter, σ, with time awake (Meisel et al., 2013). The significance of this result was shown while comparing the mean over the first 9 h awake to the mean over 30–39 h awake. Nonetheless, the σ trace (with 3 h resolution) clearly indicates that the substantial variation is associated with approaching sustained 24 h awake. Accordingly, here, the potential difference between patients and control subjects for up to a few additional hours of wakefulness should amount to minute differences, if any. This is in agreement with our findings.

### Neuronal avalanches in epilepsy: individual differences

Overall, our findings suggest that IEA periods, as opposed to non-IEA periods, deviate from the critical dynamics associated with healthy cortical function and that variations in the neuronal avalanche metrics are correlated with epileptiform activity. Despite substantial across-subject variability in neuronal avalanche metrics, we saw a consistent, without exception, intrasubject deviation in neuronal avalanche metrics between the IEA and non-IEA periods. Accordingly, within each individual, IEA periods deviated toward the supercritical regimen in all measures. Thus, it may be possible to monitor the instantaneous distance from a critical state and deviations in EIB relative to the patient-specific scale. Adopting this framework may open up new prospects for diagnosis and treatment of the disease through personalized medicine.

Additionally, we were able to demonstrate that the collected avalanches during IEA have stereotypical spatial patterns, in accordance with the activation of specific epileptic networks of the individual. Epileptic networks give rise to propagating synchronized discharges, which seem to be successfully captured by the spatiotemporal nature of neuronal avalanches. If the links between each pair of sensors are defined as their probability to coparticipate in an avalanche, the connectivity matrices of IEA periods versus non-IEA periods dissociate between abnormal connectivity (of disperse degree distribution and with nodes of significantly higher degrees) to a relatively normal connectivity. Resting state activity reflects functional networks that are correlated with the underlying anatomical connectivity (Hagmann et al., 2008; Greicius et al., 2009) and may differ between healthy subjects and patients with epilepsy (Luo et al., 2011; Wang Z. et al., 2011). In previous work, we showed that the brain dynamics of healthy subjects maintain an individual-specific operation point in a phase-space characterized by the (α, κ, σ) metrics, as each individual revealed a high consistency across evoked and resting-state activity (Arviv et al., 2015). Nonetheless, the spatiotemporal patterns of neuronal avalanches that were manifested under such similar dynamical constraints were substantially different between the two cognitive states (Arviv et al., 2015). This implies that the spatiotemporal patterns of neuronal avalanches carry additional information, which could be harnessed in future studies of abnormal connectivity in the epileptic brain.

### Neuronal avalanches in epilepsy: integration into the state of the art

The present results suggest that, in patients with drug-resistant epilepsy, there is a deviation from the critical branching process description, particularly at temporal and spatial proximity to epileptiform discharges. These results accord with those of previous studies exploring various metrics within the criticality framework. For instance, Hobbs et al. (2010) analyzed local field potential activity recorded by microelectrode arrays from cortical tissues removed from patients with epilepsy. They demonstrated both a bump in the avalanche size distribution over relatively large-size avalanches and a correlation between the branching parameter and the firing rate that was restricted to periods of elevated firing rate (Hobbs et al., 2010). Previously, the scaling exponents of neuronal avalanches and long-range temporal correlations were shown to be highly correlated in healthy brain dynamics (Palva et al., 2013). In patients with epilepsy, intracranial EEG recordings provided evidence of robust power law scaling of long-range temporal correlations, with higher scaling exponents in the seizure onset zone, even during interictal periods; these findings suggest more persistent temporal correlations, with slower decay, in epileptogenic regions (Parish et al., 2004; Monto et al., 2007). Hence, the framework of critical dynamics offers novel viewpoints that may provide mechanistic insights into this neurological disease (Monto et al., 2007).

Another illuminating aspect is the relationship between epilepsy and sleep. Sleep, particularly deep NREM sleep, increases IEA and certain seizure types, yet, REM sleep seems to suppress seizures (Bazil, 2000; Méndez and Radtke, 2001; Kotagal and Yardi, 2008). Using neuronal avalanche analysis from intracranial depth electrodes in epilepsy patients during seizure-free nights, significantly larger and longer avalanches were detected during NREM sleep than while awake or during REM sleep (Priesemann et al., 2013). Sleep deprivation also increases the occurrence of epileptiform discharges as well as seizures (Pillai and Sperling, 2006; Kotagal and Yardi, 2008). Correspondingly, the probability for larger cascade sizes and branching parameter was shown to increase during sustained wakefulness of over 24 h awake (Meisel et al., 2013). The current study complements this important perspective by showing that IEA is related to deviation from criticality toward supercritical dynamics and that avalanches metrics are correlated with the propensity for IEDs.

Moreover, recently, the hypercortical excitability of sleep-deprived healthy subjects, particularly when selectively deprived of rapid eye movement sleep, was demonstrated by using transcranial magnetic stimulation (Huber et al., 2013; Placidi et al., 2013). Here, we were able to demonstrate that the hyperexcitability of IEA is related to deviations from critical dynamics, thus providing additional support to the association between critical dynamics and EIB (Shew et al., 2009; Poil et al., 2012). Consequently, these findings may provide means to quantitatively assess changes in EIB based on noninvasive recordings from humans, which does not rely on external stimulation. Characterizing and monitoring these changes are of high importance for diagnosis and treatment of the disease.

## Footnotes

This work was supported by Israel Science Foundation Grants 754/14 and 51/11, and the Israeli Centers for Research Excellence Program of the Planning and Budgeting Committee. O.A. was supported by the Israeli Ministry of Science, Technology, and Space. This research was performed in partial fulfillment of the PhD requirements of O.A.

The authors declare no competing financial interests.

- Correspondence should be addressed to Dr. Oshrit Arviv, Brain Research Center, Bar-Ilan University, 5290002 Ramat-Gan, Israel. Oshrit.Arviv{at}biu.ac.il