Wave Turbulence and Energy Cascade in the Hippocampus

1. Introduction

Hebb's (1958) hypothesis that no psychological function can be attributed uniquely to any segment of cortex has the profound implication that cognition emerges from the coordination of activity across all spatial and temporal scales of the brain (Lashley, 1958; Allen and Collins, 2013), and that understanding the brain begins with studying the nature and role of different scales of brain activity.

Temporal scales, defined based on the frequency structure of extracellular recordings (local-field potential, LFP), are typically more accessible to observations. Their relation to spatial scales is not exactly known, but they are assumed to be in a monotonic relation to spatial scales (lower frequencies correspond to larger populations; Buzsáki and Draguhn, 2004). The Fourier spectra of hippocampal LFP recordings cover a frequency range approximately from 0.05 to 500 Hz, with spectra generally decaying as a power law. Two scales ranges are readily identified. The high end of the spectrum, say, f > 200 Hz¹ represents the microscale, mainly occupied by action potential activity of single neurons, fast synaptic time constants, ion channel opening and closing, and heat dissipation (energy sinks)². These processes are generated by microscopic neural units such as individual neurons, pairs of excitatory/inhibitory neurons, or small neuronal sequences (e.g., Lorente de No, 1938). The low end of the frequency spectrum (say, f < 60 Hz) is assumed to represent macroscopic processes encompassing several segments of the brain. For example, the theta rhythm is observed across hippocampus, entorhinal cortex, hypothalamus, prefrontal cortex, and others; e.g., Vertes and Kocsis, 1997; Buzsaki, 2002; Siapas et al., 2005). It is associated with active exploration and REM sleep, and is assumed to provide the temporal structure for the organization of local networks (Green and Arduini, 1954; Green and Machne, 1955; Vanderwolf, 1969; Lisman and Idiart, 1995; Buzsaki, 2002).

In addition to these two scales, cortical activity exhibits an intermediate (meso-) scale (e.g., Freeman, 2000a; Lubenov and Siapas, 2009; Patel et al., 2012; Muller et al., 2018; Zhang et al., 2018), which spans the frequency band 60 < f < 200 Hz, a range of oscillations collectively referred to as the gamma rhythm (Buzsáki et al., 1983; Bragin et al., 1995; Belluscio et al., 2012; Lasztóczi and Klausberger, 2014; Schomburg et al., 2014). Mesoscopic processes have been observed within brain regions associated with higher cognition, such as the neocortex and the hippocampus, that exhibit anisotropic and homogeneous mesoscopic structure³. The monotonic relation between temporal and spatial scales places mesoscopic spatial scales in the order of a few millimeters in rats, suggesting that they represent collective neuronal activity (i.e., synchronized neuronal firing; Buzsáki and Draguhn, 2004)⁴. Rhythms in the 8–200 Hz range have the spatial structure of propagating perturbations, i.e., waves (Petsche and Stumpf, 1960; Lubenov and Siapas, 2009; Patel et al., 2012, 2013; Muller et al., 2018).

Collective action responds to task behavior and intensity of behavior activity. During spatial exploration, theta and gamma increase their power and develop measurable phase coupling in response to intensity of activity (Whishaw and Vanderwolf, 1973; Morris and Hagan, 1983; Bragin et al., 1995; Chorbak and Buzsaki, 1998; Chen et al., 2011; Ahmed and Mehta, 2012; Kemere et al., 2013; Sheremet et al., 2016; Zheng et al., 2016). This suggests a relation to cognition that is particularly intriguing, because isotropic and homogeneous mesoscopic neuronal structures are the type of environment that would favor collective action over that of microscopic elements. While mesoscopic collective activity could be explained away as a marginally-significant synchronization effect of strongly non-linear microscopic units (Buzsaki, 2006), its strong relation to isotropic and homogeneous networks and its connection to behavior suggest that it might be in fact the function of this type of network. Freeman and Vitiello (2010) hypothesize, citing Lashley (1942), that mesoscale processes reflect the essential cognition step of abstraction and generalization of a particular stimulus to a category of equivalent inputs, “because they require the formation of non-local, very large-scale statistical ensembles (our emphasis)”⁵.

We therefore conjecture that mesoscopic action plays a crucial role in the integration of brain activity across scales, and that its stochastic character is essential for the fulfillment of this role. This conjecture changes significantly the interpretation of LFP recordings and possibly other observations of of hippocampal activity. For example, rather than attributing special significance to the activity of microscopic elements (like musical scores played by separate instruments), LFP recordings should be interpreted as local observations of a stochastic process representing collective action by mesoscale ensembles involving a large number of microscopic units, in which the precise state of a single unit is irrelevant.

Despite a few brilliant insights into collective action dynamics (e.g., Wilson and Cowan, 1973; Wright and Liley, 1995; Freeman, 2000a, 2006, 2007; Freeman and Vitiello, 2010; Cowan et al., 2016), a consistent theoretical approach to mesoscopic dynamics is missing. It is important to note that, because of the scale separation, mesoscopic and microscopic dynamics are different, therefore the wealth of knowledge accumulated about microscopic physics cannot be directly extended to mesoscopic processes. The notion that theory is scale-dependent is common in physics: macroscopic systems are characterized by state variables and laws that are typically inaccessible directly to microscale theories. An example is Boltzmann's celebrated H-theorem, which introduces the entropy (H) as a new state variable, and elucidates the process through which time-reversible microscopic dynamics begets the macroscopic irreversibility of evolution toward equilibrium (Gibbs, 1902; Khinchin, 1949; Toda et al., 1983; Pathria and Beale, 2011)⁶.

Here, we propose the theory of weak turbulence as a framework for the description of mesoscopic collective action and its energy balance in the hippocampus. The paper is organized as follows:

(i) In section 2, we discuss briefly data collection and analysis procedures and provide some guidance on reading the bispectral maps used to estimate non-linear coupling between Fourier components of LFP recordings.

(ii) In section 3.1, we examine stochastic features LFP recordings, using second and third order correlators (spectra and bispectra). The goal of this section is to highlight important characteristics of collective action and their dependency on energy input into the hippocampus. This perspective is important because it captures the transformation (ultimately, the time evolution) of the LFP statistical measures, thus providing a dynamical view of mesoscopic activity. In particular, the evolution of energy distribution over spatial or temporal scales is relevant for describing the energy balance of collective action under non-linear interactions.

(iii) Section 3.5 introduces and discusses the theory of weak turbulence as a framework for mesoscopic collective action. Originally formulated for hydrodynamics, the theory of turbulence has expanded in scope through the work of Richardson, Kolmogorov, and Zakharov Richardson (1922), Zakharov et al. (1992), and Kolmogorov (1941) to become the theoretical foundation of physical disciplines ranging from plasma physics, to non-linear optics, Bose-Einstein condensation, water waves, coagulation-fragmentation processes, and many others (Zakharov et al., 1992; L'vov, 1998; Nazarenko, 2011). At the center of the turbulence theory is the study of internal energy processes in non-linear, multi-scale systems with a large number of components. The non-linear character of the system allows for scales to interact, creating the conditions for cross-scale flows of energy and other conserved quantities. Non-linearity implies interaction connects evolution across scales, allowing for a cross-scale flux of energy called the turbulent cascade (Richardson, 1922; Kolmogorov, 1941). The turbulent cascade is the fundamental property defining turbulence.

(iv) In section demonstrate the capabilities of the theory by using the three-wave equations (a simplified, universal non-linear interaction modeling framework) a to re-evaluate the significance of linear and non-linear structure of hippocampal LFP recordings.

(v) The significance of the turbulence theoretical framework for understanding collective action, and its implications for cognition, are discussed in section 4.

Further details of the weak turbulence formalism are given in the Supplementary Material.

2. Materials and Methods

2.1. Subjects and Behavioral Training

Rat r539♂-maurer used in this study belongs to a cohort of 4–9 months old Fisher344-Brown Norway Rats (Taconic; see e.g, Zhou et al. 2018 for additional information about the cohort). The methods detailing the collection for these rats have been described in detail elsewhere (Zhou et al., 2018). Briefly, following acclimation to the University of Florida colony, animals were trained to traverse a circular track for food reward (45 mg, unflavored dustless precision pellets; BioServ, New Jersey; Product #F0021). During this time, their body weight was slowly reduced to 85% to their ad libitum baseline. Rats were surgically implanted with a custom single shank silicon probe from NeuroNexus (Ann Arbor, MI), designed such that thirty-two recording sites spanned across multiple hippocampal lamina.For probe preparation instructions and surgical methods, see (Vandecasteele et al., 2012; Zhou et al., 2018).

All behavioral procedures were performed in accordance with the National Institutes of Health guidelines for rodents and with protocols approved by the University of Florida Institutional Animal Care and Use Committee.

2.2. Neurophysiology

Following recovery from surgery, the rat was retrained to run unidirectionally on a circle track (outer diameter: 115 cm, inner diameter: 88 cm), receiving food reward at a single location. Following a few days of circle track running, the rats were trained to run on a digital-8 maze (121 × 101 cm, L × W). During these sessions, the local-field potential was record on a Tucker-Davis Neurophysiology System (Alachua, FL) at ~ 24 kHz (PZ2 and RZ2, Tucker-Davis Technologies). The animals position was recorded at 30 frames/s with a spatial resolution <0.5 cm/pixel.

2.3. Spectral Analysis

The spectral analysis of the LFP in the current study was based on standard techniques used for stationary signals Priestley, 1981; Papoulis and Pillai, 2002. Descriptions of the stochastic estimators, meaning, normalization procedures, as well as how to interpret the bispectral maps, are given in Hasselmann et al. (1963), Rosenblatt and Van Ness (1965), Swami et al. (2001), Elgar (1987), Elgar and Guza (1985), and Sheremet et al. (2016), and many others. Here, we remind the reader the main definitions and terminology.

Assume the LFP recordings g(t) and h(t) are realizations of zero-mean stochastic processes, stationary in the relevant statistics, with Fourier transforms G(f_n) and H(f_n), n = 1, …, N. The second and third order spectral statistics are estimated using cross-spectrum and bispectrum, defined as

where the angular brackets denote the ensemble average, the asterisk denotes complex conjugation. We will generally omit the superscript when a single time series is involved. The diagonal $S_n^{g,g}$ of the cross-spectrum matrix are power spectra. The coherence and phase lag of time series g and h are the normalized modulus and phase of the cross-spectrum,

The cross-spectrum matrix provides information about the degree of correlation and phase lags for between different time series; spectra describe the frequency distribution of the variance of processes g and h, i.e., a complete characterization of the average linear structure of the Fourier representation.

The bispectrum provides information about the phase correlations between different frequency components of the same time series (e.g., Sheremet et al., 2016; Kovach et al., 2018). The bispectrum is statistically zero if the Fourier coefficients are mutually independent, i.e., for a linear system, and will exhibit peaks at triads (f_n, f_m, f_n+m) that are phase correlated. The real and imaginary part of the bispectrum are related to the skewness ${S}$ (e.g., positive skewness corresponds to sharp barrow peaks and flat troughs) and asymmetry ${A}$ (e.g., positive asymmetry corresponds to the front of the wave being steeper than the back, similar to a saw-tooth wave) of the time series g through

where σ is its standard deviation (Haubrich and MacKenzie, 1965; Masuda and Kuo, 1981). The bicoherence b_mn and biphase Φ_mn are defined in way similar to the coherence and phase lag as the normalized modulus and the argument of the bispectrum, that is

The bispectrum definition 2 has the following symmetries: 1) B_mn = B_nm; 2) $B_{mn}=B_{-m,-n}^{*}$; and 3) $B_{m,-n}=B_{nq}^{*}$, where q = m − n, showing that the information contained in the entire discretized plane (f_n, f_m) is redundant. Symmetry 1) implies that the bispectral distribution in octants 1 is 2 is symmetric with respect to the first diagonal. Symmetry 2) shows that quadrants 1 and 3 contain equivalent, complex conjugate, information, and the same is true for quadrants 2 and 4. Finally, symmetry 3) shows that octants 1 and 8 contain equivalent information. Therefore, the smallest domain of non-redundant information is octant 1, i.e., the area bounded by the positive frequency axis and the first diagonal (pink in Figure 1). In addition, if the number of frequencies is finite n = 1, … , N, as is the case in all numerical applications (e.g., N = 12 in Figure 1), the bispectrum is only defined in the area below the second diagonal, because f_m + f_n ≤ f_N, hence the triangle appearance of the bispectral plots.

2.4. LFP Power

Because the quantity measured by LFP observations is the potential of the electromagnetic field generated by the synaptic pulses, the spectral distribution of LFP variance (i.e., LFP spectral density integrated over some frequency interval) is proportional to the energy of the electromagnetic field per unit volume. Indeed, if the LFP time series g(t) is a zero-mean, weakly non-linear stochastic process, stationary in the second order statistics (Priestley, 1981; Percival and Walden, 2009), one can show that the discrete Fourier representation has the property

where ${\sigma }_n^2=\frac{1}{2}\langle |a_n{|}^2\rangle =S_n^{gg}\Delta f$ is the variance of Fourier mode n, a_n is its amplitude, $S_n^{gg}$ is the spectral density (Equation 1), and Δf is the frequency band width. The mean energy per unit volume of an electromagnetic wave with the amplitude of the electric potential g is w = εg², where ε is the dielectric constant of the medium (e.g., Pollack and Stump, 2002; Nolting, 2016). The energy per unit volume of a stochastic electromagnetic field is therefore

In other words, the amplitude squared of the Fourier components of the LFP are proportional to the energy stored in the unit volume by that particular Fourier component. A large number of algorithms have been developed for estimating for second- and higher-order statistics (power spectral density, bispectra, skewness, asymmetry, etc) of such processes. The analysis of the LFP in the current study was based on standard techniques used for variance-stationary signals (Priestley, 1981; Papoulis and Pillai, 2002) as previously described in Sheremet et al. (2016).

2.5. Numerical Implementation

All data analysis was performed in Matlab^® (MathWorks, Natick, MA, USA) using in-house developed code, as well as code imported from the HOSAtoolbox (Swami et al., 2001) for higher order spectral analysis.

Hippocampal layers were determined from the location of current sources and sinks derived on ripple and theta events, gamma power, and the polarity of the sharp-wave (Buzsáki, 1986; Buzsáki et al., 1986; Bragin et al., 1995; Ylinen et al., 1995; Lubenov and Siapas, 2009; Fernández-Ruiz et al., 2017). Details are given in Zhou et al. 2018.

The rat speed was calculated as the smoothed derivative of position. Raw LFP records sampled at 24 kHz (Tucker-Davis system) were pre-processed by applying a 2-kHz low-pass filter and divided into segments of 2048 time samples (approx. 1 s). Spectra and bispectra were classified by speed by averaging the speed over each of the 1-s LFP segments.

The dynamical and kinetic three-wave equations were integrated using the ODE solvers provided by Matlab^®. The implementation is trivial, therefore, to save space, the codes are not provided. The authors will, however, gladly share them upon request.

3. Results

3.1. Observations: Collective-Action Response to Behavior

The relationship between rat speed and theta and gamma power (Whishaw and Vanderwolf, 1973; Morris and Hagan, 1983; Chen et al., 2011; Ahmed and Mehta, 2012; Kemere et al., 2013; Zheng et al., 2015) is easily observable for speeds roughly above 5 cm/s (e.g., Figure 2). While low speeds, say < 3 cm/s may represent behavior uncorrelated to movement (Zhou et al., 2018), higher speeds exhibit a monotonic correlation with the power of hippocampal activity (firing rates of hippocampal and entorhinal neurons increase with speed; McNaughton et al., 1983; Rivas et al., 1996; Shen et al., 1997; Hirase et al., 1999; Maurer et al., 2005; Kropff et al., 2015). This observation provides a monotonic ordering of LFP statistics, with speed playing the role of an order parameter. Because of the relation is monotonic, the expression “evolution with speed” may be used unambiguously for “change as speed increases” (similar to the standard expression “time evolution” for change as time the time parameter increases). A relation to time exists, since speed itself is fundamentally a function of time: for example, an nominal increase of the speed parameter is in fact associated with time instances of acceleration.

For simplicity, we limit the data used here for illustration to a single representative source, rat r539♂-maurer. A discussion of the consistency of theses trends across the available, small (but growing) population of rats is presented elsewhere (Sheremet et al., 2018; Zhou et al., 2018).

3.2. Variance Spectra as a Function of Speed

Here, we examine estimates of spectral density of LFP traces recorded in the str. pyramidale (CA1.pyr), radiatum (CA1.rad), lacunosum moleculare (LM) and upper blade of the dentate gyrus (DG), indexed by rat speed.

Spectral densities show a weak, but significant variability as a function of speed and layer (Figure 3). At lowest discernible speeds (v ≲ 3 cm/s), all hippocampal layers exhibit non-trivial baseline, lowest-variance spectrum, that can be characterized in general as having a power-law shape f^−α over the entire range of approximately 6–300 Hz, with the exception of the DG layer, which exhibits a two-slope shape with a break point in the neighborhood of 50 Hz. The absolute value α of the power-law exponent is loosely referred to as “spectral slope”. In the CA1 layers the spectral slope is α ≈ 2 with slight variations (Figure 3). In the DG, the slope of the lower frequency range f ≲ 50 Hz visibly smaller, α ≈ 1. At high speeds (v > 15 cm/s), the theta rhythm and its harmonics dominate the lower frequency band of the spectrum, and low-power peak appears in the gamma range between 50 Hz and 120 Hz. Overall, the slope of the spectrum decreases. The CA1 layers have similar slopes 1.4 < α < 1.7, with the DG layer again standing out at α ≈ 0.7. The power in the high-frequency tail of the spectrum (f > 200 Hz) also increases significantly. In agreement with prior research (Bragin et al., 1995), the CA1.pyr layer shows the least energetic gamma range of all the layers examined.

Details of the variability of hippocampal LFP spectra with speed are shown in Figure 3 for seven speed intervals. According to Figure 2, speed ordering is statistically equivalent to ordering by total LFP variance. The spectra are normalized by dividing them by f^−α, using the high-speed slopes 1.4, 1.7, 1.7, and 0.7 for CA1. Pyr, CA1.Rad, LM, and DG, respectively. The normalization reduces the spectra to ≈ 1 in the frequency range where it agrees to the power law, and highlights spectral peaks. Although the ordering parameter is total variance, the spectra show a remarkable monotonic ordering in all frequency bands, with the exception of the highest two speed intervals, where the evolution stagnates and perhaps reverses slightly. The normalization re-scaling highlights several features of the evolution as a function of speed. The spectra in all layers tilt as energy increases (e.g., the LM slope changes from 2.1 to 1.7). Theta and its harmonics become dominant in the low-frequency range: four peaks at multiples of 8 Hz may be seen clearly in the CA1.Pyr and DG spectrum, perhaps three in the CA1.Rad and LM spectra. Also remarkable is the peculiar way the gamma range evolves: rather than developing some broad peak in 50–120 Hz range, the gamma range growth seems to result from the s ≈ 1 domain (where s is the normalized spectrum) progressively extending into the higher frequency range, while a “bump” develops in the neighborhood of f ≈ 100 Hz. This evolution is suggestive of a “front” of energy that propagates “against some resistance” toward higher frequency. This type of evolution may be seen in all layers. Although one might expect slight differences between acceleration and deceleration states for the same speed, these are negligible (data not shown). Therefore, the classification by speed does not differentiate between acceleration sign: the transition, say, between 10 cm/s and 20 cm/s may includes both acceleration and deceleration. Thus, the transition process described by Figures 2, 4 (and also below) is reversible, therefore quasi-stationary: all states described by these spectral may be assumed quasi-equilibrium states (in other words, imagining this evolution as a number of time steps, the transition between different steps is slow enough to allow the system to reach equilibrium and every step). This is consistent with the ability of the rat to control its speed. We will return to these ideas below.

3.3. Higher-Order Spectra as a Function of Speed

Higher-order spectra provide information about cross-frequency coupling. The bispectrum, the lowest order (and hence, the most “accessible” such estimator) has been used for a long time in wave dynamics (Hasselmann et al., 1963; Rosenblatt and Van Ness, 1965; Coppi et al., 1969). The relationship between its structure and third order statistics of the time series is well-understood (e.g., Haubrich and MacKenzie; Masuda and Kuo; Elgar; for bispectral definitions and terminology, see section 3.3). Although bispectral analysis is not common in neuroscience, recent work (Kovach et al., 2018) has shown that similar, widely used estimators for phase-amplitude and amplitude-amplitude coupling are, in fact, particular implementations of bispectral estimators (containing additional restrictive assumptions that make them susceptible to misinterpretation; e.g., Hyafil, 2015). To save space, we discuss only low (v < 10 cm/s) and high (v > 35 cm/s) speed levels; intermediate levels (not shown) represent a relatively smooth transition between these two limits. Because of the substantial difference in dynamics between CA1 and DG, we confine our discussion of bispectra to the LM layer.

As previously reported (Sheremet et al., 2016), bispectral estimates also show significant variability with rat speed. Low-speed bispectra are statistically zero (Gaussian) overall (Figure 5). The bicoherence exhibits a weak peak at (f_θ, f_θ, 2f_θ), corresponding to the phase coupling between theta (f_θ ≈ 8 Hz) and a relatively broad band around its 2f_θ harmonic (see arrow, Figure 5A, bicoherence). The phase of the harmonic band is in part in opposition and part in quadrature with theta (see transition red to blue in Figure 5, biphase), thus contributing overall to negative skewness of the LFP (peak is negative in 5a, skewness), but positive asymmetry (5a, asymmetry; see definitions in section).

In contrast, high speed activity show rich phase-coupling structures, involving theta and gamma (Figure 5B). One may separate two frequency regions corresponding to the coupling between theta and its harmonics, and theta and gamma. The picture of the coupling between theta and its harmonics agrees with the spectral evolution (Figure 3) and previously reported results (Sheremet et al., 2016). The bicoherence exhibits significant levels of phase coupling, reaching as high as (5f_θ, f_θ, 6f_θ), (3f_θ, 2f_θ, 5f_θ), and (3f_θ, 3f_θ, 6f_θ). The relationship between harmonics and theta is quite diverse, with some harmonics contributing to LFP skewness, and others only to LFP asymmetry. For example (Figure 5B, skewness and asymmetry), the coupling (f_θ, f_θ, 2f_θ) generates both negative skewness and positive asymmetry, while (2f_θ, 2f_θ, 4f_θ) results in negative asymmetry only. Theta-gamma coupling is prominent at high speed, engaging theta harmonics, and contributing strongly to negative skewness (Figure 5B).

3.4. Summary of Observations

Our observations show that sorting LFP epochs by speed provides an efficient classification device that produces remarkably well-ordered spectra and bispectra. As summarized in Figure 6, with increased speed the total LFP power increases, as well as the power in the theta and gamma bands. The power in the theta band grows overall by a factor of 4 at a relatively steady rate as a function of speed, while gamma grows by about a factor of 2 and seems to plateau. Phase coupling (as measured by the bicoherence integrated over the theta/harmonics and theta/gamma frequency domains) shows a steady growth as a function of speed, accelerating at high speeds. Spectral slopes (α) decrease, consistent with an accumulation of energy in the gamma range; and phase coupling involving theta and gamma increases.

The process of evolution with speed is reversible. Recall that we refer to evolution as change with increased speed, therefore the “reverse” process is the transformation corresponding to decreasing speed. The reverse process appears to converge as v → 0 (v is speed) toward a limiting, “background” state. Ignoring for now other forms of activity not taken into account here, the background state represents “inactive” behavior, whereas the active state would be associated with significant speed levels. This hypothetical background has some remarkable properties. Although its spectrum decays faster than any active-state spectrum (higher slope α), it still contains significant power (in Figure 2, about $\frac{1}{4}$ of the most active state). In addition, the LFP is nearly Gaussian, i.e., exhibits overall no phase coupling (ignoring the weak theta signal).

The evolution with speed of hippocampal power distribution over scales is strongly suggestive of cross-scale energy exchanges, statistically directed from low to high frequencies (from large to small scales). Cross-scale energy transfers can only result from non-linear interaction between scales; however, there is no obvious reason for a statistically preferential direction of transfer, unless the physical system is perturbed in a way that forces a certain internal exchange. The apparent direction of the cross-scale energy flow is consistent with the system receiving energy at the large scales (low frequency end) and losing energy at the small scales (high frequency end). The collapse, in the range f < 12 Hz, of the background spectrum as theta itself progressively concentrates all spectral power suggests that the background spectrum is not associated with the energy input and represents a completely different process, replaced at large scale by theta during task behavior. This shift and the considerable increase in theta power are consistent with theta playing an important part in the energy input process (please see next section for further discussion).

In our observations, the evolution of the scale distribution of LFP power and non-linearity exhibits some obvious, yet unexplained features. Why do spectral slopes change? Why does phase coupling develop? Why does theta develop harmonics, but gamma does not, coupling instead to theta? How is the collective action energy balance related to the background spectrum?

3.5. The Weak Turbulence Framework

Figure 7 compares side by side the standard conceptual model of turbulence with a schematic of observed hippocampal LFP evolution with energy input (speed).

The main focus of the turbulence theory is the internal energy⁷ balance in a non-linear, multi-scale physical system whose scales interact and exchange energy (Figure 7, left). If the energy source and sink of the system are well-separated in scale, for example, the source is located at large scales, and the sink at small scales, there exists an intermediate domain, called inertial range, where the only energy process is the cross-scale energy flow. Assume that the input rate at the source is constant, and the system has initially zero energy. If the scale interaction is local (i.e., the strongest interaction is between similar scales) input energy will accumulate initially in a band of scales adjacent to the source. As the energy in that band increases, local non-linear exchanges intensify (non-linearity increases with energy), and the energy introduced into the system flows downscale. Eventually, the energy flow (non-linear cascade) reaches the dissipation (sink) scales, where it is removed from the system. The stationary state will occur when the dissipation rate matches the input rate, and will be characterized by a constant flux of energy across the inertial range, from source to sink⁸. The cross-scale energy flux is turbulent cascade (Richardson, 1922). One of the most celebrated results of hydrodynamic turbulence is the Kolmogorov's (1941) argument that stationary spectrum follows the power law k^−5/3, where scales are represented by the wavenumber k.

But for the presence of a background spectrum, the similarity between the standard turbulence model and LFP spectral evolution are striking. The analogy suggests that theta is the source of energy for collective action, and microscopic processes are as the main energy sink, while the mesoscale acts largely as an inertial window that allows for a cross-scale energy flow from source to sink. The development of the gamma peak, similar to the bottleneck effect in hydrodynamics, may be interpreted as the existence of a transitional scale right above microscopic, that has a limited energy-flux ability, thus causing an accumulation of energy at intermediate scales. The observation that spectral evolution with speed is reversible (in other words, that spectra at each speed represent an equilibrium state) is also remarkable, because it is consistent with Kolmogorov stationary spectra of turbulence.

3.5.1. Mesoscopic Turbulence on an Active Network

Therefore, turbulence framework could be useful for collective action dynamics, provided that the mesoscale may be identified with the inertial range; or, equivalently, if non-linearity dominates dissipation at mesoscale. However, while observations of propagating theta waves suggest that collective action is weakly dissipative, this seems to contradict the well-known strongly-dissipative character of microscopic neuronal dynamics⁹. Furthermore, the presence and role of the background spectrum deserves some discussion.

From the perspective of collective action, the hippocampus behaves as an “active” network, i.e., a physical system of active elements that contain a certain amount of energy which they release in explosive bursts when activated by a threshold type of trigger. The statistical state of the network may be described by a probability distribution of internal energy around a mean level, sustained by energy (pulses) channeled through network connections. Even if the mean internal energy level is below the threshold, a small number of elements will have energy exceeding the threshold, and will fire. A fraction of the burst energy is recaptured and redistributed through the network connections to maintain a mean level of internal energy. After bursts, network elements go through a “recharge” (refractory) period. Note that, in the absence of external input, bursts are the only energy source for maintaining the internal energy level. If no bursts occur, the system collapses. In an equilibrium state, the probability distribution of the internal energy should have a standard deviation large enough to maintain the mean.

If the state described above is perturbed by adding energy over a mesoscopic area, the internal energy distribution shifts toward the threshold, increasing the number of bursts. The perturbation decays or grows in time depending on whether the ratio of energy recaptured from bursts to the energy of the initial perturbation smaller or larger than 1. If this fraction is ≈ 1, the perturbation is self sustained.

We conjecture that self-sustained collective action is the behaviorally meaningful hippocampal activity. Thus, the relevant type of collective action is only weakly dissipative. Remarkably, it also likely requires a non-zero background activity (background spectrum, section 3.1). Indeed, for a self sustained perturbation, both the energy recaptured from bursts and the energy of the perturbation depend on the mean energy level of the system. A likely mechanism to maintain an equilibrium level of internal energy, and its required standard deviation, is to randomly trigger multi-scale burst patterns. This idea circles back to the concept of self-organized criticality (SOC) (SOC; Bak et al., 1988; Beggs and Plenz, 2003; Beggs and Timma, 2012; Pruessner, 2012) and metastability (Tognoli and Kelso, 2014). The featureless, power-law, Gaussian background spectrum we observe (Figure 3) is consistent with the “edge of chaos”, self-organized critical background state for optimal transmission of collective action. In this context, SOC becomes an important element of brain dynamics.

The relation between background SOC state (rest) and collective-action turbulence (intense task behavior) might also explain the behavior of the LFP power spectra in the low frequency range f < 12 Hz (see discussion at the end of section 3.1). If the excess energy for sustaining the collective action is provided (as observations suggest) by increasing theta power, much of the large scale activity is taken over by organized theta action, thus restricting the large scale extent of SOC spectra.

3.5.2. The Weak Turbulence Model

The turbulence framework may now be formalized. The turbulence theory investigates the internal balance of a multi-scale system. It studies the interaction efficiency as a function of scale (e.g., resonance conditions); the development of cross-scale coupling correlations between scales; characteristics of the long term evolution of the system; the existence of equilibrium spectra associated with cascades of conserved quantities (e.g., power law spectra that maintain a non-zero cross-scale energy flux); non-stationary evolution patterns (e.g., bottleneck formation); and other characteristics of the non-linear, cross-scale exchange mechanisms related to conservation laws.

These are precisely the features that characterized the evolution with speed of our LFP observations (e.g., Figure 4).

Here we give only the a brief elementary discussion of the basic equations that govern weak turbulence. The principles of the weak turbulence model, some derivation algebra, and some results are further discussed in the Supplementary Material. We caution, however, that presentation is an oversimplified, retold version, and that the full theory is much richer and complex. We encourage the interested reader to consult the original sources, written by the fathers of the theory: the comprehensive monograph (Zakharov et al., 1992), excellent review papers by Zakharov (1999), Newell et al. (2001), and Newell and Rumpf (2010), and the account by Nazarenko (2011), that includes more recent results.

3.5.3. Dynamical Equation

The focus of the turbulence theory is to describe the internal energy balance of a multi-scale system, i.e., the evolution of the distribution of energy over scales (the spectrum), the cross-scale energy flux (the energy cascade), associated non-linear cross-scale coupling, and other quantities. The general way to achieve this is to define the state of the system using a “state variable”, for example internal energy, and write its evolution using conservation laws. For example, in thermodynamics the rate of change of internal energy is equal to the heat and mechanical work exchanged by the system with the environment.

For a spatially-distributed system such as collective action, the state variable becomes a function of space and time (a field), and the evolution equation becomes a partial differential equation. The exact form of these equations for collective action will be presented elsewhere. Here, it suffices to note that the internal energy field for hippocampal collective action may be defined proportional to the mean neuron potential per unit area, and the conservation of the internal energy involves a balance of electric pulses sent down the network connection and the energy recaptured by the network from bursts. This approach leads to equations similar to the, Wilson and Cowan (1972,1973); Cowan et al. (2016) model, although its form and derivation methodology is different and somewhat less consistent.

To avoid having to reference the specifics of the system, for illustration purposes only we will tentatively identify the internal energy with the Hamiltonian, which will conveniently allow us to introduce the equations of weak turbulence in a general form, without having to reference the specifics of the system. The Fourier representation of the collective-action field provides a scale decomposition of the internal energy field, so in principle, the Fourier transform of the conservation law (Hamilton's equations) yields an equation for the Fourier amplitude a(k, t) of the collective action component with wave number k. The general form of this equation is:

where Fourier collective-action mode (or simply, mode) is identified by its wavenumber k and frequency f(k), (we assume that the dispersion relation has a single root). In Equation (8), the notations are: ω(k) = 2 π f(k) is the radian frequency; $\stackrel{•}{a}=\frac{da}{dt}$; V_{k; 12} = V(k; k₁, k₂) is the interaction coefficient; and δ is the Dirac delta function. We use the shorthand notation σ_1,2 = σ (k_1,2), σ_k = σ (k), and ${\delta }_{k;12}^{\sigma }=\delta ({\sigma }_k-{\sigma }_1-{\sigma }_2)$, where σ is some quantity depending on k; and also dk₁₂ = dk₁dk₂. For simplicity, we limit the discussion quadratic non-linearity (right-hand side terms of the form a²) terms, but a full description might require including cubic (a³), and possibly higher-order non-linearity (e.g., the non-linear Schrodinger equation, Newell, 1985; or the Ginzburg-Landau equation Ermentrout, 1981; Cross and Hohenberg, 1993; Passot and Newell, 1994; Ermentrout et al., 1997).

Two equations for the modulus b ≥ 0 and phase θ of a(k, t) are obtain substituting a(k, t) = b(k, t)e^{iθ(k, t)} into Equation 8

where we used the notation ${\Delta }_{k;23}^{\theta }={\theta }_k-{\theta }_1-{\theta }_2$.

Meaning of the dynamical Equation 8.

(i) Temporal (frequencies) and spatial scales (wavenumbers) are related through the dispersion relation (see the Supplementary Material), therefore only one of the parameters k and f is independent. The choice of the independent parameter is arbitrary, because the dispersion relation is invertible. The f(k) representation resolves the spatial structures and yields a time-evolution equation. The k(f) representation resolves the time structure and is appropriate for time-series analysis (e.g., LFP recordings). Therefore, equations are written here in f(k), but observational data is discussed in the k(f) representation. Below, the concepts of “frequency”, “wavenumber”, “mode” and “scale” will be treated as equivalent and interchangeable.

(ii) Equation 8 is called dynamical equation. If the Hamiltonian is identified with energy, the quantity |a|² has the physical dimensions of action (energy × time). The form of Equation (8) in universal in the sense that the details of the physics of the system are contained the coefficients ω_k and V_{k; 12} only.

(ii) Equation (8) describes the non-linear interaction of mode k with the pair of modes (k₁, k₂). A triplet of interacting modes (k, k₁, k₂) is called a “triad”. The strength of the interaction depends on the interaction coefficient V_{k; 12}. The factor ${\delta }_{k;12}^k$, resulting from the orthogonality of the Fourier representation, is a selection criterion: interacting modes satisfy the equation

It is useful to think of Equation (8) in a discretized form, e.g., replacing the integrals by sums. A schematic representation is shown in Figure 8.

(iv) Equations (9) show that non-linear interaction result in both amplitude and phase evolution. The effectiveness of non-linearity depends on modal amplitudes, on the interaction coefficient, and (importantly) on the phase mismatch ${\Delta }_{k;12}^{\theta }$. If ${\Delta }_{k;12}^{\theta }$ is large, the non-linear term oscillates fast and the effect small; if ${\Delta }_{k;12}^{\theta }$ is small, the non-linear term preserves sign over longer periods of time and the effect is significant. If ${\Delta }_{k;12}^{\theta }$=0, the contribution is maximal.

(v) The non-linear contributions to modal phase evolution (Equation 9b) have same properties as those in the amplitude Equation 9a, with the important difference of b_k appearing at the denominator. If b_k → 0 (mode k has small amplitude) the non-linear phase component becomes arbitrarily large and dominates the total phase. Therefore, the phase of low-amplitude waves is “dictated” by non-linear forcing.

(vi) The dynamical equation 8 is deterministic: may be integrated exactly to obtain the state of the system at time t if the initial value of amplitudes a(k, t₀) are known (t > t₀). To understand the general energy balance (e.g., spectral evolution), however, one is not interested in particular realizations, but in the predictions the equation provides for averaged quantities such as the spectrum or bispectrum. In practice, the averaging operator 〈·〉 used in Equations (1–2) may be replaced by averaging over initial phases (see the discussion of random phase average in Nazarenko 2011). Therefore, estimates of spectra and bispectra may be obtained by integrating Equation (8) many times with different sets of initial phases and, for example for the spectrum, averaging $|a_k{|}^2$ over all integrations.

3.5.4. The Kinetic Equation

The last remark above suggests averaging the equation itself, rather than averaging the solutions. Briefly (see Supplementary Material), averaging introduces a hierarchy of averages of amplitude products (correlators), such as $\langle a_1^{*}{a}_2\rangle$, $\langle a_1^{*}{a}_2{a}_3\rangle$, $\langle a_1^{*}{a}_2^{*}{a}_3{a}_4\rangle$, and so on, where the angular brackets denote the ensemble average. Assuming spatial homogeneity implies that

The quantity n(k) represents the action density, or “occupancy number”, or “number of particles” (by analogy with quantum mechanics). In statistics n and ${B}$ are called “spectrum” and “bispectrum”, respectively. The quasi-Gaussian assumption results in a system of two coupled equations for the spectrum and bispectrum

Assuming additional long-time regularity conditions reduces the system to single equation called the kinetic equation

where ${\delta }_{k;12}^{\omega }=\delta ({\Delta }_{k;12}^{\omega })=\delta ({\omega }_k-{\omega }_1-{\omega }_2)$.

Meaning of kinetic Equation 14.

(i) Equation (13a) highlights the dynamical significance of the bispectrum as the non-linear forcing in the evolution of the spectrum. If the bispectrum cancels, non-linear interactions cancel and the system is effectively, on average, linear.

(ii) The kinetic Equation (14) or the more general system of Equations 13, represent a stochastic, ensemble-averaged description of the system 8. Kinetic equations of the type of Equation 14 were introduced in statistical mechanics by Boltzmann (e.g., Boltzmann, 1872, 2003; Alexeev, 2004), and are powerful tools in the study of multiple-scale system.

(iii) Equation (14) states that in the long-time limit the only interactions that are effective are due to triads that satisfy the resonance conditions imposed by the factors ${\delta }_{k;12}^{\omega }{\delta }_{k;12}^k$, i.e., satisfying the conditions,

equivalent to the “maximal” effectiveness of non-linear interaction (see discussion of Equations 9). Whether or not Equation (8) has resonant triads depends on the linear properties of the physical system. The resonance conditions play an important role in the stochastic theory (e.g., Zakharov et al., 1992; Nazarenko, 2011; Anenkov and Shrira, 2018).

3.5.5. Stationary Spectra

Stationary spectra ($\stackrel{•}{n}(k)=0$) are of importance for systems whose evolution is a quasi-equilibrium process. Equation (14) has two classes of stationary solutions (e.g., Nazarenko, 2011).

The Rayleigh-Jeans (RJ) class of spectra comprises the stationary solutions k⁻¹, ω⁻¹, and ${(c_{\omega }\omega +c_{k}k)}^{-1}$, that obviously cancel the integrand in Equation (14), and correspond to zero exchange across scales. Therefore RJ spectra correspond to thermodynamic equilibrium and equipartition of momentum (nk) and energy (nω). Realistic, non-isolated systems do not typically reach thermodynamic equilibrium, therefore RJ spectra are not important for the weak turbulence framework, which includes sources and sinks as essential elements.

The Kolmogorov-Zakharov (KZ) spectra are different class of stationary spectra, that correspond to a constant spectral flux ∂_kF_q = 0, F_q(k) ≠ 0 across the inertial range. They were derived for Equation (14) by Zakharov and Filonenko (1967a,b). They are important for systems in which the dissipation sink can absorb arbitrary rates of energy. Remarkably, they are realized as a non-trivial power law spectrum n^KZ ∝ k^ν, with ν < 0, ν ≠ − 1. The slope ν of the spectrum is a value that reflects the dimensionality of the system, as well as its linear non-linear properties (homogeneity degrees of the interaction coefficient and dispersion relation).

3.6. A Demonstration: Dynamics and Kinetics of a Single Triad

Some features of the high-speed bispectra are consistent with stationary solutions of the dynamical equations 8. We illustrate this using a simplified, universal toy model derived from the framework dynamical equation by considering a single triad of modes k = κ₁, k₁ = κ₂, k₂ = κ₃, with κ₁ + κ₂ = κ₃ (see Equation 10). If interactions with all other modes are ignored (equivalent to a single blue dot in Figure 8), Equations (8) to reduce to three equations, called the “three-wave system” (see Supplementary Material; Craick 1985; Rabinovich and Trubetskov 1989). Written in amplitude/phase these are

where b_j and θ_j are amplitudes of modes κ_j, j = 1, 2, 3, and ${\Delta }_{3;12}^{\theta }={\theta }_3-{\theta }_1-{\theta }_2$. Averaged over realizations, the quantity ${\Delta }_{3;12}^{\theta }={\theta }_3-{\theta }_1-{\theta }_2$ may be identified with the biphase (see section 3.3). If the triad is resonant, the kinetic version of the three wave equation is (e.g., Rabinovich and Trubetskov, 1989)

A large body of literature is available that investigates the relevance and dynamics of single-triad interactions in many physical situations, including plasma physics (Coppi et al., 1969; Weiland and Wilhelmsson, 1977; Craick, 1985), non-linear optics (Ablowitz and Segur, 1981; Boyd, 2003), internal oceanic waves (Phillips, 1977; Craick, 1985), and other fields.

One may readily check that a stationary solution of Equations (17) is given by the conditions: ${\Delta }_{3;12}^{\Theta }=0$ or ${\Delta }_{3;12}^{\Theta }=\pi$; ${\Delta }_{3;12}^{\omega }=0$ (the triad is resonant, Equations 16), and that

The last equation has, for example, the trivial solution $b_j^2={(c_{\omega }\omega +c_{k}k)}^{-1}$(RJ spectrum, thermodynamic equilibrium). In other words, stationary state solutions of the framework equations exhibit naturally 8 a biphase of 0 or π.

Applied to the triad formed by theta and its first harmonic (θ, θ, 2θ), i.e., κ₁ = k_θ, κ₂ = k_θ, and κ₃ = k_2θ (red dot in Figure 8) this result is consistent with observations that show theta and in first harmonic are in phase (see Figure 5B, upper-right panel). The effect of this type of coupling is to sharpen the crests and flatten the troughs of the time series, generating positive skewness.

Applied to a theta-gamma triad (γ, θ, γ + θ), κ₁ = k_θ, κ₂ = k_γ, and κ₃ = k_θ + k_γ (yellow dot in Figure 8), this result is consistent with the biphase value of π in Figure 5B, upper-right panel. It is easy to check that effect of this type of coupling is that gamma envelope is maximal in theta troughs. Indeed, let φ_j = a_j cos [κ_jx − ω_jt + ψ_j], where κ_j and ω_j satisfy the resonance conditions 16. Elementary trigonometric manipulations yield the gamma envelope $\propto sin(\frac{{\omega }_1}{2}t).$ In other words, the gamma envelope is in quadrature with theta and its period is twice that of theta. The kinetic evolution of a triad is illustrated in Figure 9. The interactions described by the three-wave model result in time-reversible cyclic transfers of energy (the total energy is conserved if the system is at resonance). In contrast, the long-time, average behavior (Equation 18) shows a slow irreversible trend toward the stationary solution. If a full set of triads is taken into account, because all triads interact, a weak energy flow (turbulent cascade) develops that has the effect of driving the system of modes toward stationary distribution of energy over scales (RJ spectrum). If a sink is introduced in the small scales (e.g., high frequencies), the flow will naturally be directed toward the leaking mode, where it is taken out of the system. In this case a stationary state can be maintained only if energy is pumped into the system at the leakage rate, and a stationary state is realized if the energy is injected by the source at the rate it is lost to the sink. The cross-scale flow of energy is in this case constant at all scales. This type of stationarity corresponds to the KZ spectra.

4. Discussion

While a wealth of knowledge has accumulated in recent years about brain activity at microscopic and macroscopic scales, the need for a consistent theory of the dynamics of intermediate-scale (mesoscopic) processes has received comparatively little attention, possibly because, in some parts of the nervous system, they are the expression of activity of complex microscopic structures (e.g., for example central pattern generators can maintain specific, yet re-configurable rhythms; e.g., Rabinovich et al., 2012; Gutierrez et al., 2013; Marder et al., 2016). Although some microcircuits, with definite structure, may be able to impose an oscillation on the mesoscale, the cortex, with an isotropic and homogeneous structure, is rhythmically organized in a different manner, with a different source and might play a different role.

A number of previous studies (Lashley, 1942; Hebb, 1958; Freeman, 2000a, 2007; Freeman and Vitiello, 2010) hypothesize that mesoscopic processes in the cortex represent the essential cognition step of abstraction and generalization, and therefore provide an essential mechanism for integration of brain activity at all scales. This observation is particularly intriguing, because the isotropic and homogeneous structure of the cortex at the anatomical mesoscale (e.g., any highly recurrent region in which activity is projected back into the same region or dense inter-connectivity mediated by interneurons; Lorente de No, 1938; Freund and Buzsáki G., 1996; Buzsáki et al., 2004; Mante et al., 2013) suggests that the material support of activity in the temporal mesoscale (e.g., gamma frequency band) is collective neural activity, Freemans's (2000a) “mass action”. These ideas suggest that the physics of mesoscopic collective action in the cortex is intimately related to cognition; that the physics of collective action is, in fact, the physics of cognition.

The focus of this study is the dynamics of mesoscopic collective action and its role in the general energy balance in the brain. Despite the potentially paramount importance of the topic, studies dedicated to its physics are few (but of outstanding quality; e.g., Wilson and Cowan, 1972, 1973; Wright and Liley, 1995; Troy, 2008; Cowan et al., 2016). Here, we attempt to lay the foundation of a systematic theory of collective action. A few recent studies show that collective action in the hippocampus takes the form of propagating waves (Lubenov and Siapas, 2009; Patel et al., 2012, 2013; Muller et al., 2018; Zhang et al., 2018), but in general, information about its spatio-temporal organization is scarce.

The discovery of the strong correlation between rat speed during active exploration and hippocampal activity provides a parametrization of the evolution of hippocampal activity with behavior. Our observations of scale distribution of LFP power (spectra) and the leading order estimators of cross-phase coupling of LFP oscillations (bispectra) in the lacunosum moleculare layer show a strong and remarkably ordered evolution (representative of both CA1 and dentate gyrus layers). The LFP power and phase coupling in the theta and gamma frequency bands increases consistently with speed. The lowest levels of LFP power are associated with a featureless power law distribution. The high-power spectra exhibit distinctive spectral peaks at theta frequency and harmonics, and significantly increased levels of gamma activity. In the transition, at intermediate states, the spectrum tilts to a smaller-slope shape that extends progressively into the gamma range generating the appearance of a spectral front. The existence of a non-zero energy, lowest-level spectrum in the absence of coherent collective action is strongly suggestive of a self-organized critical background state (Buzsaki, 2006). As collective action energy increases with exploration, the decreasing spectral slope, the appearance of a spectral front, and the development of a broad gamma peak, are strongly suggestive of a transfer of energy from the low frequencies (large scales) toward high frequencies (microscopic scales; Buzsáki and Draguhn, 2004). The similarity of collective action with the general turbulence theory is striking. This motivated us to propose weak turbulence as an framework for collective action dynamics.

In summary, we propose that the mesoscopic hippocampus may be described as an isotropic and homogeneous, active network containing a macroscopic number of randomly and densely connected neural units (neurons). Collective action represents a perturbation of a background state of the active network, that may be represented as a self-organized critical state. Because collective action is macroscopic with respect to neural units, we postulate therefore that it is fundamentally stochastic (the precise state of a single unit does not matter), weakly non-linear, and weakly dissipative. These form a minimal set of features for the development of a turbulence theory of collective action. We summarize the principles of the turbulence framework and demonstrate its applicability to observations by showing that the reduced (universal) three-wave interaction toy model provides results consistent with observations.

Our proposed description provides a unified view of the physics of active networks that reconciles the theories of self organized criticality and turbulence in the hippocampus and perhaps other regions of the cortex. The central idea of the turbulence theory is the energy cascade through the inertial range of scales. As a theory of the internal balance of a multi-scale system, a full turbulence formulation of hippocampal collective action should provide answers to questions raised by the evolution the scale distribution of LFP power. Collective action physics appears to be consistent with the energy cascade concept (Richardson, 1922; Kolmogorov, 1941; Zakharov et al., 1992; Newell et al., 2001; Nazarenko, 2011). In particular, the spectral tilting (slope change) observed in the evolution from the background state to high-speed states suggests a transition from self-organized criticality to some type of stationary turbulent state of the Kolmogorov-Zakharov kind (Zakharov and Filonenko, 1967a,b; Zakharov et al., 1992; Newell et al., 2001). Do spectral slopes change because the system transitions from a SOC background state to a turbulent, stationary state supporting non-zero cross-scale energy fluxes? Are non-linear resonances band-limited, for example, forcing theta to develop self interaction, while allowing theta-gamma resonances? Is the growth of gamma power related to a spectral bottleneck (a break in cross-scale non-linear energy exchange at high frequencies)? Further research is needed to understand the structure of this evolution; a systematic exploration of hippocampal dynamics starting from first-principles governing equations is ongoing. The scope of this study is limited to investigating the effect of running speed as a proxy for energy into the system. Future experimental applications should explore transitions between sleep and wakefulness, the effects of learning and memory recall as well as the system is compromised in aging and disease. As cognition is the consequence of activity moving through the circuits of the brain, any experimental approach which alters the forcing or how the local network capture the activity (Bottleneck, Figure 7) will have the ability to test the predictions of the turbulence theory.

Remarkably, the isotropic and homogeneous mesoscale structure appears to be a general theme that is used across species and brain regions (Lorente de No, 1938; Parent and Hazrati, 1995; Marder and Bucher, 2001; Garamszegi and Eens, 2004; Apps and Garwicz, 2005; Mante et al., 2013), which suggests a “universal computational principle” with a comprehensive reconfiguration potential, especially under a priori unknown conditions (Sussillo and Abbott, 2009). The nature of this computation process is not well-understood. As often argued (e.g., Freeman, 2000b; Edelman and Gally, 2001; Frisch, 2014), cognition processes cannot resemble a human-engineered system, built based on principles such as maximum simplicity, well-defined internal interactions, explicit assignment of function, no irrelevancy, and no adventitious compensation for error. Rather, they are expected to resemble biological systems: no design, no a priori function, and for which irrelevance has no meaning (Edelman and Gally, 2001). As Frisch (2014) states it, “biological systems have an intrinsic ability to maintain functions in the course of structural changes”, such that “specific functions can obviously be constituted on the basis of structurally different elements, a biological property that is referred to under the term degeneracy (Edelman and Gally, 2001)”.

This begs the question, if collective action is fundamental for cognition, what is its role? We conclude this study by suggesting a possible answer. An intriguing paradigm of the computational function of mesoscale turbulence is offered by Liquid State Machines (LSM) models (Jager, 2002; Maass et al., 2002). LSMs are online neural network models that process a time-windowed signal in real time. Their basic function is to perform a non-linear transformation of the input, e.g., expand it into a wave field, and hold this information for a short duration, while output neurons extract local information from the field. Learning is achieved at the readout stage. Fading memory and input separability imply that LSMs are universal function approximators, and can serve as effective online classification pre-processor for the readout neurons. If the brain is a prediction machine, scrambling to assign meaning to streams of data¹⁰ in real time, processing of fragmentary information (short-time windows) is crucial. A LSM may perform fast, online, short-term memory pre-processing; learning is performed by long-term memory readouts, that can record optimal responses. It is conceivable that the hippocampus uses collective action in a way similar to a LSM, perhaps as a fast online classification machine, or as a dynamical system simulator. It seems plausible to imagine the cortex as a network of LSMs. In the least, the concept seems to agree with most natural systems.

Author Contributions

AS and AM collaborated closely in developing the turbulence framework for brain activity. AS developed part of the numerical methods and the Matlab^® tools used fort the data analysis. YZ and YQ developed part of the numerical methods, and performed the numerical analysis. JK managed data collection. AM provided the impetus for the research; designed, organized, and supervised the data collection, and provided the guidance and advice regarding all things neuroscience. AS wrote the manuscript with input from all authors.

Funding

This work was supported by the McKnight Brain Research Foundation, and NIH grants-Grant Sponsor: National Institute on Aging; Grant number: AG055544 and Grant Sponsor: National Institute on Mental Health; Grant Number: MH109548.

Conflict of Interest Statement

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.

Acknowledgments

Special thanks to S. Burke for multiple insightful discussions, advice on the manuscript redaction and for patiently reviewing manuscript drafts; and to S. D. Lovett for technical support. AS is grateful to Prof. Victor Shrira for helpful discussions and suggestions. We are grateful to the two reviewers for the thoughtful advice on improving the quality of the paper.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnsys.2018.00062/full#supplementary-material

Footnotes

References