Abstract
The classic brain criticality hypothesis postulates that the brain benefits from operating near a continuous second-order phase transition. Slow feedback regulation of neuronal activity could, however, lead to a discontinuous first-order transition and thereby bistable activity. Observations of bistability in awake brain activity have nonetheless remained scarce and its functional significance unclear. Moreover, there is no empirical evidence to support the hypothesis that the human brain could flexibly operate near either a first- or second-order phase transition despite such a continuum being common in models. Here, using computational modeling, we found bistable synchronization dynamics to emerge through elevated positive feedback and occur exclusively in a regimen of critical-like dynamics. We then assessed bistability in vivo with resting-state MEG in healthy adults (7 females, 11 males) and stereo-electroencephalography in epilepsy patients (28 females, 36 males). This analysis revealed that a large fraction of the neocortices exhibited varying degrees of bistability in neuronal oscillations from 3 to 200 Hz. In line with our modeling results, the neuronal bistability was positively correlated with classic assessment of brain criticality across narrow-band frequencies. Excessive bistability was predictive of epileptic pathophysiology in the patients, whereas moderate bistability was positively correlated with task performance in the healthy subjects. These empirical findings thus reveal the human brain as a one-of-a-kind complex system that exhibits critical-like dynamics in a continuum between continuous and discontinuous phase transitions.
SIGNIFICANCE STATEMENT In the model, while synchrony per se was controlled by connectivity, increasing positive local feedback led to gradually emerging bistable synchrony with scale-free dynamics, suggesting a continuum between second- and first-order phase transitions in synchrony dynamics inside a critical-like regimen. In resting-state MEG and SEEG, bistability of ongoing neuronal oscillations was pervasive across brain areas and frequency bands and was observed only with concurring critical-like dynamics as the modeling predicted. As evidence for functional relevance, moderate bistability was positively correlated with executive functioning in the healthy subjects, and excessive bistability was associated with epileptic pathophysiology. These findings show that critical-like neuronal dynamics in vivo involves both continuous and discontinuous phase transitions in a frequency-, neuroanatomy-, and state-dependent manner.
Introduction
The activity of neuronal populations self-organizes into long-range synchronized oscillations that regulate neuronal communication (Fries, 2005; Uhlhaas et al., 2009). The Brain Criticality Hypothesis posits that brains, like numerous other complex systems, operate near a critical point of a continuous, second-order transition (Beggs, 2008; Chialvo, 2010) between asynchronous and fully synchronous phases (Beggs and Timme, 2012; Muñoz, 2018; Palva and Palva, 2018; Wilting and Priesemann, 2019). This would lead to moderate mean synchronization and emergent power-law spatiotemporal correlations, and yield many functional benefits, such as maximal dynamic range (Kinouchi and Copelli, 2006), efficient communication (Shew et al., 2011), processing (Boedecker et al., 2012), and representational capacity (Bertschinger and Natschläger, 2004; Shriki and Yellin, 2016). Conversely, operation far from the phase transition would be associated with inadequate or excessive synchronization of neuronal activity, both of which are incompatible with healthy brain functions (Zimmern, 2020) and, for example, representative of comatose (Liu et al., 2014) and epileptic pathology (Hobbs et al., 2010), respectively.
The balance between excitation and inhibition has been proposed as a key control parameter for regulating neuronal systems to operate in the vicinity of this phase transition. As we will illustrate with a generative model, slow positive feedback regulation of neuronal synchrony could, however, lead to a discontinuous first-order phase transition and thus to the emergence of bistable dynamics (Wilson and Cowan, 1972; Freyer et al., 2012). Electrophysiological studies have already established that neuronal populations may show bistability per se For example, during sleep, slow (0.1-1 Hz) oscillations are associated with discrete up- and down-states (Steriade et al., 1993; Holcman and Tsodyks, 2005); in awake humans, an EEG study has revealed bistability in alpha band (8-14 Hz) oscillation amplitudes (Freyer et al., 2009). Nevertheless, whether these findings revealed neuronal systems operating near a first-order phase transition has remained unresolved while several converging lines of in vivo evidence support the notion of critical brain dynamics specifically near a second-order phase transition (Beggs and Timme, 2012; Palva and Palva, 2018; Wilting and Priesemann, 2019).
In this study, we asked whether awake resting-state human brain exhibits critical bistable dynamics indicative of neuronal systems operating near a first- rather than second-order phase transition. Moreover, as strong bistability is universally predictive of a potential for catastrophic shifts (Thom, 1972; Zeeman, 1976; Agu and Teramachi, 1978; Villa Martín et al., 2015), we hypothesize that high bistability neuronal synchronization dynamics would be indicative of a pathologic regimen where neuronal populations may abruptly switch from asynchronous to momentarily hypersynchronized, seizure-predisposing states.
We first used generative modeling to establish how varying degrees of bistable dynamics emerge as a consequence of introducing slow positive local feedback (Freyer et al., 2012) that is conceptually equivalent to increasing demands for limited resources (di Santo et al., 2018). We then analyzed a large body of source-reconstructed MEG and intracerebral stereo-EEG (SEEG) recordings of resting-state human brain activity. In both MEG and SEEG, we found that anatomically and spectrally widespread bistability characterized neuronal oscillations from δ (3-7 Hz) to high-γ (100-225 Hz) frequencies. In SEEG, conversely, excessive resting-state bistability was colocalized with the epileptogenic zone (EZ) and thereby associated with the pathophysiology underlying epilepsy. Bistable criticality thus constitutes a pervasive and functionally significant feature of awake resting-state brain dynamics.
Materials and Methods
The canonical model for bistability
Freyer et al. (2012) proposed the canonical Hopf bifurcation as a dynamical mechanism for bistability. They also presented evidence for bistability in a kinetic model of gene regulation, and suggested the universality of the canonical model for bistability across a broad range of biological systems. The canonical Hopf bifurcation is defined as follows:
Experimental design and statistical analysis
Bistability in the Kuramoto model
We studied first- and second-order phase transitions in a Kuramoto model with a modified noise term. The Kuramoto model is a generative model that can be used for studying the collective behaviors of a number of interconnected phase oscillators because of weak interactions (Breakspear et al., 2010; Rodrigues et al., 2016). In a Kuramoto model, the dynamics of each oscillator i is a scalar phase time series θi (θ∈ 0:2π), coupled into a population ensemble θ as follows:
In our model, the natural frequency of the oscillators has small variance; therefore, the order is correlated with coupling positively and with noise negatively. The total noise weight is fixed with a constant parameter η, and positive feedback is represented by the state-dependent noise term. As the model is tuned to operate near criticality by a specific κ value and a large ρ, the influence of the additive noise term becomes negligible, whereas the multiplicative noise exerts a strong effect on R because of the positive feedback from R. For example, when R is in the “down” mode, the multiplicative noise amplitude becomes large; therefore, R tends to remain in the down mode; when R is in the “up” mode, the multiplicative noise amplitudes are tuned down, thus causing the model to remain stuck in the up mode. Because the noise here is a Gaussian process, occasionally a very large amplitude noise sample can “kick” the model from up to down mode; vice versa, a small amplitude noise sample can kick the model from down to up mode. This results in bistability.
MEG resting-state recording and subjects
We recorded 10 min resting-state MEG data from 18 subjects (11 males, 31.7 ± 10.5, mean ± SD, yeas of age) at the BioMag Laboratory, HUS Medical Imaging Center (Helsinki, Finland). Subjects were seated in a dimly lit room and instructed to focus on a cross on the screen in front of them. Recordings were conducted at Meilahti Hospital (Helsinki, Finland). All subjects were screened for neurologic conditions. The study protocol for MEG and MRI data obtained at the University of Helsinki was approved by the Coordinating Ethical Committee of Helsinki University Central Hospital (ID 290/13/03/2013) and was performed according to the Declaration of Helsinki.
We also assessed working memory, attention, and executive functions in these subjects with a battery of neuropsychological tests. These included (see Fig. 4, x axis): Zoo Map Time, Toulouse-Pieron test, Digit Symbol Coding test, Zoo Map Plan, Digit Span forward and backward (BackDigits and ForwDigits, respectively), Letter-Number Sequencing, and Trail Making Test parts A and B. Some subjects had missing/invalid behavioral scores, but for each score, at least 16 valid scores were available for estimating neuro-behavioral correlations.
SEEG resting-state recording and subjects
We acquired 10 min of uninterrupted, seizure-free resting-state brain activity with eyes closed from 64 drug-resistant focal epilepsy patients (28 females, mean ± SD age, 30.1 ± 9.1 years; see Table 1) undergoing presurgical assessments. The subjects gave written informed consent for participation in research studies. The study protocol for SEEG, CT, and MRI data obtained in the La Niguarda Hospital were approved by the ethical committee of the Niguarda “Ca Granda” Hospital, Milan (ID 939), Italy, and was performed according to the Declaration of Helsinki.
Demographic data for the SEEG patient cohorta
Before surgery, medical doctors identified epileptogenic and seizure propagation zone by visual analysis of the SEEG traces (Cossu et al., 2005; Cardinale et al., 2013). Epileptogenic areas (generators) are the ROI that are necessary and sufficient for the origin and early organization of the epileptic activities (Luders et al., 1993). SEEG contacts recorded from such generators often show low-voltage fast discharge or spike and wave events at seizure onset. Seizure propagation areas (receivers) are recruited during seizure propagation, but they do not initialize seizures (Bartolomei et al., 2013; Jirsa et al., 2017). Contacts recorded from receivers show delayed, rhythmic modifications after seizure initiation in the generators. It is common to see regions demonstrating both generator and receiver dynamics; thus, they were identified as generator-receiver. In this study, we refer to generator, receiver, and generator-receiver collectively as EZ to distinguish them from those that were tentatively identified as healthy non-EZ regions (nEZ).
SEEG cortical sampling statistics
Preprocessing of the SEEG data yielded 7019 SEEG contacts in various cortical and subcortical gray matter locations. For investigating the LRTCs and bistability of cortical dynamics that were tentatively considered as normal, contacts recorded from subcortical structures and EZs were excluded. Contacts with >2.5% samples identified as “spiky” were also excluded. Thereby, ∼3 of 7 available contacts were excluded. From the resulting 4142 contacts, a small fraction of the contacts cannot be reliably assigned to a parcel by the segmentation software (Arnulfo et al., 2015b), these therefore were also excluded. This resulted in 4122 contacts (66.8 ± 24.5 per patient, range: 4-123) for analyses. Although the cortical sampling was heterogeneous across patients, with 4124 cortical nEZ contacts, we were able to cover 90 of the 100 Schaefer parcels with each parcel sampled by at least 3 subjects and at least 10 contacts.
Constructing surrogate data for statistical tests
To determine the chance level observations of detrend fluctuation analysis (DFA) and bistability index (BiS) in null hypothesis data (see Fig. 2D,E) (i.e., without embedded nonlinear critical-like structures) but with the same power spectrum as the real physiological signals (Lancaster et al., 2018), phase randomized Fourier transform surrogates of the broad-band time series were constructed for all MEG parcels and SEEG contacts (NMEG = 6800 and NSEEG = 4142). The surrogate broad-band data were filtered into narrow-band data, and their DFA and BiS estimates were subsequently computed. Real values were then compared against the distribution of the surrogate data for each frequency, and those samples with estimates >95th percentile of the surrogate data were considered as significant (see Fig. 2).
To examine whether any Yeo systems showed different BiS and DFA estimates from the null hypothesis (H0) that all the estimates were randomly distributed across the brain regions (see Fig. 3C–E), we constructed 105 label-shuffled surrogates to derive the 95th percentile of the H0 CIs for each of the MEG and SEEG band-collapsed θ-α and γ-band DFA and BiS brain maps.
MEG data acquisition
A 306-channel MEG system (204 planar gradiometers and 102 magnetometers) with a Vectorview-Triux (Elekta-Neuromag) was used to record 10 min eyes-open resting-state brain activity from 18 healthy adult subjects at the BioMag Laboratory (HUS Medical Imaging Center). For cortical surface reconstruction, T1-weighted anatomic MRI scans were obtained at a resolution of 1 × 1 × 1 mm with a 1.5 T MRI scanner (Siemens). This study was approved by the ethical committee of Helsinki University Central Hospital and was performed according to the Declaration of Helsinki. Written informed consent was obtained from each subject before the experiment.
MEG data preprocessing and filtering
The Maxfilter software with temporal signal space separation (Elekta Neuromag) was used to suppress extracranial noise in sensors and to interpolate bad channels (Taulu and Simola, 2006). Independent component analysis (MATLAB Fieldtrip toolbox, http://fieldtrip.fcdonders.nl) was next used to identify and remove components that were correlated with ocular (identified using the EOG signal), heartbeat (identified using the magnetometer signal as a reference), or muscle artifacts (Oostenveld et al., 2011). The FreeSurfer software (http://surfer.nmr.mgh.harvard.edu/) was used for subject MEG sources reconstruction, volumetric segmentation of MRI data, surface reconstruction, flattening, cortical parcellation, and neuroanatomical labeling with the Schaefer-2018 atlas (Schaefer et al., 2018) in which each cortical parcel is assigned to a functional system (Yeo et al., 2011), which informed later systems-level analysis. We here used resolutions of 400 and 100 parcels.
The MNE software package was used to create head conductivity models and cortically constrained source models with 5000-7500 sources per subject and for the MEG-MRI coregistration and for the preparation of the forward and inverse operators (Hamalainen and Sarvas, 1989; Hämäläinen and Ilmoniemi, 1994). For each MEG subject, data were source-reconstructed and collapsed to 400 parcels using a reconstruction-accuracy optimizing source-vertex-to-parcel collapsing method (Korhonen et al., 2014). The broadband time series of these parcels were then filtered into narrow-band time series using a bank of 20 Morlet filters with m = 5 and log-linearly spaced center frequencies ranging from 2 to 225 Hz. For group-level analyses, subject DFA and BiS estimates were morphed from 400 to 100 parcels.
SEEG data acquisition
Resting-state brain activity from 64 drug-resistant focal epilepsy patients (28 females, 30.1 ± 9.1 years, mean ± SD; see Table 1) was acquired as monopolar local field potentials (LFPs) from brain tissue with platinum-iridium, multilead electrodes using a 192-channel SEEG amplifier system (Nihon-Kohden Neurofax-110) at 1 kHz sampling rate. Each penetrating electrode shaft has 8-15 contacts, and the contacts were 2 mm long, 0.8 mm thick, and had an intercontact border-to-border distance of 1.5 mm (Dixi medical). The anatomic positions and amounts of electrodes varied according to surgical requirements (Cardinale et al., 2013). On average, each subject had 17 ± 3 (mean ± SD) shafts (range 9-23) with a total of 153 ± 20 electrode contacts (range 122-184, left hemisphere: 66 ± 54, right hemisphere: 47 ± 55 contacts, gray-matter contacts: 110 ± 25). The subjects gave written informed consent for participation in research studies and for publication of their data. This study was approved by the ethical committee (ID 939) of the Niguarda “Ca' Granda” Hospital, Milan, and was performed according to the Declaration of Helsinki.
SEEG preprocessing and filtering
Cortical parcels were extracted from presurgically scanned T1 MRI 3D-FFE (used for surgical planning) using the Freesurfer package (Fischl et al., 2002). A novel nearest-white-matter referencing scheme (its merits discussed in Arnulfo et al., 2015a) was used for referencing the monopolar SEEG LFP signals. An automated SEEG-electrode localization method was next used to assign each SEEG contact to a cortical parcel of Schaefer 100-parcel atlas with submillimeter accuracy (Arnulfo et al., 2015b). The SEEG electrodes were implanted to probe the suspected EZs while inevitably passing through healthy cortical structures. Contacts located at EZ are known to pick up frequent interictal spikes and would bias DFA estimates (Monto et al., 2007). Therefore, EZ contacts and contacts recorded from subcortical regions, such as thalamus, hippocampus, and basal ganglia were excluded from analysis.
Nevertheless, interictal events (IIEs), such as spikes, can be occasionally observed at nEZ locations in some subjects during rests. These IIE are characterized by high-amplitude fast temporal dynamics as well as by widespread spatial diffusion, which need to be excluded to avoid bias to DFA and BiS estimates. We followed approach used in to identify such IIEs. Briefly, each SEEG contact broad-band signal was partitioned into nonoverlapping windows of 500 ms in length; a window was tagged as “spiky” and discarded from LRTCs and bistability analyses when at least three consecutive samples exceeding 7 times the SD above the channel mean amplitude. Last, narrow-band frequency amplitude time series was obtained by convoluting the broad-band SEEG contact time series with Morlet wavelets that were identical to that of MEG data.
Estimating LRTCs using DFA
LRTCs in 1D time series can be assessed with several metrics (Eke et al., 2000); and in this study, linear DFA was used to assess specifically how fast the overall root mean square (RMS) of local fluctuations grows with increasing sampling period (Linkenkaer-Hansen et al., 2001). An estimated DFA exponent reflects the finite-size power-law scaling in narrow-band amplitude fluctuations based on the assumption that the gradual evolution of a mono-fractal process time series would result in a normal distribution where the fluctuation can be captured by the second-order statistical moments, such as variance. The computation of DFA can be described briefly as follows (for rationales and technical details of the algorithm, see Hardstone et al., 2012):
The signal profile (X) of a signal was computed by computing the cumulative sum of a demeaned narrow-band amplitude of a MEG parcel or SEEG contact time series.
A vector of window widths (T) was defined in which the widths were linearly spaced on log10 scale between 10 and 90 s. The same scaling range was used across frequencies and for both MEG and SEEG (i.e., identical vector of T). The lower boundary of 10 s was set to safely avoid high nonstationarity and the filter artifacts (i.e., 20 cycles of the slowest rhythm of 2 Hz); the upper bound of 90 s was 15% of total sample of the resting-state recording.
For each window width t∈T, X was partitioned into an array of temporal windows, in which each window was of length t, with 25% overlap between windows W(t).
For each window w∈W(t), a detrended signal wdetrend was obtained by removing the linear trend (i.e., subtracting the least-square fit of the samples from the samples of w, and then taking the RMS of wdetrend).
Finally, F(t), the detrended fluctuation of window size t, was obtained by computing the mean of RMS (wdetrend).
By repeating Step 3 for all window lengths of T defined in Step 2, F, a vector of F(t), t∈T, was obtained. The DFA exponent is the slope of the trend line of F as a function of T on log-log scale (see Fig. 1I).
Different value ranges of the DFA exponent are commonly thought to indicate different temporal dynamics in a 1-dimensional time series (Hardstone et al., 2012). A DFA exponent of 0.5 indicates that the time series is indistinguishable from a random walk without memory; an exponent between 0.5 and 1 indicates that the time series has a memory with positive correlations; an exponent between 0 and 0.5 indicates that the time series has a memory with negative (anti-) correlations; and an exponent between 1 and 2 indicates a nonstationary time series. When obtaining the signal profile in Step 1, the narrow-band oscillation amplitudes are treated as the “rate of change.” Hence, the scaling relationship between a temporal window size t (in seconds) and the detrended fluctuation F(t) indicates how large the fluctuations in the “work” of a neuronal oscillation could demonstrate over time. A sublinear DFA exponent (i.e., 0.5 < DFA < 1) can be interpreted as the neuronal oscillations showing a tendency to conserve energy. For example, a DFA exponent of 0.75 means that, as the window size wi doubles, the RMS in signal profile is expected to increase by 68.2%. On the other hand, a superlinear exponent indicates the growth of fluctuations is unbounded. For example, a DFA exponent of 1.2 means that as the window size doubles, the RMS is expected to increase by 115%.
Estimating BiS
The BiS of a power time series (R2) is derived from model comparison between a bimodal or unimodal fit of its probability distribution function (pdf); a large BiS means that the observed pdf is better described as bimodal, and when BiS → 0 the pdf is better described as unimodal. We followed the approach used previously (Freyer et al., 2012; Roberts et al., 2015) to compute BiS. First, to find the pdf of power time series R2, the empirical R2 was partitioned into 200 equal-distance bins and the number of observations in each bin was tallied. Next, maximum likelihood estimate was used to fit a single-exponential function (i.e., the square of a Gaussian process follows an exponential pdf) as follows:
Next, the Bayesian information criterion (BIC) was computed for the single- and bi-exponential fitting as follows:
Last, the BiS estimate is computed as the log10 transform of difference between the two BIC estimates as follows:
Morphing MEG and SEEG data into a standard atlas
The MEG and SEEG group-level analyses were conducted in a 100-parcel standardized Schaefer atlas (Schaefer et al., 2018). Although the SEEG electrode locations were variable across subjects, we have previously confirmed with random cohort-split analyses that only using half of these subjects could readily produce the same group level observations (Arnulfo et al., 2020). When morphing narrow-band DFA and BiS estimates of individual SEEG contacts into the Schafer atlas, the resulting group-level parcel estimates were heterogeneous in terms of sampling. For example, one parcel may contain observations from a varying number of electrodes and/or subjects. Hence, the group-level estimate of each parcel was the median of all observations, and the estimate for each of the seven Yeo subsystems was the median of its constituent parcels. Furthermore, 10 parcels sampled by ≤3 subjects and 10 contacts were excluded from group-level analysis. The group mean parcel metrics were identical to that of raw data (see Fig. 2D,E, overlay curves).
Supervised classification for EZs
The group-level frequency clustering analysis revealed that much of the narrow-band data were topologically correlated (see Fig. 3). Hence, for the classification task, 20 narrow-band metrics were also collapsed into four frequency clusters as δ, θ−α, β, and γ band (see Fig. 5A). As subjects varied greatly in their DFA and BiS estimates, band-collapsed data were normalized within subjects as [X-median(X))./max(X-median(X)], and thereby the differences between EZ and nEZ within subjects remained. The effect size of differences between band-collapsed and normalized DFA and BiS estimates were assessed with Cohen's d and compared with the 99th percentile of Cohen's d observed from 1000 EZ-nEZ label-shuffled surrogate data (see Fig. 5C).
The feature importance of these neuronal estimates was assessed with the SHapley Additive exPlanations (SHAP) values (Lundberg et al., 2017). In addition to the neuronal scores, the contact location in Yeo systems was also included as an additional feature (see Fig. 4D). The SHAP values is a generic metric to explain any tree-based model by explicating the local and global interpretability of features, which advances the transparency that conventional classification approaches lack. For solving the EZ-classification problem, the nonparametric random-forest method (Hastie et al., 2009) was used. The random-forest algorithm is a machine learning method that uses bootstrapped training dataset and combines the simplicity of decision trees with extended flexibility to handle new data.
Data availability
The data to support the main results can be downloaded from https://figshare.com/ndownloader/files/41233434, which includes the following: extended discussions on the theory for bistability and criticality, relevant generative models, candidate mechanisms for resource loading and state-dependent noise; the description of the MEG and SEEG analyses pipeline; SEEG sampling statistics; additional traces showing bistable critical dynamics with exemplary bistability and DFA fitting; interim results to support the MEG and SEEG band-clustering by narrow-band topological similarity; interim results to support the correlations between LRTCs and bistability on group-level and within Yeo functional systems; anatomic specificity of the θ−α and γ band DFA and BiS estimates and corresponding statistical tests (see Fig. 3); statistical information to support the MEG behavioral correlations of DFA and BiS estimates (see Fig. 4); and technical details and interim results to support the supervised and unsupervised EZ classification (see Fig. 5).
Raw data and patient details cannot be shared because of Italian governing laws and Ethical Committee restrictions. Intermediate data, final processed data, and all code that support the findings of this study are available from the corresponding authors on reasonable request.
Results
Positive local feedback is a control parameter for bistable criticality in silico
To assess the emergence of bistability in synchronization dynamics and its relationship with criticality, we built a variant of the classic Kuramoto model that is a universal generative model of synchronization dynamics (see Materials and Methods). The classic Kuramoto model has a single control parameter, κ, that defines the coupling strength between neuronal oscillators (Breakspear et al., 2010). Synchrony among the oscillators is quantified by the “order” parameter R (Eq. 5; Fig. 1A), and a gradual increase of κ leads to different level of order and associated critical dynamics (Fig. 1B, top).
Bistability is caused by elevated positive local feedback. A, The order parameter R(t) increases (top) as phases (Ph.) of the oscillators in the model become increasingly synchronized (bottom). Re/Im, Real/imaginary part of the complex R. B, Exemplary segments of order R time series when the model is in the sub-, super-, critical, and bistable phases indicated in F. C–E, Criticality estimates as a function of the local positive feedback (ρ) and neuronal coupling (κ). C, The mean order. Two arrows indicate the amount of coupling required for the model to transition from asynchrony (R = 0.1) to hypersynchrony (R = 0.9). D, The DFA exponent as a measure of LRTCs in the fluctuations of R(t). E, The BiS assessed from the R2 time series. C–E, Each pixel is the mean of 50 simulations. F, An overlaid regimen map based on observation from C–E; classic criticality (with a second-order phase transition) is associated with small positive feedback ρ (black dashed); bistable criticality is seen at mid-to-high degree of ρ (enclave inside red line). G, Probability distribution of R in bistable (top) and classic (bottom) criticality. The peak DFA (black line) coincided with the phase transition (i.e., moderate R). H, Probability density (pdf) of R and (I) the detrended fluctuations (DF) and the DFA exponents (the slope) of the time series from B, color-coded.
To study how bistability emerges in the dynamics of order, we introduced a second parameter ρ to scale local positive feedback. Such feedback is a generic mechanism for bistability, and here it was implemented as state-dependent noise (Eqs. 2–4) that is conceptually equivalent to state dependency in canonical models (Izhikevich, 2007; Freyer et al., 2012; Breakspear, 2017), resource-loading mechanisms in conceptual models (di Santo et al., 2016; Buendía et al., 2020), and various synaptic strengths in firing rate models (Wilson and Cowan, 1972).
When the model is controlled by weak feedback (i.e., small values of ρ), the model behaved identically to the classic critical-like ensemble with a second-order phase transition where order monotonically increased as κ increased (Fig. 1C). At moderate order (i.e., at the phase transition between low and high order), power-law long-range temporal correlations (LRTCs) emerged and delineated a critical region (Fig. 1D). Here, LRTCs were assessed using the detrended fluctuation analysis exponent (DFA), and DFA > 0.6 demarcates the critical region.
With strengthening feedback, the order parameter became progressively bistable (Fig. 1B, bottom) as evidenced by increasing values of the BiS (Fig. 1E), an index of the relative fit of a bimodal versus a unimodal pdf to the order power time series (R2). The bistable dynamics were seen exclusively within the critical region (Fig. 1F). High bistability (Fig. 1B,F,G, red triangles) was characterized by (1) a sudden increase in the order in response to a relative small increase in coupling compared with classic criticality (Fig. 1C, arrows); (2) a sharp DFA peak indicating a narrowing of the critical region (Fig. 1G); (3) a clear bimodal distribution of order (Fig. 1H); and (4) abnormal super-linear LRTCs (DFA > 1, Fig. 1I), indicating unbounded growth of oscillation amplitudes (see Materials and Methods). These in silico findings thus showed that synchronization of oscillators may exhibit a continuum between classic second-order critical and bistable first-order critical dynamics when under the influence of positive local feedback.
Bistable criticality characterized brain dynamics in vivo
We next assessed neuronal bistability and critical dynamics in 10 min human resting-state brain activity in SEEG (N = 64) and source-reconstructed MEG (N = 18). For SEEG data, we first restricted analysis to neocortical gray matter contacts outside of the EZ (Figs. 2–4). Although the anatomic sampling with SEEG is heterogeneous across patients, the present cohort size yielded essentially a full coverage of the cerebral cortex. We estimated LRTCs using DFA and bistability with BiS for narrow-band (Morlet) SEEG and MEG source amplitude time series that predominantly reflect local cortical synchronization dynamics.
Bistability and LRTCs were robust, large-scale phenomena in the resting-state brain. A, B, Five minutes of broad-band and narrow-band power (R2) time series from (A) a MEG parcel in visual area (Vis) in 1 subject and (B) an SEEG contact in middle frontal gyrus in one patient. Insets, Bistability as narrow-band traces switching between “up” and “down” states. C, Group-level probability (z axis) distribution of narrow-band (y axis) mean amplitude (R). D, DFA exponents. E, BiS estimates. Data were pooled over all nEZ SEEG contacts or MEG parcels; subject and contact/parcel number indicated in C. Black lines indicate mean of real data. Red dashed lines indicate 99th percentile of surrogate observation. F–H, Examples of narrow-band DFA and BiS probability distribution as indicated by colored arrows in D, E.
Bistability and LRTCs were coexisting, correlated phenomena in resting-state MEG and SEEG. Neuroanatomical similarity (Spearman's correlation) between group-average narrow-band BiS and DFA estimates of (A) MEG and (B) SEEG in Schaefer 100-parcel atlas. Red boxes represent frequency clusters showing high similarity. C, Narrow-band group-averaged estimates were collapsed into θ−α (5.4–11 Hz) and γ (40–225 Hz) band based on similarity shown in A. White-out columns in SEEG data represent excluded parcels because of insufficient sampling. D, E, Parcel-wise group-average θ−α band BiS maps for (D) MEG and (E) SEEG. F, Kruskal–Wallis one-way ANOVA for group-level differences in DFA and BiS estimates between Yeo systems. Dashed line indicates –log10(p value) > 1.3 (i.e., p < 0.05). Correlations between group-average parcel BiS and DFA estimates in θ−α (cross) and γ band (circles) in (G) MEG and (H) SEEG, −log10(p) > 5 (FDR-corrected). I, Spearman's correlations between within-subject-average BiS and DFA estimates in Yeo systems (subject NMEG per system =18; NSEEG per system = 50 ± 9.4, range: 36–60, variable SEEG subject N per system because of heterogamous spatial sampling).
Executive functions were correlated with θ−α band DFA and BiS estimates in MEG subjects. A, Spearman's correlation between subject neuropsychological test scores and within subject mean parcel θ−α band DFA. B, BiS estimates collapsed over parcels. Dashed lines indicate 5th and 95th percentile of correlations for surrogate data (Nsurrogate = 105, FDR-corrected). C, D, Subject Zoo map time rank and θ−α band parcel-collapsed (C) DFA and (D) BiS estimates. Each marker represents 1 subject. E, Fraction of parcels that showed significant correlation between neuropsychological test scores and individual parcel θ−α band DFA and (G) BiS estimates (p < 0.05, FDR-corrected). F, Parcels showing significant correlations between Zoom map time scores and θ−α band DFA and (H) BiS estimates.
Bistability was anatomically widespread and spectrally prevalent
Visual inspection of the narrow-band MEG and SEEG amplitude time series revealed abundant examples of bistability as intermittent switching between low- and high-amplitude oscillations (Fig. 2A,B). Comparison with surrogate data showed that both MEG source signals and SEEG electrode contact LFP signals exhibited significant (>95th percentile CI of the surrogate data) bistability and LRTCs across broad frequencies (Fig. 2C–E). MEG showed peak DFA and BiS estimates in the α (11 Hz) frequency band, whereas in SEEG, the BiS peak extended over δ (2-4 Hz), θ (4-7.6 Hz), and α (10-13 Hz) bands (Fig. 2F–H). Across frequencies, the narrow-band DFA and BiS estimates in SEEG were overall greater than in MEG (Fig. 2E–G).
Neuroanatomical structures of bistability and LRTCs were correlated
We next studied the neuroanatomical structure of bistability and inspected its anatomic relationship with LRTCs across frequencies. We first collapsed narrow-band BiS and DFA estimates for MEG parcels (400) and SEEG contacts into a standard atlas of 100 cortical parcels. Next, the neuroanatomical similarity within bistability and LRTCs and between them were assessed by computing Spearman's correlations between narrow-band parcel BiS and DFA estimates.
Both MEG and SEEG showed high anatomic similarity between neighboring frequencies (Fig. 3A,B). Correlations between slow and fast rhythms were negative in MEG and weak in SEEG. This indicated that regions tended to show bistability and criticality in a cluster of high or low frequencies, but not both. Based on the neuroanatomical similarities (Fig. 3A,B), red boxes), we collapsed narrow-band estimates into θ−α (5.4-11 Hz) and γ-band (45-225 Hz) for further analyses (Fig. 3C). The partitioning of β (15-30 Hz) band was not consistent and thus was not included.
The MEG and SEEG θ−α band BiS maps showed distinct neuroanatomical patterns. We therefore next asked whether DFA and BiS metrics differed between the Yeo functional systems by testing the null hypothesis (H0) that all observed values were randomly distributed across brain regions. For each of the MEG and SEEG band-collapsed DFA and BiS brain maps, we constructed 105 label-shuffled surrogates to derive the H0 CIs. In MEG, visual (VIS), somatomotor (SM), and dorsal attention (DAN) networks (Fig. 3C,D) exhibited greater BiS than expected by chance (p < 0.05, two-tailed permutation test). SEEG showed high BiS in frontoparietal (FP), ventral attention (VAN), default network (DEF), and limbic (LIM) systems (Fig. 3C,E). VIS showed the lowest BiS in SEEG (p < 0.05, two-tailed permutation test), although these values were similar to those observed in MEG.
A Kruskal–Wallis test for variance among subjects' BiS and DFA estimates revealed that, in SEEG, individuals showed different levels of BiS and DFA estimates between systems (Fig. 3F) with bistability greater in DEF, FP, and LIM than in VIS and SM (unpaired t test, p < 0.05, FDR-corrected). There was no statistically significant regional variation in MEG data.
In both MEG and SEEG, group-average parcel bistability was correlated with LRTCs (Fig. 3G,H). We validated this analysis in narrow-band frequencies and found the results to converge well. To further validate this relationship, we averaged parcel BiS and DFA within subjects for each Yeo system and found that the subject BiS were indeed correlated with their DFA estimates on systems level (Fig. 3I).
Bistability was functionally significant in healthy MEG subjects
We next asked whether bistability and LRTCs would predict individual differences in cognition. We assessed working memory, attention, and executive functions with eight neuropsychological tests (see Materials and Methods). The BiS and DFA estimates of the 100 parcels were averaged to four neurophysiological scores for each subject: DFAθ−α, DFAγ, BiSθ−α, and BiSγ, and we then tested their correlation with the neuropsychological scores. The θ−alpha band BiS and DFA were negatively correlated (Fig. 4A,B) with execution time in the “Zoo Map” flexible planning task (p < 0.05, FDR-corrected for the eight neuropsychological tests using the Benjamini–Hochberg procedure). The negative correlation indicated that subjects with greater θ−α band bistability and larger LRTCs completed the task faster (Fig. 4C,D), which is in line with prior observations linking LRTCs with cognitive flexibility (Simola et al., 2017).
To inspect the neuro-behavioral correlations in greater anatomic detail, we computed Spearman's correlations between neuropsychological scores and individual parcel BiS and DFA estimates. A large fraction of cortical parcels showed significant neuro-behavioral correlations of θ−α band BiS and DFA estimates with Zoo Map time, but not with other neuropsychological scores (Fig. 4E,G, p < 0.05, FDR-corrected for 100 parcels using the Benjamini–Hochberg procedure for each neuropsychological test). The correlations of Zoo Map time with DFA estimates were most pronounced in fronto-parietal, limbic, somatosensory, and visual areas (Fig. 4F), whereas the correlations with BiS estimates were widespread (Fig. 4H).
Excessive bistability characterized the EZ
Excessive bistability may predispose complex systems to catastrophic events. Under the influence of strong state-dependent noise, our model demonstrated increased sensitivity to coupling strength (Fig. 1E), which suggests that strong bistability could be an early sign of shift toward supercritical hypersynchrony events (i.e., the clonic phase of epileptic seizures) (Jiruska et al., 2013). We thus asked whether bistability estimated from seizure-free, interictal activity-free resting-state SEEG recording could be informative about epileptic pathophysiology. In particular, we addressed whether bistability could delineate the EZ and dissociate EZ signals from signals in nEZ contacts that reflect more healthy forms of brain activity.
Representative time series (Fig. 5A,B) showed that the EZ contacts did not show conspicuous epileptic IIE, and the sparse IIEs were removed from analysis where found (see Materials and Methods). Interestingly, elevated >80 Hz bistability of the EZ contact was already a visually salient characteristic and stronger in EZ than in a nearby nEZ contact from the same region. We assessed bistability and LRTCs in narrow-band frequencies at the group level for nEZ- versus EZ-electrode contacts (Fig. 5C,D). Collapsing narrow-band DFA and BiS estimates into broader frequency bands revealed significant differences between nEZ- and EZ-electrode contacts in β- and γ-band BiS estimates with effect sizes (Cohen's d) of 0.5 and 0.65, respectively (Fig. 5E). There was also a difference between nEZ- and EZ-electrode contacts in the δ-band DFA exponent with a Cohen's d of 0.2.
Bistability showed strong predictive power for epileptogenic pathophysiology. A, B, Five minutes of broad-band traces and narrow-band power (R2) time series of an EZ (A) and an nEZ (B) cortical location recorded with two distinct electrode shafts in 1 subject. The two contacts were 19.7 mm apart within the supervisor frontal gyrus, and they were referenced with the same nearest white matter contact (Arnulfo et al., 2015a). C, Average normalized narrow-band BiS and (D) DFA estimates for all EZ (pink) and nEZ contacts (green) of the patient cohort. Shades represent 25th and 75th percentiles. E, The effect size of differences between EZ and nEZ contacts in frequency-collapsed BiS (red) and DFA (black). Dashed line indicates 99th percentile observation from surrogate data (Nsurrogate =1000). F, Feature importance estimated using SHAP values. G, The AUC of receiver operating characteristics averaged across subjects (black) and the AUC of pooled within-subject classification results (blue) when using (1) DFA alone, (2) BiS alone, (3) D&B, and (4) D&B plus contact loci in Yeo systems (D&B(Y)). Dashed lines indicate 99th percentile of AUC observed from 1000 surrogates created independently for each of the four feature sets. H–J, Post hoc inspection of results derived using D&B(Y) feature set (black marker in G). H, Spearman's correlation (p < 10−6, n = 55) between individual AUC and within-subject mean Cohen's d between EZ and nEZ in band-collapsed DFA and BiS. I, Receiver operating characteristics of classification within subjects (thin lines) and mean receiver operating characteristic (thick). J, Within-patient prediction precision as a function of TPR indicated by the magenta box from I. Red marker represents the population mean. Precision = true positive ÷ reported positive.
These group-level findings suggest that both bistability and LRTCs could constitute informative features for classifying nEZ- and EZ-electrode contacts. We thus conducted an EZ-vs-nEZ classification analysis using random forest algorithm (Breiman, 2001) and with frequency-collapsed BiS and DFA estimates as neuronal features, with the electrode contact location in Yeo systems as an additional feature. The cross-validation for the classification was a partition of 80%:20% (training:test) with 500 iterations. This revealed a reliable outcome with the area under curve (AUC) for the receiver operating characteristic reaching AUC = 0.8 ± 0.002 (mean ± SD). To identify the most informative components for the classifier, we assessed global and within-subject feature importance with the SHAP values (Lundberg et al., 2017). The SHAP values corroborated that γ- and beta-band BiS estimates were indeed the most important features, followed by contact location (within a Yeo system), and δ-band DFA (Fig. 5F).
Given these encouraging group-level and contact-classification results, we quantified the within-subject accuracy of neuronal bistability in localizing the epileptogenic area. We used leave-one-out validation so that the EZ-vs-nEZ contact classification was performed for each patient with the rest of the patients serving as training data. Additionally, to independently evaluate the contributions of BiS and DFA estimates to classification accuracy, we implemented the classification with four feature sets: DFA alone, BiS alone, combining DFA and BiS (D&B), and combining D&B and SEEG contact location in Yeo systems.
Overall, the within-subject classification accuracy for EZ contacts was higher than chance level across all feature combinations (Fig. 5G). Classification with all features yielded the best performance at an average AUC of 0.7 ± 0.14 (Fig. 5G, black marker). BiS alone yielded a greater AUC than DFA alone. Including the contact-brain system as an additional feature to D&B increased the AUC by 0.06. The subject AUC values were correlated with the subject-specific differences in DFA and BiS estimates between EZ and nEZ (r = −0.53, p < 10−6) (Fig. 5H) and were not affected by the total numbers of contacts, EZ contacts, nor the ratio of EZ and nEZ contacts (Pearson's correlation coefficient, r = −0.06, p = 0.66; r = −0.07, p = 0.61; r = −0.09, p = 0.50; respectively). Finally, the classifier yielded an average precision of 0.74 ± 0.30 (mean ± SD). While the true positive rate was 0.24 ± 0.17, the false-positive rate was only 0.03 ± 0.03 (Fig. 5G,H), which shows that most EZ contacts identified with the bistability-based classification were correct, although the classifier did not identify all true EZ contacts.
Discussion
We investigated whether the awake resting-state human brain exhibits critical bistability indicative of neurons operating near a first- rather than second-order phase transition. We advance here the first comprehensive modeling and experimental evidence for critical-like neuronal oscillation dynamics in a continuum between a second- and a first-order phase transition. While classical models of criticality have considered one control parameter (e.g., the excitation-inhibition balance) to operate near a second-order phase transition, our findings suggested that an additional control parameter, positive local feedback, is required to operate in such a continuum. We showed in MEG and SEEG data that bistability in neuronal oscillations, likely because of an underlying first-order phase transition, is likely an important characteristic of brain dynamics and functionally relevant to both the healthy and diseased human brain functions.
Both second-order (Beggs, 2008; Chialvo, 2010) and first-order (Millman et al., 2010; Scarpetta et al., 2018) phase transitions, as well as a hybrid type (Buendía et al., 2020, 2022), have been proposed to account for the critical phenomena observed in models and in actual brains across species in recent studies. While a large body of brain criticality literature has been dedicated to the study of phase transitions and associated bifurcation mechanisms for avalanche dynamics, fewer efforts have considered the same subject matter in the long-range temporal correlations of neuronal oscillations.
Feedback mechanisms (Buendía et al., 2020) and local hysteresis (Buendía et al., 2021, 2022) are known to invoke rich dynamics near criticality. In line with the canonical bistability models (Thom, 1972; Freyer et al., 2012), we show that state-dependent feedback could lead to bistable oscillations in a model conventionally used for generating unimodal oscillations. Notably, bistability occurred within a regimen of scale-free LRTCs dynamics, which suggested an underlying first-order phase transition near criticality. The positive feedback in our model hence represented a necessary second control parameter for shifting the model in a continuum between a second- and a first-order phase transition. While this conceptual model is suitable for generating narrow-band oscillations in mesoscopic ensembles with uniform connectivity, it is limited and does not produce the large-scale critical phenomena with anatomic and frequency specificity that we observed in MEG and SEEG data.
A positive feedback loop is thought to be a generic mechanism (Sornette and Ouillon, 2012; Hugo et al., 2013) for bistability in a wide range of modeled and real-world complex systems, including the canonical sand-pile model (di Santo et al., 2016) and its variations (Buendía et al., 2020), ecosystems (Kéfi et al., 2007; Villa Martín et al., 2015), gene regulatory networks (Dubnau and Losick, 2006; Freyer et al., 2012; Kuwahara and Soyer, 2012), intracellular compartments (Mitrophanov and Groisman, 2008; Bednarz et al., 2014), and network models of spiking neurons (Cowan et al., 2016; di Santo et al., 2018).
The positive feedback in our model was represented by a state-dependent control parameter. In more detailed models for ensemble activity similar to that measured in SEEG and MEG, bistability can be introduced through multiple state-dependent mechanisms. Several theoretical studies have suggested that state dependency could be a slowly fluctuating physiological parameter that reflects excitability (Jirsa et al., 2014; Cowan et al., 2016) and resource demand (di Santo et al., 2018). Unlike the constant control parameters in our model, various factors in vivo, such as arousal levels and tasks (Fontenele et al., 2019; Fosque et al., 2021), and individual variability (Fosque et al., 2022; Fuscà et al., 2023), may be characterized by dynamically changing control parameters affecting the neuronal operating point in the state space. For microscopic neuronal dynamics, three mechanisms have been proposed to account for feedback and state dependency (Doiron et al., 2016), whereas the mechanisms for mesoscopic and macroscopic state dependency remain unclear.
In MEG and SEEG, we observed significant bistability and LRTCs in spontaneous amplitude fluctuations of neuronal oscillations from 2 to 225 Hz across the neocortex. In line with our modeling results, MEG and SEEG bistability and LRTCs were positively correlated with within but not between θ−α, and γ bands on whole-brain level (Fig. 3G,H) and within functional systems across individuals (Fig. 3I). These cohort-level, frequency-specific correlations suggested that: (1) bistable dynamics are likely dependent on distinct frequency bands of neuronal oscillations as previously reported for LRTCs and synchrony (Palva et al., 2013; Zhigalov et al., 2015; Fuscà et al., 2023); and (2) brain dynamics follow a gradient from “low LRTCs and weakly bistable” to “high LRTCs and strongly bistable.” As the key evidence for functional significance, higher θ−α band bistability in healthy MEG subjects was a trait-like predictor (Palva et al., 2013) for better cognitive performance, whereas in the SEEG subjects, excessive β and γ band bistability was associated with epileptogenic pathophysiology, implying a “catastrophic shift” in the epileptic brain (Thom, 1972).
Bistability, in general, could be associated with a dichotomy of both beneficial and detrimental outcomes. Organisms can operate in a bistable mode that is thought to reflect a dynamic motif favorable to adaptation and survival (Dubnau and Losick, 2006; Freyer et al., 2012; Bednarz et al., 2014). However, excessive bistability characterizes catastrophic shifts in many complex systems (Sornette and Ouillon, 2012; Hugo et al., 2013), such as irrevocable environmental changes (Barnosky et al., 2012; Boerlijst et al., 2013; Villa Martín et al., 2015), wars and conflicts (Díaz, 2017), and sudden violent vibrations in aerodynamic systems (Xin and Shi, 2015). In healthy MEG subjects, bistability and LRTCs were correlated, and higher θ−α band BiS and DFA estimates predicted better cognitive performance. In epileptic patients, excessive β- and γ-band BiS, but not DFA, best characterized EZ. The findings from these two datasets suggest a functional gradient, wherein moderate bistability reflect functional advantages, while excessive bistability be a sign of pathologic hyperexcitability. The pathologic bistability could be associated with hyperexcitability, excessive synchrony, high resource demands, and likely subsequent oxidative stress and tissue damage (Salim, 2017). This speculation is in accordance with biophysical models of seizures that suggest a crucial role of a discontinuous transition (Freyer et al., 2012; Breakspear, 2017) in generalized seizures (Breakspear et al., 2006; Jirsa et al., 2014).
With invasive SEEG, we found consistent and accurate performance of the BiS estimates in EZ localization, which suggests a great potential for broader clinical utility (e.g., using noninvasive MEG or EEG). Future work could exploit the presence of widespread bistability to large-scale biophysical models of neural dynamics, building on the analytic link between the simplified model used here and physiologically derived neural mass and mean field models. Whereas the simple model yields dynamical insights, the large-scale biophysical models are crucial for understanding biological mechanisms (Breakspear, 2017), including those that describe seizure propagation in individual patient brain networks (Proix et al., 2017; Jirsa et al., 2023).
In theory, power-law scaling (Beggs and Timme, 2012; Beggs, 2022) and bistability (Freyer et al., 2012) may be seen in noncritical systems. Hence, to claim criticality, strong criteria must be satisfied (Sethna et al., 2001). In both human MEG and SEEG data, we have shown that local LRTCs are correlated with large-scale phase synchrony (Fuscà et al., 2023) and avalanche dynamics (Palva et al., 2013; Zhigalov et al., 2015); synchrony and scale-free networks are also correlated with distinct modular structures (Zhigalov et al., 2017). Here, we present novel evidence for colocalization of LRTCs and bistability that is relevant for cognitive function. These correlations between categorically different estimates of local and global critical dynamics offer a strong support to the brain criticality hypothesis. Building on these novel results, we propose that human brain dynamics operates in a continuum between the classic second- and the first-order phase transition near criticality, which is best described by at least two control parameters: the E/I balance and positive local feedback mechanism.
Footnotes
This work was supported by Academy of Finland SA 266402, 303933, and SA 325404 to S.P.; Academy of Finland SA 253130 and 296304 to J.M.P.; Ella and Georg Ehrnrooth Foundation 14-10553-2 to S.H.W.; and Finnish Cultural Foundation 00220071 to S.H.W. We thank our Italian collaborators from the Center of Epilepsy Surgery “C. Munari” at Niguarda Hospital, in particular Annalisa Rubino and Dr. Francesco Cardinale, for support; Dr. James Roberts for sharing code; and Drs. Serena di Santo, Stewart Heitmann, and Alexander Zhigalov for insightful discussions.
The authors declare no competing financial interests.
- Correspondence should be addressed to Sheng H. Wang at sheng.wang{at}helsinki.fi or J. Matias Palva at matias.palva{at}helsinki.fi
