Critical brain wave dynamics of neuronal avalanches

1 Introduction

The complexity of oscillatory and wave patterns across a wide range of spatial and temporal scales of brain activity results in multiple independent models for these activity patterns. 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]. Crucial events, as a manifestation of criticality, have been discussed in [21] using Diffusion Entropy Analysis [22]. It was also hypothesized that an existence of crucial events facilitates information transmission in various states of brain functioning [23, 24].

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 that 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] that can be phenomenologically described using the ideas of the self-organized criticality [25, 26], and extended to the mean-field theory of the self-organized branching processes (SOBP) [27–29].

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 [30–32], 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 [33–36] 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 [36–38], but the usual agreement is that avalanches belong to the mean-field directed percolation universality class, which 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 directed percolation ones [20].

In this paper we show that the above important observational phenomenon, the so-called “neuronal avalanches”, which have been noted and studied for almost 2 decades, can be naturally explained by the wave activity model. Our recently described theory of weakly evanescent brain waves (WETCOW) originally developed in [1, 2] and then reformulated in a general Hamiltonian framework [39] 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 non-linear equations for both the temporal and spatial evolution of brain wave modes that includes all possible non-linear interaction between propagating modes at multiple spatial and temporal scales and degrees of non-linearity. 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.

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. Both temporal and spatial scaling expressions analytically derived from non-linear amplitude/phase evolutionary equations show excellent agreement with the experimental neural avalanche probability spectra reproducing not only 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 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 and 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 non-linear wave dynamics that is responsible for the wide range of activities in the brain tissue, ranging from the linear coherently propagating waves to the non-linear incoherent asynchronous spiking, including in between the peculiar power law-like coherence of the neuronal avalanches.

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 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 non-linearity that characterizes anharmonicity of the non-linear wave modes due to non-resonant 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 [36], sometimes even including inter-mode coupling [38], or are generated from some ad hoc artificial distributions [40], and require the addition of stochastic noise properties [41], etc. This emphasizes generality and importance of our derivation.

2 Weakly evanescent brain waves

Beginning from our non-linear Hamiltonian formulation of the WETCOW theory [39], we have for an anharmonic wave mode

where a is a complex wave amplitude and a^† is its conjugate. The amplitude a denotes either temporal a_k(t) or spatial a_ω(x) wave mode amplitudes that are related to the spatiotemporal wave field ψ(x, t) through a Fourier integral expansions.

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 spatiotemporal wave field ψ(x, t) is a superposition of multiple waves, that may include neuronal firings, membrane sub-threshold oscillations, LFPs, etc. The frequency ω and the wave number k of the wave modes satisfy dispersion relation D(ω, k) = 0, and ω_k and k_ω 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 Γaa^† in Eq.1 denotes the harmonic (quadratic) part of the Hamiltonian with either the complex valued frequency Γ = iω + γ or the wave number Γ = ik + λ that both include a pure oscillatory parts (ω or k) and possible weakly excitation or damping rates, either temporal γ or spatial λ. The second anharmonic term is cubic in the lowest order of non-linearity and describes the interactions between various propagating and non-propagating wave modes, where α, β_a and ${\beta }_{a^{†}}$ are the complex valued strengths of those different non-linear processes. As it was shown in [1, 2] the inverse proportionality of frequency and wave number in the dispersion relation 6) results in the third order expressions for the non-resonant coupling between multiple waves.

Distribution of various charges in brain tissue, including free ionic charges (sodium, potassium, etc), bonded macromolecular volume charges, membrane polarization and/or surface charges, etc., determines brain electrodynamics. The voltages and currents measured in real brains are produced by those electrodynamic processes that in the most general form can be represented by the Maxwell equations together with state or motion equations for the brain matter, particularly by a generalized Ohm’s law, that describes electro-diffusive flow of charged particles through inhomogeneous media (that may include both concentration and voltage gradients). The neuron action potential itself is nothing more than propagating non-linear electrostatic wave described by the same electrodynamics formalism. A set of derivations that lead to this description was presented in details in [1, 2, 39] 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

Using the electrostatic potential E = −∇ψ, Ohm’s law J = σ ⋅E (where σ ≡ {σ_ij} is an anisotropic conductivity tensor), a linear electrostatic property for brain tissue D = ɛE, assuming that the scalar permittivity ɛ is a “good” function (i.e. it does not go to zero or infinity everywhere) and taking the change of variables ∂x → ɛ∂x′, the charge continuity equation for the spatial-temporal evolution of the potential ψ can be written in terms of a permittivity scaled conductivity tensor Σ = {σ_ij/ɛ} as

where we have included a possible external source (or forcing) term $\mathcal{F}$. For brain fiber tissues the conductivity tensor Σ might have significantly larger values along the fiber direction than across them. The charge continuity without forcing i.e., $(\mathcal{F}=0)$ can be written in tensor notation as

where repeating indices denote summation. Simple linear wave analysis, i.e. substitution of ψ ∼ exp(−i (k ⋅ r − Ωt)), where k is the wavenumber, r is the coordinate, Ω = ω + iγ is the frequency and t is the time, gives the following complex dispersion relation:

which is composed of the real and imaginary components:

Although in this general form the electrostatic potential ψ, as well as the dispersion relation D(Ω, 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_ω(x) wave mode amplitudes can be used to define the time dependent wave number energy spectral density Π_k(t) or the position dependent frequency energy spectral density Π_ω(x) for the spatiotemporal wave field ψ(x, t) as

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 [39] as

where the single mode amplitude a_n again denotes either a_k or a_ω, a ≡ {a_n} and $r_{nm}=R_{nm}{e}^{i{\mathrm{\Delta }}_{nm}}$ is the complex network adjacency matrix with R_nm providing the coupling power and Δ_nm taking into account any possible differences in phase between network nodes. This description includes both amplitude $\mathfrak{R}(a)$ and phase $\mathfrak{I}(a)$ mode coupling and as shown in [39] 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 provide 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 [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 used for expansion of the spatiotemporal wave field ψ into a set of wave modes imply 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 [39]. Thus we will consider an ensemble of non-interacting 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 (Γ, β_a, ${\beta }_{a^{†}}$, α) 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₀ ≤ k ≤ k₁ or an arbitrary frequency ω from a range of wave frequencies ω₀ ≤ ω ≤ ω₁.

3 Single anharmonic mode criticality

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

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. As frequencies and wave numbers for linear waves satisfy the dispersion relation 6), they are related and the same Hamiltonian expression 1) can be used either for temporal a_k(t) or spatial a_ω(x) wave amplitudes. Therefore, we note that although (10) is an equation for the temporal evolution, the spatial evolution of the mode amplitudes a_ω(x) can be described by a similar equation substituting temporal variables by their spatial counterparts, i.e., (t, ω, γ) → (x, k, λ).

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

gives the set of equations.

where Ψ_a ≡ ϕ − δ_a, ${\mathrm{\Psi }}_{a^{†}}\equiv \varphi -{\delta }_{a^{†}}$.

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

so that.

where.

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

where Φ = Φ_a − Φ_ϕ.

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

as ϕ_e = ϕ₀ ≡ const and ${\tilde{A}}_e=\tilde{\omega }/R_{\varphi }\cos {\varphi }_0\equiv$ const. This shows that for α > R_a there exist critical values of ω̃ and A_e, where the stationary solution disappears and is replaced by non-linear oscillations, such that.

which can also be expressed in terms of critical value of one of the unscaled variables, either ω or γ

This stationary solution provides the locus of the saddle node on an invariant circle bifurcation point at where the non-linear spiking oscillations occur (as was shown both in [1, 2] and in [39]).

4 Effective spiking rate

The effective period ${\mathcal{T}}_s$ of spiking solutions of (21) and (22) (or its inverse–either the firing rate 1/${\mathcal{T}}_s$ or the effective firing frequency $2\pi /{\mathcal{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 ϕ) as

giving the unscaled effective spiking period T_s and the effective firing frequency ω_s

with the periodic amplitude $\tilde{A}$ reaching the maximum ${\tilde{A}}_{\mathit{max}}=1/(\alpha -R_a)$ and the minimum ${\tilde{A}}_{\mathit{min}}=1/(\alpha +R_a)$ for $d\tilde{A}/d\tau =0$ when ϕ = Φ and ϕ = Φ + π respectively.

The expressions (28) and (29) are more general than typically used expressions for the scaling exponent in the close vicinity |γ − γ_c|≪ γ_c of the critical point [42–44]. They allow recovery of the correct T limits both at γ → γ_c with the familiar $T\sim 1/\sqrt{{\gamma }_c-\gamma }$ scaling and at γ ∼ 0 with the period T approaching T₀ as T ∼ T₀ + O (γ²) ≡ 2π/ω + O(γ²), where T₀ is the period of linear wave oscillations with the frequency ω. In the intermediate range 0 < γ < γ_c the expressions (28) and (29) show reasonable agreement (Figure 1) with peak–to–peak period/frequency estimates from direct simulations of system (12) and (13).

FIGURE 1. [Left] Comparison of the analytical expression (29) for the effective spiking frequency ω_s =2π/T_s (red) and the frequency estimated from numerical solution of Eq. 21 and Eq. 22 (blue) as a function of the criticality parameter γ/γ_c. In the numerical solution only γ was varied and the remaining parameters were the same as parameters reported in [39] [Right] Spiking solutions for typical parameters producing temporal ((21) and (22), red) and spatial ((38) and (39), blue) spiking profiles where some signal of width δt_s or δl_s was detected and surrounded by quiet area with the total effective period T_s or L_s.

5 Temporal probability of single spike detection

As the periodic solution of Eq. 21 and Eq. 22 in the 0 ≪ γ ≪ γ_c range looks like linear waves at γ close to zero, but transforms to spike as γ increases, we can approximate the probability of detecting a single spike by a ratio of a spike peak duration (recorded above some threshold) to a total peak-to-peak time. Taking into account that the initial phase of spiking solutions of Eq. 21 and Eq. 22 is a random variable uniformly distributed on [0, 2π] interval, the probability that a spike (either positive or the more frequently experimentally reported negative [12, 13]) with duration width δt_s and with the total period between the spikes (T_s) will be detected is simply δt_s/T_s–where the distance between spikes is determined as the time interval needed for 2π radian phase change, that is the effective spiking period T_s. Assuming initially that the spike width δt_s does not change when approaching the critical point ω_c, δt_s can be approximated by some fixed fraction of the linear wave period, i.e., δt_s ∼ π/ω, that gives for the probability density

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 Π_ω(x) (or with the power spectral density).

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

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

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 [12, 13].

6 Multi-mode avalanche probability

The above single wave mode analysis shows that the probability density $P_k^{\{T\}}$ for any arbitrary selected wave mode k with arbitrary chosen threshold follows a power law distribution with -2 exponent, therefore, a mixture of multiple wave modes that enters into the spatiotemporal wave field ψ(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. Figure 2 compares the avalanche distribution for all wave modes with the -2 exponent.

FIGURE 2. (Top) The avalanche durations distribution for all wave modes compared with the -2 exponent (Bottom) WETCOW modes randomly distributed and propagated on a 1,000 by 1000 grid. Two examples of temporal signal snapshots with different values of signal threshold are shown (color pallet encodes the change of frequencies from the smallest (blue) to the largest (red). Localized regions of wave activity in the spatial domain are clearly evident.

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^{\prime }\le T+\mathrm{\Delta }T}^a$ that an avalanche happens at any time between 0 and T + ΔT, where ΔT is some small binning interval used by the above experimental method, can be expressed as

Since the probability

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

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

An expansion of Eq. 35 in the leading order of ΔT gives for the avalanche probability density p^a(T)

that is the avalanche probability density p^a(T) shows the same T⁻² 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

7 Spatial spike detection probability

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

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

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

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 non-linear equations suggests that the non-linear 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 non-linear waves was presented as an 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 Eq. 38 and Eq. 39 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 ω as.

where L is the linear spatial scale related to the wave number as k = 2π/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

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 non-linear anharmonic oscillations described by the (21) and (22) only exists for frequencies and wave numbers that are above the critical frequency ω_c or the critical wave number k_c values that define maximal possible temporal T_c or spatial L_c scales of the non-linear 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 Eq. 1, and hence in the equations Eq. 21 and Eq. 22 or (38) and (39), 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 [36, 38]. 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 levels 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 do not require to know the details of the coupling between nodes, essentially viewing all nodes as completely non-interacting. 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 [38]. 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 δt_s used in Eq. 30 and 32 (or the spike length for spatial spiking in Eq. 43 and Eq. 44) 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 Eq. 21 but the spatial analysis with equation Eq. 38 is exactly the same).

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

where τ_± = τ_max ± δτ, and τ_max is the location of spiking peak. Neglecting the spike shape asymmetries, i.e., assuming that τ_± correspond to symmetric changes in both the amplitudes ${\tilde{A}}_{\pm }=\tilde{A}({\tau }_{\pm })={\tilde{A}}_{\mathit{max}}-\delta \tilde{A}$, and the phases Φ_± = Φ(τ_±) = Φ ± δΦ, we can then estimate the spike duration δt_s ≡ (τ₊ − τ₋)/γ as

where, similar to the spiking period estimation in Eq. 28, 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 δΦ some fixed value that is smaller or around a quarter of the period (i.e., δΦ ≲ π/2) can be chosen, and R = ˜ω_cR_a/R_ϕ.

An expression (48) 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 δl_s/L_s again substituting $L=\sqrt{S}$ and $dL=dS/(2\sqrt{S})$, for several values of the S_c parameter (Figure 3), as well as for several values of the phase shift Φ (Figure 4). The spectra clearly show again the same power law dependence with -3/2 exponent as was reported in [12, 13] followed by a steep falloff sufficiently close to the critical point. What is interesting, however, is that the spectra for Φ = π/2 (and this is the phase shift value used for spiking solutions reported in [1, 2, 39]) 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 Figure 3 as well.

FIGURE 3. [Left] Analytical probability density spectra as a function of brain waves criticality parameter S/S_c show excellent agreement with the experimental avalanche data [Right, from [12, 13]] reproducing not only the overall shape of the spectra with the -3/2 power exponent at the initial scale free part of the spectra and the steep falling edge in the vicinity of the critical point, but also reproduce the fine details such as the small bump-like flattening of the spectra at the transition from -3/2 leg to the steep falling edge that is clearly evident in experimental spectra.

FIGURE 4. Analytical probability density spectra multiplied by a ${(S/S_c)}^{3/2}$ as a function of brain waves criticality parameter S/S_c plotted for several values of the phase shift Φ.

9 Conclusion

Brain activity in general and neuronal avalanches in particular show an abundance of very complex and strangely organized activity patterns. Understanding the nature and the origin of cascades of synchronized activity in the cortex has multiple implications to understanding of organization of cortical functioning. Although originally neuronal avalanches were detected in vitro using multi-electrode arrays in 2D slices of cultivated cortex cultures [12, 13], there are now multiple experimental data of in vivo avalanche recordings [45–47] involving optical recordings as well [48].

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, Figure 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. Figure 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 Figure 6.

FIGURE 5. Examples of complete wave mode trajectory snapshots for two randomly selected parameters and initial conditions. The trajectories was randomly selected from an ensemble of 10⁶ WETCOW modes used for generation of probability distributions of Figure 6.

FIGURE 6. Plots of spatial (A) and temporal (B) probability density spectra obtained by binning oscillatory signal of ensemble of 10⁶ WETCOW modes randomly distributed and propagated through cortical tissue. Two examples of temporal signal (dots or “spikes”) snapshots with different values of signal threshold are shown in (C) and (D) (color pallet encodes the change of frequencies from the smallest (blue) to the largest (red).

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 non-linear weakly evanescent cortical wave (WETCOW) modes [1, 2]. The Hamiltonian framework developed for these WETCOW modes in [39] is advantageous in that it explicitly uncovers the reciprocity of the temporal and the spatial dynamics of the evolutionary equations. Therefore, in the non-linear regime sufficiently close to the critical point the spatial behavior of the wave modes displays features similar to the properties of their non-linear 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 non-linear interactions of weakly evanescent cortical waves.

Both temporal and spatial scaling expressions analytically estimated from the non-linear 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 non-linear terms described by a general wave Hamiltonian necessarily produces anharmonic wave modes with temporal and spatial scaling properties that follow scale free power laws.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

VG and LF developed the theoretical formalism, performed the analytic calculations and performed the numerical simulations. Both VG and LF contributed to the final version of the manuscript.

Funding

LF and VG were supported by NIH grant R01 AG054049.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References