Skip to main content
Advertisement
  • Loading metrics

Assessing criticality in pre-seizure single-neuron activity of human epileptic cortex

  • Annika Hagemann,

    Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing

    Affiliation Max-Planck-Institute for Dynamics and Self-Organization, Göttingen, Germany

  • Jens Wilting,

    Roles Conceptualization, Writing – review & editing

    Affiliation Max-Planck-Institute for Dynamics and Self-Organization, Göttingen, Germany

  • Bita Samimizad,

    Roles Data curation, Writing – review & editing

    Affiliation Department of Epileptology, University of Bonn Medical Centre, Bonn, Germany

  • Florian Mormann,

    Roles Data curation, Supervision, Writing – review & editing

    Affiliation Department of Epileptology, University of Bonn Medical Centre, Bonn, Germany

  • Viola Priesemann

    Roles Conceptualization, Investigation, Methodology, Supervision, Writing – original draft, Writing – review & editing

    viola.priesemann@ds.mpg.de

    Affiliations Max-Planck-Institute for Dynamics and Self-Organization, Göttingen, Germany, Bernstein Center for Computational Neuroscience (BCCN) Göttingen, Germany

Abstract

Epileptic seizures are characterized by abnormal and excessive neural activity, where cortical network dynamics seem to become unstable. However, most of the time, during seizure-free periods, cortex of epilepsy patients shows perfectly stable dynamics. This raises the question of how recurring instability can arise in the light of this stable default state. In this work, we examine two potential scenarios of seizure generation: (i) epileptic cortical areas might generally operate closer to instability, which would make epilepsy patients generally more susceptible to seizures, or (ii) epileptic cortical areas might drift systematically towards instability before seizure onset. We analyzed single-unit spike recordings from both the epileptogenic (focal) and the nonfocal cortical hemispheres of 20 epilepsy patients. We quantified the distance to instability in the framework of criticality, using a novel estimator, which enables an unbiased inference from a small set of recorded neurons. Surprisingly, we found no evidence for either scenario: Neither did focal areas generally operate closer to instability, nor were seizures preceded by a drift towards instability. In fact, our results from both pre-seizure and seizure-free intervals suggest that despite epilepsy, human cortex operates in the stable, slightly subcritical regime, just like cortex of other healthy mammalians.

Author summary

In epilepsy patients, the brain regularly fails to control its activity, resulting in epileptic seizures. So far, it is not fully understood why the brains of epilepsy patients are susceptible to seizures and what the mechanism behind seizure generation is. We investigated epilepsy from the perspective of collective neural dynamics in the brain. It has been hypothesized that epileptic seizures might be a tipping over from stable, so-called subcritical, dynamics (which are commonly found in healthy brains) to unstable, so-called supercritical dynamics. We therefore examined two potential scenarios of seizure generation: (i) epileptic brain areas might generally operate closer to instability, which would make epilepsy patients generally susceptible to seizures, or (ii) epileptic brain areas might slowly drift towards instability before seizure onset. To test these two hypotheses, we analyzed activity of single neurons recorded with micro-electrodes in epilepsy patients. Contrary to widespread expectation, we found no evidence for either scenario, thus no evidence that epilepsy involves a transition to supercritical collective neural dynamics. In fact, our results from both seizure-free and pre-seizure intervals suggest that the human epileptic brain operates in the stable regime, just like the brains of other healthy mammalians.

Introduction

The existence of epileptic seizures suggests that cortical networks self-organize to a state that is prone to instability. Interestingly, from a computational perspective, a working point at the border to instability has advantages for information processing, because it renders a network highly sensitive to input and fosters non-stereotypical responses to stimuli [14]. It is therefore hypothesized that cortical networks have to strike a balance between maximizing their sensitivity and variability, while maintaining a safety margin to instability [5, 6]. Indeed, evidence is accumulating that cortical networks in vitro as well as in healthy, awake animals operate close to a critical point which marks a transition between stable (subcritical) and unstable (supercritical) dynamics [1, 2, 623].

Epileptic seizures have been hypothesized to reflect a transition to supercritical, unstable dynamics [2431], thus presenting a failure in keeping a safety-margin from instability. Multiple computational models suggest that entering a seizure could correspond to a change in the brain’s dynamical state, via a bifurcation or critical transition [3133]. Furthermore, experimental evidence suggests that neural activity during epileptic seizures and epileptiform activity deviates from healthy activity in indicators of criticality [26, 27], and neural activity during seizure termination has shown signatures of a critical transition across different recording levels [34]. Recently, long-term changes in signatures of critical slowing down, together with epileptiform discharges, were shown to be a reliable indicator for seizure risk [35]. In contrast to these findings, recent studies on human iEEG found no consistent warning signals of a critical transition prior to seizures [36, 37], raising the question of whether, and under which conditions, the framework of criticality can capture the mechanisms that underlie seizure generation.

Assuming that seizure onset does indeed correspond to loosing the safety margin to supercriticality, it remains open whether the safety margin is in general smaller in certain areas of epileptic brains, or whether the safety margin is lost prior to seizure onset [32]. In the first scenario, one expects brain areas, in which seizures emerge, to generally operate closer to instability than the same areas in a healthy human brain, so that small fluctuations can push the system into a seizure (Fig 1b). In the second scenario, one expects that the distance to criticality diminishes systematically before seizure onset, making the onset predictable (Fig 1c).

thumbnail
Fig 1. Comparison of the distance to criticality between the hemisphere containing the epileptic focus and the contralateral hemisphere.

a Branching process approximation of activity propagation in the brain. Both excitatory and inhibitory connections are summarized to an effective excitation. Depending on the branching parameter m, dynamics are stable (m < 1), unstable (m > 1) or critical (m = 1). b, c Hypothetical scenarios of seizure generation. The epileptic focus might, in general, operate closer to criticality and consequently enter the supercritical regime in an unpredictable manner because of small fluctuations (scenario A). Alternatively, dynamics might systematically drift towards instability before seizure onset. If that drift is sufficiently slow, seizures might be predictable (scenario B). d The branching parameter of neuronal activity in MTL across recordings and patients shows no significant difference between ipsilateral and contralateral hemisphere (mixed effect ANOVA with patientID as random effect). e Estimated in ipsilateral and contralateral hemisphere for multiple recordings of two exemplary patients. While several patients showed a consistent difference between hemispheres across recordings, this difference is not predicted by the location of the epileptic focus (orange: reference recording, red: pre-seizure recording, see S2 Fig for all patients). f Comparison of the branching parameter in ipsilateral and contralateral hemisphere for different subregions of the MTL, (hippocampus (H), amygdala (A), parahippocampal cortex (PHC) and entorhinal cortex (EC)). The hippocampus shows a just significant difference between ipsilateral and contralateral hemisphere, with a smaller m in the ipsilateral hemisphere. The other subregions show no significant difference (p < 0.05/4 required after Bonferroni correction). Box plots show median and quartiles, while whiskers extend to the rest of the distribution, except for outliers (points beyond 1.5 x interquartile range).

https://doi.org/10.1371/journal.pcbi.1008773.g001

We investigated these two hypotheses using extracellular spike recordings from patients with focal epilepsy. In particular, we addressed two questions: (i) Do the brain regions in which the seizures emerge generally operate closer to criticality? (ii) Does the distance to criticality change systematically prior to seizure onset? To that end, we assessed the distance to criticality based on the branching parameter m (Fig 1a). The branching parameter m characterizes the spreading of neural activity. If m is smaller than one, one action potential (spike) on average triggers less than one action potential in the subsequent time step, and the neural network is stable (subcritical regime); if m is larger than one, runaway activity may emerge and the system is unstable (supercritical regime), and m = 1 marks the transition between the two (critical state) [38, 39].

By characterizing activity spread, the branching parameter incorporates multiple network properties, whose alterations may be associated with epilepsy, including the excitability of neurons, the excitation-inhibition ratio, and the connectivity [40, 41]. More generally, an increasing branching parameter indicates more temporally correlated neural activity [12] and is therefore a signature of critical slowing, which is expected when a system approaches a critical transition [25].

To quantify m, we made use of a novel, unbiased estimator that only requires knowing the number of spikes sampled from a small set of neurons to return a reliable estimate of the branching parameter [39]. Importantly, the estimator is invariant under subsampling [42], thus it can infer the propagation of activity in a network even when recording only a small fraction of all neurons [12, 43].

To assess whether epileptic areas generally operate closer to instability, we compared the medial temporal lobe (MTL) containing the epileptic focus to the same area in the contralateral hemisphere, as an approximation of a healthy control. To assess whether the distance to criticality diminishes systematically before seizure onset, we monitored m during the last 10 minutes before seizure onset.

Results

We estimated the distance to criticality in human MTL based on single unit activity from n = 20 patients with focal epilepsy. A precise quantification of the distance to criticality has become possible with a novel, unbiased estimator [39]. The estimator had returned a branching parameter of m ≈ 0.98 for various different mammalian species, which corresponds to stable, slightly subcritical dynamics [12, 39]. We applied the same estimator to single unit activity in human MTL, both from the hemisphere containing the epileptic focus (ipsilateral hemisphere), and the contralateral one.

Consistent with recent studies in awake mammalians [5, 6, 12, 13, 44], we found single unit activity in human MTL to reflect dynamics close to criticality, but with a small safety margin to instability (see Fig 1). The recorded activity At was vastly consistent with a branching process, as most recordings clearly showed the exponential decay that is expected for the autocorrelation functions (83 out of 91 recordings in the contralateral hemisphere, 93 out of 105 recordings in the ipsilateral hemisphere, see S1 Fig). Across subjects and recordings, the branching parameter indicated dynamics in a slightly subcritical regime, with most recordings being close to criticality (0.9 < m < 1, see Fig 1). Hence, our results for human MTL are consistent with those in other mammalian species, suggesting that mammalian cortex in general self-organizes to a slightly subcritical regime.

We tested the hypothesis that the ipsilateral hemisphere generally operates closer to the critical point, which would make it prone to tipping over to unstable dynamics (Fig 1b). However, we found no significant difference between ipsilateral and contralateral hemisphere for across recordings and patients (mixed effect ANOVA with patientID as random effect, p ≈ 0.84). Within some of the patients, however, there were consistent differences in the distance to criticality between hemispheres (Fig 1e and S2 Fig). Out of the 8 patients for whom there were sufficient recordings for a significance test, we found in 2 patients, in 2 patients and no significant difference in the remaining 4 patients. Our result thus suggests that there can be consistent differences in the distance to criticality between hemispheres but that the location of the epileptic focus does not predict that difference consistently across patients.

To test whether the distance to criticality in the ipsilateral hemisphere is only altered in specific subregions of the MTL, we analyzed the recorded activity At separately in each of the subregions, including amygdala, entorhinal cortex, parahippocampal cortex and hippocampus. Consistent with the above results, we found that single unit activity in all subregions reflects a subcritical regime (Fig 1f). To test for differences between ipsilateral and contralateral hemisphere, we jointly took into account the effects of brain area, location of the epileptic focus and patient-ID. We found a significant interaction effect of brain area and focus (mixed effect ANOVA with patientID as random effect, p ≈ 0.021). The post-hoc comparisons of the distance to criticality within each subregion revealed a significant difference between ipsilateral and contralateral hemisphere in hippocampus (mixed effect ANOVA with patientID as random effect, p ≈ 0.012). However, in contrast to our expectation, the results indicate that the hippocampus of the ipsilateral hemisphere operated further away from criticality than the contralateral hippocampus. In all other subregions, we found no significant differences between ipsilateral and contralateral hemisphere (Fig 1f).

The autocorrelation structure of neuronal activity, and therefore the branching parameter m, can change depending on vigilance state and circadian rhythm [35, 37, 45, 46]. Furthermore, anti-epileptic drugs modulate cortical excitability and were shown to affect the branching parameter [47, 48]. To make sure that temporal differences do not overshadow differences between the ipsilateral and the contralateral hemisphere, or wrongfully give rise to apparent differences, we performed the same analyses using only paired recordings (recordings obtained from the same patients during the same time intervals). Consistent with the results shown in Fig 1d and 1f, we found a significant difference between the ipsilateral and the contralateral hemisphere in hippocampus, and no significant differences in the remaining subregions and the full MTL (S3 Fig).

Previous studies have reported changes in different characteristics of neuronal activity from seconds up to a time scale of days prior to seizure onset [30, 35, 49, 50]. To investigate whether epileptic seizures are caused by systematically loosing the safety margin to supercriticality (scenario B, Fig 1c), we estimated the branching parameter m(t) in a time-resolved manner during pre-seizure recordings (Fig 2a, see also S6 Fig for all patients). While the branching parameter showed variations over time, none of the patients consistently showed a systematic trend prior to seizure onset.

thumbnail
Fig 2. No consistent change in the distance to criticality before seizure onset.

a Time-resolved estimation of for the full MTL within the last 10 min prior to seizure onset showed variability, but no consistent behavior during the last minutes before seizure onset. Plots show example traces of multiple seizures from one patient. Each trace corresponds to one pre-seizure period, shown both for the ipsilateral and the contralateral hemisphere. Traces of all patients and seizures shown in S6 Fig. b No significant difference in the distance to criticality between the last 5 min prior to seizure onset and the previous 5 min—neither in the MTL containing the epileptic focus, nor in the contralateral one (p-values of Wilcoxon signed rank test). c Separate analysis of the individual subregions of MTL. Within none of the subregions there were significant differences between the last 5 min before seizure onset and the previous 5 min. Box plots in c show results from the hemisphere containing the epileptic focus. d Variability of the time-resolved branching parameter in the last 5 min prior to seizure onset and the previous 5 min. Variability is quantified using the median absolute deviation from the median logarithmic distance to criticality (p-values of Wilcoxon signed rank test).

https://doi.org/10.1371/journal.pcbi.1008773.g002

As a more coarse measure, we compared the branching parameter of the first half to the second half of pre-seizure recordings (Fig 2b). Again, we found no significant difference when approaching seizure onset, neither in the ipsilateral, nor in the contralateral hemisphere (Wilcoxon signed rank test on the logarithmic distance to criticality, p ≈ 0.71, and p ≈ 0.96, respectively). Finally, we applied the same analysis to the individual subregions of MTL. None of the individual subregions showed a significant difference between the first and the second half of pre-seizure recordings (Fig 2c). The results were consistent when additionally accounting for the patientID as a random effect. Thus, across patients and subregions of MTL, we found no evidence for dynamics approaching supercriticality during the last 10 minutes before seizure onset.

Instead of a systematic trend in the distance to criticality before seizure onset, it is conceivable that the variability in the distance to criticality increases. If the variability is high compared to the safety margin, dynamics could enter the supercritical regime spontaneously, driven random fluctuations. To investigate whether fluctuations in the distance to criticality increase when approaching seizure onset, we compared the variability in the time-resolved distance to criticality in the first half and the second half of pre-seizure recordings (Fig 2d). However, we found no significant change when approaching seizure onset across patients and recordings.

So far, we have used the electrographic focus for our analyses, i.e. the brain area in which seizures emerged, determined by clinicians prior to surgery. However, surgery outcomes of the patients differ (see S1 Table). We therefore ran the main analyses again, including only patients that were seizure-free after the surgery (Engel scale 1a, 10 out of 20 patients). Consistent with our previous findings, we found no significant change before seizure onset in these patients (S8 Fig). Furthermore, we found no significant differences between ipsilateral and contralateral hemispheres neither in the full MTL, nor in the individual subregions (S7 Fig). In particular, the previously observed increased distance to criticality in ipsilateral hippocampus was not present in this reduced dataset. Note, however, that the number of recordings was smaller compared to the full dataset, such that the required effect sizes for significant results were correspondingly larger.

Discussion

We started off with the hypotheses that a brain area affected by epilepsy might (i) generally operate closer to criticality (and instability), or (ii) systematically move from the stable, subcritical to the unstable, supercritical regime before seizure onset. However, we found no evidence for either hypothesis when evaluating single unit activity of human medial temporal lobe. Instead, we found that both hemispheres, the one containing the epileptic focus, and the contralateral one operate in a slightly subcritical regime. Thereby our results in human medial temporal lobe are in line with those on single unit activity from multiple other mammalian species that all show reverberating, slightly subcritical dynamics [12].

We found no evidence that the hemisphere, in which seizures emerge, generally operates closer to instability. Although there were systematic differences between ipsilateral and contralateral hemisphere in some patients, these differences were not consistent across patients. Therefore, the observed differences may only reflect patient-specific local differences in the distance to criticality and cannot be attributed to the epileptogenicity of the hemisphere based on our data. One potential explanation for the lack of systematic difference in our data is that the epileptic focus might only consist in a small, localized population of neurons that was not always captured by the recording electrodes. This hypothesis is supported by evidence from epileptic seizures indicating a core territory of recruited neurons that is surrounded by a so-called ictal penumbra [51]. On the other hand, a variety of studies suggest that epilepsy is a network disease, and pre-seizure changes have been observed in neurons far away from the epileptic focus [29, 49]. Thus, it is conceivable that the entire brain operates in a different dynamical regime, which would explain the lack of a systematic difference between the hemispheres. Investigating this hypothesis would ideally require single neuron data from healthy humans, which are, however, practically not accessible with existing recording techniques.

Interestingly, we found that the branching parameter in the hippocampus of the ipsilateral hemisphere was smaller than in the contralateral hippocampus. This contrasts our initial hypothesis of a reduced distance to criticality in the ipsilateral hemisphere. A potential explanation for the increased distance to criticality could be the destruction of neuronal tissue (e.g. mesial temporal sclerosis), which is frequently observed in MTL epilepsy [52]. Neuronal cell loss can result in a more segregated network, reducing the effective coupling between neurons and therefore the branching parameter. Another potential reason are compensatory mechanisms that reduce the coupling strength between neurons because of repeated excessive activity in that area (e.g. via homeostatic plasticity).

We found no evidence for systematic transitions from the sub- to the supercritical regime before seizure onset. As all our recordings have been obtained during seizure-free intervals (reference or pre-seizure), this finding does not contradict the experimental evidence that epileptiform activity and seizures proper might reflect a transition to the supercritical regime [26, 27]. In the literature, the reported time window of putative pre-seizure changes varies considerably from seconds to minutes and hours [30, 35, 49, 53, 54]. Even on longer time scales, recent studies suggest fluctuations in seizure susceptibility, including circadian, multidien and multimonth rhythms [55, 56]. Thus potentially, the transition to supercriticality happens on a time scale that is not detectable in the analyzed 10-minute interval. The transition could either be so rapid that our time-resolved analysis cannot capture it (faster than the 80 seconds analysis window). Or the transition could occur on a much slower time scale, such that changes during the last 10 minutes are too small to be detectable. Once the system is sufficiently close to criticality, it could spontaneously be driven into a seizure by noise. This hypothesis is supported by a recent study that found signatures of critical slowing on a timescale of hours to days before seizure onset in human iEEG [35]. A former study investigating putative pre-seizure periods of 30 minutes to 4 hours, on the other hand, did not find such signatures in human iEEG [37]. A recent study on ECoG data of human patients with focal epilepsy found clear systematic trends in signatures of critical slowing during the last 30 min before seizure onset, however these trends were not consistent across patients (increase in the signal’s autocorrelation coefficient (ACC) in 4/12 patients, decrease in 4/12 and no significant change in the ACC in 4/12 patients) [30]. High frequency activity in hippocampal slices of rats, on the other hand, showed a consistent increase in the ACC as well as other signatures of critical slowing prior to seizures [30]. Thus it remains to be investigated precisely, under which conditions (e.g. types of epilepsy, types of recordings, considered time window) the framework of criticality can serve to predict an impending epileptic seizure.

Our analyses might miss signatures of a transition to supercritical dynamics because of not capturing the recruited set of neurons, or because single unit activity presents too fine a spatial scale. Single neuron activity before and during an epileptic seizure was shown to be highly heterogeneous, with some neurons increasing their firing rates, others decreasing, and a considerable fraction not changing at all [49]. In fact, evidence from spatially extended spike recordings from epilepsy patients suggests that there is a sharp boundary between areas with increased, hypersynchronous spiking and adjacent areas with low-level, unstructured activity during focal seizures [51]. This implies that recruitment of neurons to seizures occurs only in small areas and that single unit recordings potentially do not capture the recruited neurons [51]. In line with this idea, the study that found signatures of critical slowing during seizure termination found these signatures only on the more coarse recording levels (EEG, ECoG, LFP), not in multi unit activity [34]. Therefore, the sparsely recorded single unit activity may not be the ideal type of signal to identify seizure onset dynamics.

Our estimation of the distance to criticality relied on the branching process approximation. Mathematically, the branching process represents a generic model of how activity propagates in a network, but clearly, it is quite simplistic, and does not account for all the biological complexity in the brain. In particular, we do not distinguish excitatory and inhibitory connections. Instead, we use the branching parameter m to describe the effective spreading of activity in the network. Any alteration in m may thus reflect one or a combination of mechanisms, like altered synaptic strength, excitability, or excitation-inhibition ratio. Despite this simplification, a previous study has shown that the branching process model does indeed reproduce statistical properties of networks with inhibition, if one assumes that the excitatory and inhibitory contributions can be described by an effective excitation [12].

The branching process formalism has proven powerful, because in contrast to classical approaches to estimate criticality, it is (i) invariant under subsampling, returning a reliable estimate of activity spread and stability even if only a tiny fraction of neurons is sampled [5, 39, 42, 57]; (ii) it returns a precise, quantitative estimate of the distance to criticality, which enables comparison across studies, recording conditions, brain areas, and species; and (iii) it requires comparably little data, thereby enabling our time-resolved analysis. Importantly, one finds a good match between the branching process and many experimentally observed features of cortical dynamics, like spectra, Fano factors, inter-spike interval distributions, and a clear exponential decay of the autocorrelation function [6, 12, 47, 58]. Together, these aspects clearly support the validity of the branching process approximation, making it a powerful tool to assess the stability of network dynamics.

In the literature, a variety of alternative mechanisms for seizure generation have been proposed. While we focused on criticality using the framework of a branching process, several other models describe seizure generation by a critical transition (bifurcation), but with different order and control parameters [31, 32, 35, 59]. While the specific dynamics of these models differ, they all involve a deterministic change in a control parameter that could, in principle, be detectable. In fact, it is expected that all mechanisms, in which the transition to seizures involves a second order phase transition, will show signatures of critical slowing [25]. As the branching parameter m is an indicator for critical slowing, it should increase if cortical dynamics approach a second order phase transition, independent of the specific underlying path. Thus, our analysis is not restricted to a branching process model, but measures the proximity of cortical dynamics to a second order phase transition.

On the other hand, there are models that do not describe seizures via bifurcations, but with multi-stable models, in which seizure onset is not driven by deterministic parameter changes, but by random fluctuations [32, 60]. In these models, seizures would be inherently unpredictable, as random fluctuations in the dynamics drive the system into the seizure state. Such a model could, in principle, explain our negative finding. However, there is accumulating evidence that seizures are, to a certain extent, predictable, which indicates that seizures are likely not only driven by noise [35, 61, 62].

Finally, a combination of mechanisms is possible: deterministic changes in the stability could lead to an increased seizure risk, accompanied by random fluctuations that drive the system into the seizure [63]. In this case, it is conceivable that the deterministic changes happened on a timescale longer than detectable in our 10-minute analysis, and entering the seizure was then driven by random fluctuations.

In addition to the variety of proposed seizure generating mechanisms, there is strong evidence that seizure onset, progression, and termination can follow different dynamical pathways [33, 56, 63]. Such differences can exist not only between patients, but also within individual patients. This renders finding a common mechanism or warning signal inherently difficult and may provide a further explanation of why we find no consistent effect within and across patients.

We would like to conclude with a general remark on the approach of this study. We started off with two clear hypotheses about how changes in the distance to criticality might be related to seizure dynamics (Fig 1b and 1c). With the intracranial spike recordings, we could test these hypotheses, but we found no evidence for either of them. Faced with these negative results, one can either revise the hypothesis, investigate a different set of parameters or observables, e.g. by switching from a branching process instability to a framework of metastability or oscillation, and eventually one might find a hypothesis for which the data provides significant results. However, this approach can easily lead to false-positive reports. Therefore, we decided to stay with the original hypotheses, which had been formulated in all detail in a grant proposal, and then report the negative results. We believe that this is scientifically the most rigorous path to follow, though not necessarily the most rewarding. As a consequence, we have to report negative results on these specific hypotheses. However, as pointed out above, it does not exclude that other perspectives on critical phenomena might still play a role in seizure dynamics.

In summary, we found no evidence for a transition towards supercritical dynamics in pre-seizure single neuron activity of human cortex. This finding is in line with previous studies finding no warning signals of a critical transition prior to epileptic seizures [36, 37], but stands in contrast with a recent study demonstrating the effectiveness of signatures of critical slowing for seizure prediction [35]. It remains to be investigated precisely whether single-unit recordings, despite providing the most direct information on neuronal activity, cannot resolve pre-seizure changes in criticality. Furthermore, it is conceivable that seizure activity proper is indeed supercritical, with dynamics crossing the critical thresholds only seconds before seizure onset. Alternatively, seizure generation might be a qualitatively different process that cannot be captured by a linear stability parameter m.

Materials and methods

Ethics statement

All patients had given written informed consent to participate in this study, which was approved by the Medical Institutional Review Board in Bonn.

Acquisition and pre-processing of intracranial recordings

We analyzed intracranial recordings from n = 20 patients with medically intractable focal epilepsy. The data was recorded at the Department of Epileptology at the University of Bonn Medical Center. For pre-surgical evaluation, patients were implanted with depth electrodes in different regions of the medial temporal lobe, including hippocampus (H), amygdala (A), parahippocampal cortex (PHC) and entorhinal cortex (EC). The location of microwires was verified using post-implantational CT scans co-registered to pre-implantational MRI scans. All wire bundles were confirmed to be located in the designated target regions in each patient. Each electrode contained 8 microwires, with which single-unit recordings could be performed. Recordings were performed continuously for pre-surgical monitoring for a typical duration of 7-14 days. Data was sampled at 32 kHz and filtered between 0.1 and 9000 Hz. Spike sorting was performed using the Combinato package [64], using the standard parameters proposed by the authors. Sorted units were classified manually as single units, multi-units, or artifacts, using the Combinato GUI [64]. The main classification criterion for putative units was the signal’s shape, but also the amplitude and distribution of inter-spike-intervals, which allows excluding artifacts due to e.g. supply voltage. For further analysis, we only used spikes of identified single units and excluded artifacts and multi-units. Thereby, we minimized the number of artifacts, which can potentially impact subsequent analyses (multi-units typically contain more artifacts than identified single units).

For each patient, we analyzed one 10-minute reference recording, obtained in a seizure-free interval after the surgery, as well as several pre-ictal recordings, spanning 10 minutes prior to seizure onset. Only clinical seizures were included in the analysis. Pre-ictal recordings end at seizure onset, which was determined by two board-certified EEG readers. Patients in which the epileptic focus could not clearly be assigned to one hemisphere were excluded from the analysis. 17 out of the 20 patients performed the surgery, and 10 patients were seizure-free afterwards (Engel scale 1a). S1 Table summarizes the analyzed patients and recordings. In total, we started with 20 patients and 116 recording periods (20 reference, 96 pre-seizure). After spike-sorting and excluding recordings that did not return any single unit, 107 recording periods (20 reference, 87 pre-seizure) remained. For each of the recording periods, we obtained data from multiple subregions, which included hippocampus, amygdala, parahippocampal cortex and entorhinal cortex, but did not always cover all subregions in each patient. For the number of recordings in each hemisphere and subregion, see S1 Text.

Branching process approximation

We use the branching process as a minimal model for spike propagation in the brain. The branching process is a stochastic model describing the number of active neurons At in discrete time bins of length Δt. Each active unit i at time t activates a random number Yt,i of units in the subsequent time step. In addition, there is an external input ht into the system, accounting for input from sensory modalities, from other brain areas, or spontaneous activation of individual neurons. The total number of active neurons is then given by (1)

Taking the conditional expectation value yields the autoregressive representation (2) where m = 〈Yt,i〉 is the branching parameter, and h = 〈ht〉 is the average input.

A branching process is stable for m < 1 (subcritical regime) and unstable for m > 1 (supercritical regime). The critical point (m = 1) separates the two regimes and marks the critical point of a second-order phase transition. For more details, see [6, 65].

Definition of the activity At

The activity At is defined as the number of active neurons in discrete time bins Δt. Implanted depth electrodes can, however, only record a tiny fraction of all neurons, and hence one only observes a subset of the activity At. Such spatial subsampling can lead to strong biases in inferred system properties [5, 11, 18, 42, 57]. However, for estimating the branching parameter m, we have developed a method that overcomes the systematic bias [39, 65]. In brief, it shows that the autocorrelation strength, which is central when inferring m, is biased by a factor B; however, the bias factor B can be partialled out, so that we can obtain an unbiased estimate.

In the following, we thus also denote the sampled activity by At. It is defined as the number of sampled active neurons at time t. To obtain At from recorded spike times, all spikes recorded in one brain area are pooled and binned to Δt = 4 ms time bins. The time step Δt was chosen to reflect the propagation time of spikes between neurons.

Definition of the autocorrelation

To estimate m, the autocorrelation function C(k) at time lags k has to be estimated from the recorded activity At: (3) where and denote the mean activity of the original and the delayed time series, respectively, and T the duration of the recording. This definition of the autocorrelation function C(k) is equivalent to the standard definition of the Pearson correlation coefficient , with standard deviations , , as long as is a stationary process and thereby . If the activity At is consistent with a stationary process with autoregressive representation (PAR), C(k) decays exponentially [39].

Estimation of the distance to criticality

We used our open source toolbox [43], which implements the Multistep-Regression (MR) estimator to infer the branching parameter m from spike recordings [39]. The estimator is invariant under subsampling, i.e. it yields consistent estimates for the branching parameter of the whole system even if only a small fraction of all neurons is recorded.

For a stationary branching process, it can be shown that the autocorrelation of At follows C(k) ∝ mk. The branching parameter m can therefore be estimated by fitting an exponential decay (4) to the measured autocorrelation C(k). In the last term of 4, we rewrote the autocorrelation function in terms of the intrinsic timescale τ = −Δt/log(m) [39, 58]. The additional offset D in the fit function f(k) accounts for contributions with long timescales that do not decay substantially within the recording time, and it compensates for small non-stationarities in At [58]. The factor B is the bias factor of the autocorrelation and depends on the subsampling.

In short, given the activity At, estimation of m is performed in two steps [39]:

  1. Compute the autocorrelation function C(k) for different time delays k (Eq 3).
  2. Fit an exponential decay to the autocorrelation function to obtain an estimate for the intrinsic timescale and the branching parameter .

All analyses were performed using the python toolbox of the MR estimator [43]. We used time delays k ∈ [4 ms, 1600 ms], which is on the order of several autocorrelation times of our data.

Confidence intervals of estimates for single recordings were obtained via a block bootstrap procedure: recordings were divided into segments of 20 s length and estimation was performed on random subsets of segments (see stationarymean method in [43]).

Exclusion criteria of MR estimation

For reliable estimation of the branching parameter m, the data must be consistent with a stationary PAR which implies that the autocorrelation C(k) must be consistent with an exponential decay. Furthermore, reliable estimation requires a minimum number of recorded spikes, because the variance in estimates increases with decreasing number of non-zero activity entries [39]. To ensure that a given time series fulfills the requirements for reliable estimation, we implemented two consistency checks:

  1. Number of non-zero time bins must be at least .
  2. The exponential decay must fit the autocorrelation better than a simple linear function. (Otherwise, the data cannot be considered consistent with a stationary PAR).

If a recording was not consistent with either of these requirements, it was excluded from the analysis. Note that at critiality, the intrinsic timescale is expected to diverge. The autocorrelation function of dynamics very close-to-critical on a limited number of time steps can therefore be mistaken with a linear decay. We nevertheless implemented the second exclusion criterion, as it captured erroneous fits due to non-stationarities in the recording or high signal-to-noise ratios in our dataset. In general, however, the use of the second criterion must be assessed carefully for any given dataset, to ensure that it does not introduce a bias towards subcritical dynamics.

Comparison of ipsilateral and contralateral hemisphere

For each recording, the spikes from the hemisphere containing the epileptic focus were combined and binned to a single activity time series . Equally, all recorded spikes from the contralateral hemisphere were binned to . Estimated branching parameters of ipsilateral and contralateral hemisphere were compared by performing a mixed effects ANOVA, taking into account the patientID as a random effect. As the distribution of was not consistent with a normal distribution, we performed the mixed effects ANOVA on the logarithmic distance to criticality . Specifically, we used the python interface pymer4 [66] to the R-function lmer [67], with model specification . As differences in between patients are possible, we additionally compared ipsilateral and contralateral hemisphere separately within each of the patients.

The same approach was applied individually to each subregion of the MTL. In this case, spikes were binned separately for each subregion in each hemisphere. To disentangle potential effects of patient-ID, subregion and location of the epileptic focus, we conducted a mixed effects ANOVA, where patient-ID was modeled as a random effect. As above, we used the python interface pymer4 [66] to the R-function lmer [67], with model specification . In addition, we performed post-hoc pairwise comparisons of of ipsilateral versus contralateral hemisphere in all subregions using the model specification .

To test whether significant differences are only present in patients with successful surgery, the comparisons of ipsilateral and contralateral hemisphere were conducted again considering only patients with ideal surgery outcome (Engel scale 1a). Furthermore, to make sure that temporal differences in the distance to criticality (e.g. due to circadian rhythm or different level of medication) do not overshadow an effect, the analyses were conducted again only taking into account paired recordings, i.e. recordings obtained from the same patients during the same period of time.

Time-resolved estimation before seizure onset

To analyze changes in the distance to criticality prior to seizure onset, we applied the MR estimator with a sliding window approach. The crucial parameter to choose is the window size Lw, which represents a trade-off between sufficiently short windows for high temporal resolution, and sufficiently long windows that provide enough data for a consistent estimation. Based on the average activity and the average number of non-zero activity entries in our data, we chose a window size of Lw = 80 s (see S4 and S5 Figs). We used overlapping windows with a fixed window-step Lstep = 2 s. As a more coarse measure of pre-seizure changes, we splitted the pre-seizure recordings into two parts and estimated separately for both parts. Estimates of the first and the second part were then compared using the Wilcoxon signed rank test, where parts of the same recording were considered pairs. Similar to the ipsilateral/contralateral analysis, we performed the statistical test on the logarithmic distance to criticality . As an additional sanity check, we furthermore ran a mixed effect ANOVA, accounting for patientID as a random effect (model specification for each brain area. The results were consistent across the different testing procedures.

To assess whether the variability in the distance to criticality changes before seizure onset, we determined the variability in the time-resolved estimates of m in the first and the second part of the recordings. Specifically, we computed the median absolute deviation (MAD) of the time-resolved distance to criticality from the median distance to criticality in each recording part. We then compared the MADs of first and second part using the Wilcoxon signed rank test.

Supporting information

S1 Fig. Autocorrelation functions (ACF) of single unit activity in human medial temporal lobe for two example patients.

ACFs were widely consistent with exponential decays with offset (lines show the fitted exponential f(k) = Bmk+ D, where m is the estimated branching parameter). The upper row of each patient corresponds to activity in the ipsilateral hemisphere, the lower row to the contralateral hemisphere. The different plots in each row show different recordings, both reference (gray) and pre-seizure recordings (blue).

https://doi.org/10.1371/journal.pcbi.1008773.s001

(TIFF)

S2 Fig. Patient-wise comparison of the branching parameter m between the hemisphere containing the seizure onset zone and the contralateral hemisphere for both pre-seizure recordings (red) and reference recordings (orange).

While there was no consistent difference between ipsilateral and contralateral hemispheres across patients, within some of the patients there was a consistent trend in either direction (p-values of two-sided Mann Whitney-U test).

https://doi.org/10.1371/journal.pcbi.1008773.s002

(TIFF)

S3 Fig. Comparison of ipsilateral versus contralateral hemisphere using only paired data, i.e. pairs of data that were recorded simultaneously in the ipsilateral and the contralateral hemisphere.

a The branching parameter across recordings and patients shows no significant difference between ipsilateral and contralateral hemisphere. b Comparison of the branching parameter in ipsilateral and contralateral hemisphere for different subregions of the MTL, (hippocampus (H), amygdala (A), parahippocampal cortex (PHC) and entorhinal cortex (EC)). The hippocampus shows a just significant difference between ipsilateral and contralateral hemisphere, with a smaller m in the ipsilateral hemisphere. The other subregions show no significant difference (p < 0.05/4 required after Bonferroni correction).

https://doi.org/10.1371/journal.pcbi.1008773.s003

(TIFF)

S4 Fig. Average population firing rate R = 〈Attt, and number non-zero activity entries , for our dataset of pre-seizure recordings.

Averages 〈Att are computed over all time steps of the respective recording. Across recordings, the median population firing rate is q50 = 25.8 Hz and the 25% quantile is q25 = 10.4 Hz. The median fraction of non-zero activity entries is q50 = 9.9%, i.e. on average, more than 90% of the time bins contain no spike.

https://doi.org/10.1371/journal.pcbi.1008773.s004

(TIFF)

S5 Fig. Choice of the window size for the time-resolved estimation of the branching parameter.

a Estimated as a function of the window size for simulated branching processes. Error bars are 95% confidence intervals of estimates on 500 trials. Simulation parameters were adjusted to approximately match the median values of 〈At〉 and of our data set (m = 0.95, α = 0.004, Δt = 4 ms, h = 1). b, c for different window sizes in two example recordings. Error bars are 95% confidence intervals of estimates from n = 1100 partially overlapping segments of the recording (window step size of of 400 ms). The variance of estimates increases with decreasing window size. For the results shown in Fig 2c and S6 Fig, we chose a window size of 80 seconds (black arrow), representing a compromise between sufficiently high temporal resolution and sufficiently low variance.

https://doi.org/10.1371/journal.pcbi.1008773.s005

(TIFF)

S6 Fig. Time-resolved estimation of within the last 10 min prior to seizure onset for all patients and recordings.

Each trace corresponds to one pre-seizure period, shown both in the ipsilateral and contralateral hemisphere. Seizure onset of each trace was at t = 600 s and estimated was assigned to the middle of the respective window (window size 80 s). Recordings, in which more than 5% of segments were not consistent with the requirements for MR estimation were excluded. Note that activity enters the supercritical regime in a few segments of patient 18, but otherwise remains in the subcritical regime across patients and recordings.

https://doi.org/10.1371/journal.pcbi.1008773.s006

(TIFF)

S7 Fig. Comparison of ipsilateral versus contralateral hemisphere using only recordings from the 10 patients with 1a outcome after surgery.

a The branching parameter across recordings and patients shows no significant difference between ipsilateral and contralateral hemisphere. b Comparison of the branching parameter in ipsilateral and contralateral hemisphere for different subregions of the MTL, (hippocampus (H), amygdala (A), parahippocampal cortex (PHC) and entorhinal cortex (EC)). None of the regions shows a significant difference between ipsilateral and contralateral hemisphere (p < 0.05/4 required after Bonferroni correction).

https://doi.org/10.1371/journal.pcbi.1008773.s007

(TIFF)

S8 Fig. Analysis of pre-seizure changes using only recordings from the 10 patients with 1a outcome after surgery.

a No significant difference in the distance to criticality between the last 5 min prior to seizure onset and the previous 5 min—neither in the MTL containing the epileptic focus, nor in the contralateral one (p-values of Wilcoxon signed rank test). b Separate analysis of the individual subregions of MTL in the hemisphere containing the epileptic focus. Due to the small number of rescordings, we report the individual estimates without performing statistical tests.

https://doi.org/10.1371/journal.pcbi.1008773.s008

(TIFF)

S1 Text. List of recording numbers that were excluded based on the exclusion criteria for MR estimation.

https://doi.org/10.1371/journal.pcbi.1008773.s009

(PDF)

S1 Table. Intracranial recordings from epilepsy patients that were analyzed in terms of the distance to criticality.

Recordings span both hemispheres (left and right) and different subregions of MTL, including hippocampus (H), amygdala (A), parahippocampal cortex (PHC) and entorhinal cortex (EC). Individual recordings of the same patient can span different subsets of the listed brain areas. Numbers of recordings after spike sorting, before applying exclusion criteria of MR estimation. The surgery outcomes were evaluated according to the Engel scale. No entry means that no surgery has been performed. Dataset provided by the Department of Epileptology in Bonn.

https://doi.org/10.1371/journal.pcbi.1008773.s010

(PDF)

References

  1. 1. Shew WL, Plenz D. The functional benefits of criticality in the cortex. The neuroscientist. 2013;19(1):88–100.
  2. 2. Haldeman C, Beggs JM. Critical Branching Captures Activity in Living Neural Networks and Maximizes the Number of Metastable States. Physical Review Letters. 2005;94(5):058101.
  3. 3. Kinouchi O, Copelli M. Optimal dynamical range of excitable networks at criticality. Nature physics. 2006;2(5):348.
  4. 4. Bertschinger N, Natschläger T. Real-time computation at the edge of chaos in recurrent neural networks. Neural computation. 2004;16(7):1413–1436.
  5. 5. Priesemann V, Wibral M, Valderrama M, Pröpper R, Le Van Quyen M, Geisel T, et al. Spike avalanches in vivo suggest a driven, slightly subcritical brain state. Frontiers in Systems Neuroscience. 2014;8. pmid:25009473
  6. 6. Wilting J, Dehning J, Pinheiro Neto J, Rudelt L, Wibral M, Zierenberg J, et al. Operating in a Reverberating Regime Enables Rapid Tuning of Network States to Task Requirements. Frontiers in Systems Neuroscience. 2018;12. pmid:30459567
  7. 7. Beggs JM, Plenz D. Neuronal avalanches in neocortical circuits. The Journal of Neuroscience: The Official Journal of the Society for Neuroscience. 2003;23(35):11167–11177.
  8. 8. Friedman N, Ito S, Brinkman BA, Shimono M, DeVille RL, Dahmen KA, et al. Universal critical dynamics in high resolution neuronal avalanche data. Physical review letters. 2012;108(20):208102. pmid:23003192
  9. 9. Zierenberg J, Wilting J, Priesemann V. Homeostatic Plasticity and External Input Shape Neural Network Dynamics. Physical Review X. 2018;8(3):031018.
  10. 10. Neto JP, Spitzner FP, Priesemann V. A unified picture of neuronal avalanches arises from the understanding of sampling effects. arXiv preprint arXiv:191009984. 2019;.
  11. 11. Wilting J, Priesemann V. 25 years of criticality in neuroscience—established results, open controversies, novel concepts. Current opinion in neurobiology. 2019;58:105–111.
  12. 12. Wilting J, Priesemann V. Between Perfectly Critical and Fully Irregular: A Reverberating Model Captures and Predicts Cortical Spike Propagation. Cerebral Cortex (New York, NY: 1991). 2019;29(6):2759–2770.
  13. 13. Dahmen D, Grün S, Diesmann M, Helias M. Second type of criticality in the brain uncovers rich multiple-neuron dynamics. Proceedings of the National Academy of Sciences. 2019;116(26):13051–13060.
  14. 14. Munoz MA. Colloquium: Criticality and dynamical scaling in living systems. Reviews of Modern Physics. 2018;90(3):031001.
  15. 15. Hahn G, Petermann T, Havenith MN, Yu S, Singer W, Plenz D, et al. Neuronal avalanches in spontaneous activity in vivo. Journal of neurophysiology. 2010;104(6):3312–3322. pmid:20631221
  16. 16. Cocchi L, Gollo LL, Zalesky A, Breakspear M. Criticality in the brain: A synthesis of neurobiology, models and cognition. Progress in neurobiology. 2017;158:132–152.
  17. 17. Aburn MJ, Holmes C, Roberts JA, Boonstra TW, Breakspear M. Critical fluctuations in cortical models near instability. Frontiers in physiology. 2012;3:331.
  18. 18. Ribeiro TL, Copelli M, Caixeta F, Belchior H, Chialvo DR, Nicolelis MA, et al. Spike avalanches exhibit universal dynamics across the sleep-wake cycle. PloS one. 2010;5(11):e14129. pmid:21152422
  19. 19. Ribeiro TL, Ribeiro S, Belchior H, Caixeta F, Copelli M. Undersampled critical branching processes on small-world and random networks fail to reproduce the statistics of spike avalanches. PloS one. 2014;9(4):e94992.
  20. 20. Fagerholm ED, Scott G, Shew WL, Song C, Leech R, Knöpfel T, et al. Cortical entropy, mutual information and scale-free dynamics in waking mice. Cerebral cortex. 2016;26(10):3945–3952. pmid:27384059
  21. 21. Ponce-Alvarez A, Jouary A, Privat M, Deco G, Sumbre G. Whole-brain neuronal activity displays crackling noise dynamics. Neuron. 2018;100(6):1446–1459.
  22. 22. Ma Z, Turrigiano GG, Wessel R, Hengen KB. Cortical circuit dynamics are homeostatically tuned to criticality in vivo. Neuron. 2019;104(4):655–664.
  23. 23. Marshel JH, Kim YS, Machado TA, Quirin S, Benson B, Kadmon J, et al. Cortical layer–specific critical dynamics triggering perception. Science. 2019;365(6453):eaaw5202. pmid:31320556
  24. 24. Hsu D, Chen W, Hsu M, Beggs JM. An open hypothesis: is epilepsy learned, and can it be unlearned? Epilepsy & behavior: E&B. 2008;13(3):511–522.
  25. 25. Scheffer M, Bascompte J, Brock WA, Brovkin V, Carpenter SR, Dakos V, et al. Early-warning signals for critical transitions. Nature. 2009;461(7260):53–59. pmid:19727193
  26. 26. Meisel C, Storch A, Hallmeyer-Elgner S, Bullmore E, Gross T. Failure of Adaptive Self-Organized Criticality during Epileptic Seizure Attacks. PLOS Computational Biology. 2012;8(1):e1002312.
  27. 27. Arviv O, Medvedovsky M, Sheintuch L, Goldstein A, Shriki O. Deviations from Critical Dynamics in Interictal Epileptiform Activity. Journal of Neuroscience. 2016;36(48):12276–12292.
  28. 28. Hobbs JP, Smith JL, Beggs JM. Aberrant Neuronal Avalanches in Cortical Tissue Removed From Juvenile Epilepsy Patients. Journal of Clinical Neurophysiology. 2010;27(6):380.
  29. 29. Yan J, Wang Y, Ouyang G, Yu T, Li Y, Sik A, et al. Analysis of electrocorticogram in epilepsy patients in terms of criticality. Nonlinear Dynamics. 2016;83(4):1909–1917.
  30. 30. Chang WC, Kudlacek J, Hlinka J, Chvojka J, Hadrava M, Kumpost V, et al. Loss of neuronal network resilience precedes seizures and determines the ictogenic nature of interictal synaptic perturbations. Nature neuroscience. 2018;21(12):1742–1752. pmid:30482946
  31. 31. Meisel C, Kuehn C. Scaling effects and spatio-temporal multilevel dynamics in epileptic seizures. PLoS One. 2012;7(2):e30371.
  32. 32. Da Silva FL, Blanes W, Kalitzin SN, Parra J, Suffczynski P, Velis DN. Epilepsies as dynamical diseases of brain systems: basic models of the transition between normal and epileptic activity. Epilepsia. 2003;44:72–83.
  33. 33. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014;137(8):2210–2230.
  34. 34. Kramer MA, Cash SS. Epilepsy as a Disorder of Cortical Network Organization, Epilepsy as a Disorder of Cortical Network Organization. The Neuroscientist. 2012;18(4):360–372.
  35. 35. Maturana MI, Meisel C, Dell K, Karoly PJ, D’Souza W, Grayden DB, et al. Critical slowing down as a biomarker for seizure susceptibility. Nature Communications. 2020;11(1):1–12. pmid:32358560
  36. 36. Milanowski P, Suffczynski P. Seizures Start without Common Signatures of Critical Transition. International Journal of Neural Systems. 2016;26(08):1650053.
  37. 37. Wilkat T, Rings T, Lehnertz K. No evidence for critical slowing down prior to human epileptic seizures. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2019;29(9):091104.
  38. 38. Harris TE. The Theory of Branching Processes. Grundlehren der mathematischen Wissenschaften. Berlin Heidelberg: Springer-Verlag; 1963. Available from: https://www.springer.com/de/book/9783642518683.
  39. 39. Wilting J, Priesemann V. Inferring collective dynamical states from widely unobserved systems. Nature Communications. 2018;9(1):2325.
  40. 40. Hsu D, Beggs JM. Neuronal avalanches and criticality: A dynamical model for homeostasis. Neurocomputing. 2006;69(10):1134–1136.
  41. 41. Badawy R, Freestone D, Lai A, Cook M. Epilepsy: ever-changing states of cortical excitability. Neuroscience. 2012;222:89–99.
  42. 42. Priesemann V, Munk MH, Wibral M. Subsampling effects in neuronal avalanche distributions recorded in vivo. BMC neuroscience. 2009;10(1):40.
  43. 43. Spitzner F, Dehning J, Wilting J, Hagemann A, Neto J, Zierenberg J, et al. MR. Estimator, a toolbox to determine intrinsic timescales from subsampled spiking activity. arXiv preprint arXiv:200703367. 2020;.
  44. 44. Tomen N, Rotermund D, Ernst U. Marginally subcritical dynamics explain enhanced stimulus discriminability under attention. Frontiers in systems neuroscience. 2014;8:151.
  45. 45. Meisel C, Klaus A, Vyazovskiy VV, Plenz D. The interplay between long-and short-range temporal correlations shapes cortex dynamics across vigilance states. Journal of neuroscience. 2017;37(42):10114–10124.
  46. 46. Priesemann V, Valderrama M, Wibral M, Le Van Quyen M. Neuronal avalanches differ from wakefulness to deep sleep–evidence from intracranial depth recordings in humans. PLoS Comput Biol. 2013;9(3):e1002985.
  47. 47. Meisel C. Antiepileptic drugs induce subcritical dynamics in human cortical networks. Proceedings of the National Academy of Sciences. 2020;117(20):11118–11125.
  48. 48. Meisel C, Schulze-Bonhage A, Freestone D, Cook MJ, Achermann P, Plenz D. Intrinsic excitability measures track antiepileptic drug action and uncover increasing/decreasing excitability over the wake/sleep cycle. Proceedings of the National Academy of Sciences. 2015;112(47):14694–14699.
  49. 49. Truccolo W, Donoghue JA, Hochberg LR, Eskandar EN, Madsen JR, Anderson WS, et al. Single-neuron dynamics in human focal epilepsy. Nature Neuroscience. 2011;14(5):635–641. pmid:21441925
  50. 50. Gast H, Niediek J, Schindler K, Boström J, Coenen VA, Beck H, et al. Burst firing of single neurons in the human medial temporal lobe changes before epileptic seizures. Clinical Neurophysiology. 2016;127(10):3329–3334. pmid:27592159
  51. 51. Schevon CA, Weiss SA, McKhann G Jr, Goodman RR, Yuste R, Emerson RG, et al. Evidence of an inhibitory restraint of seizure activity in humans. Nature communications. 2012;3:1060. pmid:22968706
  52. 52. Falconer M. Mesial temporal (Ammon’s horn) sclerosis as a common cause of epilepsy: etiology, treatment, and prevention. The Lancet. 1974;304(7883):767–770.
  53. 53. Badawy R, Macdonell R, Jackson G, Berkovic S. The peri-ictal state: cortical excitability changes within 24 h of a seizure. Brain. 2009;132(4):1013–1021.
  54. 54. Litt B, Lehnertz K. Seizure prediction and the preseizure period. Current opinion in neurology. 2002;15(2):173–177.
  55. 55. Baud MO, Kleen JK, Mirro EA, Andrechak JC, King-Stephens D, Chang EF, et al. Multi-day rhythms modulate seizure risk in epilepsy. Nature communications. 2018;9(1):1–10. pmid:29311566
  56. 56. Karoly PJ, Goldenholz DM, Freestone DR, Moss RE, Grayden DB, Theodore WH, et al. Circadian and circaseptan rhythms in human epilepsy: a retrospective cohort study. The Lancet Neurology. 2018;17(11):977–985. pmid:30219655
  57. 57. Levina A, Priesemann V. Subsampling scaling. Nature Communications. 2017;8:15140.
  58. 58. Murray JD, Bernacchia A, Freedman DJ, Romo R, Wallis JD, Cai X, et al. A hierarchy of intrinsic timescales across primate cortex. Nature Neuroscience. 2014;17(12):1661–1663. pmid:25383900
  59. 59. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014;137(8):2210–2230.
  60. 60. Suffczynski P, Kalitzin S, Da Silva FL. Dynamics of non-convulsive epileptic phenomena modeled by a bistable neuronal network. Neuroscience. 2004;126(2):467–484.
  61. 61. Brinkmann BH, Wagenaar J, Abbot D, Adkins P, Bosshard SC, Chen M, et al. Crowdsourcing reproducible seizure forecasting in human and canine epilepsy. Brain. 2016;139(6):1713–1722. pmid:27034258
  62. 62. Cook MJ, O’Brien TJ, Berkovic SF, Murphy M, Morokoff A, Fabinyi G, et al. Prediction of seizure likelihood with a long-term, implanted seizure advisory system in patients with drug-resistant epilepsy: a first-in-man study. The Lancet Neurology. 2013;12(6):563–571. pmid:23642342
  63. 63. Freestone DR, Karoly PJ, Cook MJ. A forward-looking review of seizure prediction. Current opinion in neurology. 2017;30(2):167–173.
  64. 64. Niediek J, Boström J, Elger CE, Mormann F. Reliable Analysis of Single-Unit Recordings from the Human Brain under Noisy Conditions: Tracking Neurons over Hours. PLoS ONE. 2016;11(12). pmid:27930664
  65. 65. Priesemann V, Levina A, Wilting J. Assessing Criticality in Experiments. In: The Functional Role of Critical Dynamics in Neural Systems. Springer; 2019. p. 199–232.
  66. 66. Jolly E. Pymer4: connecting R and Python for linear mixed modeling. Journal of Open Source Software. 2018;3(31):862.
  67. 67. Bates D, Mächler M, Bolker B, Walker S. Fitting linear mixed-effects models using lme4. arXiv preprint arXiv:14065823. 2014;.