Neuronal avalanches: Sandpiles of self-organized criticality or critical dynamics of brain waves?

Abstract

Analytical expressions for scaling of brain wave spectra derived from the general nonlinear wave Hamiltonian form show excellent agreement with experimental “neuronal avalanche” data. The theory of the weakly evanescent nonlinear brain wave dynamics [Phys. Rev. Research 2, 023061 (2020); J. Cognitive Neurosci. 32, 2178 (2020)] reveals the underlying collective processes hidden behind the phenomenological statistical description of the neuronal avalanches and connects together the whole range of brain activity states, from oscillatory wave-like modes, to neuronal avalanches, to incoherent spiking, showing that the neuronal avalanches are just the manifestation of the different nonlinear side of wave processes abundant in cortical tissue. In a more broad way these results show that a system of wave modes interacting through all possible combinations of the third order nonlinear terms described by a general wave Hamiltonian necessarily produces anharmonic wave modes with temporal and spatial scaling properties that follow scale free power laws. To the best of our knowledge this has never been reported in the physical literature and may be applicable to many physical systems that involve wave processes and not just to neuronal avalanches.

Graphical abstract

Keywords

nonlinear waves / critical exponent / Hamiltonian system / neuronal avalanches / critical dynamics

1 Introduction

In this paper, we show that an important observational phenomenon, the so-called “neuronal avalanches”, which have been noted and studied for almost two decades, can be naturally explained by the wave activity model. Both temporal and spatial scaling expressions analytically derived from nonlinear amplitude/phase evolutionary equations show excellent agreement with the experimental neuronal avalanche probability spectra. The model is not only able to reproduce general average power law exponent values and falloffs in the vicinity of the critical point, but also finds some very subtle but nevertheless clearly experimentally evident fine details, like bumps in the transition region at the edge of the power leg of the spatial probability spectra. Overall, the quantitative theoretical analysis presented in the paper clearly shows the relevance of the wave Hamiltonian framework to the neuronal avalanche dynamics. The analysis also suggests that instead of relying on clever but ad hoc analogies between live brain tissues and lifeless loose sand piles often used to construct a phenomenological statistical description of the scaling exponents, both the origin of the critical phenomena and the physics behind the neuronal avalanches can be understood from the same nonlinear wave dynamics that is responsible for the wide range of activities in the brain tissue. Those different activities are ranging from the linear coherently propagating waves to the nonlinear incoherent asynchronous spiking, including in between the peculiar power law-like coherence of the neuronal avalanches.

The standard view of brain electromagnetic activity classifies this activity into two significant but essentially independent classes. The first class includes a variety of the oscillatory and wave-like patterns that show relatively high level of coherence across a wide range of spatial and temporal scales [3]. The second class focusses on the asynchronous, seemingly incoherent spiking activity at scales of a single neuron and often uses various ad hoc neuron models [4–8] to describe this activity. Linking these two seemingly disparate classes to explain the emergence of oscillatory rhythms from incoherent activity is essential to understanding brain function and is typically posed in the form using the construct of networks of incoherently spiking neurons [9–11].

Coherent macroscopic behavior arising from seemingly incoherent microscopic processes naturally suggests the influence of critical phenomena, a potential model from brain activity that was bolstered by the experimental discovery of the “neuronal avalanches” [12, 13] where both spatial and temporal distributions of spontaneous propagating neuronal activity in 2D cortex slices were shown to follow scale-free power laws. This discovery has generated significant interest in the role and the importance of criticality in brain activity [14–20], especially for transmitting or processing information [21].

Although the precise neuronal mechanisms leading to the observed scale-free avalanche behavior is still uncertain after almost 20 years since their discovery, the commonly agreed upon paradigm is that this collective neuronal avalanche activity represents a unique and specialized pattern of brain activity. This pattern exists somewhere between the oscillatory, wave-like coherent activity and the asynchronous and incoherent spiking. Central to this claim of neuronal avalanches as a unique brain phenomena is that they do not show either wave-like propagation or synchrony at short scales, and thus constitute a new mode of network activity [12, 13]. This new activity can be phenomenologically described using the ideas of the self-organized criticality [22, 23], and extended to the mean-field theory of the self-organized branching processes (SOBP) [24–26].

However, despite the success of the SOBP theory in describing neuronal avalanche statistical properties, i.e., replicating the power law exponents based on the criticality considerations, the SOBP theory provides no explanation about the physical mechanisms of the critical behavior and its relationship to the development of the observed collective neuronal “avalanche” behavior. Because similar statistics can result from several mechanisms other than critical dynamics [27–29], it is essential to have a physical model that explains the relationship between the statistical properties and the existence, if any, of critical neural phenomena arising from the actual collective behavior of neuronal populations. While it is generally accepted in that some form of critical phenomena is at work, this has led to the presupposition of ad hoc descriptive models [30–33] that exhibit critical behavior but provide no insight into the actual physical mechanisms that might produce such critical dynamics. It has been suggested that the brain can be at the edge of a synchronization phase transition [33–35]. Another universally accepted view is that avalanches belong to the mean-field-directed percolation universality class [20]. But the percolation does not seem to be compatible with a synchronization transition, as synchronization transitions do not fulfill spatial correlations observed in experiments, and the exponents tend to differ from the exponent obtained by directed percolation [20].

In this paper, we show that our recently described theory of weakly evanescent brain waves (WETCOW), originally developed in Refs. [1, 2] and then reformulated in a general Hamiltonian framework [36] provides a physical theory, based on the propagation of electromagnetic fields through the highly complex geometry of inhomogeneous and anisotropic domain of real brain tissues, that explains the broad range of observed seemingly disparate brain wave characteristics. This theory produces a set of nonlinear equations for both the temporal and spatial evolution of brain wave modes that includes all possible nonlinear interactions between propagating modes at multiple spatial and temporal scales and degrees of nonlinearity. This theory bridges the gap between the two seemingly unrelated spiking and wave “camps” as the generated wave dynamics includes the complete spectra of brain activity ranging from incoherent asynchronous spatial or temporal spiking events, to coherent wave-like propagating modes in either temporal or spatial domains, to collectively synchronized spiking of multiple temporal or spatial modes. In this paper, we further demonstrate that the origin of these “avalanche” properties emerges directly from the same theory that produces this wide range of activity and does not require one to posit the existence of either new brain activity states, nor construct analogies between brain activity and ad hoc generic “sandpile” models.

We emphasize that although the general WETCOW theory describes complex interactions between modes, the explanation for neuronal avalanches and their attendant scaling properties presented in this paper are based on the analysis of a single wave mode with a completely arbitrary set of mode parameters. This includes arbitrary amplitude, phase, frequency, and criticality. No interaction between modes, except a general form of the third order nonlinearity that characterizes anharmonicity of the nonlinear wave modes due to nonresonant interaction of the linear modes, is needed to derive the scaling result. Thus a key result of this paper is the demonstration that neuronal avalanches and their attendant scaling properties are obtained within the simplest form of the WETCOW theory where mode coupling is ignored, but significantly without the ad-hoc and physically implausible assumptions typically made that the parameters of all network nodes are either constant and the same for all nodes [33], sometimes even including inter-mode coupling [35], or are generated from some ad-hoc artificial distributions [37], and require the addition of stochastic noise properties [38], etc. This emphasizes the generality and the importance of our derivation.

The generalization of our approach to multiple wave modes is very important for understanding the details of the functioning of the actual neuronal networks. In particular, WETCOW theory can be useful to relate biological learning to wave dynamics, demonstrating how wave dynamics can not only augment our understanding of cognition but provide the basis for a novel class of engineering analogs for both software and hardware learning systems that can operate with the extreme energy and data efficiency characteristics of biological systems that facilitate adaptive resilience in dynamic environments [39]. Additionally, multiple practical applications of the WETCOW model to analyses of complex clinical problems are possible. One example is an analysis of a complex problem of significant scientific and social importance, autism spectrum disorder (ASD), where the WETCOW-based method was able to detect spatially resolved frequency-dependent cognitive network anomalies with a high level of statistical evidence [40].

2 Weakly evanescent brain waves

Beginning from our nonlinear Hamiltonian formulation of the WETCOW theory [36], we have for an anharmonic wave mode

(1)\begin{array}{l}{H}^s(a,a^{†})=\mathrm{\Gamma }a{a}^{†}+a{a}^{†}\left[{\beta }_{a}a+{\beta }_{a^{†}}{a}^{†}-2\alpha {\left(a{a}^{†}\right)}^{1/2}\right],\end{array}

where a is a complex wave amplitude and a^{†} is its conjugate. H^s(a,a^{†}) denotes the Hamiltonian form for a single wave mode. The amplitude a denotes either temporal a_k(t) or spatial a_{\omega }(x) wave mode amplitudes that are related to the spatiotemporal wave field (electrostatic potential \psi (x,t)) through a Fourier integral expansions

(2)\begin{array}{l}{a}_k(t)={\displaystyle \frac{1}{2\pi }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\psi (x,t){\mathrm{e}}^{-\mathrm{i}\left(kx+{\omega }_{k}t\right)}\mathrm{d}x,}\end{array}

(3)\begin{array}{l}{a}_{\omega }(x)={\displaystyle \frac{1}{2\pi }\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\psi (x,t){\mathrm{e}}^{-\mathrm{i}\left(k_{\omega }x+\omega t\right)}\mathrm{d}t,}\end{array}

where for the sake of clarity we use one dimensional scalar expressions for spatial variables x and k, but it can be easily generalized for the multi dimensional wave propagation as well. The frequency \omega and the wave number k of the wave modes satisfy dispersion relation D(\omega ,k)=0, and {\omega }_k and k_{\omega } denote the frequency and the wave number roots of the dispersion relation (the structure of the dispersion relation and its connection to the brain tissue properties has been discussed in [1]).

The first term \mathrm{\Gamma }a{a}^{†} in (1) denotes the harmonic (quadratic) part of the Hamiltonian with either the complex valued frequency \mathrm{\Gamma }=\mathrm{i}\omega +\gamma or the wave number \mathrm{\Gamma }=\mathrm{i}k+\lambda that both include a pure oscillatory parts (\omega or k) and possible weakly excitation or damping rates, either temporal \gamma or spatial \lambda. The second anharmonic term is cubic in the lowest order of nonlinearity and describes the interactions between various propagating and nonpropagating wave modes, where \alpha, {\beta }_a and {\beta }_{a^{†}} are the complex valued strengths of those different nonlinear processes.

A set of derivations that lead to this description was presented in details in [1, 2, 36] and is based on considerations that follow from the most general form of brain electromagnetic activity expressed by Maxwell equations in inhomogeneous and anisotropic medium

\begin{array}{l}\mathrm{\nabla }\cdot \mathit{D}=\rho ,\mathrm{\nabla }\times \mathit{H}=\mathit{J}+{\displaystyle \frac{\mathrm{\partial }\mathit{D}}{\mathrm{\partial }t}}\Rightarrow {\displaystyle \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}}+\mathbf{\nabla }\cdot \mathit{J}=0.\end{array}

Using the electric field \mathit{E}=-\mathrm{\nabla }\psi, Ohm’s law \mathit{J}=\mathit{\sigma }\cdot \mathit{E} (where \mathit{\sigma }\equiv \{{\sigma }_{ij}\} is an anisotropic conductivity tensor), a linear electrostatic property for brain tissue \mathit{D}=\epsilon \mathit{E}, assuming that the scalar permittivity \epsilon is a “good” function (i.e., it does not go to zero or infinity everywhere) and taking the change of variables \mathrm{\partial }x\to \epsilon \mathrm{\partial }{x}^{\mathrm{\prime }}, the charge continuity equation for the spatial-temporal evolution of the potential \psi can be written in terms of a permittivity scaled conductivity tensor \mathbf{\Sigma }=\{{\sigma }_{ij}/\epsilon \} as

(4)\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}\left({\mathrm{\nabla }}^2\psi \right)=-\mathbf{\nabla }\cdot \mathbf{\Sigma }\cdot \mathrm{\nabla }\psi +\text{embassy}\text{F},\end{array}

where we have included a possible external source (or forcing) term \text{embassy}\text{F}. An anisotropic material has different properties in different directions. For example, porous media formed with a set of aligned conductive sheets and filled with some non-conductive dielectric between them will allow current flows very easily through each sheet, but won’t allow any current from one sheet to the adjacent one. This means that the component of the conductivity tensor in the direction normal to the sheet’s plane will be zero. Similarly, for brain fiber tissues the conductivity tensor \mathbf{\Sigma } might have significantly larger values along the fiber direction than across them. The charge continuity without forcing, i.e., (\text{embassy}\text{F}=0) can be written in tensor notation as

(5)\begin{array}{l}{\mathrm{\partial }}_t{\mathrm{\partial }}_i^2\psi +{\mathrm{\Sigma }}_{ij}{\mathrm{\partial }}_i{\mathrm{\partial }}_j\psi +\left({\mathrm{\partial }}_i{\mathrm{\Sigma }}_{ij}\right)\left({\mathrm{\partial }}_j\psi \right)=0,\end{array}

where repeating indices denote summation and \mathbf{\nabla }\cdot \mathbf{\Sigma }\cdot \mathrm{\nabla }\psi \equiv {\mathrm{\partial }}_i({\mathrm{\Sigma }}_{ij}{\mathrm{\partial }}_j\psi )={\mathrm{\Sigma }}_{ij}{\mathrm{\partial }}_i{\mathrm{\partial }}_j\psi +\left({\mathrm{\partial }}_i{\mathrm{\Sigma }}_{ij}\right)\left({\mathrm{\partial }}_j\psi \right). Simple linear wave analysis, i.e. substitution of \psi \sim \exp (-\mathrm{i}(\mathit{k}\cdot \mathit{r}-\mathrm{\Omega }t)), where \mathit{k} is the wavenumber, \mathit{r} is the coordinate, \mathrm{\Omega } is the frequency and t is the time, gives the following complex dispersion relation:

(6)\begin{array}{l}D(\mathrm{\Omega },\mathit{k})=-\mathrm{i}\mathrm{\Omega }{k}_i^2-{\mathrm{\Sigma }}_{ij}{k}_i{k}_j-\mathrm{i}{\mathrm{\partial }}_i{\mathrm{\Sigma }}_{ij}{k}_j=0,\end{array}

which is composed of the real and imaginary components:

(7)\begin{array}{l}\gamma \equiv \mathrm{\Im }\mathrm{\Omega }={\mathrm{\Sigma }}_{ij}{\displaystyle \frac{k_i{k}_j}{k^2}},\omega \equiv \mathrm{\Re }\mathrm{\Omega }=-{\displaystyle \frac{{\mathrm{\partial }}_i{\mathrm{\Sigma }}_{ij}{k}_j}{k^2}}.\end{array}

Although in this general form the electrostatic potential \psi, as well as the dispersion relation D(\mathrm{\Omega },\mathit{k}), describe three dimensional wave propagation, we have shown [1, 2] that in anisotropic and inhomogeneous media some directions of wave propagation are more equal than others with preferred directions determined by the complex interplay of the anisotropy tensor and the inhomogeneity gradient. While this is of significant practical importance, in particular because the anisotropy and inhomogeneity can be directly estimated from non-invasive methods, for the sake of clarity we focus here on the one dimensional scalar expressions for spatial variables x and k that can be easily generalized for the multi dimensional wave propagation as well.

The multiple temporal a_k(t) or spatial a_{\omega }(x) wave mode amplitudes can be used to define the time dependent wave number energy spectral density a_k(t) or the position dependent frequency energy spectral density a_{\omega }(x) for the spatiotemporal wave field \psi (x,t) as

(8)\begin{array}{l}{a}_k(t)=|a_k(t){|}^2,a_{\omega }(x)=|a_{\omega }(x){|}^2,\end{array}

or alternatively we can add additional length or time normalizations to convert those quantities to power spectral densities instead.

The network Hamiltonian form that describes discrete spectrum of those multiple wave modes was presented in Ref. [36] as

(9)\begin{array}{rl}H(\mathit{a},{\mathit{a}}^{†})=& {\displaystyle \sum _n[H^s(a_n,a_n^{†})}\\ & +\sum _{m\ne n}\left(a_n{r}_{nm}{a}_m^{†}+a_n^{†}{r}_{nm}^{\ast }{a}_m\right)],\end{array}

where H^s(a_n,a_n^{†}) denotes the Hamiltonian form (1) for a single wave mode n, the single mode amplitude a_n again denotes either a_k or a_{\omega }, \mathit{a}\equiv \{a_n\} and r_{nm}=R_{nm}{\mathrm{e}}^{\mathrm{i}{\mathrm{\Delta }}_{nm}} is the complex network adjacency matrix with R_{nm} providing the coupling power and {\mathrm{\Delta }}_{nm} taking into account any possible differences in phase between network nodes. This description includes both amplitude \mathrm{\Re }(a) and phase \mathrm{\Im }(a) mode coupling and as shown in Ref. [36] allows for significantly unique synchronization behavior different from both phase coupled Kuramoto oscillator networks and from networks of amplitude coupled integrate-and-fire neuronal units.

The third order nature of the theory is of particular interest and provides the theory with a broad range of applicability. It has distinctly different characteristics than the harmonic oscillator. Of particular importance is the fact that the third order terms become important when wave amplitudes are high enough but only if or until higher order terms are absent or suppressed by some physical mechanism. This suppression becomes significant in incorporating the anisotropic inhomogeneous and resistive nature of brain tissues. An important consequence derived in Ref. [1] is that the inverse frequency–wave number proportionality of the linear wave dispersion guarantees that the resonant terms higher than the third order are not important and can be neglected and, at the same time, the non-resonant third order terms (that are typically excluded when compared to the resonant terms) should now be retained resulting in the third order form of Hamiltonian (1). It is our contention, and the subject of future studies, that the anharmonic third order forms (1) and (9) are not brain specific and can be used to describe oscillations and waves in active media abundant in various areas of physics.

Although the Fourier integrals (2,3) used for expansion of the spatiotemporal wave field \psi into a set of wave modes imply the presence of a large (actually infinite) number of modes in the network Hamiltonian (9) the derivation of neuronal avalanches is evident even without this generality of this coupling between modes expressed by the coupling parameters r_{nm} as it was done in Ref. [36]. Thus we will consider an ensemble of noninteracting modes, effectively setting r_{nm}=0, for the analysis of this paper. But importantly we will not make any assumptions about parameters of the single mode Hamiltonian form (1), assuming that all parameters (\mathrm{\Gamma }, {\beta }_a, {\beta }_{a^{†}}, \alpha) are arbitrary and do not carry any mode dependence. This is a non-trivial point worth emphasizing, as it is a departure from the extant literature wherein the ad-hoc, and physically implausible, assumption of the equivalence of network nodes is made. Therefore, we will proceed with our analysis of a single mode amplitude a suppressing all subscripts and indices, and assuming that a denotes a_n where n may represent either an arbitrary wave number k from a range of wave numbers k_0\le k\le k_1 or an arbitrary frequency \omega from a range of wave frequencies {\omega }_0\le \omega \le {\omega }_1.

3 Single anharmonic mode criticality

An equation for the nonlinear oscillatory amplitude a then can be expressed as a derivative of the Hamiltonian form

(10)\begin{array}{l}{\displaystyle \frac{\mathrm{d}a}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{\partial }{H}^s}{\mathrm{\partial }{a}^{†}}}\equiv \mathrm{\Gamma }a+{\beta }_{a^{†}}a{a}^{†}+{\beta }_a{a}^2-\alpha a(a{a}^{†}{)}^{1/2},\end{array}

after removing the constants with a substitution of {\beta }_{a^{†}}=1/2{\tilde{\beta }}_{a^{†}} and \alpha =1/3\tilde{\alpha } and dropping the tilde. We note that although (10) is an equation for the temporal evolution, the spatial evolution of the mode amplitudes a_{\omega }(x) can be described by a similar equation substituting temporal variables by their spatial counterparts, i.e., (t,\omega ,\gamma )\to (x,k,\lambda ).

Splitting (10) into an amplitude/phase pair of equations using a=A{e}^{i\varphi }, assuming {\beta }_a=\widetilde{{\beta }_a}{\mathrm{e}}^{-\mathrm{i}{\delta }_a}, {\beta }_{a^{†}}={\tilde{\beta }}_{a^{†}}{\mathrm{e}}^{\mathrm{i}{\delta }_{a^{†}}}, and scaling the variables as

(11)\begin{array}{l}A=\gamma \tilde{A},t={\displaystyle \frac{\tau }{\gamma }},\omega =\tilde{\omega }\gamma ,\end{array}

gives the set of equations

(12)\begin{array}{l}{\displaystyle \frac{\mathrm{d}\tilde{A}}{\mathrm{d}\tau }}=\tilde{A}+{\tilde{A}}^2\left({\beta }_{a^{†}}\cos {\mathrm{\Omega }}_{a^{†}}+{\beta }_a\cos {\mathrm{\Omega }}_a-\alpha \right),\end{array}

(13)\begin{array}{l}{\displaystyle \frac{\mathrm{d}\varphi }{\mathrm{d}\tau }}=\tilde{\omega }+\tilde{A}\left(-{\beta }_{a^{†}}\sin {\mathrm{\Omega }}_{a^{†}}+{\beta }_a\sin {\mathrm{\Omega }}_a\right),\end{array}

where {\mathrm{\Omega }}_a\equiv \varphi -{\delta }_a, {\mathrm{\Omega }}_{a^{†}}\equiv \varphi -{\delta }_{a^{†}} (please see [36] for more details).

These equations can further be cast into a more compact form by defining

(14)\begin{array}{l}\mathit{\beta }=(\begin{array}{c}{\beta }_a\\ {\beta }_{a^{†}}\end{array}),\mathit{u}=(\begin{array}{c}{\mathrm{e}}^{\mathrm{i}{\delta }_a}\\ {\mathrm{e}}^{\mathrm{i}{\delta }_{a^{†}}}\end{array}),\mathit{v}=(\begin{array}{c}\phantom{-}\mathrm{i}{\mathrm{e}}^{\mathrm{i}{\delta }_a}\\ -\mathrm{i}{\mathrm{e}}^{\mathrm{i}{\delta }_{a^{†}}}\end{array}),\end{array}

so that

(15)\begin{array}{l}{z}_a=\mathit{\beta }\cdot \mathit{u}=X_a+\mathrm{i}{Y}_a,\end{array}

(16)\begin{array}{l}{z}_{\varphi }=\mathit{\beta }\cdot \mathit{v}=X_{\varphi }+\mathrm{i}{Y}_{\varphi },\end{array}

where

(17)\begin{array}{l}{R}_a=\left|z_a\right|=\sqrt{X_a^2+Y_a^2},\end{array}

(18)\begin{array}{l}{R}_{\varphi }=\left|z_{\varphi }\right|=\sqrt{X_{\varphi }^2+Y_{\varphi }^2},\end{array}

(19)\begin{array}{l}{\mathrm{\Phi }}_a=\arg (z_a)=\arctan {\displaystyle \frac{Y_a}{X_a}},\end{array}

(20)\begin{array}{l}{\mathrm{\Phi }}_{\varphi }=\arg (z_{\varphi })=\arctan {\displaystyle \frac{Y_{\varphi }}{X_{\varphi }}},\end{array}

whereupon (12) and (13) can be rewritten as

(21)\begin{array}{l}{\displaystyle \frac{\mathrm{d}\tilde{A}}{\mathrm{d}\tau }}=\tilde{A}+{\tilde{A}}^2\left[R_a\cos (\varphi -\mathrm{\Phi })-\alpha \right],\end{array}

(22)\begin{array}{l}{\displaystyle \frac{\mathrm{d}\varphi }{\mathrm{d}\tau }}=\tilde{\omega }+\tilde{A}{R}_{\varphi }\cos \varphi ,\end{array}

where \mathrm{\Phi }={\mathrm{\Phi }}_a-{\mathrm{\Phi }}_{\varphi }.

An equilibrium (i.e., \mathrm{d}\tilde{A}/\mathrm{d}\tau =\mathrm{d}\varphi /\mathrm{d}\tau =0) solution of (21) and (22) can be found from

(23)\begin{array}{l}-{\displaystyle \frac{R_{\varphi }}{\tilde{\omega }}}\cos \varphi +R_a\cos (\varphi -\mathrm{\Phi })-\alpha =0,\end{array}

as {\varphi }_e={\varphi }_0\equiv const and {\tilde{A}}_e=\tilde{\omega }/R_{\varphi }\cos {\varphi }_0\equiv const. This shows that for \alpha >R_a there exist critical values of \tilde{\omega } and A_e, such that

(24)\begin{array}{l}{\tilde{\omega }}_c={\displaystyle \frac{R_{\varphi }\cos {\varphi }_c}{\alpha +R_a\cos ({\varphi }_c+\mathrm{\Phi })}},{\tilde{A}}_c={\tilde{\omega }}_c/R_{\varphi },\end{array}

(25)\begin{array}{l}{\varphi }_c=\arctan \left[{\displaystyle \frac{R_a\sin \mathrm{\Phi }}{\sqrt{{\alpha }^2-(R_a\sin \mathrm{\Phi }{)}^2}}}\right],\end{array}

which can also be expressed in terms of critical value of one of the unscaled variables, either \omega or \gamma

(26)\begin{array}{l}{\omega }_c=\gamma {\tilde{\omega }}_c,\mathrm{o}\mathrm{r}{\gamma }_c={\displaystyle \frac{\omega }{{\tilde{\omega }}_c}}.\end{array}

The linearized system of equations at (A_c,{\varphi }_c) then can be obtained as

(27)\begin{array}{rl}\frac{\mathrm{d}\tilde{A}}{\mathrm{d}\tau }=& \left(1+2{\tilde{A}}_c\left[R_a\cos ({\varphi }_c-\mathrm{\Phi })-\alpha \right]\right)\tilde{A}\\ & -{\tilde{A}}_c^2{R}_a\sin ({\varphi }_c-\mathrm{\Phi })\varphi ,\end{array}

(28)\begin{array}{l}{\displaystyle \frac{\mathrm{d}\varphi }{\mathrm{d}\tau }}=R_{\varphi }\cos {\varphi }_c\tilde{A}-{\tilde{A}}_c{R}_{\varphi }\sin {\varphi }_c\varphi .\end{array}

Taking for example a symmetric case \mathrm{\Phi }=\pi the eigenvalues of the Jacobian matrix are

(29)

Thus, the equilibrium solution provides the locus of the saddle node on an invariant circle bifurcation point where the nonlinear spiking oscillations occur (as was shown both in Refs. [1, 2] and in Ref. [36]).

4 Effective spiking rate

The effective period of spiking {\text{embassy}\text{T}}_s (or its inverse – either the firing rate 1/{\text{embassy}\text{T}}_s or the effective firing frequency 2\pi /{\text{embassy}\text{T}}_s) can be estimated from (22) by substituting {\tilde{A}}_c for \tilde{A} (assuming that the change of amplitude \tilde{A} is slower than the change of the phase \varphi) as

(30)\begin{array}{l}{\text{embassy}\text{T}}_s={\displaystyle \underset{0}{\overset{2\pi }{\int }}\frac{\mathrm{d}\varphi }{\tilde{\omega }+{\tilde{\omega }}_c\cos \varphi }=\frac{2\pi }{\sqrt{{\tilde{\omega }}^2-{\tilde{\omega }}_c^2}},}\end{array}

giving the unscaled effective spiking period T_s and the effective firing frequency {\omega }_s

(31)\begin{array}{l}{T}_s=\frac{{\text{embassy}\text{T}}_s}{\gamma }={\displaystyle \frac{2\pi }{\omega \sqrt{1-{\gamma }^2/{\gamma }_c^2}}=\frac{2\pi }{\omega \sqrt{1-{\omega }_c^2/{\omega }^2}},}\end{array}

(32)\begin{array}{l}{\omega }_s={\displaystyle \frac{2\pi }{T_s}=\omega \sqrt{1-{\omega }_c^2/{\omega }^2},}\end{array}

with the periodic amplitude \tilde{A} reaching the maximum {\tilde{A}}_{max}=1/(\alpha -R_a) and the minimum {\tilde{A}}_{min}=1/(\alpha +R_a) for \mathrm{d}\tilde{A}/\mathrm{d}\tau =0 when \varphi =\mathrm{\Phi } and \varphi =\mathrm{\Phi }+\pi respectively. For a harmonic oscillator the amplitude is constant (does not change at all) and the phase is changing rapidly. In our model, we assume that the wave modes are anharmonic, but the anharmonicity is not very large, therefore, we can assume that the phase variable is changing faster than the amplitude variable and substiture {\tilde{A}}_c for \tilde{A}.

The expressions (31) and (32) are more general than typically used expressions for the scaling exponent in the close vicinity |\gamma -{\gamma }_c|\ll {\gamma }_c of the critical point [41–43]. They allow recovery of the correct T limits both at \gamma \to {\gamma }_c with the familiar T\sim 1/\sqrt{{\gamma }_c-\gamma } scaling and at \gamma \sim 0 with the period T approaching T_0 as T\sim T_0+O({\gamma }^2)\equiv 2\pi /\omega +O({\gamma }^2), where T_0 is the period of linear wave oscillations with the frequency \omega. In the intermediate range the expressions (31) and (32) show reasonable agreement (Fig.1) with peak–to–peak period/frequency estimates from direct simulations of the system (12) and (13).

5 Temporal probability of single spike detection

Taking into account that the initial phase of spiking solutions of (21) and (22) is a random variable uniformly distributed on [0,2\pi ] interval, the probability that a spike (either positive or the more frequently experimentally reported negative [12, 13]) with duration width \delta t_s and with the total period between the spikes (T_s) will be detected is simply \delta t_s/T_s – where the distance between spikes is determined as the time interval needed for 2\pi radian phase change, that is the effective spiking period T_s. Assuming initially that the spike width \delta t_s does not change when approaching the critical point {\omega }_c, \delta t_s can be approximated by some fixed fraction of the linear wave period, i.e., \delta t_s\sim \pi /\omega, that gives for the probability density

(33)\begin{array}{l}{P}_k^{\{\omega \}}(\omega ,{\omega }_c)\sim {\omega }^{-1}\sqrt{{\omega }^2/{\omega }_c^2-1},\end{array}

for every wave mode with the wavenumber k. It should be noted that the probability density P_k^{\{\omega \}} has no relation to, and should not be confused, with the frequency energy spectral density {\mathrm{\Pi }}_{\omega }(x) (or with the power spectral density).

Transforming the frequency dependence of the wavenumber spectra P_k^{\{\omega \}} to the temporal domain (T=2\pi /\omega, T_c=2\pi /{\omega }_c)

(34)\begin{array}{rl}\underset{{\omega }_c}{\overset{\mathrm{\infty }}{\int }}{P}_k^{\{\omega \}}(\omega ,{\omega }_c)\mathrm{d}\omega & =\underset{0}{\overset{T_c}{\int }}{P}_k^{\{\omega \}}\left(\frac{2\pi }{T},\frac{2\pi }{T_c}\right)\frac{2\pi }{T^2}\mathrm{d}T\\ & =\underset{0}{\overset{T_c}{\int }}{P}_k^{\{T\}}(T,T_c)\mathrm{d}T,\end{array}

gives for the temporal probability density P_k^{\{T\}}

(35)\begin{array}{l}{P}_k^{\{T\}}(T,T_c)\sim T^{-2}\sqrt{1-T^2/T_c^2},\end{array}

hence the scaling of the temporal probability density P_k^{\{T\}} follows the power law with −2 exponent with additional \sqrt{1-T/T_c} falloff in close vicinity of the critical point in agreement with temporal scaling of neuronal avalanches reported in Refs. [12, 13].

6 Multi-mode avalanche probability

The above single wave mode analysis shows that the probability density P_k^{\{T\}} for any arbitrarily selected wave mode k with an arbitrarily chosen threshold follows a power law distribution with −2 exponent, therefore, a mixture of multiple independent wave modes that enters into the spatiotemporal wave field \psi (x,t) with different amplitudes and different thresholds will again show nothing more than the same power law distribution.

To clearly demonstrate that the probability density function P_k^{\{T\}} of finding a spike reflects the avalanche duration distribution we conducted a simple numerical experiment using a procedure that replicates the original experimental method of computing neuronal avalanches employed in the original papers by Beggs and Plenz [12, 13]. We used 10⁶ wave modes with arbitrary parameters and computed avalanches by cutting the temporal series with a threshold (converting the event to a single dot or “spike”), then binning the signal using a time equal to the average inter-spike interval. After that, an avalanche duration is given by the time between two empty bins. Fig.2 compares the avalanche distribution for all wave modes with the −2 exponent.

Similar proof of equivalence of the single mode probability density function P_k^{\{T\}} and a probability density of a multi-mode avalanche event obtained by the method replicating the experimental demarcation of the quiescence, that we will denote as p^a(T), can be also derived using simple analytical considerations. The probability P_{0\le T^{\mathrm{\prime }}\le T+\mathrm{\Delta }T}^a that an avalanche happens at any time between 0 and T+\mathrm{\Delta }T, where \mathrm{\Delta }T is some small binning interval used by the above experimental method, can be expressed as

(36)\begin{array}{rl}{P}_{0\le T^{\mathrm{\prime }}\le T+\mathrm{\Delta }T}^a=& \underset{0}{\overset{T+\mathrm{\Delta }T}{\int }}{p}^a(T^{\mathrm{\prime }})\mathrm{d}{T}^{\mathrm{\prime }}=\underset{0}{\overset{T}{\int }}{p}^a(T^{\mathrm{\prime }})\mathrm{d}{T}^{\mathrm{\prime }}\\ & +\underset{T}{\overset{T+\mathrm{\Delta }T}{\int }}{p}^a(T^{\mathrm{\prime }})\mathrm{d}{T}^{\mathrm{\prime }}\approx P_{0\le T^{\mathrm{\prime }}\le T}^a+p^a(T)\mathrm{\Delta }T.\end{array}

Since the probability

(37)\begin{array}{l}{P}_{k_j}^{\{T\}}(T,T_{c_j})\mathrm{\Delta }T={\nu }_j{T}^{-2}\sqrt{1-T^2/T_{c_j}^2}\mathrm{\Delta }T,\end{array}

(where {\nu }_j is an arbitrary mode specific proportionality constant) describes the probability of finding a signal for a single mode j (j=1,\dots ,N) in a time interval between T and T+\mathrm{\Delta }T, the probability that the condition for detection of a multi mode avalanche is recorded in the same interval can be expressed as

(38)\begin{array}{r}{p}^a(T)\mathrm{\Delta }T=\left[1-P_0(T-\mathrm{\Delta }T)\right]\times P_0(T),\end{array}

\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{P}_0(T)=\prod _{j=1}^N\left(1-P_{k_j}^{\{T\}}(T,T_{c_j})\mathrm{\Delta }T\right),

where all wave modes are assumed to be independent. The second factor P_0(T) represents the probability that there is no signal for any of the modes detected between T and T+\mathrm{\Delta }T. The first factor (1-P_0(T-\mathrm{\Delta }T)) makes sure that no avalanche was recorded in the previous \mathrm{\Delta }T bin, that is a signal for at least one mode has been found in the interval between T-\mathrm{\Delta }T and T.

An expansion of (38) in the leading order of \mathrm{\Delta }T gives for the avalanche probability density p^a(T)

(39)\begin{array}{l}{p}^a(T)\approx T^{-2}{\displaystyle \sum _{j=1}^N{\nu }_j\sqrt{1-T^2/T_{c_j}^2},}\end{array}

that is the avalanche probability density p^a(T) shows the same T^{-2} scaling as the probability density of finding signal for a single mode.

If additionally the criticality parameters T_{c_j} for all wave modes k_j are assumed to be the same (T_{c_j}\equiv T_c) then the avalanche probability density scaling takes exactly the same form as the single mode probability density

(40)\begin{array}{l}{p}^a(T)\sim P_k^{\{T\}}(T)\sim T^{-2}\sqrt{1-T^2/T_c^2}.\end{array}

7 Spatial spike detection probability

Due to the reciprocity of the temporal and spatial representations of the Hamiltonian form (1) equations for the spatial wave amplitude have the same form as the temporal equations (21) and (22)

(41)\begin{array}{l}{\displaystyle \frac{\mathrm{d}\tilde{A}}{\mathrm{d}\xi }=\tilde{A}+{\tilde{A}}^2\left[R_a\cos (\varphi -\mathrm{\Phi })-\alpha \right],}\end{array}

(42)\begin{array}{l}{\displaystyle \frac{\mathrm{d}\varphi }{\mathrm{d}\xi }=\tilde{k}+\tilde{A}{R}_{\varphi }\cos \varphi ,}\end{array}

under similar scaling (the spatial equivalent of (11)) of the wave amplitude, the coordinate, and the wave number

(43)\begin{array}{l}A=\lambda \tilde{A},x={\displaystyle \frac{\xi }{\lambda },k=\tilde{k}\lambda .}\end{array}

In the spatial domain, this leads to the critical parameters {\tilde{A}}_c and {\tilde{k}}_c

(44)\begin{array}{rl}{\tilde{k}}_c& =\frac{R_{\varphi }\cos {\varphi }_c}{\alpha +R_a\cos ({\varphi }_c+\Phi )},{\tilde{A}}_c={\tilde{k}}_c/R_{\varphi },\\ {\varphi }_c& =\arctan \left[\frac{R_a\sin \Phi }{\sqrt{{\alpha }^2-(R_a\sin \Phi {)}^2}}\right].\end{array}

Although our simple one dimensional scaling estimates do not take into account the intrinsic spatial scales of the brain, e.g., cortex radius of curvature, cortical thickness, etc., nevertheless, even in this simplified form the similarity between spatial and temporal nonlinear equations suggests that the nonlinear spatial wave behavior will generally look like spiking in the spatial domain where some localized regions of activity are separated by areas that are relatively signal free and this separation will increase near the critical point. Exactly this behavior was reported in the original experimental studies of the neuronal avalanches [12, 13], where it was stated that the analysis of the contiguity index revealed that activity detected at one electrode is most often skipped over the nearest neighbors. Interestingly, this experimental observation of near critical nonlinear waves was presented as the indicator that the activity propagation is not wave-like. But we see here that they are directly explained within the context of the WETCOW wave model. Of significant practical importance will be the effects of the intrinsic spatial scales of the brain that will certainly affect the details of the spatial critical wave dynamics and so their inclusion will be important for more completely characterizing the details of brain criticality and will be the focus of future investigations.

Using the spatial equations (41) and (42) similar scaling results can be obtained for the wave number k and the linear spatial dimension L probabilities for every wave mode with the frequency \omega as

(45)\begin{array}{l}{P}_{\omega }^{\{k\}}(k,k_c)\sim k^{-1}\sqrt{k^2/k_c^2-1},\end{array}

(46)\begin{array}{l}{P}_{\omega }^{\{L\}}(L,L_c)\sim L^{-2}\sqrt{1-L^2/L_c^2},\end{array}

where L is the linear spatial scale related to the wave number as k=2\pi /L.

The linear spatial dimension of the avalanche L is related to its area S on a 2 dimensional surface as L=\sqrt{S}, hence

(47)\begin{array}{rl}\underset{0}{\overset{L_c}{\int }}{P}_{\omega }^{\{L\}}(L,L_c)\mathrm{d}L& =\underset{0}{\overset{S_c}{\int }}\frac{P_{\omega }^{\{L\}}\left(\sqrt{S},\sqrt{S_c}\right)}{2\sqrt{S}}\mathrm{d}S\\ & =\underset{0}{\overset{S_c}{\int }}{P}_{\omega }^{\{S\}}(S,S_c)\mathrm{d}S,\end{array}

(48)\begin{array}{l}{P}_{\omega }^{\{S\}}(S,S_c)\sim S^{-3/2}\sqrt{1-S/S_c},\end{array}

hence the spatial probability scaling for the size S follows the power law with −3/2 exponent again with additional \sqrt{1-S/S_c} falloff in close vicinity of the critical point, that is also in agreement with experimentally reported spatial scaling of neuronal avalanches [12, 13]. We would like to mention that the nonlinear anharmonic oscillations described by (21) and (22) only exists for frequencies and wave numbers that are above the critical frequency {\omega }_c or the critical wave number k_c values that define maximal possible temporal T_c or spatial L_c scales of the nonlinear oscillations. If the finite system sizes are below those maximal values the cutoffs will be defined by the system scales.

We would like emphasize again the generality of our analysis that makes no assumptions about parameters used in Hamiltonian form (1), and hence in the equations (21) and (22) or (41) and (42), analytically deriving scaling valid for a wide (and arbitrary) range of those parameters. This is in striking difference from analyses and results based on oversimplified ad-hoc numerical studies of synchronization in networks [33, 35]. Those typical numerical analysis studies consider networks of completely identical individual nodes sometimes even globally connected with completely identical weights. Therefore, all these studies require artificial (and significantly high level of) noise added to each node just to be able to impose some range of scales into the system. This is an artificial and, as demonstrated here, unnecessary complication. The consequence of such models is that they are capable of obtaining something that resembles scale free behavior with exponent values that are in general rather vague and strongly noise dependent. Without this sufficiently strong noise those studies of course are not capable to show any scale free behavior. It is essential to realize that such models are thus highly dependent on the noise properties, and less so on the actual properties of the brain tissue itself as in the WETCOW theory, which is the critical link to practical applications of any brain activity theory. By contrast, no externally induced stochasticity in the form of additional noise term is required for our analysis.

Another important point is that for deriving scale free exponents in our approach we don’t require to know the details of the coupling between nodes, essentially viewing all nodes as completely noninteracting. Presence of interactions in the form of (9) will not modify our analysis, and will not require any of the common ad-hoc assumptions of identical global coupling between nodes [35]. When coupling between some of the nodes in (9) is sufficiently strong and these nodes are completely synchronized, we can always replace this subset of completely synchronized nodes by a single node and continue again with the same presented in this paper “coupling-independent” node analysis.

8 Effects of criticality on spike length

The assumption of the fixed spike duration \delta t_s used in (33) and (35) (or the spike length for spatial spiking in (45) and (46)) can be improved by estimating the scaling of the spike width as a function of the criticality parameter from the amplitude equation (we will use the temporal form of the equation (21) but the spatial analysis with equation (41) is exactly the same).

Dividing equation (21) by \tilde{A} and taking an integral around some area in the vicinity of the amplitude peak {\tilde{A}}_{max} we can write

(49)\begin{array}{l}{\displaystyle \underset{{\tilde{A}}_{-}}{\overset{{\tilde{A}}_{+}}{\int }}\frac{1}{\tilde{A}}\mathrm{d}\tilde{A}=\underset{{\tau }_{-}}{\overset{{\tau }_{+}}{\int }}\mathrm{d}\tau +\underset{{\mathrm{\Phi }}_{-}}{\overset{{\mathrm{\Phi }}_{+}}{\int }}\frac{{\tilde{\omega }}_c}{R_{\varphi }}\frac{R_a\cos (\varphi -\mathrm{\Phi })-\alpha }{\tilde{\omega }+{\tilde{\omega }}_c\cos \varphi }\mathrm{d}\varphi ,}\end{array}

where {\tau }_{\pm }={\tau }_{max}\pm \delta \tau, and {\tau }_{max} is the location of spiking peak. Neglecting the spike shape asymmetries, i.e., assuming that {\tau }_{\pm } correspond to symmetric changes in both the amplitudes {\tilde{A}}_{\pm }=\tilde{A}({\tau }_{\pm })={\tilde{A}}_{max}-\delta \tilde{A}, and the phases {\mathrm{\Phi }}_{\pm }=\mathrm{\Phi }({\tau }_{\pm })=\mathrm{\Phi }\pm \delta \mathrm{\Phi }, we can then estimate the spike duration \delta t_s\equiv ({\tau }_{+}-{\tau }_{-})/\gamma as

(50)\begin{array}{l}\delta t_s={\displaystyle \frac{1}{\gamma }\underset{\mathrm{\Phi }-\delta \mathrm{\Phi }}{\overset{\mathrm{\Phi }+\delta \mathrm{\Phi }}{\int }}\frac{1-R[\cos \mathrm{\Phi }+\cos (\varphi -\mathrm{\Phi })]}{\tilde{\omega }+{\tilde{\omega }}_c\cos \varphi }\mathrm{d}\varphi ,}\end{array}

where, similar to the spiking period estimation in (31), we again assume that the main contribution comes from the change of the oscillation phase, hence {\tilde{A}}_c can be substituted for \tilde{A}. For \delta \mathrm{\Phi } some fixed value that is smaller or around a quarter of the period (i.e., \delta \mathrm{\Phi }\lesssim \pi /2) can be chosen, and R={\tilde{\omega }}_c{R}_a/R_{\varphi }.

An expression (50) can be evaluated in closed form but we do not include it here and instead plotted the final spatial probability density spectra P(S/S_c), similarly obtained from the expression for \delta l_s/L_s again substituting L=\sqrt{S} and \mathrm{d}L=\mathrm{d}S/(2\sqrt{S}), for several values of the S_c parameter (Fig.3), as well as for several values of the phase shift \mathrm{\Phi } (Fig.4). The spectra clearly show again the same power law dependence with −3/2 exponent as was reported in Refs. [12, 13] followed by a steep falloff sufficiently close to the critical point. What is interesting, however, is that the spectra for \mathrm{\Phi }=\pi /2 (and this is the phase shift value used for spiking solutions reported in Refs. [1, 2, 36]) recover even the fine structure of the scaling and clearly show the small bump at the end of the scale free part of the spectra where the local probability deflects from the initial −3/2 power exponent and flattens first before turning in to the steep falloff. These small bumps are evident in all experimental spectra [12, 13] shown in Fig.3 as well.

9 Conclusion

One of the properties of the WETCOW wave modes is that the anisotropy structure of brain conductivity as well as the structure of brain inhomogeneity favors their propagation in the outer regions of the cortex (see, for example, Fig.2 of [1, 2]). Neuronal avalanches are measured in the most external layer of the cortex and, usually, introducing the electrodes deeper in the cortical columns will eliminate the scale-free distributions. Therefore, it seems to be an interesting problem to check the whole-brain scale free distribution in the region of typical propagation of WETCOW wave modes. To do this numerical experiment we generated an ensemble of 10⁶ WETCOW modes distributed and randomly propagating through inhomogeneous and anisotropic cortical tissue. Fig.5 shows two randomly selected snapshots of wave mode trajectories that were generated using the procedure described in details in [1] and propagate in the surface-like 2D manner in the external layer of the cortex. Using the same procedure, that replicates the original experimental neuronal avalanche detection method, that is thresholding and then binning the wave signal into dots or “spikes”, we again see that the WETCOW modes show scale free behavior as shown in Fig.6.

In summary, in this paper we have presented an analysis of temporal and spatial probability density spectra that are generated due to the critical dynamics of the nonlinear weakly evanescent cortical wave (WETCOW) modes [1, 2]. The Hamiltonian framework developed for these WETCOW modes in Ref. [36] is advantageous in that it explicitly uncovers the reciprocity of the temporal and the spatial dynamics of the evolutionary equations. Therefore, in the nonlinear regime sufficiently close to the critical point the spatial behavior of the wave modes displays features similar to the properties of their nonlinear temporal dynamics that can be described as spatial domain spiking, with localized regions of wave activity separated by quiescent areas, with this spatial spiking intermittence increasing near the critical point. Similar spatial behavior was observed experimentally in neuronal avalanches, when activity detected at one electrode was typically skipped over the nearest neighbors. This was interpreted as evidence that avalanche spatial intermittency is not wave-like in nature [12, 13]. Our results demonstrate the contrary, however: the spatial patterns are the direct result of nonlinear interactions of weakly evanescent cortical waves.

Both temporal and spatial scaling expressions analytically estimated from the nonlinear amplitude/phase evolutionary equations show excellent agreement with the experimental neuronal avalanche probability spectra reproducing not only the general average power law exponent values and falloffs in the vicinity of the critical point, but also finding some very subtle but nevertheless clearly experimentally evident fine details, like bumps in the transition region at the edge of the scale free part of the probability spectra.

In a more general way these results may be applicable not only to neuronal avalanches but to many other physical systems that involve wave processes as they show that a system of wave modes interacting through all possible combinations of the third order nonlinear terms described by a general wave Hamiltonian necessarily produces anharmonic wave modes with temporal and spatial scaling properties that follow scale free power laws.

References

[1] V. L. Galinsky , L. R. Frank . Universal theory of brain waves: From linear loops to nonlinear synchronized spiking and collective brain rhythms. Phys. Rev. Research, 2020, 2: 023061

[2] V. L. Galinsky , L. R. Frank . Brain waves: Emergence of localized, persistent, weakly evanescent cortical loops. J. Cogn. Neurosci., 2020, 32(11): 2178

[3] G.Buzsaki, Rhythms of the Brain, Oxford University Press, 2006

[4] A. L. Hodgkin , A. F. Huxley . A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol., 1952, 117(4): 500

[5] R. FitzHugh . Impulses and physiological states in theoretical models of nerve membrane. Biophys. J., 1961, 1(6): 445

[6] J. Nagumo , S. Arimoto , S. Yoshizawa . An active pulse transmission line simulating nerve axon. Proceedings of the IRE, 1962, 50(10): 2061

[7] C. Morris , H. Lecar . Voltage oscillations in the barnacle giant muscle fiber. Biophys. J., 1981, 35(1): 193

[8] E. M. Izhikevich . Simple model of spiking neurons. IEEE Trans. Neural Netw., 2003, 14(6): 1569

[9] W.GerstnerW.M. KistlerR.NaudL.Paninski, Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition, Cambridge University Press, New York, NY, USA, 2014

[10] A. Kulkarni , J. Ranft , V. Hakim . Synchronization, stochasticity, and phase waves in neuronal networks with spatially-structured connectivity. Front. Comput. Neurosci., 2020, 14: 569644

[11] R. Kim , T. J. Sejnowski . Strong inhibitory signaling underlies stable temporal dynamics and working memory in spiking neural networks. Nat. Neurosci., 2021, 24(1): 129

[12] J. M. Beggs , D. Plenz . Neuronal avalanches in neocortical circuits. J. Neurosci., 2003, 23(35): 11167

[13] J. M. Beggs , D. Plenz . Neuronal avalanches are diverse and precise activity patterns that are stable for many hours in cortical slice cultures. J. Neurosci., 2004, 24(22): 5216

[14] D. Plenz , T. L. Ribeiro , S. R. Miller , P. A. Kells , A. Vakili , E. L. Capek . Self-organized criticality in the brain. Front. Phys. (Lausanne), 2021, 9: 639389

[15] N. Friedman , S. Ito , B. A. Brinkman , M. Shimono , R. E. DeVille , K. A. Dahmen , J. M. Beggs , T. C. Butler . Universal critical dynamics in high resolution neuronal avalanche data. Phys. Rev. Lett., 2012, 108(20): 208102

[16] D. R. Chialvo . Emergent complex neural dynamics. Nat. Phys., 2010, 6(10): 744

[17] J. M. Beggs , N. Timme . Being critical of criticality in the brain. Front. Physiol., 2012, 3: 163

[18] V. Priesemann , M. Wibral , M. Valderrama , R. Pröpper , M. Le Van Quyen , T. Geisel , J. Triesch , D. Nikolić , M. H. Munk . Spike avalanches in vivo suggest a driven, slightly subcritical brain state. Front. Syst. Neurosci., 2014, 8: 108

[19] B. Cramer , D. Stöckel , M. Kreft , M. Wibral , J. Schemmel , K. Meier , V. Priesemann . Control of criticality and computation in spiking neuromorphic networks with plasticity. Nat. Commun., 2020, 11(1): 2853

[20] A. J. Fontenele , N. A. P. de Vasconcelos , T. Feliciano , L. A. A. Aguiar , C. Soares-Cunha , B. Coimbra , L. Dalla Porta , S. Ribeiro , A. J. Rodrigues , N. Sousa , P. V. Carelli , M. Copelli . Criticality between cortical states. Phys. Rev. Lett., 2019, 122(20): 208101

[21] L. J. Fosque , R. V. Williams-García , J. M. Beggs , G. Ortiz . Evidence for quasicritical brain dynamics. Phys. Rev. Lett., 2021, 126(9): 098101

[22] P. Bak , C. Tang , K. Wiesenfeld . Self-organized criticality: An explanation of the 1/f noise. Phys. Rev. Lett., 1987, 59(4): 381

[23] P. Bak , C. Tang , K. Wiesenfeld . Self-organized criticality. Phys. Rev. A, 1988, 38(1): 364

[24] S. Zapperi , K. B. Lauritsen , H. E. Stanley . Self-organized branching processes: Mean-field theory for avalanches. Phys. Rev. Lett., 1995, 75(22): 4071

[25] K. Bækgaard Lauritsen , S. Zapperi , H. E. Stanley . Self-organized branching processes: Avalanche models with dissipation. Phys. Rev. E, 1996, 54(3): 2483

[26] C. W. Eurich , J. M. Herrmann , U. A. Ernst . Finite-size effects of avalanche dynamics. Phys. Rev. E, 2002, 66(6): 066137

[27] C. Bédard , H. Kröger , A. Destexhe . Does the 1/f frequency scaling of brain signals reflect self-organized critical states. Phys. Rev. Lett., 2006, 97(11): 118102

[28] J. Touboul , A. Destexhe . Can power-law scaling and neuronal avalanches arise from stochastic dynamics. PLoS One, 2010, 5(2): e8982

[29] J. Touboul , A. Destexhe . Power-law statistics and universal scaling in the absence of criticality. Phys. Rev. E, 2017, 95(1): 012413

[30] P. A. Robinson , C. J. Rennie , J. J. Wright . Propagation and stability of waves of electrical activity in the cerebral cortex. Phys. Rev. E, 1997, 56(1): 826

[31] D. P. Yang , P. A. Robinson . Critical dynamics of Hopf bifurcations in the corticothalamic system: Transitions from normal arousal states to epileptic seizures. Phys. Rev. E, 2017, 95(4): 042410

[32] P. A. Robinson , C. J. Rennie , D. L. Rowe . Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys. Rev. E, 2002, 65(4): 041924

[33] S.di SantoP.VillegasR.BurioniM.A. Muñoz, Landau–Ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization, Proc. Natl. Acad. Sci. USA 115(7), E1356 (2018), Available at:

[34] E. D. Gireesh , D. Plenz . Neuronal avalanches organize as nested theta- and beta/gamma-oscillations during development of cortical layer 2/3. Proc. Natl. Acad. Sci. USA, 2008, 105(21): 7576

[35] V. Buendía , P. Villegas , R. Burioni , M. A. Muñoz . Hybrid-type synchronization transitions: Where incipient oscillations, scale-free avalanches, and bistability live together. Phys. Rev. Research, 2021, 3(2): 023224

[36] V. L. Galinsky , L. R. Frank . Collective synchronous spiking in a brain network of coupled nonlinear oscillators. Phys. Rev. Lett., 2021, 126(15): 158102

[37] E. Ott , T. M. Antonsen . Long time evolution of phase oscillator systems. Chaos, 2009, 19(2): 023117

[38] I. V. Tyulkina , D. S. Goldobin , L. S. Klimenko , A. Pikovsky . Dynamics of noisy oscillator populations beyond the Ott−Antonsen ansatz. Phys. Rev. Lett., 2018, 120(26): 264101

[39] V.L. GalinskyL.R. Frank, Critically synchronized brain waves form an effective, robust and flexible basis for human memory and learning, 2022)

[40] L.FrankV.GalinskyJ.TownsendR.A. MuellerB.Keehn, Imaging of brain electric field networks, 2023)

[41] Y.Kuramoto, Chemical Oscillations, Waves, and Turbulence, Dover Books on Chemistry Series, Dover Publications, Incorporated, 2013

[42] H. Daido . Generic scaling at the onset of macroscopic mutual entrainment in limit-cycle oscillators with uniform all-to-all coupling. Phys. Rev. Lett., 1994, 73(5): 760

[43] J. D. Crawford . Scaling and singularities in the entrainment of globally coupled oscillators. Phys. Rev. Lett., 1995, 74(21): 4341

RIGHTS & PERMISSIONS

Higher Education Press