Adaptive proximity to criticality underlies amplification of ultra-slow fluctuations during free recall

TL;DR

The brain may use a delicate balance near a critical state to generate slow neural waves and transition into spontaneous thought.

Contribution

A simple computational model shows how cortical networks near criticality can amplify slow fluctuations and switch to generative modes.

Findings

A random recurrent network near criticality reproduces ultra-slow fluctuations seen in human visual cortex.

Slight increases in network gain selectively amplify slow waves without destabilizing the system.

The model matches intracranial EEG data during rest and free recall, supporting the role of criticality in spontaneous brain activity.

Abstract

Ultra-slow fluctuations are a hallmark of spontaneous cortical activity. We examine the hypothesis that these dynamics arise from recurrent neuronal networks operating near a phase-transition point, a state marked by “critical slowing down”. In such networks, a subtle shift toward criticality should selectively amplify slow fluctuations, providing a lever that can switch the cortex from quiet rest into self-generated behavior. Using a simple random recurrent network, we reproduce this amplification effect. The resulting spectra closely match intracranial electroencephalography from human visual cortex recorded during rest and during category-specific visual free recall. In particular, the model captures the experimentally observed enhancement of slow fluctuations during recall. These simulations reveal a parsimonious mechanism that explains spontaneous ultra-slow activity and enables rapid transitions between spontaneous states, suggesting that dynamic tuning toward criticality may be a general strategy by which cortical networks enter a generative mode.

Our brains never stand completely still: even at rest, neural activity drifts in very slow waves that last seconds to minutes. Where do these sluggish rhythms come from, and how do they help us shift from quiet rest into free, creative thought? We explored these questions using a minimalist computational model. The model is a network of simple “neurons” whose connections can be strengthened or weakened by a single gain knob. When the gain is set just below a critical tipping point, the network exhibits “critical slowing down,” creating ultra-slow fluctuations like those measured in real brains. Nudging the gain only slightly higher selectively boosts the slowest waves without destabilizing the network. This two-step behavior reproduces human intracranial recordings: the power at slow frequencies rises markedly in visual cortex when people freely recall images compared with passive rest. Our results suggest that cortical circuits normally hover near criticality and can edge even closer on demand, using critical slowing to shift rapidly from baseline activity to a generative mode that supports recall and other spontaneous, free behavior.

Figures

The impact of gain change on the dynamics of a recurrent system.The figure provides an illustrative and simplifying explanation of the critical slowing down (CSD) phenomenon in the dynamics of a recurrent system. CSD takes place as the system crosses between its sub-critical domain (G < 1, here depicted in blue) and its super-critical domain (G > 1, here in green). In this near-critical state (G ≲ 1, here depicted in yellow), when the system’s gain (also termed ‘control’) parameter is increased to the vicinity of (but still below) the critical point, it displays cardinal slowing down and amplification of its dynamics. Panel a delineates the expected activity in a random recurrent network at varying proximity to its critical point. For a transient stimulus, activity in the network will (a) dissipate in the sub-critical domain (G < 1, in the bottom sub-panel) (b) diverge to the maximal neuronal rate (black dashed line) in the super-critical domain (G > 1, top sub-panel) and (c) tend to become persistent as the gain approaches G ≲ 1 (middle sub-panel). For a continuous signal, the power spectrum will indicate rise in the amplitude of slow-wave activity only at the domain of G ≲ 1. Panels b and c provide further illustration of the CSD phenomenon exploring the dynamics of a simple one-dimensional non-linear system, driven by the quadratic equation of f(x)= x2+r. Panel b depicts the system and its solutions as the value of r, the control parameter, is modified between super-critical (r > 0, in green), sub-critical (r < 0, blue) and near-critical (yellow, r ≲ 0) states. Panel c shows how the system dynamics slows down as r approaches 0. The curve corresponding to r > 0 diverges exponentially. See a more detailed discussion in the text.

Ultra-slow spontaneous oscillations of neuronal activity in the brain and in a recurrent network simulation.Empirical iEEG results during wakeful-rest (adapted from [7]), presented here in panels a and b, are compared and contrasted with corresponding fitted simulation results in panels d and e, respectively, from an artificial neural network model. Panel a presents the HFB signal from homological left and right auditory cortices. Panel b presents the autocorrelation (black line) and the cross-correlation (orange line) of the HFB PSD obtained across the two auditory cortex hemispheres. Panel c presents simulation analysis results placing computational bounds on expected values of the control parameter G, as it grows from the sub-critical to the critical range. The upper part of the panel shows the simulation results for maximum (blue line) and average (green line) pair correlation, placing a bound (red interval) on the G (at the range 0.91 – 0.96) for the simulated correlations to match empirical results in panel a. The lower part of panel c shows the log ratio of simulated auto- to cross-correlation mean power at low frequency (orange line) indicating best match to empirical results at G* ≈0.95~ (red cross). Panel d demonstrates a corresponding signal from two highly correlated units in the artificial network at G*. Panel e depicts, for comparison, our simulation results of the activity fluctuations’ mean auto- and cross-correlation at G*. The location of the knee frequency in the autocorrelation spectra (at ~1 Hz) of the empirical and simulated data is further compared in Appendix B, Fig B in S1 File.

Amplification of slow-fluctuations in the brain and in the recurrent network.iEEG HFB signal from targeted free recall target (red) and non-targeted (green) neuron populations (adapted from [13]) in panel a, is compared with the simulated spontaneous activity in “recall” target (with additive noise) vs non-target “rest” blocks in panel b. The power spectra obtained from all intra-cranially recorded free-recall target sites (red) relative to the non-targeted population (blue) is shown in panel c. Simulated results are overlaid on the empirical PSD results in panel d, showing near fit, in expectation, over multiple realizations. The mean and standard deviation are shown for the free-recall target (yellow) and non-target (green).

Effect of network connectivity strength and additive noise on the spontaneous fluctuations’ PSD.Panels a and b demonstrate the CSD of low-range frequencies in an arbitrary network setting as the control parameter approaches the value of 1 (colors indicate different values of G). In panel a, the cases of 0-connectivity (i.e., γ = 0) and of G* are included (blue and green bold, respectively). In the inset (panel b), the gradual verge-of-criticality slowing down is highlighted. Panel c shows a similar amplification of slow frequencies, here as a function of the amount of additive noise (see Methods Iadd) overlaid on PSD results from [13]. Importantly, the network is set at the vicinity of criticality (G*). Note that the different scales used on the ‘Power’ axis in the different panels are aligned with the relevant references (panels a and b are scaled to match [30], whereas panel c follows the scaling used in [13]).

Impact of changes in input noise spectral profile on the network PSD.Panels a-b demonstrate the simulated near-critical network PSD when the injected noise is high-pass filtered (panel a) relative to the initial unfiltered white noise (panel b). Insets c – d present the noise spectra in use for each of the main panels, respectively.

The impact of network size and network connectivity strength on the amplification of the slow spontaneous activity fluctuations.The heatmap illustrates that amplification of normalized low-frequency power is driven mainly by the network’s proximity to criticality, as quantified by the control parameter, G.

Equations

Full text

Introduction

Ultra-slow spontaneous activity fluctuations, also known as “resting-state” activity [1], have been consistently observed across all studied cortical regions [2–4]. The power spectra of these fluctuations have been characterized as obeying a power-law profile [5–7] and they emerge in well-defined networks across the cortex, uncovering unique patterns of functional connectivity [3,8,9]. These spontaneous fluctuations have also been linked to ramping anticipatory signals, as well as other spontaneous behavioral manifestations, such as pupil diameter fluctuations [10] and eye-drifts [11]. We previously proposed [12,13] and recently demonstrated, a direct link between the spontaneous activity fluctuations and free and creative verbal behaviors [14,15]. However, it is still not clear how the characteristic power spectrum of the spontaneous fluctuations - dominated by ultra-slow frequencies - emerges out of the typically fast neuronal elements. Several non-recurrent mechanisms were proposed, including protracted synaptic currents that may integrate inputs over slow timescales [16], and power-law spike-frequency adaptation introducing slow dynamics [17]. We and others have previously suggested that these slow dynamics may emerge from recurrent neuronal connections [12,18–22].

Near-critical dynamics as a potential mechanism. Here, we extend this notion and examine the hypothesis that the cortex operates in the vicinity of a critical transition point [23–26]. In such a system, the slow resting state fluctuations can be explained in the context of a fundamental property of near-critical recurrent systems - the slowing down of stochastic fluctuations as these systems approach criticality. Termed “critical slowing down” [CSD; 27–29], the phenomenon describes the universal tendency of systems near a transition point in their dynamics to exhibit slower time scales and take longer time to relax to a steady state equilibrium following perturbations. Hence, the slowing down of spontaneous fluctuations near criticality may explain their unique power spectra as observed in the resting state networks. Indeed, some prior modeling work has tentatively proposed such a link [30,31].

Recent experimental and theoretical evidence supports the idea that the healthy cortex operates near a critical point, exhibiting slow time scales and scale-free dynamics [31–37]. Scale-free spatiotemporal cascades of neuronal activity, termed neuronal avalanches, have been found in in-vitro and in-vivo animal studies as well as in large-scale human fMRI, EEG and MEG recordings [32,38–40]. These avalanches can be well described by the framework of critical branching processes [41] and they can coexist with oscillations [42]. Other studies indicate edge-of-chaos criticality, in which the system hovers not on the brink of explosive growth but rather in a regime of richly chaotic yet bounded dynamics [43–45]. Another type of critical transitions discussed in the literature is pattern formation criticality, which refers to a transition into self-organized spatial patterns of activity, such as stationary bumps or traveling waves [22,46–48].

Importantly, near-critical dynamics optimize information processing [26,49–51] and deviations from critical dynamics are associated with disorders in information processing [52–55]. We have previously demonstrated the phenomenon of CSD in a recurrent network model of early visual processing, where the effective time scales grow longer as the network approaches a critical point during its learning process [22]. For a recent review of critical brain dynamics, see [23,56].

A straightforward prediction of the near-critical dynamics is that bringing the system closer to criticality should result in biasing the power spectra of the spontaneous fluctuations towards even slower frequencies. Indeed, recent experimental findings derived from direct recordings of brain activity in patients are compatible with this notion. Thus, in a study exploring brain mechanisms of visual recall, Norman et al. [13] found that resting state fluctuations displayed an increase in baseline neuronal firing rate and power of the slow spontaneous fluctuations during free recall. This local baseline shift - i.e., increase in average firing rate - was predicted to bring the targeted neurons closer to the decision threshold and enhance the probability of free recall during such periods.

It is important to emphasize that this amplification of the fluctuation amplitude can serve a significant cognitive and functional role. Free behavior relies on the ability to rapidly and flexibly switch between active and rest conditions. For instance, musicians can easily switch into free improvisation mode from a state of rest. If the slow spontaneous fluctuations indeed underlie free behaviors, a mechanism must be present in the brain to regulate these fluctuations. In the study of Norman et al [13], exploring this idea in the domain of free visual imagery and recall, it was found that resting state fluctuations displayed an increase in baseline and power of the slow spontaneous fluctuations specifically when subjects switched from rest to an active recall state. Here, we propose that this local dynamical shift, which increases the probability of targeted neurons to cross the decision threshold, can be elegantly explained by the CSD phenomenon.

Specifically, we propose that the related amplification of the spontaneous fluctuations may be due to bringing the network closer to its critical point, thereby leading to CSD. This provides a mechanism for the efficient implementation of a state-change from rest to active, yet spontaneous, behaviors. We further speculate that the proximity to criticality is controlled through neuromodulation, allowing dynamic regulation of cortical states

Network model and tuning of proximity to criticality. To explore the suggested mechanism, we simulated a random recurrent network of rate-based units at varying proximities to criticality. These rate models can be rigorously derived from large-scale conductance-based networks that model cortical dynamics [57]. Our network architecture is adapted from the one proposed by Chaudhuri et al. [30] to simulate intra-cranial EEG (iEEG) dynamics. Proximity to criticality is controlled by tuning relevant biophysically meaningful parameters. The precise tuning is determined by analyzing the recurrent interaction matrix using considerations from random matrix theory assuming a large enough network size, typically in the order of a few hundred neurons [58–60]. Near-criticality is required to reproduce the observed scale-free power spectrum of iEEG recordings [30]. Setting the model closer to its critical point may account for the unique dynamics observed during free recall [13].

To provide further insight, Fig 1 schematically illustrates CSD in a random recurrent neural network and in a one-dimensional non-linear system near its critical point. Throughout this work, we refer to “criticality” in the following dynamical-systems sense: in the subcritical domain, perturbations decay and spontaneous activity remains bounded, whereas in the super-critical regime even small fluctuations grow exponentially, driving the mean firing rate toward saturation. This mean-field, non-oscillatory bifurcation yields avalanche-like dynamics rather than the edge-of-chaos or pattern formation transitions discussed above. Panel 1a shows the expected activity in a random recurrent network at varying proximities to its critical point, as set by the network gain parameter, $[eqn]$ . The critical point is defined when $[eqn]$ (see Methods for further mathematical substantiation). For a transient stimulus, activity will dissipate when the network gain is set to $[eqn]$ (blue arrow), diverge to the maximal neuronal rate (in black dashed line) when $[eqn]$ (green arrow), and tend to persist at enhanced slow fluctuation frequencies when $[eqn]$ , namely slightly below the critical point (yellow arrow). For a continuous signal, the steady-state power spectrum indicates an increase in the amplitude of slow activity only for $[eqn]$ .

The impact of gain change on the dynamics of a recurrent system.The figure provides an illustrative and simplifying explanation of the critical slowing down (CSD) phenomenon in the dynamics of a recurrent system. CSD takes place as the system crosses between its sub-critical domain (G < 1, here depicted in blue) and its super-critical domain (G > 1, here in green). In this near-critical state (G ≲ 1, here depicted in yellow), when the system’s gain (also termed ‘control’) parameter is increased to the vicinity of (but still below) the critical point, it displays cardinal slowing down and amplification of its dynamics. Panel a delineates the expected activity in a random recurrent network at varying proximity to its critical point. For a transient stimulus, activity in the network will (a) dissipate in the sub-critical domain (G < 1, in the bottom sub-panel) (b) diverge to the maximal neuronal rate (black dashed line) in the super-critical domain (G > 1, top sub-panel) and (c) tend to become persistent as the gain approaches G ≲ 1 (middle sub-panel). For a continuous signal, the power spectrum will indicate rise in the amplitude of slow-wave activity only at the domain of G ≲ 1. Panels b and c provide further illustration of the CSD phenomenon exploring the dynamics of a simple one-dimensional non-linear system, driven by the quadratic equation of f(x)= x2+r. Panel b depicts the system and its solutions as the value of r, the control parameter, is modified between super-critical (r > 0, in green), sub-critical (r < 0, blue) and near-critical (yellow, r ≲ 0) states. Panel c shows how the system dynamics slows down as r approaches 0. The curve corresponding to r > 0 diverges exponentially. See a more detailed discussion in the text.

Panels 1b & 1c offer a didactic illustration of CSD in a simple one-dimensional system. The dynamics are characterized by a simple quadratic equation: $[eqn]$ , where $[eqn]$ is a time-dependent variable satisfying $[eqn]$ , and $[eqn]$ is the system’s control parameter. For $[eqn]$ the system is super-critical. At this domain, there are no fixed points and $[eqn]$ always diverges to infinity. At $[eqn]$ the system is sub-critical and two fixed points arise at $[eqn]$ (panel 1b; filled circles correspond to stable fixed points and hollow circles to unstable fixed points). Close to a given fixed point, the dynamics can be linearized and the deviation from the fixed point can be approximated by an exponential temporal profile, $[eqn]$ , where $[eqn]$ (note that for stable fixed-points the derivative is negative and the dynamics converge). The magnitude of the time constant, $[eqn]$ , that controls the dynamics approaches infinity as the system approaches its critical point (where $[eqn]$ ). Thus, the dynamics slow down as the system approaches the critical point, as illustrated in panel 1c. Furthermore, at the bifurcation, the linear term vanishes, and the dynamics become dominated by higher-order terms, leading to power-law scaling behavior. This transition to power-law dynamics is a hallmark of critical systems, reflecting the absence of a characteristic timescale and the emergence of scale-invariant fluctuations

In summary, we propose that the unique dynamics observed by Norman et al. [13] may be accounted for by a closer proximity of the local cortical network to a critical point during free recall. This enhances the gain of the local network, making it more sensitive to noise and driving it closer to the decision bound. We propose a model of a randomly connected recurrent network subject to a white noise input for this putative mechanism. Our key hypothesis is that this amplification in the slow activity amplitude stems from the CSD phenomenon, and can be obtained by a relatively small shift towards criticality. It is likely that any model that manifests such critical phase transition behavior and power law spectra would be equivalently sufficient. Our findings support the hypothesized slowing down modulation near-criticality, showing good fit with the empirical data. They support the hypothesis that the phenomenon of ultra-slow, scale-free resting state fluctuations and their modulation during intrinsically motivated tasks reflect a basic and ubiquitous behavior common to all suitably configured recurrent systems.

Results

To model the iEEG signal, we simulated a sparse random recurrent network, realizing a firing rate model (Methods). The summed activity of a sample of units represented the corresponding iEEG [30]. The proximity to the critical transition between stable and unstable dynamics is quantified by the control parameter

$[eqn]$ ,

where $[eqn]$ is the neural gain, $[eqn]$ is network sparseness and $[eqn]$ the mean synaptic strength. This relationship is obtained using tools from random-matrix theory to evaluate the spectrum of eigenvalues of the matrix that governs the linearized network dynamics [Methods; 30,59]. At $[eqn]$ , the dominant eigenvalue crosses from negative to positive, leading to a transition into runaway (supercritical) activity. To satisfy the large-network assumption underlying this analysis, the model size was set to $[eqn]$ units (Methods; Appendix B, Fig A in S1 File).

Throughout the main simulations $[eqn]$ was fixed at 0.2 and $[eqn]$ at 49.881 $[eqn]$ (following [30]), while $[eqn]$ was varied to move the system toward or away from the critical point. The time constant ( $[eqn]$ ) in these simulations was set to 20 ms, which is close to well-known estimates for cortical neurons’ membrane time constant [61–63]. A sampling parameter, $[eqn]$ , defined the fraction of network nodes summed from the model to produce the simulated iEEG, and background drive was provided by a white-noise input, $[eqn]$ (Methods).

To study the model dynamics, we also used an additive white-noise current, $[eqn]$ . This was motivated by $[eqn]$ . Tuning just these two parameters, independently, the neural gain, $[eqn]$ , and the additive white-noise, $[eqn]$ , enabled us to match the slow-frequency power spectrum and reproduce the observed shift from rest to free-recall (see Methods and Fig 3).

Ultra-slow spontaneous oscillations of neuronal activity in the brain and in a recurrent network simulation.Empirical iEEG results during wakeful-rest (adapted from [7]), presented here in panels a and b, are compared and contrasted with corresponding fitted simulation results in panels d and e, respectively, from an artificial neural network model. Panel a presents the HFB signal from homological left and right auditory cortices. Panel b presents the autocorrelation (black line) and the cross-correlation (orange line) of the HFB PSD obtained across the two auditory cortex hemispheres. Panel c presents simulation analysis results placing computational bounds on expected values of the control parameter G, as it grows from the sub-critical to the critical range. The upper part of the panel shows the simulation results for maximum (blue line) and average (green line) pair correlation, placing a bound (red interval) on the G (at the range 0.91 – 0.96) for the simulated correlations to match empirical results in panel a. The lower part of panel c shows the log ratio of simulated auto- to cross-correlation mean power at low frequency (orange line) indicating best match to empirical results at G ≈0.95~ (red cross). Panel d demonstrates a corresponding signal from two highly correlated units in the artificial network at G*. Panel e depicts, for comparison, our simulation results of the activity fluctuations’ mean auto- and cross-correlation at G*. The location of the knee frequency in the autocorrelation spectra (at ~1 Hz) of the empirical and simulated data is further compared in Appendix B, Fig B in S1 File.*

The analysis unfolds in two stages. First, the model is calibrated to the resting-state HFB spectra measured by [7], which establishes quantitative bounds on how close the network can lie to the critical point. Working within that window, we then apply the same model to the free-recall data of [13]. We demonstrate that a modest gain increase, well inside the limits set by the data of [7], reproduces the baseline amplification and ultra-slow power boost observed in target populations. This progression from rest to task links the two empirical datasets highlighting how proximity to criticality can enable rapid, goal-directed reconfiguration of network dynamics.

To that end, we first examined whether our network model could reproduce the dynamical characteristics of the spontaneous resting state activity fluctuations as recorded in the human cortex, and found that it did when set close to its critical point. This is demonstrated in Fig 2, which compares real iEEG recordings and simulated network dynamics. It has been previously shown [64–66] that the power of the LFP in the high frequency broadband (HFB) provide a reliable index of the average firing rates of neurons within the recording site. Note that this HFB measure should be distinguished from the direct recordings of the LFP signal itself. Panel 2a shows high correlation levels (r = 0.56 Pearson) in spontaneous firing rate fluctuations based on HFB power from two auditory cortex electrical contacts recorded across hemispheres during rest (obtained from [7]). The autocorrelation (black) and cross-correlation (orange) PSD (often referred simply as PSD and cross-PSD) of the HFB obtained from multiple auditory cortex sites across the two hemispheres during rest in Nir et al are shown in panel 2b.

To investigate the influence of the proximity to criticality on the typical correlation levels between pairs of units in the network, we simulated multiple network realizations, gradually increasing the control parameter in each realization from the subcritical ( $[eqn]$ = 0.8) to the near-critical ( $[eqn]$ = 0.99) range. The upper part of Panel 2c presents the simulation’s maximum and average pair correlation in blue and green, respectively. Results indicate that the network first reaches the correlation level observed by Nir et al only at $[eqn]$ = 0.91 (blue curve), and that this becomes the mean correlation level at $[eqn]$ = 0.96 (green curve), setting a bound on the operating range of $[eqn]$ (see red interval).

Further “fine-tuning” of the actual value of the control parameter $[eqn]$ was established by computing the ratio of the cross- and auto-correlation power at low frequencies as the system approaches criticality. Previous work demonstrated that spatial cross-correlations arise more markedly than temporal autocorrelation during CSD [29,67]. The lower part of panel 2c shows the simulation results of this computation, in the range $[eqn]$ = 0.8 to 0.99, over multiple realizations of the model network. Our findings show that the autocorrelation is less influenced by the growth of $[eqn]$ relative to the cross-correlation, and thus, their exponent ratio changes rather steeply as a function of $[eqn]$ (orange line), providing the approximation $[eqn]$ * $[eqn]$ (see red cross) for the log ratio observed in the empirical findings (0.84). Panel 2d demonstrates correlated activity fluctuations of two simulated units in the recurrent network (corresponding to Fig 2a), whereas panel 2e depicts the simulated auto- and cross-correlation PSD results for this estimated value $[eqn]$ . For later reference, and given the inherent difficulty in determining an exact value for $[eqn]$ , we define a typical resting-state range, $[eqn]$ , within which $[eqn]$ is expected to lie. Note that results for other values of τ were explored to optimize the fit of the “knee” bend in the PSD curve as presented in Appendix B, Fig B in S1 File, showing that a τ value of 20 ms provides better fit than alternative settings at 2 or 200 ms. Note also that in the very near vicinity of the critical point ( $[eqn]$ ) the slow mode takes hold of neuronal correlations in the network, leading to mean correlation $[eqn]$ .

Our simulation focused on the low-frequency range, where it faithfully reproduced the auto- to cross-correlation power ratio found in Nir et al. (as demonstrated in panel 2c). We did not attempt to optimize the mid- to high-frequency fit (panel 2e relative to 2b). At present, the model captures autocorrelation power more accurately than cross-correlation; refining this aspect remains for future work. The sole purpose of the simulation based on Nir et al. was to bound the admissible values of the global gain parameter, $[eqn]$ , establishing the range later used to model the rest-vs-recall shifts observed by Norman et al. [13].

Next, we examined whether the recurrent network simulation using the same network size (240 units) and a similar parameter setting (as presented in Fig 2) could also fit the baseline shift and amplitude amplification of the slow spontaneous fluctuations observed during targeted free recall behavior, as reported in [13]. In this regard, the parameters $[eqn]$ , $[eqn]$ and $[eqn]$ remained unchanged from the previous simulation.

Fig 3 compares the intra-cranially recorded HFB baseline shift amplification with the simulated neuronal activity. To that end, few remaining parameters of the simulated network were again fitted for optimal power spectrum match with results obtained by the intracranial HFB signal PSD. As for the setting of the decay constant, in addition to the optimal knee-location fit obtained in Fig 2 (indicating τ = 20 ms), another approximation method was employed, a Lorentzian approach, following [30]. The result of this approximation is shown in Appendix B, Fig C in S1 File. Thus, we establish once again, now based on two distinct methods, that setting this parameter near the biophysically-grounded-estimate of the membrane time constant optimizes the fit of the knee frequency location, and here also of the PSD curve for mid-range frequencies (> 0.1 Hz and < 10 Hz). To adjust the curve slope in the higher frequencies (> 10 Hz), the sampling fraction was set to α = 0.01. The simulation was run for 1800s, with the first 600s omitted as a stabilization period (by means of signal stationarity testing) before analysis of the remaining 1200s.

Panel 3a depicts an example of the experimentally derived spontaneous HFB fluctuations in two neuronal populations - one that was targeted to elicit free recall (red) and the other that was not targeted (green). The baseline shift of the fluctuations for the targeted area is readily visible. It should be noted that in [13], a generalized enhancement of slow fluctuations was found for the free recall conditions; however, the study did not find clear evidence for a selective enhancement in specifically targeted vs non targeted categories, likely due to the mixed category selectivity and low SNR of the clinical recording sites. By contrast, in the model simulation, we assume a fully selective neuronal group - and the simulation results should be taken as a prediction of such category-selective effects in future studies. For comparison, Panel 3b shows the simulated activity fluctuations for a network when a uniform baseline (fixed value) shift, $[eqn]$ (see Methods), was added to the white noise inputs (red) and without this added fixed value shift (green). Note that the gap between the two curves can easily be fine-tuned via setting the relative value and ratio of $[eqn]$ and $[eqn]$ , without compromising other measurable properties. Importantly, similar amplification in the baseline firing rate was observed also by increasing the network gain, γ, towards the critical point, implying a minor 1% change of the control parameter $[eqn]$ .

Panel 3c provides a more comprehensive analysis by showing the mean PSD across all intra-cranially recorded sites in the targeted neuronal population (red PSD) vs. the mean PSD across all sites in the non-targeted population (blue PSD). To make the comparison of these neuronal recordings to the simulations easier, we overlay the simulated results on top of the brain recordings in panel 3d. In this panel, the mean and standard deviation for multiple realizations (specifically, 8 arbitrary random-seed network connectivity initializations) of simulated target (yellow) vs non-target (green) blocks are presented. Here, the fit in PSD was reached by adding a fixed value shift $[eqn]$ of 12.5 pA to the white noise during simulated targeted (orange thick curve) but not to the non-targeted blocks (green thick curve). The control parameter was set within the $[eqn]$ range, as in Fig 2. As in [13], the increase in the power of the ultra-slow modulation in the simulated PSD amplitude during recall relative to rest blocks over the different realizations was statistically significant (for the range < 0.2 Hz; p < 0.05, FDR corrected, signed-rank test). Goodness of fit analysis indicated that the simulated network (mean PSD, for both rest and recall conditions) captures a relevant segment of the empirical spectra (in the range 0.05 - 10 Hz) rather well (log-domain $[eqn]$ = 0.95); after correcting for the autocorrelation between neighboring frequency bins ( $[eqn]$ = 14 vs. 256 raw points; see Methods), the corresponding two-tailed t-test indicated significance ( $[eqn]$ , $[eqn]$ ).

Importantly, a comparable amplification of the targeted ultra-slow frequency PSD was also achieved through direct modulation of the synaptic connectivity matrix. Notably, increasing the network gain, γ, by 1% within the scope of $[eqn]$ * was sufficient to replicate the empirically observed enhancement of slow frequencies during free recall. In Appendix B, Fig D in S1 File, we show a distinct single realization PSD result from these simulations using the same color scheme mentioned above. Note the similarity in modulations of the PSD, and, specifically, the amplification in the ultra-slow-frequency domain of the PSD, evident in the targeted relative to non-targeted simulated blocks.

In the following subsections, we systematically examine the phenomenon of CSD in our recurrent random network model, focusing on its fundamental properties rather than fitting specific data. Accordingly, some parameter settings are adjusted freely to enhance clarity in the illustrations.

To examine the impact of network connectivity on the PSD over a broader modeling range and greater detail we have systematically examined the influence of changing network and connectivity parameters over a wide range starting from well below the dynamical critical point. This is depicted in Fig 4. In panels 4a-b, a simplified network setting with 100 units was used to explore the influence of the control parameter $[eqn]$ on the spectra, via gradual wide range adjustments of γ, while keeping all other parameters fixed. The simulation lasted 800s, with the first 200s omitted in the analysis to allow signal stabilization; τ was set again to 20 ms, and α was set to 1.0.

Panel 4a shows the effect of setting $[eqn]$ over its full range from 0.0 to 1.0, whereas panel 4b focuses on the narrow range near the critical point. Three key findings emerge from this analysis: first, a gradual increase can be observed in the amplitude of slow frequencies that accelerates near the critical point of $[eqn]$ (when $[eqn]$ ); second, the medium to high frequency range of the PSD remains unchanged even when all network connections are set to zero, effectively isolating the neurons from each other; and third, the slow-frequency PSD slope for $[eqn]$ * (green line), which is, interestingly, close to the scale-free exponent of negative 2 [5]. In panel 4c, the PSD results for the same simulation as in Fig 3 are shown, with the control parameter in vicinity of $[eqn]$ *, but with gradually increasing levels of baseline shift of the noisy input, $[eqn]$ . Similarly to the modulation of connectivity strength towards criticality in $[eqn]$ , when the network is already near-critical, additive noise leads, as well, to a slowdown in the ultra-slow frequencies. To further investigate, we compared the influence of additive noise on the slowing down phenomenon in networks at varying proximities to the critical point. Fig E in Appendix B in S1 File shows that additive noise significantly increases the amplitude of low frequencies when the network is near the critical point but has practically no effect when the network is far below the critical point.

The emergence of a scale-free power-law in the fluctuations of simulated activity (even in individual isolated neurons, for mid-to high frequencies) raises the concern that the network behavior merely reflected the dynamical structure of the injected input noise. To examine this possibility, we modified the injected noise profile by high-pass filtering and observed the effect on the network PSD for parameters near criticality. The simulation was 800s long, with the first 200s excluded for signal stabilization. The time constant τ was set here to 20 ms, as before. The control parameter, $[eqn]$ , was set in vicinity to criticality and α was set to 1.0. Fig 5 displays the simulation results. We also used the same noise profiles to simulate isolated neurons (zero-connectivity network), with the resulting PSD shown in Appendix B, Fig F in S1 File.

Amplification of slow-fluctuations in the brain and in the recurrent network.iEEG HFB signal from targeted free recall target (red) and non-targeted (green) neuron populations (adapted from [13]) in panel a, is compared with the simulated spontaneous activity in “recall” target (with additive noise) vs non-target “rest” blocks in panel b. The power spectra obtained from all intra-cranially recorded free-recall target sites (red) relative to the non-targeted population (blue) is shown in panel c. Simulated results are overlaid on the empirical PSD results in panel d, showing near fit, in expectation, over multiple realizations. The mean and standard deviation are shown for the free-recall target (yellow) and non-target (green).

Effect of network connectivity strength and additive noise on the spontaneous fluctuations’ PSD.Panels a and b demonstrate the CSD of low-range frequencies in an arbitrary network setting as the control parameter approaches the value of 1 (colors indicate different values of G). In panel a, the cases of 0-connectivity (i.e., γ = 0) and of G are included (blue and green bold, respectively). In the inset (panel b), the gradual verge-of-criticality slowing down is highlighted. Panel c shows a similar amplification of slow frequencies, here as a function of the amount of additive noise (see Methods Iadd) overlaid on PSD results from [13]. Importantly, the network is set at the vicinity of criticality (G*). Note that the different scales used on the ‘Power’ axis in the different panels are aligned with the relevant references (panels a and b are scaled to match [30], whereas panel c follows the scaling used in [13]).*

Impact of changes in input noise spectral profile on the network PSD.Panels a-b demonstrate the simulated near-critical network PSD when the injected noise is high-pass filtered (panel a) relative to the initial unfiltered white noise (panel b). Insets c – d present the noise spectra in use for each of the main panels, respectively.

As can be seen, even drastic changes in the spectral power profile of the input noise had only minor impact on the PSD low frequency power-law profile of the network, leaving the mid-to-high profile nearly unchanged. This is evident in the striking discrepancy between the PSD profile of the input noise (Fig 5 insets c-d) and the PSD of the entire network (panels 5a-b). Note a bump at ~3 Hz, a consequence of setting the high-pass filter’s cutoff at that frequency. Of particular interest is the insensitivity of the network even when all slow input frequencies were drastically quenched to a tenth of their original power (see panel 5a relative to 5c). Similar behavior for the PSD of mid-to-high range frequencies is preserved also in the isolated neuron network (Appendix B, Fig E in S1 File).

Finally, to analyze the influence of network size on the amplification of the spontaneous fluctuations’ slow frequencies as the network approaches criticality, we scanned the parameter space of both network connectivity and size. In the analysis, we iteratively generated network realizations of larger and larger scale (150–350 units), and increased the control parameter G (via γ) from subcritical to near-critical range, while collecting the slow (normalized by fast) frequency power values.

The results are shown as a two-dimensional heat map (Fig 6). The x and y axes in this map depict network size and control parameter settings, respectively, whereas the color per pixel denotes the degree of slow frequency amplification (more details in Methods). The heatmap clearly demonstrates that modulation of the control parameter towards criticality is the dominant factor driving the slow-frequency power amplification effect. A mild trend for larger amplification in slow frequency power is shown for network size > 200. This trend appears to become negligible as networks grow larger.

The impact of network size and network connectivity strength on the amplification of the slow spontaneous activity fluctuations.The heatmap illustrates that amplification of normalized low-frequency power is driven mainly by the network’s proximity to criticality, as quantified by the control parameter, G.

Discussion

Ultra-slow spontaneous fluctuations are one of the most widely studied phenomena of the human brain. Numerous studies have revealed that these fluctuations are organized along highly complex but informative functional connectivity patterns [3,68–70] suggested to be related to cognitive habits, traits and pathologies [71–73]. It is also proposed that they drive self-generated, voluntary, and creative behaviors [13–15,74].

Despite extensive research on the properties and dynamics of these slow fluctuations [1,3,7,75], the underlying neuronal mechanism remains unknown. The slow dynamics and intricate organization of these fluctuations may suggest a complex underlying mechanism, but we argue that the opposite is true. We hypothesize that critical slowing down (CSD), a fundamental principle in the dynamics of recurrent networks, underlies the behaviors of the ultra-slow spontaneous fluctuations during both rest and task conditions.

Using a simple random recurrent network near criticality and unstructured white noise inputs, we successfully simulated the CSD dynamics and reproduced the power spectrum of the spontaneous (resting state) activity fluctuations as reflected in the intracranially measured HFB. Furthermore, a small uniform increase in synaptic efficacy, or alternatively, a small DC shift in this noise input, near the critical transition point of the network dynamics, were sufficient to simulate the amplitude enhancement documented using intra-cranial (iEEG) recordings during targeted free recall in the visual cortex. The idea that changes in the efficacy of recurrent interactions can modulate local intrinsic timescales was also explored in the context of selective attention [76]. Specifically, their computational modeling work proposes that increased intrinsic timescales during selective attention are associated with a shift of network dynamics closer to a critical point.

Recently, modeling of stochastic noise coupled with slow synapses in a recurrent spiking network successfully simulated the slow buildup anticipating voluntary movements [77]. Wider-scale whole-brain modelling indicates that the richness of resting-state dynamics may arise from a delicate balance between neuromodulatory gain and stochastic drive. In a primate connectome model, small amounts of noise let the network roam a repertoire of metastable configurations beyond those accessible in a purely deterministic regime [78]. A later study demonstrated that the same gain-noise interplay breaks symmetry and carves out a low-dimensional manifold matching empirical co-activation patterns [79]. Our minimal network formalizes this mechanism by demonstrating analytically how tuning a single gain parameter near the critical point both elevates noise sensitivity and stretches time scales into the ultra-slow domain. Complementary studies reveal that nudging gain toward criticality [31] or allowing inhibitory plasticity to self-tune this balance [80] likewise promotes realistic, noise-driven spontaneous activity, underscoring the generality of the framework we provide.

It is important to clarify that our model was chosen not by virtue of its likely optimal fit to the experimental data - but rather because it was the most basic model of a recurrent network we could conceive of. The model was previously shown to account for the 1/f envelope of the power spectrum of the raw iEEG data [30]. Furthermore, a major advantage of the simplicity of this model is that it allows for explicit control of the proximity to the critical point. By controlling the network gain, we can move a sub-critical model to a model that is near the critical threshold and explore the changes in the dynamics. In this respect, we effectively employ a spectrum of models, namely various random recurrent networks, of different levels of connectivity and size that generate iEEG like activity, at varying proximity to the critical transition. The crucial point is that only under parameters that enable CSD, namely near criticality, were we able to replicate the experimental data, and specifically the amplification in ultra-slow fluctuations.

Based on these simulations, we propose that the primary driver of slow spontaneous fluctuations in neuronal activity is stochastic, likely, but not necessarily, balanced noise that critically accumulates through a highly specific, recurrent, connectivity structure. The exact source of this noise is unknown and could originate from various processes. One possibility involves intrinsic neuronal mechanisms, such as thermal, synaptic, and ion channel fluctuations [81,82]. These processes could collectively contribute to the random membrane potential fluctuations consistently observed across all types of cortical neurons [1,82,83]. Another possibility is that the power spectrum of incoming spiking activity from subcortical sites exhibits characteristics of white noise. A recent study on activity in the subthalamic nucleus demonstrated that while the spectrum of synaptic activity follows a power law with an exponent around -2, the spectrum of the resulting spike rate is nearly flat [84]. It is important to clarify, that our model equations capture the dynamics of synaptic activity, which is the primary contributor to iEEG signals [85].

In this work, we simulated two key features of spontaneous fluctuations: their power spectra on the one hand and modulation during an intrinsically motivated task, on the other. However, we did not address how complex functional connectivity patterns may emerge from such noise fluctuations. We have previously suggested that these patterns are shaped by synaptic structures and reflect prior learning [73,86]. This hypothesis predicts that a recurrent network with updatable connectivity based on suitable learning rules will spontaneously reproduce a repeatedly activated specific pattern when network neurons, having learned the pattern, are injected with white noise near criticality. While we focus on the power spectra and modulation during task-induced targeting of categorical boundaries, future simulations could investigate the emergence of such functional connectivity patterns.

Our simulation focused specifically on ultra-slow fluctuations and their targeted amplification. While other aspects of the [13] study are intriguing, they may require more sophisticated network models and additional analyses. The dynamics of the recurrent network when crossing the critical point, and its relationship with the concept of ‘ignition’ [12], is another fascinating area of study. However, these aspects are beyond the scope of our current work and are left for future exploration.

As noted above, our model does not attempt to capture the full extent of the task and the consequent behavior. Rather, here we focus on the slowing down and amplification of the overall activity fluctuations. Our assumption is that shifting the network closer to criticality or enhancing the external noise may transition the network into a generative mode, where stored activity patterns are more likely to cross a relevant threshold. The advantage of the model in our context, is that the proximity to criticality can be directly controlled, allowing for a systematic examination of its effects on the resulting power spectrum.

Earlier studies have demonstrated the emergence of a scale-free 1/f power-law spectrum using a rate model simulation near criticality [30,61]. However, they focused primarily on the tipping point of criticality ( $[eqn]$ ) of the network parameter setting. In contrast, we aimed to investigate the minimal conditions required for the scale-free PSD to arise, examining network connectivity weights far below the critical point, as well as the impact of the input noise profile on the network’s neuronal activity.

Intracranial recordings of single neurons and the high frequency broadband (HFB) amplitudes of the LFP, which has been shown to index the average firing rate of cortical neurons [64,66], reveal the power spectra of the spontaneous fluctuations in the human cerebral cortex during rest (panels 2c and 3c). These spectra consist of two power-law exponents, with slower frequencies following a shallower exponent than higher frequencies. Remarkably, a similar profile emerges in a random recurrent network with fully stochastic connection strengths and input noise under the appropriate parameters (compare to panels 2d and 3d).

These results were obtained using minimal, biophysically-oriented, parameter settings, consisting of a physiologically plausible membrane time constant (~20 ms), sampling of only a small part of the network nodes (α), which is compatible with the localized nature of the empirical recording sites, and control of the system’s slowing down by a single gain parameter ( $[eqn]$ ).

We first discuss the high frequency slope of the HFB PSD. Analyzing the behavior of the network under zero connectivity (panel 4a) revealed that the resultant fluctuations in membrane potential approximated a power-law beyond a certain “knee” frequency, even when the spectra of the stochastic noise injected into each neuron was flat (white noise) or biased towards higher frequencies (see Methods and Appendix B, Fig E in S1 File). Thus, each unit, when decoupled from the rest of the network, maintains the slow frequency band and attenuates the higher frequencies in a scale-free manner, following a Lorentzian power spectrum ([30] and Appendix B, Fig C in S1 File). A mathematical analysis presented in SI section 1 shows that the dynamics of a scale-free power spectrum with an exponent of negative 2 is indeed the expected result for the fast-mode HFB PSD. Recent studies of activity fluctuations in isolated cortical and hippocampal neurons also reveal that individual neurons are capable of generating stochastic fluctuations obeying scale-free power-law profiles as well as long time-scales, even in the absence of network connections [87–89].

Examining the impact of network connectivity on the spontaneous fluctuations reveals a specific amplification of slow frequencies’ amplitude in the spontaneous fluctuations spectra as the network approaches its dynamical phase transition point. This is shown when separately examining the impact of changing the connectivity strength, i.e., the efficacy of intra-neuronal communication weight, and the richness of connectivity as reflected in network size (Fig 6). A successful amplification of slow frequencies can only be achieved within a regime in which these two parameter spaces overlap. However, only random network configurations were explored in this study, and a better structured connectome arrangement, whether shaped by innate prearrangement or as learned by prior experiences, may present other near-critical behaviors. Nevertheless, it should be noted that even in extremely small, properly configured networks possessing fewer than 20 neurons, a powerful amplification of ultra-slow fluctuations can be obtained, suggesting that ultra-slow fluctuations lasting many seconds do not necessarily depend on large-scale cortical networks for their generation, as also suggested by previous work [90].

Our simulation results are compatible with the notion that the spontaneous fluctuations can be selectively amplified when behavior switches into a free recall state. Furthermore, they predict that when participants attempt to visualize different items belonging to a certain category of visual stimuli, their slow spontaneous fluctuations should be selectively amplified, specifically in the corresponding neuronal representations. This can happen through a modulation towards criticality via connectivity gain or via additive external noise (see panels 3a and 3c and more details in [13]). Such amplification enhances the probability of these fluctuations to cross the decision bound, specifically in the category-selective area. This modulation toward a critical threshold thus provides an elegant mechanism for rapid and flexible control over the expression of any spontaneous behavior within an intentionally selected categorical boundary. For example, when we attempt to freely imagine different types of faces, we rarely, if ever, mistakenly picture a building. As we demonstrate here, modulation of proximity to criticality and the resulting targeted CSD offer a plausible mechanism for such precise cognitive distinctions. However, it should be emphasized that such specific, category selective, enhancement of slow frequencies has not been demonstrated experimentally yet, likely due to low selectivity and poor SNR of the recordings available in [13] study. Thus, our simulation of this result should be taken as a prediction for future studies with a higher spatial resolution.

When local populations’ activity fluctuations transiently cross the decision bound, this may lead to rapid self-amplification termed “local ignition”. Such ignitions have been experimentally observed during crossing the perceptual awareness threshold using iEEG recordings [12,91]. Additionally, near a critical point, network sensitivity increases, enabling better encoding and recall of information through balanced response to inputs and adaptation to noise (M. [22,92]).

Our simulations of a random recurrent network demonstrate that a simple shift in noise level or in gain modulation accurately reproduces the intra-cranially recorded slow frequency amplification observed in the empirical data during category-selective visual recall (panel 3d and Appendix B, Fig D in S1 File, respectively). The key factor inducing the amplification effect is the network’s proximity to its critical transition point, regardless of the specific mechanism controlling its excitation.

It is important to note that empirical brain data and recurrent network simulations specifically relate to fluctuations in neuronal firing rate or, in cases where such measures are unavailable, to the broadband high frequencies (gamma) envelope of the LFP. Previous studies have shown that the gamma envelope reflects fluctuations in the average firing rates of local neuronal populations [64,66]. These measures of firing rate fluctuations should be distinguished from fluctuations in the local field potential, i.e., the LFP, itself. Thus, it should be emphasized that the power-spectrum of the LFP should not be confused with the power spectrum of the firing rate or HFB fluctuations, as the two types of fluctuations can even be antagonistic to each other under different experimental paradigms [93].

Moreover, previous studies have suggested that the LFP power spectrum is optimally fitted by setting the dynamics control parameter precisely at the tipping point of criticality (i.e., at $[eqn]$ ), which could potentially imply fundamental instability and a network easily crossing between subcritical to supercritical states [30,61]. However, our optimal fit for the HFB power spectra implies that relevant cortical networks operate below this critical transition point (i.e., at $[eqn]$ ), where stability is still maintained but close enough to critical dynamics to allow the system to computationally benefit from near-criticality effects.

Temporally structured drive can push cortical circuits towards criticality. Cultured networks, for instance, shift toward avalanche-like dynamics when exposed to patterned inputs [94]. Recent whole-brain studies likewise report brief approaches to the critical boundary during executive control and creative ideation [95,96]. Our model captures these observations with a single control knob: a modest gain increase (G ≈ 0.94 → 0.96) that modulates the system from a stable, reverberating sub-critical regime towards the critical transition point, where noise-driven excursions enable rapid, metastable switching. This mechanism provides a plausible neural substrate for the explore–exploit cycle implicated in creative cognition. Importantly, the operating point remains sub-critical, consistent with the “edge-of-criticality” framework proposed by [97], while avoiding the sampling artefacts that may falsely suggest full transition into critical mode [98].

Wider implications of the CSD phenomenon should be noted. For example, it is well known that sleep deprivation can lead to over excitability and even to epileptic seizures. Indeed, it is proposed that sleep deprivation influences the proximity to criticality [99,100]. It is also known that sleep deprivation slows down reaction times and increases the presence of slow frequencies [28,99,101]. These phenomena may all be related. It may be that extensive recall tasks, which bring the system closer to criticality, will not only manifest in CSD but even trigger seizures in those who are highly susceptible.

A key hypothesis of our study is that cortical activity is modulated closer to a critical point in order to shift to a more generative mode, which is required under a wide range of conditions, such as creative and associative thinking [13,15,102]. The prediction is that behaviors requiring a generative mode would lead to slowing down of the dynamics compared to resting state. Note that this may potentially manifest in other signatures of criticality, such as measures related to edge-of-chaos dynamics and neuronal avalanches.

Finally, it is worth noting that the relationship between the LFP and firing rates is more appropriately captured by measuring the changes in the exponent of the LFP power-law PSD, as observed in our previous intra-cranial analysis of visual responses [103]. While these issues could potentially be examined in our recurrent network simulations, they have not been pursued in this work.

Overall, the findings from our simulations suggest that the ultra-slow spontaneous fluctuations are an emergent network phenomenon and support the idea that the cortex operates near a critical point, allowing for rapid and flexible control over spontaneous behaviors while maintaining stability. Furthermore, they provide a mechanistic explanation for the ability to flexibly target specific categories for free behavior. This work makes specific predictions and opens new avenues for investigating the role of criticality in brain function and could potentially lead to the development of new computational models for understanding the neural basis of cognition and behavior.

Methods

Random recurrent network as a rate model

Rate models can be derived from the well-studied 102 [104]. They were explored within the context of iEEG data in recent literature [30,61]. In previous work [57], we have shown that such rate models, describing the neuronal firing rates, can be faithfully derived from large-scale conductance based networks resembling cortical organization. The simplicity of rate models makes them highly useful in studying properties of large networks.

Here, we considered a random recurrent network model consisting of N rate-based units representing single neurons or groups of neurons (see Figs 2-6). We closely follow the model proposed in [30] and adapted its implementation by ([105]; see: github.com/gyyang/KimYang2017_Cell) to explain the power spectrum of iEEG recordings. The rate dynamics of the j-th unit over time in this model is given by:

$$[eqn]$$

where rj is the j-th unit activity level, and τ is a characteristic integration time. The [+] brackets denote a rectified linear (ReLU) activation function. The γ variable denotes the slope of the firing rate – current curve (f-I curve). The weights $[eqn]$ describe the strength of the connection from the presynaptic unit k to the postsynaptic unit j, and $[eqn]$ denotes the external input to the j-th unit (drawn here randomly from a “standard normal” distribution).

The dynamical variables, rj, can be interpreted as representing synaptic variables, which effectively low-pass filter the neuronal firing rates [57].

Criticality in the random recurrent network

We consider a network architecture comprising a non-hierarchical cluster of recurrently connected neurons characterized by a semi-linear f-I curve with gain γ. Following [30], we choose the connections in this network (the entries of matrix W) to be sparse and random. More specifically, each entry is nonzero with probability p, and the nonzero entries are drawn from a normal distribution, with positive mean μ_conn_ and variance σ_conn_, normalized by the size of the network, N.

When set as a rate model (assumed to be represented by the iEEG high-frequency broadband power modulation), we hypothesize it depicts the summed recurrent synaptic current into a fraction of neurons in this network, along with an intrinsic complementary membrane noise or some form of external noise.

Previous work [30] demonstrated how the product of three main parameters controls the proximity to the network’s critical point. These are: $[eqn]$ , $[eqn]$ , and $[eqn]$ . At the tipping point of criticality, they satisfy the relation $[eqn]$ , see [30]. Thus, we define a control parameter $[eqn]$ , and require that it lies in the vicinity of the critical point - G $[eqn]$ 1. In a near-critical state, either one of these arguments composing the control parameter could be varied to yield similar results towards the crossing point into a supercritical state.

To provide additional insight in the present context, we note that near a steady-state, an algebraic linearized approximation of the spontaneous fluctuations dynamics can be obtained. We briefly summarize the details given in [30,58]:

The eigenvalue spectrum of such random matrices is composed of a cluster of complex eigenvalues with negative real part, and a single isolated real eigenvalue, which corresponds to the uniform eigenvector (namely, to the mean network activity).The matrix that controls the linearized dynamics is a random matrix whose dominant eigenvalue is given by $[eqn]$ . As the overall connection strength $[eqn]$ approaches 1, the network’s dominant feedback loop weakens its damping; activity therefore decays more and more slowly, until at $[eqn]$ it no longer decays at all and it experiences a transition from stable to unstable dynamics.

Importantly, the underlying analysis relies on random matrix theory and consequently assumes that the network has a large scale [59,106], specifically, a minimum of a few hundred units (see Appendix B, Fig A in S1 File). This scale is fundamental because random matrix theory draws on statistical properties that accurately reflect the system’s behavior only when applied to large matrices. Adequate network size is thus vital for capturing the delicate dynamics and stability features essential for close analysis of proximity to criticality.

Critical slowing down is manifested in the system’s temporal dynamics, as they are linked to the inverse of the abovementioned dominant eigenvalue (see SI for additional detail). As the system approaches criticality, signals reverberate in the network for increasingly longer time, producing a slowly decaying autocorrelation function and increasing the power of low frequencies.

For the sake of convenience, the following parameters were kept fixed across all simulation runs presented in the paper: μ_conn_ = 49.881 pA/Hz, σ_conn_ = 4.988 pA/Hz, *p *= 0.2.

Unless otherwise noted, the parameter that was used to control the proximity to criticality was the slope of the f-I curve, $[eqn]$ (for $[eqn]$ , we obtain $[eqn]$ ). We did however verify that qualitatively equivalent results can be reached by modifying G using other parameters. The mean and std of $[eqn]$ , the noisy input, was set to $[eqn]$

The characteristic integration time and PSD profile’s pattern

As for the characteristic integration time τ, we explored values in the range of 2 – 200 ms. Using a simple linearized mathematical approximation (see Appendix A; Appendix B, Fig C in S1 File and [30]) it can be shown that the commonly observed “knee”-bend in in the power spectral density (PSD) profile appears near $[eqn]$ , and that modulation of slow-frequencies is also significantly influenced by this the integration time constant, implying that the fit tuning process should carefully consider its setting.

Our findings here indicate that using the previously reported timescale of the membrane time constant in the proximity of τ = 20 ms [61–63] led to an optimal fit regarding the location of this knee in the activity rate PSD curve as well as the related slow frequencies slope. That is, for both of the cases explored in [7,13], as demonstrated in Figs 2 and 3 respectively.

The sampling parameter α, defined as the fraction (i.e., 0 < α < 1) of network nodes selected for summation, plays a significant role in determining the slope of medium range frequencies (> 0.1 Hz and < 10 Hz) in the PSD curve, since it influences the mixture of slow and fast dynamical modes exposed in the system. To best adjust the fit in the iEEG PSD mid-to high frequencies (> 1 Hz), α was set between 0.01 <= α <= 0.1 [30]. When constructing the PSD, we sum the activities of αN random nodes and compute the spectrum of the summed activity. The larger the value of α, the larger the contribution of the dominant slow mode to the mid-range PSD slope, since this mode is shared among the network nodes.

Simulation and parameter fitting efforts

Simulation code, including full implementation of the rate model and plotting of figures, was developed in Python, taking advantage of suitable open-sourced libraries’ methods (including, but not limited to, SciPy and Matplotlib), facilitating tuning efforts for fitting the empirical results (details below).

To determine a stationary signal (and the period to be omitted before signal stabilization), the ‘Augmented Dickey-Fuller’ method was iteratively applied, using the Python statsmodel library’s adfuller method, extending the delay, until the remaining segment’s stability reached p-value < 0.01. For the results depicted in Fig 2, the signal was low-pass filtered (for < 0.1 Hz) as in [7], using a Butterworth method. Correlation, auto- and cross-correlation coefficients were computed using time-lagged correlation methods. In panel 2d, the activity of two arbitrarily selected highly-correlated neurons is presented. In panel 2e, the power spectrum was computed using Matplotlib mlab’s PSD, implementing the Welch method.

A simple amplitude scaling factor k was applied to match the Welch method results with those presented by direct use of FFT power in [7]. A regression line curve was fitted using a curve-fitting method in the range of 0.1 – 10 Hz. Following [7], the autocorrelation was computed using signals from the first sampling only (resembling a single cluster of neurons in the network), and cross-correlation was computed iteratively between unit pairs of both samples. The ratio of the power law exponents of these cross- and auto-correlation functions at low frequencies. Finally, the power spectrum in panel 2e was obtained by applying the PSD method to the mean of these autocorrelation (black line) and cross-correlation (orange line) curves, which is then plotted in log scale.

For Fig 3, the simulation was arranged as consecutive “rest” vs “free-recall” 1800 s blocks. Additive noise $[eqn]$ was defined as constant value within each simulation run, so that $[eqn]$ $[eqn]$ . For the purpose of sampling a satisfactory distribution for statistical analysis, the simulation was repeated using different random seed values. The multiple realizations generated were independently compared with the empirical data. Other than that, similar simulation and computational methods were applied as in Fig 2.

Testing for various possible confounds, the simulation configuration for the results presented in Figs 4 – 6 was suitably adapted, but remained similar in concept to the settings in Fig 3. In data preparation for Fig 5, white noise was sampled in advance, applying suitable filters (using scipy signal.butter methods) over the random (white noise) normal distribution based temporal signal. For the map presented in Fig 6, we iterated through the parameter space of network connectivity and size, iteratively generating random network realizations of larger and larger scale. Thus, as the control parameter G was modified (via γ) from subcritical to near-critical range, values of the slow (normalized by fast) frequency power modulation were collected.

Goodness-of-fit and degrees-of-freedom correction

To quantify how closely the model reproduced the empirical PSD we computed LMSE, as the coefficient of determination in log-power space, $[eqn]$ .

Adjacent frequency bins are not independent because Welch segmentation and tapering introduce spectral leakage. We therefore replaced the nominal sample size $[eqn]$ by an effective sample size $[eqn]$ [107] obtained

from the residual series $[eqn]$ :

$$[eqn]$$

where $[eqn]$ is the lag-k autocorrelation of $[eqn]$ ; summation stops once $[eqn]$ . The corrected $[eqn]$ replaced $[eqn]$ , yielding conservative, assumption-compatible p-values that account for spectral redundancy.

The human iEEG data and its analysis

The empirical data used in this study has been previously reported in [7,13]. In both studies, intracranial recordings were obtained from patients with pharmacologically resistant epilepsy who were undergoing a surgical evaluation process for the treatment of their condition. In the first study [7], data were recorded from five participants (over multiple sessions) under the conditions of wakeful rest. During the wakeful rest, subjects’ eyes were closed and the room was dark and quiet, however, alertness was occasionally monitored. The second study [13] included twelve participants who performed a free recall task, as follows. Immediately after a 200 s period of closed-eyes resting state, participants were presented with 14 different pictures of famous faces and popular landmarks in a pseudorandom order (each picture was repeated four times, and shown for 1500 ms, with 750 ms inter-stimulus interval). After viewing the pictures, participants put on a blindfold (or closed their eyes) and began a short interference task of counting back from 150 in steps of 5 for 40 s. Upon completion, recall instructions were presented, guiding the patients to recall items from only one category at a time, starting with faces in one run and with places in a second run. Participants were asked to verbally describe each picture they recall, as soon as it comes to mind, with 2–3 prominent visual features. The duration of this free recall phase was 2.5 min per each category. In both experiments analyses of iEEG data were focused on changes in high-frequency broadband amplitude (HFB, also known as high-gamma, defined as the range 40–100 Hz in the first study and 60–160 Hz in the latter one), shown to be a reliable marker of local neuronal population activity [66].

Supporting information

S1 FileFig A. Relationship between network size and the effective control parameter ( $[eqn]$ ), illustrating the variability in proximity to criticality for small networks. The real part of the largest eigenvalue in the linearized matrix representation determines the network gain. To achieve critical dynamics, the effective gain should be 1. The figure shows the distribution of the effective gain as a function of network size when the parameters are set to yield a theoretical value of 1. Quantiles represent 100 random realizations of the network connectivity for each distinct network size. The analysis clearly demonstrates how the variance diminishes with increasing network size. Note (in the inset) that when the network size exceeds 200 units, the effective gain is very close to 1. For small networks, the network can be effectively sub- or super-critical, even when theoretically set to be near-critical. These results stem from the fact that the relevant analytical findings regarding random matrices hold only in the limit of very large networks. Fig B. Different settings of τ in fitting the PSD curve knee frequency. Parameter fitting of [7] empirical iEEG HFB signal results, using an artificial neural network model with 400 neurons, explored the optimal setting for τ. The autocorrelation (black line) and cross-correlation (orange line) PSD curve data, from two highly correlated neurons in the artificial network, are presented in panels a - c for different values of τ, 2, 20 and 200 ms, respectively. As indicated by arrows pointing to the knee position within the empirical autocorrelation (black) and cross-correlation (orange) PSD curve data, only panel b (reflecting τ = 20 ms) demonstrates an accurate match. Fig C. A Lorentzian fitting of free recall vs. rest HFB PSD curves. The power spectrum of the Lorentzian function is near-flat at low-frequencies and shows $[eqn]$ scaling at high frequencies, with a transition point set by the time-constant of the exponential. The fitting process maximized correlation between respective Lorentzian and HFB PSD curves. The different fitting for the sub- and near-critical curve fits lies in the setting of fslow which was slightly adapted. Fig D. Slow-fluctuations power increase in iEEG and in the recurrent network simulation. Single realization simulation results from the fitting of Norman et al 2018 free-recall vs. rest PSD curves are demonstrated using additive noise and gain control. In both panels, the iEEG power spectra obtained from empirically recorded free-recall target sites (red), relative to the non-targeted population (blue), are compared with respective simulated results. In the left panel, additive noise is used to modulate slow fluctuation power increase whereas the F-I curve gain parameter, γ, is modulated for this purpose in the right panel. Note that Panel 4c in the main text extends the fit shown here by demonstrating the influence of higher levels of additive noise on the amplification in power of the slow fluctuations. Fig E. The influence of additive noise on the ultra-slow frequency power. Simulation results from our recurrent random network model far below the critical point (left panel) and at its very near vicinity (right panel) indicate that the significant amplification of the ultra-slow frequency power due to increase in additive noise is dependent on the proximity to criticality. In both panels, the x-axis indicates additive noise and y-axis denotes the normalized slow-frequency power amplification. Results are based on 10 random realizations for each level of additive noise. Normalization is two-staged, first obtaining the ratio of slow over fast frequency power per each realization, and then dividing the obtained value by the 0-additive-noise mean. Statistical analysis indicates a significant effect only for the case of near critical simulation results (R^2^ = 0.862, p < 10^-20^). Fig F. Impact of different input noise spectral profiles on the PSD of a zero-connectivity network. Panels a-b demonstrate a simulated zero-connectivity network PSD when the injected noise is high-pass filtered (panel a) relative to the initial unfiltered white noise (panel b). The cutoff frequency was 2 Hz for the high-pass filter. As in Fig 5 of the main text, the scale-free spectral profile in the mid-to-high frequency range is qualitatively preserved.(DOCX)