## Abstract

An inhomogeneous anisotropic physical model of the brain cortex is presented that predicts the emergence of nonevanescent (weakly damped) wave-like modes propagating in the thin cortex layers transverse to both the mean neural fiber direction and the cortex spatial gradient. Although the amplitude of these modes stays below the typically observed axon spiking potential, the lifetime of these modes may significantly exceed the spiking potential inverse decay constant. Full-brain numerical simulations based on parameters extracted from diffusion and structural MRI confirm the existence and extended duration of these wave modes. Contrary to the commonly agreed paradigm that the neural fibers determine the pathways for signal propagation in the brain, the signal propagation because of the cortex wave modes in the highly folded areas will exhibit no apparent correlation with the fiber directions. Nonlinear coupling of those linear weakly evanescent wave modes then provides a universal mechanism for the emergence of synchronized brain wave field activity. The resonant and nonresonant terms of nonlinear coupling between multiple modes produce both synchronous spiking-like high-frequency wave activity as well as low-frequency wave rhythms. Numerical simulation of forced multiple-mode dynamics shows that, as forcing increases, there is a transition from damped to oscillatory regime that can then transition quickly to a nonoscillatory state when a critical excitation threshold is reached. The resonant nonlinear coupling results in the emergence of low-frequency rhythms with frequencies that are several orders of magnitude below the linear frequencies of modes taking part in the coupling. The localization and persistence of these weakly evanescent cortical wave modes have significant implications in particular for neuroimaging methods that detect electromagnetic physiological activity, such as EEG and magnetoencephalography, and for the understanding of brain activity in general, including mechanisms of memory.

## INTRODUCTION

Most approaches to characterizing brain dynamical behavior are based on the assumption that signal propagation along well-known anatomically defined pathways, such as major neural fiber bundles, tracts, or groups of axons (down to a single-axon connectivity), should be sufficient to deduce the dynamical characteristics of brain activity at different spatio-temporal scales. Experimentally, data on which this assumption is employed range from high-temporal-resolution neural oscillations detected at low spatial resolution by EEG/magnetoencephalography (MEG; Monai, Inoue, Miyakawa, & Aonishi, 2012; Ingber & Nunez, 2011) to high-spatial-resolution fMRI resting-state mode oscillations detected at low temporal resolution (Li et al., 2015; Zalesky, Fornito, Cocchi, Gollo, & Breakspear, 2014). As a consequence, a great deal of research activity is directed at the construction of connectivity maps between different brain regions (e.g., the Human Connectome Project; Van Essen et al., 2013) and using those maps to study dynamical network properties with the help of different models of signal communication through this network along structurally aligned pathways (Hallez et al., 2005; Haueisen et al., 2002; Tuch, Wedeen, Dale, George, & Belliveau, 1999), for example, by allowing input from different scales or introducing axon propagation delays (Nunez & Srinivasan, 2014).

However, recent detection (Muller et al., 2016) of cortical wave activity spatio-temporally organized into circular wave-like patterns on the cortical surface, spanning the area not directly related to any of the structurally aligned pathways but nevertheless persistent over hours of sleep with millisecond temporal precision, presents a formidable challenge for network theories to explain such a remarkable synchronization across a multitude of different local networks. In addition, many studies show evidence that electrostatic field activity in animal or human hippocampus (as well as cortex) are traveling waves (Muller, Chavane, Reynolds, & Sejnowski, 2018; Zhang, Watrous, Patel, & Jacobs, 2018; Lubenov & Siapas, 2009) that can affect neuronal activity by modulating the firing rates (Czurkó, Huxter, Li, Hangya, & Muller, 2011; Stewart, Quirk, Barry, & Fox, 1992; Fox, Wolfson, & Ranck, 1986) and may possibly play an important functional role in diverse brain structures (Weiss & Faber, 2010). Perhaps more importantly, it has been experimentally shown (Chiang, Shivacharan, Wei, Gonzalez-Reyes, & Durand, 2019; Shivacharan, Chiang, Zhang, Gonzalez-Reyes, & Durand, 2019; Qiu, Shivacharan, Zhang, & Durand, 2015; Zhang et al., 2014) that periodic activity can self-propagate by endogenous electric fields through a physical cut in vitro that destroys all mechanisms of neuron-to-neuron communication.

In this paper, we investigate a more general physical wave mechanism that allows cortical surface wave propagation in the cross-fiber directions because of the interplay between tissue inhomogeneity and anisotropy in the thin surface cortex layer. This new mechanism has been overlooked by previous models of brain wave characterization and thus is absent from current network pathway reconstruction and analysis approaches.

The motivation behind this paper is twofold: (1) to provide an analytical framework for deriving circular wave-like patterns given the complex anisotropy and conductivity in the brain and (2) to show that the model can explain many of the observed properties of the real data using different models of tissue anisotropy and conductivity varying from the simplest fixed-position independent anisotropy to the complex anisotropy created by an ensemble of multiple fibers. The framework and the model can be used to inform neuroscientists about what is the “optimal” way to investigate wave patterns in the brain and can be used in diffusion simulations (Frank, Zahneisen, & Galinsky, 2020; Berry, Regner, Galinsky, Ward, & Frank, 2018; Baxter & Frank, 2013) to optimize diffusion MRI (dMRI) acquisitions to get a better sense of the underlying anisotropy and conductivity and thus get a more precise characterization of brain waves.

The main claim of this paper is that there is a simple and elegant physical mechanism behind the existence of these cross-fiber waves that can explain the emergence and persistence of wave loops and wave propagation along the highly folded cortical regions with a relatively slow damping. The lifetime of these wave-like cortex activity events can significantly exceed the decay time of the typical axon action potential spikes and thus can provide “memory-like” response in the cortical areas generated as a result of “along-the-axon” spiking activations. This new mechanism may provide an alternative approach for the integration of microscopic brain properties and for the development of a “physical” model for memory.

Our nonlinear theory is a significant departure from existing approaches, which are based on local nonlinear model fitting in single neuron models, typically by means of empirical selection of a set of model parameters, like firing thresholds, fixed frequencies of oscillations of slow and fast subsystems of the model, and so forth. Instead, our model uses actual tissue and geometry properties and thus presents the possibility of making quantitative hypotheses testable with actual measurements that elucidate the dependence of the wave parameters (e.g., frequency) on tissue parameters and geometry.

Because this framework for nonlinear brain wave dynamics gives rise to brain waves not determined by neuronal geometry and connectivity, but instead by tissue inhomogeneity and anisotropy, waves can traverse regions between neurons and thus the framework provides a mechanism for ephaptic coupling and related complex phenomena such as development or disruption of spiking synchronization. There is currently no physical theory for such effects despite their importance in a range of problems from normal brain functioning to disorders such as epilepsy. Our model also facilitates the practical study of these effects by inclusion of a single (or multiple) neuron-scale input to EEG/MEG/MRI estimation and analysis. The brain wave dynamics-based mechanism may also play a role as an effective biophysically motivated scale-dependent learning-like mechanism showing how the brain implements back-propagation/feed-forward approaches of artificial neural network models. The spatial and temporal properties of this brain wave synchronization mechanism may be used to estimate the importance of different spatial and temporal scales and to select the most appropriate interactions between scales for joint estimation and analysis and ultimately may allow the expansion of estimation, analysis, and stimulation prediction to the whole hierarchy of scales, from membrane (subneuron) to the entire brain.

## METHODS

### Linear Theory

We will start with the most general form of description of brain electromagnetic activity using Maxwell equations in a medium:
$∇·D=ρ,∇×H=J+∂D∂t⇒∂ρ∂t+∇·J=0$
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 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
$∂∂t∇2φ=−∇·Σ·∇φ+𝓕,$
(1)
where we have included a possible external source (or forcing) term 𝓕. 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., 𝓕 = 0) can be written in tensor notation as
$∂t∂i2φ+Σij∂i∂jφ+∂iΣij∂jφ=0,$
(2)
where repeating indices denote summation. Simple linear wave analysis, that is, substitution of φ ∼ exp(−i(k · r − 𝒪t)), where k is the wave number, r is the coordinate, 𝒪 is the complex frequency, and t is the time, gives the following complex dispersion relation:
$𝒟𝒪k=i𝒪ki2+Σijkikj+i∂iΣijkj=0,$
(3)
which is composed of the real and imaginary components:
$γ≡𝔍𝒪=Σijkikjk2,ω≡ℜ𝒪=−∂iΣijkjk2$
(4)
The condition for nonevanescence (or weak evanescence) is that the oscillatory (i.e., imaginary) component of ϕ, characterized by the frequency ω, is much larger than the decaying (i.e., real) component, characterized by the damping γ; that is, the condition |γ/ω| ≪ 1 must be satisfied. This requirement is clearly not satisfied if reported average isotropic and homogeneous parameters are used to describe brain tissues. For typical low-frequency (<10 Hz) white and gray matter conductivity and permittivity (i.e., from Gabriel, Lau, & Gabriel, 1996a, 1996b; εGM = 4.07 · 107 ε0, εWM = 2.76 · 107 ε0, σGM = 2.75 · 10−2 S/m, and σWM = 2.77 · 10−2 S/m, where ε0 = 8.854187817 · 10−12 F/m is the vacuum permittivity), the damping rate γ is in the range of 75–115 s−1, which would be expected to give strong wave damping.

To better understand the effects of brain tissue microstructure and macrostructure on the manifestation of propagating brain waves, it is instructive to consider two idealized tissue models. In the first model (Figure 1A), all brain fibers are packed in a half space aligned in z direction and their number decreases in x direction in a relatively thin layer at the boundary.

Figure 1.

(A) Schematic picture of half-plane packing of fibers. The uniform area of fibers oriented along z direction (shown in green) is bounded by a thin region (magenta) where the inhomogeneity (and the gradient of the conductivity) may be present (a sketch of one possible conductivity profile is shown at the bottom). (B) Schematic picture that can be used as a crude 2-D approximation of fold. The direction of fiber conductivity has only x and z components, and all quantities are assumed to be uniform in y direction.

Figure 1.

(A) Schematic picture of half-plane packing of fibers. The uniform area of fibers oriented along z direction (shown in green) is bounded by a thin region (magenta) where the inhomogeneity (and the gradient of the conductivity) may be present (a sketch of one possible conductivity profile is shown at the bottom). (B) Schematic picture that can be used as a crude 2-D approximation of fold. The direction of fiber conductivity has only x and z components, and all quantities are assumed to be uniform in y direction.

We assume that small cross-fiber currents can be characterized by a small parameter ε and introduce the conductivity tensor as
$Σ=ευευευευευ0ευ0υ,$
(5)
where υ ≡ υ(x). For the υ(x) dependence, we will assume that the conductivity is changing only through a relatively narrow layer at the boundary (as illustrated in Figure 1A) and the conductivity gradient is directed along the x axis. Then, we will look for a solution for the potential φ located in the thin layer of inhomogeneity (i.e., we substitute x → εx) that depends on z and y only as φ = φ(z) + εφ(y).
$ε0:∂∂t∂2φ∥∂z2+υ∂2φ∥∂z2+∂υ∂x∂φ∥∂z=0$
(6)
$ε1:∂∂t∂2φ⊥∂y2+∂υ∂x∂φ⊥∂y+Oε=0$
(7)
where ε0 and ε1 denote the zeroth and first orders of ε power, respectively.

Equation 6 describes a potential along the fiber direction and is a damped oscillator equation that has a decaying solution. However, Equation 7 describes a potential perpendicular to the fiber direction and does not include a damping term, hence it describes a pure wave-like solution that propagates in the thin layer transverse to the main fiber direction. Thus, although this wave-like solution φ has a smaller amplitude than along the fiber action potential φ, it can nevertheless have a much longer lifetime.

We would like to stress that Equations 6 and 7 are here only for illustrative purposes to reiterate a relatively obvious but often overlooked consideration that, under anisotropic inhomogeneous conditions, some direction may happen to be better suited for wave propagation than the others. Those equations should not be considered as the most important part as the complete dispersion relation in Equation 3 will be used to study wave propagation later.

To account for geometric variations, we construct a slightly more complex 2-D model that can be viewed as a very crude approximation of a cortical fold (Figure 1B). The equation for the φ(y) will again be of a wave type, similar to Equation 7, with the addition of a z component of the conductivity gradient and with a similar wave-like solution that will include an additional term induced by the inhomogeneity in z.

The sole purpose of those examples is to provide illustrations that dissipative media with complex structure may show surface wave-like solutions. Surface waves at the boundary of various elastic media have been extensively studied and used in various areas of science (including acoustics, hydrodynamics, plasma physics, etc.) since the work of Lord Rayleigh (1885). The existence of surface waves at the dissipative medium boundary is also known (Kartashov & Kuzelev, 2014).

### Cortical Fold Model

To extend the above illustrative analysis to a more realistic model of brain tissue architecture, we extracted volumetric structural brain parameters from high-resolution anatomical (HRA) MRI datasets as well as brain fiber anisotropy from diffusion-weighted MRI datasets. All anatomical MRI and dMRI datasets are from the Human Connectome Project (Van Essen et al., 2013). We first provide details about the procedures used for generating inhomogeneous and anisotropic components of the permittivity scaled conductivity tensor Σ (more details for all of the processing steps can be found in Galinsky & Frank, 2014, 2015, 2019).

The cortical fold model is represented by an anisotropy tensor σij(r) scaled by an inhomogeneous density ρ(r), such that the total conductivity tensor Σij(r) is defined by
$Σijr=ρrσijr,$
(8)
where σij(r) is an anisotropic, positive semidefinite symmetric tensor and ρ(r) is a continuous function of position vector r. The inhomogeneous density function ρ(r) is estimated from HRA MRI data by processing it with spherical wave decomposition (Galinsky & Frank, 2014; skull stripping, field of view normalization, noise filtering) and then registering to MNI152 space (Fonov et al., 2011; Fonov, Evans, McKinstry, Almli, & Collins, 2009) with SYM-REG (Galinsky & Frank, 2019). The final 1-mm3 182 × 218 × 182 volume is used as the inhomogeneous density ρ(r). For the anisotropy tensor σij(r), several different test cases of increasing complexity were employed, to more clearly elucidate the separate effects of the different conductivity anisotropy and inhomogeneity parameters.

#### Fixed Anisotropy Orientation and Value

The simplest case is that, in which three different values of fixed σ(1), σ(2), and σ(3), fixed (location-independent) anisotropic tensors were assigned to every location in the brain.
$σ1=ε000ε0001σ2=131+2ε1−ε1−ε1−ε1+2ε1−ε1−ε1−ε1+2εσ3=1000ε0001,$
(9)
where (for ε = 0) σ(1) represents the conductivity tensor that will be used in the analytical solution of Equation 22 in the next section, where currents only flow in z direction; σ(2) represents different orientations, where currents are allowed in the direction of 45° relative to all axes; and σ(3) represents more complicated current anisotropy, with crossing currents flowing in x and z directions, but with no currents in the y direction.

#### Varying Anisotropy Orientation and Fixed Anisotropy Value

In the next slightly more complicated case, we kept the anisotropy value fixed but allowed for variations in its orientation. These variations were estimated from actual data. Multiple diffusion direction and diffusion gradient strength MRI data were used to estimate diffusion tensor 𝔻ij (Niethammer, Estepar, Bouix, Shenton, & Westin, 2006). The anisotropy orientation was estimated using eigenvector d(1) of the diffusion tensor 𝔻ij with the largest eigenvalue $λd1$. The anisotropy tensor σij(r) was defined as
$σr=RTσ1R,$
(10)
where the value of anisotropy is constant across the volume (ε = 0, 0.01, and 0.1 were used for different test examples). R is a rotation matrix between directions (0, 0, 1) and d(1) ≡ (x, y, z) that can be expressed as
$R=1−x21+z−xy1+zx−xy1+z1−y21+zy−x−y1−x2+y21+z$
(11)

#### Varying Anisotropy Orientation and Value

For the next level of complexity, we allowed the anisotropy value to also vary along with orientation, and these values were also estimated from actual data. Similarly to the previous case, the diffusion tensor 𝔻ij was estimated for each voxel, and eigenvector d(1) with the largest eigenvalue $λd1$ was used to define the anisotropy tensor major axis. The position-dependent anisotropy value was defined by introducing parameter ε as either $λd2$/$λd1$ or ($λd2$ + $λd3$)/2/$λd1$ (in both cases resulting in axisymmetric form of the conductivity tensor) and also using two different values $λd2$/$λd1$ and $λd3$/$λd1$ for σ22 and σ11.

#### Varying Anisotropy Orientation and Value Estimated from Microstructure Fiber Anisotropy

The microstructure anisotropy at the level of a single cell and its extracellular vicinity may be significantly higher than detected by dMRI estimates. Recent measurements of effective impedance that involve both intercellular and extracellular electrodes (Gomes et al., 2016) show significantly lower conductivity (by orders of magnitude) than conductivity obtained by measurements between extracellular electrodes only. Although attributing this to lower than typically assumed conductivity (or higher impedance) of extracellular medium may be questionable (Barbour, 2017), this clearly confirms high anisotropy for the effective conductivity that should be used for the analysis of effects of intercellular and membrane sources (and these are the most important type of sources for the generation of surface brain waves considered in this paper). This may be important for future generation and refinement of the macroscopic conductivity estimates from microstructure data.

We assumed here that the anisotropic form of $σij1$ with ε = 0.001, 0.01, or 0.1 can be used to describe different levels of intercellular conductivity as well as microstructure extracellular conductivity in the vicinity of cell membranes (where z direction corresponds to the direction of the fiber). Using full-brain tractography results (Galinsky & Frank, 2016), we generated the anisotropic part of the conductivity tensor σ(r) in voxel at r location as an average over all fiber orientations assuming N fibers inside the voxel with Rk orientation matrix for every fiber k, that is,
$σr=1N∑kNRkTσ1Rk,$
(12)
with the complete form of the conductivity tensor Σ again given by Equation 8.

### Wave Trajectory Integration

To obtain the trajectory of a brain wave with frequency ω emitted at a point r0 with an initial wave vector k0, we substitute the inhomogeneous anisotropic conductivity tensor Σ in the wave dispersion relation (imaginary part of Equation 3)
$𝔍𝒟ωk=ωk2+∂iΣijkj=0,$
(13)
and obtain wave ray equations as
$dtdτ=∂𝔍𝒟ωk∂ω=k2,⇒t=∫0k2dτ$
(14)
$dωdτ=−∂𝔍𝒟ωk∂t=0,⇒ω=const$
(15)
$drldτ=∂𝔍𝒟ωk∂kl=2ωkl+∂iΣil,$
(16)
$dkldτ=−∂𝔍𝒟ωk∂rl=−∂l∂iΣijkj.$
(17)
We integrate these equations starting at τ = 0, t = 0, r = r0, and k = k0 to trace the wave trajectory t(τ), r(τ), and k(τ), respectively.
For numerical integration, it is beneficial to split k into magnitude k and direction $k˜$ parts (k = k$k˜$, and |$k˜$|2 = 1) and rewrite the last two equations as
$drldτ=−2∂iΣijk˜jk˜l+∂iΣil,$
(18)
$dk˜ldτ=−∂l∂iΣijk˜j+∂m∂iΣijk˜jk˜mk˜l,$
(19)
where in the last equation the right hand side is orthogonal to $k˜$ to guarantee that |$k˜$|2 = 1, and the algebraic expression in Equation 13 was substituted for the differential equation for wave vector magnitude k.

#### Spherical Cortex Shell Model and Wave Loop Emergence

Brain wave trajectories within the convoluted geometry of the cortex will necessarily be complicated, so it is useful to consider first a highly idealized geometry that provides a clearer intuitive picture, and more readily highlights the emergence, of the wave loops uniquely predicted by our theory. For this, we used a simple spherical shell geometry with shell of fixed thickness h ≈ 1.5–3 mm (based on typical human brain cortical thickness dimensions) spread over a hemisphere of radius R ≈ 75 mm, with all parameters kept constant inside the hemisphere (for r < R) and changing as a function of radius r in a cortical layer. Even without taking into account the known strong anisotropy of neural tissue, these simple geometric considerations provide for the longest waves (with the smallest amount of damping) with
$γω=Σijkikj∂iΣijkj≈Sk2S/hk≈kh∼0.02–0.04.$
Anisotropy (Σ < Σ) will reduce this estimate even further,
$γω=Σijkikj∂iΣijkj≈Σ⊥k⊥2Σ∥/hk⊥≈Σ⊥Σ∥kh
thus further strengthening the condition necessary to support stable waves.
This simple spherical shell model can be used to illustrate the natural appearance of loop-type cortical waves, which can be easily understood from simple geometrical optics arguments. Using the dispersion relation in Equation 3 with only a single component of the conductivity tensor ΣzzS(r), the wave frequency ω can be expressed as
$ω=−1dSkzzrdrk2;$
(20)
where we have neglected the term S$kz2$ because kh ≪ 1. Then, from the geometrical optics ray equations
$drldt=∂ω∂kl,dkldt=−∂ω∂rl,$
(21)
we can get
$dxdt=−2ωkxk2,dkxdt=−dωxdrr,dydt=−2ωkyk2,dkydt=−dωydrr,dzdt=−ω2kz2−k2k2kz,dkzdt=−1zdωz2drr+ω.$
(22)
These equations will generate rays inside the spherical cortex shell showing wave propagation across both the fibers and the conductivity gradient in the cortex subregion where ωdω/dr < 0. For kz = k/$2$ and z = $−ωr/dω/dr$, the wave path has the simplest form—the wave follows the same loop through the cortex over and over again. Different families of frequencies ω and wave vectors k will result in the appearance of cortical wave loops at different cortex locations.

#### Wave Dissipation and Excitation

A useful formulation of the change of the wave energy W is as a path integral along its trajectory, parameterized by t(τ), r(τ), and k(τ):
$W=W0exp−∫0γdisrτ−γexcrτdt=W0exp−∫0Σijrτk˜iτk˜jτ−γexcrτk2τdτ$
(23)
with any appropriate model form of wave excitation γexc. The generality and flexibility of the path integral formulation facilitate the investigation of a broad range of interesting questions of brain wave dynamics, that is, to identify potential active area where wave intensity may grow as a result of certain frequency and/or spatial distributions of neuronal spiky activation, to estimate the spatial extent of coherent activation area as a result of some particular point sources, or to find when and/or how these waves may potentially become important in triggering neuronal firing, that is, to study possible mechanisms of synchronization, feedback, and so forth.

### Nonlinear Theory

The above linear analysis can be extended to include nonlinear effects resulting from a more complex description of the physical processes involving the tissue properties. We assume for simplicity a 2-D symmetric form of the conductivity tensor with constant diagonal terms Σxx and Σyy (where Σyy is along the fiber conductivity, Σxx < Σyy) and position-dependent off-diagonal terms Σxy that are changing linearly with y through a relatively narrow layer at the boundary so that the conductivity gradient exists only inside this layer and is directed along the y axis. We will only be interested in a 1-D solution for the potential φ(x) located in this thin layer of inhomogeneity that can be described by the reduced equation
$∂t+γd∂x2φ+Ω∂xφ=𝓕,$
(24)
where γd = Σxx and Ω = ∂yΣxy.
The source term 𝓕 can be assumed to have a frequency-independent forcing part with a linear growth rate γe representing some averaged input from random spiking activity and an additional term that describes the nonlinear amplitude–phase coupling of the firing rate to the wave field itself (Czurkó et al., 2011; Buzsáki, 2002; Stewart et al., 1992; Fox et al., 1986),
$𝓕=−γeφ−𝒩φ.$
(25)
The solution φ can be sought as a Fourier integral expansion
$φxt=∫akteikx+ωktdk+c.c.,$
(26)
for wave modes with frequencies ωk and wave numbers k (where “c.c.” denotes complex conjugate), which results in a set of coupled equations for time-dependent complex amplitudes ak(t) ≡ a(k, t)
$dakdt=γek2−γdak+1k2𝒩k,$
(27)
where the wave mode frequencies are inversely proportional to the wave number
$ωk=Ω/k,∣k∣>k0=2π/L,$
(28)
and
$𝒩k=12π∫𝒩φe−ikx+ωktdx.$
(29)
represents the nonlinear terms. These nonlinear terms can be characterized in terms of resonant and nonresonant coupling modes, each of which gives rise to unique features of the wave fields.

#### Resonant Coupling of Wave Modes

The nonlinear terms 𝒩k will include a sum of inputs from multiple waves, that is, k = $Σin$ki where n is the order of the nonlinearity. Those resonant conditions will give rise to coupling terms that include various combinations of exp(ik$Σin$ ± ωki)t). Additional requirements for frequency resonances (ωk = $Σin$ ± ωki) produce wave turbulence-like (Nazarenko, 2011; Zakharov, L'vov, & Falkovich, 1992) selection rules for the coupling terms that are similar to phase coupling terms in a ring of connected oscillators (Kuramoto, 2002; Kuramoto & Battogtokh, 2002).

We would like to highlight the importance of the inverse relationship between the frequency and wave number contained in Equation 28. For waves having typical dispersion properties, that is, with the frequencies directly proportional to the wave numbers (ωkkm, where m > 0 is a constant), the maximum oscillatory frequency is increasing and going to infinity with increasing wave numbers. In this case, the nonlinear terms produce a direct cascade of wave energy (Nazarenko, 2011; Zakharov et al., 1992) constantly generating larger and larger frequencies (for brain oscillations' relevant examples and discussion, please see Sheremet, Qin, Kennedy, Zhou, & Maurer, 2019). However, for the inversely proportional wave dispersion, like the waves considered in this paper, the wave energy will be cascaded into smaller frequencies, thus providing a natural mechanism for synchronization of high frequency spiking input and emergence of low frequency rhythms.

This model is also able to characterize another important phenomenon whose existence is supported by an abundance of experimental data—feedback between field potential and firing rate (Czurkó et al., 2011; Buzsáki, 2002; Stewart et al., 1992; Fox et al., 1986). The feedback can be represented through nonlinear coupling. This will be demonstrated using the simplest quadratic form 𝒩(φ) = φ(x, t)2 for the coupling, which can arise through many different processes. More complex feedback can be generated by higher-order coupling terms of course, but that discussion is beyond the scope of this current paper. The quadratic form of coupling results in
$𝒩k=∫δk±k′±k″e−iωk−ωk′−ωk″tak′ak″dk′dk″.$
(30)
Using symmetry conditions ak = $ak*$ and ωk = −ωk (these conditions appear because φ(x) is a real value function, i.e., they are a consequence of φ*(x) = φ(x)), this can be rewritten as
$𝒩k=∫[e−iδω−ωk−k′tak′ak−k′+e−iδω+ωk−k′tak′ak−k′*+e−iδω−ωk+k′tak′ak+k′+e−iδω+ωk+k′tak′ak+k′*]dk′.$
(31)
where δω = ωk − ωk.
The corresponding conditions for the frequency resonances ωk = ±ωk ± ωk allow the expression of the nonlinear resonant coupling $𝒩kR$ by extraction of only the relevant terms as
$𝒩kR∼ak−2ak−1*+ak−1ak1+ak1*ak2,$
(32)
where only three of four wave number resonances appear, as the resonance k + k′ + k″ = 0 is not possible (Nazarenko, 2011; Zakharov et al., 1992), and k−2 = k(3 − $5$)/2, k−1 = k(1 − $5$)/2, k1 = k(−1 − $5$)/2, and k2 = k(3 + $5$)/2 are the real solutions of quadratic equations 1/k ± 1/k′ ± 1/|kk′| = 0.

#### Nonresonant Coupling of Wave Modes

An important addition to these coupling terms arises from the inverse proportionality of frequency and wave number in the dispersion relation in Equation 28. The difference of frequencies of nonlinear nonresonant harmonics is decreasing and going to zero with increasing wave number, thus effectively allowing a closed form expression for the limit of k → ∞, an effect that is absent for coupling of waves with directly proportional dispersion. To illustrate this, we will estimate the nonresonant nonlinear input $𝒩k0nR$ to the k0 wave mode.
$𝒩k0nR=∫[e−iδω1ktakak−k0*+e−iδω2ktakak−k0+e−iδω3ktakak+k0*+e−iδω4ktakak+k0]dk,$
(33)
where
$δω1k=δω0k+ωk−k0,δω2k=δω0k−ωk−k0,δω3k=δω0k+ωk+k0,δω4k=δω0k−ωk+k0.$
where δω0(k) = ωk0 − ωk.
The approximate expression for the forced oscillations solution can be obtained assuming that all the forcing input originates from the scales of k = k0 (the forcing term γe/k2 in Equation 27 is largest when k = k0). Therefore, we will derive the nonresonant input term $𝒩k0nR$ only for k = k0, thus neglecting nonlinear and damping terms for any k > k0 (more correctly, for any k that is not in resonance with k0 or for k > k2 = k0(3 + $5$/2) ≈ 2.618k0). At the limit k → ∞, all frequency deltas δω1–4(k) → ωk0 ≡ ω0 and kk0k, hence approximately, we can estimate the nonresonant term $𝒩k0nR$ as
$𝒩k0nR≈2e−iωk0t∫akak*+ak2dk≈e−iωk0t∫ak+ak*2dk.$
(34)
To estimate forced oscillation terms required in evaluation of the integral in Equation 34, one can write from Equations 27, 31, and 33 that
$dakdt=1k2[eiδω1ktak0ak−k0+eiδω2ktak0ak−k0*+eiδω3ktak0ak+k0+eiδω4ktak0ak+k0*].$
(35)
Looking again for an approximate large k solution (akk0ak for kk0) and keeping only terms that include ak0 (assuming that the amplitude ak0 is small and can be considered constant relative to any of the δω(k) terms), we can approximately write that
$ak≈−iak0k2Σj=14Cjkeiδωjktδωjk$
(36)
$≈−iCkak0k2ωk0eiωk0t,$
(37)
where C1…4(k) and C(k) = ΣCi(k) are complex integration constants that we assume to have random phases with the amplitude independent of k; hence, cross terms in C(k)C(k) (and their complex conjugates) at different frequencies vanish because of the random phases, and only with frequency terms involving k0 are left. Therefore, the nonresonant input term $𝒩k0nR$ in Equation 34 depends on ak0 and t as
$𝒩k0nR≈2C˜3ωk02k03e−iωk0tak0ak0*,$
(38)
where C(k)C*(k) = $C˜$.
A more accurate estimation of $𝒩k0nR$ is possible but requires evaluation of the more complicated integrals of the form
$InR=∫e−iδωi1k−δωi2k−δωi3k±k0tk2k±k02δωi2kδωi3k±k0dk=1ωk02k03∫e−iωk01−1/k±1/k±1tk2k±121−1/k±1/k±12dk,$
(39)
Nevertheless, despite the complexity of these expressions, they have the same eiωk0t asymptotic behavior for t → ∞.
Therefore, an equation for the longest wave length brain mode ak0 that integrates the nonlinear nonresonant input from smaller spatial scales can be written as
$dak0dt=ak0k02γ+βe−iωk0tΩ2∣k0∣ak0*−2α∣ak0∣ak0,$
(40)
where γ describes the excitation strength and β is the strength of nonresonant coupling. The last term (with the parameter α) was included to ensure that coupling does not produce an overall mean field excitation as well as to ensure that, in the limit of vanishing coupling (β = 0), the solution of Equation 40
$ak0t=γC0γexp−γt/k02+2αk02,$
(41)
(where C0 is a constant) has the same 1/$k02$ asymptotic behavior for t → ∞ as the solution of Equation 24 obtained with time and space scale-independent forcing.

#### Nonresonant, Longest Wavelength Brain Modes

Equation 40 can be converted to a system of equations for the amplitude A and phase B (ak0 = AeiB) as
$dAdt=Ak02γ+βAcosB+ωk0t−δAΩ2∣k0∣−2αA2,$
(42)
$dBdt=−βAΩ2k03sinB+ωk0t−δB,$
(43)
where δA and δB were added to introduce tunable phase delays (Abrams & Strogatz, 2004). The system of Equations 42 and 43 includes nonresonant input from all wave modes, but the resonant term in Equation 32 should also be added together with additional equations for resonant wave amplitudes that participate in resonant coupling.

#### Multiple Nonlinear Modes and Critical Points

Finally, we considered a model that combines a chain of inputs from nonlinear modes generated because of resonant terms in Equation 32 into a set of nonresonant mode in Equation 40, which results in a system of equations for mode amplitudes ak for k = k0kN
$dakndt=aknkn2γ+βΩ2∣kn∣e−iωknt+iδakn*−2α∣akn∣akn+λkn2akn−2akn−1*+akn−1akn+1+akn+1*akn+2$
(44)
where the parameter λ describes the strength of resonant coupling between modes and akn = 0 for n < 0 and n > N.

### Datasets and Simulations

To obtain a realistic model of brain tissue architecture, we extracted volumetric structural brain parameters from HRA MRI datasets as well as brain fiber anisotropy from diffusion-weighted MRI datasets. All anatomical MRI and dMRI datasets are from the Human Connectome Project (Van Essen et al., 2013). The details for all the processing steps can be found in Galinsky and Frank (2014, 2015, 2019). More refined procedures for constructing the conductivity tensor and anisotropy in different brain regions, cortical areas in particular, would clearly be beneficial and will be addressed in the future.

The derived inhomogeneous and anisotropic tissue model was used to construct a map of existence conditions for cortex surface waves presented in the Linear Models and Estimates section. The same tissue model complimented by different strategies for estimating the conductivity tensor anisotropy described in the Fixed Anisotropy Orientation and Value section to the Varying Anisotropy Orientation and Value Estimated from Microstructure Fiber Anisotropy section was used for wave packet simulations presented in the Wave Simulation of Cortical Loops section using the simulation procedure described in the Wave Trajectory Integration section.

To estimate actual amplitudes in the human brain for an ensemble of wave packets (without identification and separation of any individual wave packet), we used several multimodal EEG and MRI datasets acquired by the Javitt Group at the Nathan Kline Institute for Psychiatric Research (Javitt & Sweet, 2015; Javitt, 2009). The data include task and resting state EEG and FMRI recordings at 2 msec for about 30 min in total and HRA T1 volumes. We followed the procedure for generation of functional EEG modes from combined EEG and FMRI data described in detail in Galinsky, Martinez, Paulus, and Frank (2018).

Several models were used in exploring the predictions of this theory. (1) The inhomogeneous and anisotropic tissue model supplemented by EEG and FMRI datasets was used for the amplitude and spectra estimates presented in the Cortical Fold Model Estimates section. (2) The nonlinear model with only nonresonant mode coupling in the Nonresonant Coupling of Wave Modes section and Nonresonant, longest wavelength brain modes section was used to study an appearance of critical point and synchronized spiking as well as an appearance of bursts in the Critical Point and Synchronized Spiking section. (3) The nonlinear model that includes both resonant and nonresonant mode coupling in the Multiple Nonlinear Modes and Critical Points section was used to study an appearance of hypersynchronicity and low-frequency quasiperiodicity in the Hypersynchronicity and Low-Frequency Quasiperiodicity section. (4) An experimental in vivo data that used an implanted multielectrode array for measuring population activity within the visual cortex of awake, freely viewing ferrets at three different developmental stages (Fiser, Chiu, & Weliky, 2004) was used to compare with the simulated correlation properties of unusual inverse proportionality of the frequency and wave number presented in the Evidence in Experimental Data section.

## RESULTS

### Linear Models and Estimates

To provide an illustration of where the conductivity anisotropy and inhomogeneity can form appropriate conditions for cortex surface wave generation, we created plots generated from actual MRI data in a normal human brain that can be used to characterize the ratio |γ/ω| using the axisymmetric form of the conductivity tensor with varying anisotropy orientation and value (see Varying Anisotropy Orientation and Value section). To construct the ratio, we calculated two vectors, ∂iΣij and Σijki, and compared their norms. For wave vector k, we used a vector with the same direction as in ∂iΣij vector and magnitude |k| = |∇Σ|/|Σ|; that is, our intention is to compare norms of dissipative and wave-like terms at kh ≈ 1 where h is the cortical thickness. Figure 2 shows plot of the ratio of |Σijki| and |∂iΣij| at two different depths inside the brain with dissipative regime in (A) the inner cortex or the innermost layers of the cerebral cortex and subcortical areas as well as wave-like conditions in (B) the outer cortex or outermost layers of gray matter of the cerebral hemispheres. The inner cortex shown in Figure 2A clearly displays prevalence of dissipation with the ratio below 1, whereas the outer cortex shown in Figure 2B allows for wave-like cortex activity in most of the locations with the ratio above 2.

Figure 2.

A ratio of wave-like to dissipation effects shown on cortical surfaces at two different values of depth from the top of the brain to inside the cortex for inhomogeneous and anisotropic normal human brain tissue model described in then Cortical Fold Model section. The color scheme shows the range for this ratio from 0 to 4 and uses shades of green to mark regions where the ratio is less than 1, that is, where the dissipative term dominates; shades of red where the ratio is between 1 and 2; and shades of blue where the ratio is more than 2, that is, where the wave-like term is more than two times dominant.

Figure 2.

A ratio of wave-like to dissipation effects shown on cortical surfaces at two different values of depth from the top of the brain to inside the cortex for inhomogeneous and anisotropic normal human brain tissue model described in then Cortical Fold Model section. The color scheme shows the range for this ratio from 0 to 4 and uses shades of green to mark regions where the ratio is less than 1, that is, where the dissipative term dominates; shades of red where the ratio is between 1 and 2; and shades of blue where the ratio is more than 2, that is, where the wave-like term is more than two times dominant.

Both anisotropy and inhomogeneity are important for the existence of the cortex surface waves. For example, for inhomogeneous but isotropic tissue, the conductivity tensor Σij will simply be Sδij, where S is a scalar inhomogeneous conductivity. Therefore, both the phase velocity
$vph=ωkkk=−∇·Σ·kk4k=−∇S·kk4k$
and the group velocity
$vgr=∂ω∂k=−∇·Σk2+2∇·Σ·kk4k=−∇Sk2+2∇S·kk4k$
will include terms ∇S and ∇S · k, meaning that those waves are not able to propagate normally to the local conductivity gradient ∇S. This restriction is absent in cortex areas when both anisotropy and inhomogeneity are present (because of tensor products in ∇ · Σ and ∇ · Σ · k).

### Wave Simulation of Cortical Loops

A remarkable consequence of the weakly evanescent cortical wave theory—the emergence of stable wave loops described theoretically in the previous section—can be demonstrated through numerical simulation of wave propagation in a thin dissipative inhomogeneous and anisotropic cortical layer. Wave packets were simulated using Equation 21 where the general form of anisotropic dispersion relation in Equations 3 and 4 is used.

This method was used to randomly generate wave trajectories in a regime with γeffLloopgr ≤ 1 where γeff is an effective wave dissipation rate, that is, the difference between average dissipation and spiking activation rates for a wave packet propagating with group velocity υgr along some characteristic loop of Lloop length. The spatial snapshots of the dynamics of the behavior of these waves are shown for two specific geometries in Figure 3: (1) an idealized spherical shell and (2) actual brain morphology constructed directly from MRI data. The idealized spherical model illustrates in the most straightforward way two remarkable features of the cortical waves: Their initial trajectory appears highly complex and spatially extended yet converges on a stable loop pattern. This persistent localized cortical wave pattern is precisely the simple loop pattern predicted by Equation 22, despite the complex initial spatio-temporal pattern of the initial wave trajectory. As further predicted, the simulations informed by the real human data with the same cortex fold geometry as in Figure 2 (center and right) also produce stable loop structures, now embedded within the complex geometry of the cortical folds. All wave packets were initialized with random parameters assuming the presence of spiky activation sources (not shown). The patterns structure is defined by an interplay of anisotropy and inhomogeneity and may be as simple as a circular loop for a symmetric spherical geometry or may be rather intricate for a cortical fold geometry with a complex anisotropy.

Figure 3.

Examples of wave trajectories obtained in simulation of wave propagation in real-data cortical fold tissue model of the Cortical Fold Model section. A–C show the complete trajectories, and D and E show the emergent stable wave loops. The spherical cortex shell model of the Spherical Cortex Shell Model and Wave Loop Emergence section is used for A and D, and the cortical fold model of the Cortical Fold Model section is used for B, C, E, and F. The colors encode wave propagation: red, left/right; green, anterior/posterior; and blue, dorsal/ventral.

Figure 3.

Examples of wave trajectories obtained in simulation of wave propagation in real-data cortical fold tissue model of the Cortical Fold Model section. A–C show the complete trajectories, and D and E show the emergent stable wave loops. The spherical cortex shell model of the Spherical Cortex Shell Model and Wave Loop Emergence section is used for A and D, and the cortical fold model of the Cortical Fold Model section is used for B, C, E, and F. The colors encode wave propagation: red, left/right; green, anterior/posterior; and blue, dorsal/ventral.

To demonstrate the generality and consistency of this result, Figure 4 shows the final persistent wave patterns for nine simulations that were initialized with wave packets of random parameters (frequency, wave number, location, etc.) as well as with different types of anisotropy (as described in the Cortical Fold Model section). All clearly show an emergence of localized persistent closed-loop patterns at different spatial and temporal scales, from scales as large as global—the whole brain rotational wave activity experimentally detected in Muller et al. (2016)—to as small as the resolution used for the cortical layer thickness detection. One natural line of investigation, beyond the scope of the current work, would be an in-depth comparison of loop patterns presented in Figures 3 and 4 with experimental findings. The patterns in the plots are identified as self-connected closed parts of wave trajectory that are followed at least several times over and over again. As the cortical surface of the brain includes relatively loosely connected surfaces of the left and right hemispheres, it seems likely that unilateral patterns will be encountered more often than structures that symmetrically span both hemispheres. These figures show that wave packets emitted at random spatial location at random propagation directions with random wave parameters (frequency and wave number) self-organize in loop patterns that are defined by anisotropy and inhomogeneity structures of the tissue.

Figure 4.

Additional random examples of wave loops obtained in simulation of wave propagation in real-data cortical fold tissue model of the Cortical Fold Model section. The colors again encode wave propagation: red, left/right; green, anterior/posterior; and blue, dorsal/ventral. The collection of movies showing wave packet propagation dynamics and loop formation is available at https://www.dropbox.com/sh/aads9hluqqzqfi6/AAC4Jzt8a-pH4FfJEf_6JUnBa?dl=0__;!!Mih3wA!VDUnBPduahevOhOiOqYHTZZHh4H7dA_b9_L4xp0hCJOVds515-icRsBCuQp6$and https://ucsdcloud-my.sharepoint.com/:f:/g/personal/vgalinsky_ucsd_edu/EhzB6j-1WLVDmiV094llcYoBKIK8OzsVDO61o3EuWJvV1A?e=tutVY0__;!!Mih3wA!VDUnBPduahevOhOiOqYHTZZHh4H7dA_b9_L4xp0hCJOVds515-icRjonU9D1$.

Figure 4.

Additional random examples of wave loops obtained in simulation of wave propagation in real-data cortical fold tissue model of the Cortical Fold Model section. The colors again encode wave propagation: red, left/right; green, anterior/posterior; and blue, dorsal/ventral. The collection of movies showing wave packet propagation dynamics and loop formation is available at https://www.dropbox.com/sh/aads9hluqqzqfi6/AAC4Jzt8a-pH4FfJEf_6JUnBa?dl=0__;!!Mih3wA!VDUnBPduahevOhOiOqYHTZZHh4H7dA_b9_L4xp0hCJOVds515-icRsBCuQp6$and https://ucsdcloud-my.sharepoint.com/:f:/g/personal/vgalinsky_ucsd_edu/EhzB6j-1WLVDmiV094llcYoBKIK8OzsVDO61o3EuWJvV1A?e=tutVY0__;!!Mih3wA!VDUnBPduahevOhOiOqYHTZZHh4H7dA_b9_L4xp0hCJOVds515-icRjonU9D1$.

### Cortical Fold Model Estimates

To determine a possible energy distribution across different frequencies for these cortical waves, some knowledge about the forcing term 𝓕 in Equation 1 is required. A rough estimate for this distribution in the form of a power spectra scaling can be carried out using some simple assumptions. Assuming that the forcing consists of spiking input localized at random locations and times, it can be described as a sum of delta functions, 𝓕 = ∑iAiδ(tti)δ(rri), corresponding to a flat forcing frequency spectrum in the Fourier domain. Then, from the last term in Equation 2 |(∂iΣij)(∂jφ)| ∼ ω$φω2$ ∼ const, we can estimate the exponent α = 2 in a power law scaling of the potential φ frequency spectrum (i.e., for |φω|2 ∼ ω−α).

The presence of the cortical wave loops described above can modify this ω−2 dependence and thus modify the spectrum in such a way that spectrum peaks are generated that correspond to these loop wave currents. From Equation 20, we can estimate the range of frequencies where those loops can possibly be present. Taking 1/r ∼ 1/R, dS/dr ∼ Σzz/h, z = $−ωr/dω/dr$$Rh$ gives a frequency estimate as f = ω/2π ∼ Σzz/(2πk$2Rh$), hence for the largest and smallest wave numbers defined by the smallest (1.5 mm) and largest (75 mm) loop radii, the frequency f spans the range of 1.2–92 Hz (υph ≈ 0.002–7 m/s). The large-scale circular cortical waves of Muller et al. (2016; 9–18 Hz, 2–5 m/s) are clearly inside this range.

Figure 5 shows a fixed time snapshots of distribution of electric fields on the skull for two subjects (A, B) and on one of the cortical surfaces (C, D) as well as temporal evolution of the scalar potential at a single point for one of the alpha-band wave oscillation modes in the cortex (E, F). The time course estimates show possible maximum wave amplitude variations of up to 5–10 mV. However, it should be recognized that we present this time course merely to demonstrate that the predicted amplitudes fall within an experimentally observable range. However, the utility of such a plot is problematic, because it represents only the amplitude at a single spatial location of an extended, coherent, and evolving wave pattern. In fact, this has implications for highly localized in vivo measurements not uncommon in EEG studies, which might detect a “ripple” that is, in fact, the brief measurement of a passing wave. This might explain the appearance of “ripples” in memory studies, for example, Squire, Genzel, Wixted, and Morris (2015).

Figure 5.

Examples of electric fields and scalar potentials for two subjects from a set of multimodal EEG and MRI datasets (Javitt & Sweet, 2015; Javitt, 2009): (A, B) fixed time snapshots of electric field distribution on the skull, (C, D) fixed time snapshots of electric field distribution on one of the cortical surfaces, and (E, F) temporal evolution of the electric potential at a single point for alpha-band wave oscillations in the cortex. The color corresponds to the electric field direction (red, left/right; green, anterior/posterior; and blue, dorsal/ventral). The temporal evolution shows possible maximum wave amplitude variations of up to 5–10 mV. (Note that the circular isocontours of these fields are not to be confused with the wave loops that are the subject of this paper.)

Figure 5.

Examples of electric fields and scalar potentials for two subjects from a set of multimodal EEG and MRI datasets (Javitt & Sweet, 2015; Javitt, 2009): (A, B) fixed time snapshots of electric field distribution on the skull, (C, D) fixed time snapshots of electric field distribution on one of the cortical surfaces, and (E, F) temporal evolution of the electric potential at a single point for alpha-band wave oscillations in the cortex. The color corresponds to the electric field direction (red, left/right; green, anterior/posterior; and blue, dorsal/ventral). The temporal evolution shows possible maximum wave amplitude variations of up to 5–10 mV. (Note that the circular isocontours of these fields are not to be confused with the wave loops that are the subject of this paper.)

Evidence for the existence of these wave-loop-induced spectral peaks is shown in Figure 6, which shows the spectral power of the EEG signal for six subjects (Martínez et al., 2015; Woldorff et al., 2002), averaged over all sensors. The dashed lines outline the predicted f−2 in the lower (f < 1.2 Hz) and higher (f > 92 Hz) parts of the spectra. The dashed–dotted vertical lines denote the frequency range where the cortical loops may exist, agreeing very well with typically observed EEG excessive activity range, from low-frequency delta (0.5–4 Hz) to high-frequency gamma (25–100 Hz) bands.

Figure 6.

Spectral power of EEG signal collected with 64-sensor array and averaged over all sensors for six independent subjects is shown in six panels. The dashed lines outline the predicted f−2 in the lower (f < 1.2 Hz) and higher (f > 92 Hz) parts of the spectra. The dashed–dotted vertical lines denote the frequency range where the cortical wave loops may be generated. Both the slope and the range agree very well with typically observed values.

Figure 6.

Spectral power of EEG signal collected with 64-sensor array and averaged over all sensors for six independent subjects is shown in six panels. The dashed lines outline the predicted f−2 in the lower (f < 1.2 Hz) and higher (f > 92 Hz) parts of the spectra. The dashed–dotted vertical lines denote the frequency range where the cortical wave loops may be generated. Both the slope and the range agree very well with typically observed values.

### Critical Point and Synchronized Spiking

We will investigate the behavior of the nonlinear, nonresonant, longest wavelength brain modes in Equations 42 and 43 and will consider the complete system in Equation 44 later.

Figure 7 shows the results of numerical solution of the system in Equations 42 and 43 for several different sets of parameters. The time evolution of the highest frequency, longest wavelength mode exhibits a variety of types of oscillatory behavior, ranging from slightly nonlinear modified sinusoidal shapes to more nonlinear looking shapes similar to network-attributed alpha waves or μ-shaped oscillations (Buzsáki, 2006). Increase in the level of activation γ produces nonlinear signals with the spike-like shape of a single neuron firing.

Figure 7.

Examples of nonlinear waveforms (A–C) and emergence of synchronized spiking (D, E) obtained by simulation of the system in Equations 42 and 43. A–C display transformation from weakly nonlinear oscillations shown by blue dotted lines for γ = 0.75 to more strongly nonlinear regime (solid line; γ = 1.5 [A, B] and 2.25 [C]). D–F show the strongest nonlinear spiking-like time evolution of the potential ϕ (solid line; γ = 2.55 [D, E] and 2.96 [F]) and its transformation to nonoscillatory (blue dotted line) regime for γ = 3 (for all plots, the values of ωk0, k0, α, and β were set to be equal to 1, and time and amplitude units are arbitrary).

Figure 7.

Examples of nonlinear waveforms (A–C) and emergence of synchronized spiking (D, E) obtained by simulation of the system in Equations 42 and 43. A–C display transformation from weakly nonlinear oscillations shown by blue dotted lines for γ = 0.75 to more strongly nonlinear regime (solid line; γ = 1.5 [A, B] and 2.25 [C]). D–F show the strongest nonlinear spiking-like time evolution of the potential ϕ (solid line; γ = 2.55 [D, E] and 2.96 [F]) and its transformation to nonoscillatory (blue dotted line) regime for γ = 3 (for all plots, the values of ωk0, k0, α, and β were set to be equal to 1, and time and amplitude units are arbitrary).

It is interesting that this spiking-like solution of system in Equations 42 and 43 appears near the critical point; the oscillatory state (red lines in the bottom row of Figure 7) undergoes bifurcation and transitions to nonoscillatory regime (blue lines in the bottom row of Figure 7) as γ reaches the value above some critical point. To illustrate the reason for this transition, we will consider the simplest case of δA = 0 and δB = π/2 (although different δA,B values can be used for a similar analysis as well). The nonoscillatory regime can be reached if dA/dt → 0 and dB/dt → −ωk0 as t → ∞. Then, from Equations 42 and 43, one can write that, at t → ∞,
$γA−ωk0A−2αA2=0,βAcosB0=−ωk0,$
where B0 is some arbitrary constant phase. Therefore, the nonoscillatory state requires that γ satisfies to
$γ=ωk01−2αβcosB0.$
(45)
Hence, for
$ωk01−2αβ<γ<ωk01+2αβ,$
(46)
the nonoscillatory solution is not possible. The simulations shown in the bottom panels of Figure 7 confirm that the critical γ value is indeed 3 when ωk0 = α = β = 1. A similar analysis when δA = δB gives the critical γ value equals to (2 – cos(π/3))/sin(π/3) ≈ 1.732.

As a next step, we employed a more complex expression for the total input from the nonresonant terms by including a sum of all InR integrals in Equation 39 instead of a single eiωk0t exponent input. Figure 8 shows simulation results for several parameter sets with the same values as were used for plots of Figure 7k0 = β = k0 = 1, δA = 0). The numerical solution shows a more complex behavior that now includes modulation of spiking rate with a lower frequency and emergence of burst-like train of spikes, effects often observed in different types of neuronal activity (Gerstner, Kistler, Naud, & Paninski, 2014).

Figure 8.

Examples of spiking rate modulation (A, C) and emergence of burst-like train of spikes (B, D) obtained by simulation of the system in Equations 42 and 43 when the exponential term was replaced by InR integrals in Equation 39 with the region of integration set to 50k0 < k < 1000k0 (for all plots, the values of ωk0, k0, α, and β were set to be equal to 1, and time and amplitude units are arbitrary).

Figure 8.

Examples of spiking rate modulation (A, C) and emergence of burst-like train of spikes (B, D) obtained by simulation of the system in Equations 42 and 43 when the exponential term was replaced by InR integrals in Equation 39 with the region of integration set to 50k0 < k < 1000k0 (for all plots, the values of ωk0, k0, α, and β were set to be equal to 1, and time and amplitude units are arbitrary).

#### Hypersynchronicity and Low-Frequency Quasiperiodicity

The system in Equation 44 includes input from multiple nonlinear modes; therefore, it can result in a more complex solution behavior than single nonlinear mode solutions of Figures 7 and 8. Figures 9 and 10 show results of numerical simulation of the system in Equation 44, clearly indicating that weak nonlinear resonant coupling between only three modes with frequencies of ωk0, 2ωk0/(1 + $5$), and 2ωk0/(3 + $5$) is capable of explaining an emergence of periodic activity with frequencies up to 100–1000 times lower than the linear frequencies of participating modes. We would like to emphasize again that the system in Equation 44 cannot be separated into traditional “slow” and “fast” subsystems, hence the low-frequency component cannot be explained by a modulation (Rinzel, 1987) of the “fast” subsystem with oscillations of the “slow” part.

Figure 9.

Emergence of low-frequency modulation of high-frequency synchronized spiking as a result of an increase of weak resonant coupling obtained by simulation of the multimode system in Equation 44. A–C show the complete simulations, and D–F show the blowup of the region from t = 19300 to t = 19500. The total potential ϕ is plotted with black, and different colors show the oscillations of the individual modes. The value of γ is 1.535, which is sufficiently far from criticality but nevertheless large enough to modify an effective period for k0 mode to be close to that of k1 (for all plots, the values of ωk0, k0, α, and β were set to be equal to 1, and δA = δB = δ = 3π/4; time and amplitude units are arbitrary).

Figure 9.

Emergence of low-frequency modulation of high-frequency synchronized spiking as a result of an increase of weak resonant coupling obtained by simulation of the multimode system in Equation 44. A–C show the complete simulations, and D–F show the blowup of the region from t = 19300 to t = 19500. The total potential ϕ is plotted with black, and different colors show the oscillations of the individual modes. The value of γ is 1.535, which is sufficiently far from criticality but nevertheless large enough to modify an effective period for k0 mode to be close to that of k1 (for all plots, the values of ωk0, k0, α, and β were set to be equal to 1, and δA = δB = δ = 3π/4; time and amplitude units are arbitrary).

Figure 10.

Emergence of spike trains in simulation of the multimode system in Equation 44 when γ is sufficiently close to criticality. A weak coupling produces jumps from subcritical to supercritical regimes with regular low-frequency quasiperiodicity. A–C show the complete simulations, D–F show the blowup of the region from t = 16600 to t = 20000, and G–I show the blowup of the region from t = 19300 to t = 20000. Values of γ, δA, and δB are the same in each column. The total potential ϕ is plotted with black, and different colors show the oscillations of the individual modes (for all plots, the values of ωk0, k0, α, and β were set to be equal to 1 and the resonant coupling λ was 0.05; time and amplitude units are arbitrary).

Figure 10.

Emergence of spike trains in simulation of the multimode system in Equation 44 when γ is sufficiently close to criticality. A weak coupling produces jumps from subcritical to supercritical regimes with regular low-frequency quasiperiodicity. A–C show the complete simulations, D–F show the blowup of the region from t = 16600 to t = 20000, and G–I show the blowup of the region from t = 19300 to t = 20000. Values of γ, δA, and δB are the same in each column. The total potential ϕ is plotted with black, and different colors show the oscillations of the individual modes (for all plots, the values of ωk0, k0, α, and β were set to be equal to 1 and the resonant coupling λ was 0.05; time and amplitude units are arbitrary).

In Figure 9, the high-frequency spiking is generated with the level of activation γ = 1.535. This activation level is yet relatively far from criticality but produces spikes with an effective rate that is close to the next linear resonance frequency. The first, second, and third columns clearly show that a small increase of the resonant coupling (0.001, 0.01, and 0.05, respectively) results in appearance of component with a significantly lower frequency.

Figure 10 shows several simulations with the level of activation that is close to criticality for the selected set of parameters in each column. The small resonant coupling λ = 0.05 in this case results in a more profound effect of quasiperiodic shift of oscillations back and forth from subcritical to supercritical regimes effectively turning spiking on and off with low frequency. Prediction of the actual period of nonlinear low-frequency oscillations from the model parameters, for example, distances from the critical points and other resonances, phase delays, and so forth, is an interesting open question that will be considered in future work.

#### Evidence in Experimental Data

As the dispersion relation in Equation 3 (or Equation 4) shows inverse proportionality of the frequency and the wave number, it is interesting to see if signatures of this unusual wave dispersion can be found in experimental data that involve spatio-temporal correlation measurements (Smith & Kohn, 2008; Fiser et al., 2004; Destexhe, Contreras, & Steriade, 1999). The type of proportionality between the frequency and the wave number in the linear wave dispersion relation is not directly related to the temporal/spatial correlation relationships, and it is possible for waves with the inversely proportional dispersion to encounter either directly or inversely proportional relationships between temporal and spatial correlations using different single-point wave correlation properties (e.g., different slopes of power spectra, different spectral frequency ranges). Nevertheless, using a Fourier integral expansion in Equation 26 for the potential φ(x, t) of the waves with the dispersion relation in Equation 4, one can write the two-point spatial Cx(ζ, t1, t2) and temporal Ct(x1, x2, τ) correlation functions as
$Cxζt1t2=∫−∞∞φxt1φx+ζt2dx=∫k1k2ak2cosωkt1−t2+kζdk,$
(47)
$Ctx1x2τ=∫−∞∞φx1tφx2t+τdt=∫k1k2ak2cosωkτ+kx1−x2dk,$
(48)
where k1 and k2 select an appropriate spectral range and |ak|2 defines an appropriate wave power spectrum shape. The temporal and spatial correlations presented in Figure 3 of Fiser et al. (2004) are in a very good agreement with wave correlations expressed in the above form. Figure 3 of Fiser et al. (2004) includes pairs of temporal/spatial correlation functions for three different age groups. There are different trials in each age group obtained with differently correlated visual stimulus (dashed, dotted, and solid lines). The in-group differences of correlations can obviously be attributed to the stimulus difference. What is interesting with respect to our model is that a difference between groups clearly reproduces, appropriate for our wave model, dependence between temporal and spatial correlations in each group. For example, temporal correlations show more tightly grouped correlations in a narrow region close to the origin progressing from P30–P32 to P44–P45 to P83–P90. At the same time, spatial correlations clearly show an increase of the absolute value and a decrease of the slope (as expected from our model) for the same group progression. To graphically illustrate this fact, we generated simulated correlations for our wave model using the same parameters for both the temporal and spatial correlations and plotted the resulting curves on top of the original figure in red color. Figure 11 shows that the simulated correlations follow the experimental trend very well.
Figure 11.

Spatio-temporal pattern of stimulus-evoked and spontaneous visual cortical activity in awake-behaving ferrets and for brain wave simulations. The experimental data from awake ferrets at three different ages are reprinted from Fiser et al. (2004) showing the correlation functions computed for dark spontaneous activity (solid black lines) as well as evoked activity to natural scene (dashed black lines) and random noise films (dotted black lines). Thin horizontal black lines show plots of correlation functions computed at each age for each condition using randomly shuffled binned spikes. The columns P30–P32, P44–P45, and P83–P90 describe developmental changes among the three age groups (mode details are available from Fiser et al., 2004). The clear trend of the temporal correlations' shift to the smaller period oscillations is evident in the left-to-right progression of the top row. At the same time, the bottom row shows, for the spatial correlations, an overall increase of the correlation values and the spread of the correlations to larger spatial distances. The solid red lines show the simulated correlation functions for inversely proportional dispersion wave considered in this paper using the same parameters for temporal and spatial simulated correlation functions in each row. The simulated functions show good agreement and clearly follow the experimental trend.

Figure 11.

Spatio-temporal pattern of stimulus-evoked and spontaneous visual cortical activity in awake-behaving ferrets and for brain wave simulations. The experimental data from awake ferrets at three different ages are reprinted from Fiser et al. (2004) showing the correlation functions computed for dark spontaneous activity (solid black lines) as well as evoked activity to natural scene (dashed black lines) and random noise films (dotted black lines). Thin horizontal black lines show plots of correlation functions computed at each age for each condition using randomly shuffled binned spikes. The columns P30–P32, P44–P45, and P83–P90 describe developmental changes among the three age groups (mode details are available from Fiser et al., 2004). The clear trend of the temporal correlations' shift to the smaller period oscillations is evident in the left-to-right progression of the top row. At the same time, the bottom row shows, for the spatial correlations, an overall increase of the correlation values and the spread of the correlations to larger spatial distances. The solid red lines show the simulated correlation functions for inversely proportional dispersion wave considered in this paper using the same parameters for temporal and spatial simulated correlation functions in each row. The simulated functions show good agreement and clearly follow the experimental trend.

## DISCUSSION

Loop-like waves have been observed in measured human brain activity but never formally quantified, nor their existence and characteristics ever explained by a physical model. Our paper derives cortex wave dispersion relation from well-known relatively basic physical principles and provides illustrations of why and how cortical tissue inhomogeneity and anisotropy influence propagation magnitude, time scales, and directions and support extended and highly structured regions of existence in dissipative media using simple 1-D and 2-D anisotropic models as well as more realistic full-brain models based on a set of parameters extracted from real diffusion and structural MRI acquisitions. The wave properties (frequency ranges, phase and group velocities, possible spectra) are then compared with real EEG wave acquisitions. Finally, examples of wave propagation are studied analytically and numerically using a simple idealized but informative spherical shell cortex model (i.e., thin inhomogeneous layer around a sphere with homogeneous anisotropic conducting medium) as well as a more realistic, anisotropic, inhomogeneous full-brain model with actual cortical fold geometry that clearly shows the emergence of localized persistent wave loops or rotating wave patterns at various scales, including at similar scales of global rotation recently detected experimentally (Muller et al., 2016).

An important aspect of the cortical wave model we presented is that it is based on relatively simple but physically motivated averaged electrostatic properties of human neuronal tissue within realistic data-derived brain tissue distributions, geometries, and anisotropy. Although the source of these averaged tissue properties includes the extraordinarily complex network of neuronal fiber connections supporting a multitude of underlying cellular, subcellular, and extracellular processes, we demonstrate that the inclusion of such details for the activation/excitation process is not necessary to produce coherent, stable, macroscopic cortical waves. Instead, a simple and elegant physical model for wave propagation in a thin dissipative inhomogeneous and anisotropic cortical layer of these averaged properties is sufficient to predict the emergence of coherent, localized, and persistent wave loop patterns in the brain.

First, we found evidence for wave loops using simulated data with a variety of models for anisotropy and conductivity values. Then, we used real data to obtain a spatial and temporal representation of these waves. We found that the ranges of parameters for the waves produced by our model are in agreement with those presented by several studies (Muller et al., 2018; Zhang et al., 2018; Lubenov & Siapas, 2009) that present evidence that electrostatic field activity in several areas of animal and human brains are traveling waves that can affect neuronal activity by modulating the firing rates (Czurkó et al., 2011; Stewart et al., 1992; Fox et al., 1986) and may possibly play an important functional role in diverse brain structures (Weiss & Faber, 2010). The natural self-organization of these traveling waves into loop-like structures that our model produces agrees well with recently detected (Muller et al., 2016) cortical wave activity spatio-temporally organized into circular wave-like patterns on the cortical surface. Self-propagation of endogenous electric fields through a physical cut in vitro when all mechanisms of neuron-to-neuron communication has been destroyed (Chiang et al., 2019; Shivacharan et al., 2019; Qiu et al., 2015; Zhang et al., 2014) can potentially be attributed to these waves as well.

Our results are consistent with the experimental results presented in Muller et al. (2016). We note that, in the experimental detection of circular wave-like patterns shown in Muller et al. (2016), one of the primary motivating studies for our work, the actual amplitudes for individual wave packets were not reported. We therefore only used the propagation properties (e.g., wave velocities, frequencies, and wavelengths) given in that paper for consistency checks with our model estimates. Although our order of magnitude estimates for wave loop amplitudes are consistent with those of the typical local oscillating potentials, it is our view that the main challenge in the detection and identification of signatures of those wave-like loop patterns in Muller et al. (2016) is not the local amplitude measurements but rather the separation of these oscillations from the other significant components in the signals and the demonstration that their spatial and temporal coherence is in agreement with expected (i.e., predicted based on the local tissue characteristics) wave propagation properties.

The abundance of regular and bursting oscillatory patterns across a wide range of spatial and temporal scales of brain electromagnetic activity makes a question of their interaction an important issue that has been widely discussed in the literature (Gerstner et al., 2014; Buzsáki, 2006; Buzsáki & Draguhn, 2004). The standard approach involves representing the brain as a large network of coupled oscillators (Goel & Ermentrout, 2002; Frank, Daffertshofer, Peper, Beek, & Haken, 2000) and using this as a testbed for the study of network wave propagation, mechanisms of synchrony, possibly deriving some mean field equations and properties, and so forth. However, such models are necessarily descriptive, and their relationship to actual physical properties of either actual brain tissue properties or the electromagnetic waves they support is tenuous.

In this paper, we employ a different approach that uses properties of brain waves in realistic brain tissue types and architectures derived in a general form from relatively basic physical principles. We demonstrate that the peculiar inverse proportionality of the wave linear dispersion combined with nonlinear resonant and nonresonant coupling of multiple wave modes produces a remarkably feature-rich nonlinear system that is able to reproduce many seemingly unrelated regimes of brain activity that have been observed experimentally throughout a wide range of spatial and temporal scales. The different regimes include high-frequency spiking-like activity occurring near the critical point of the equation that integrates multiple nonresonant wave modes and low-frequency oscillations that emerge when weak resonant coupling is present in the vicinity of the critical point. The strongly nonlinear regime exists sufficiently close to the critical point where the solution bifurcates from oscillatory to nonoscillatory behavior. The weak resonant coupling then demonstrates a mechanism that constantly moves the system back and forth from subcritical to supercritical domains turning the spiking on and off with low-frequency quasiperiodicity.

To describe this complex behavior, we show for the first time that the inverse proportionality of frequency and wave number in the brain wave dispersion relation permits the characterization of a limiting form for the signals in terms of a large number of wave modes as a summation of nonresonant wave harmonics, thus allowing a closed analytical form of the nonlinear equations that integrates and thus combines the collective nonresonant input from multiple wave modes. Following the ideas of wave turbulence (Nazarenko, 2011; Zakharov et al., 1992), we also show that the resonant coupling between those high-frequency nonlinear wave modes can provide an effective universal mechanism for the emergence of low-frequency wave rhythms.

In turn, we know from past work that coherence is important for brain communication (Fries, 2015). Thus, we used modeling to probe the coherence between these brain waves using nonlinear wave dynamics. The nonlinear analysis of multimode wave dynamics shows that the unusual dispersion properties of those waves provide a universal physical mechanism for emergence of low frequencies from high-frequency oscillations. The simple quadratic nonlinearity introduced as a coupling source for the wave model allowed the derivation of an equation for a nonlinear form of those waves by taking a limit for a large number of nonresonantly interacting wave modes, which we emphasize is a limit that exists only because of unusual dispersion properties of the waves. The collective input from those nonresonant modes results in nonlinear spiking-like solutions of this equation and an existence of a bifurcation point from oscillatory to nonoscillatory regime. The multimode nonlinear system that includes both nonresonant and resonant coupling between multiple modes shows emergence of low-frequency modulations as well as strongly nonlinear low-frequency quasiperiodic oscillations from subcritical to supercritical regimes.

We would like to emphasize that all variety of models used for a description of action potential neuron spikes, starting from the seminal model by Hodgkin and Huxley (1952) and finishing with many dynamical integrate-and-fire models of neuron (Gerstner et al., 2014), are based on an approximation of several local neuron variables, for example, membrane currents, gate voltages, and so forth, and defining the relations between these local properties. Contrary to this and rather unexpectedly, Equation 40 is obtained through an integration of a large number of oscillatory brain wave modes nonresonantly interacting in inhomogeneous anisotropic media and shows spiking pattern solutions emerging as a result of this nonresonant multimode interaction rather than as a consequence of empirical fitting of a nonlinear model to several locally measured parameters. It is also important to note that Equation 40 cannot be separated into “fast” and “slow” parts as typically required for the functioning of “traditional” neuron models. Because of this, we would like to reiterate that this equation should not be viewed as a single neuron model and should not be considered as an alternative to any of the single neuron models (Morris & Lecar, 1981; Nagumo, Arimoto, & Yoshizawa, 1962; FitzHugh, 1961). It describes a mechanism for the generation of synchronous spiking activity as a result of a collective input from many nonresonant wave modes. The transition to the synchronous spiking activity occurs in the vicinity of the critical point where a bifurcation from oscillatory to nonoscillatory state happens, thus indirectly supporting the subcriticality hypothesis (Wilting & Priesemann, 2018) of brain operation.

We would like to mention two new, rather important, and not entirely obvious features that appear in the nonlinear system in Equation 44, but absent for phase-coupled oscillators (Kuramoto, 2002; Kuramoto & Battogtokh, 2002). First, the system in Equation 44 may show multiple critical point transitions corresponding to multiple linear resonant frequencies ωki as activation level γ increases. Second, sufficiently close to the critical point, the strong modification of an effective wave mode frequency by the nonresonant input from multiple wave modes may result in nonlinear resonances with different modes that are not possible for linear waves, thus providing a mechanism for the emergence of unexpected oscillations difficult to explain by more simplistic models.

We emphasize that the nonlinear regime of those surface wave modes and their possible effects on coupling between different frequencies and between different spatial scales depends on geometry as well as the frequencies. Thus, it is difficult to say what geometries or frequencies may be important in coupling without a proper simulation study, especially when oscillations at multiple scales propagating in different patterns are taking part in the coupling. Moreover, coherent oscillations that may appear as persistent patterns will clearly be more important for modulating firing than chaotic oscillations of the same (and probably much higher) amplitude. The numerical wave simulation construction of our theoretical approach is based on realistic tissue properties and is capable of directly incorporating both large-scale geometries of hippocampus or cortex and small-scale inhomogeneity/anisotropy of local fiber bundles that is well suited for studies of the relative importance of various frequencies or geometrical regions in providing significant neuronal coupling.

We would also like to highlight the fact that our inhomogeneity/anisotropy-related wave mechanism is not dependent on the specific source of the activity and thus is applicable to both excitation and/or propagation of synaptic-related activity as well as to the synaptic-free wave propagation shown experimentally. Nevertheless, we want to emphasize that, in this initial paper, we have avoided any computational modeling of tissue conductivity and permittivity distributions at the cellular or subcellular microscale level, as this is a complex problem well beyond the scope of this initial paper. More importantly, this is not necessary for demonstrating the significant effects generated by the average tissue properties, whatever their source may be.

It is worth noting that there has been a growing interest in recent years in developing a “physics theory of the brain” (for a review, see Lynn & Bassett, 2019) or in using physics-inspired scale-free or multiscale network models (Bassett & Bullmore, 2017; Betzel & Bassett, 2017). However, such theories typically consist of the application of known models from other fields of physics applied only by analogy (e.g., statistical mechanics of brain networks). Very few of these physics models of the brain actually include and model physics processes in the brain. In contrast, we have developed a theory of brain waves based on the physical principles of wave propagation in the known inhomogeneous and anisotropic nature of actual brain tissue, using morphological features derived directly from real imaging data.

### Limitations

We note that the wave mechanism described in this paper does not dismiss or underestimate in any way the excitation and/or propagation of synaptic-related activity (Bao & Wu, 2003; Wu, Guan, Bai, & Yang, 2001; Wu, Guan, & Tsau, 1999; Fukuda, Hata, Ohshima, & Tsumoto, 1998). It is formally possible, although quite computationally expensive, to generate tissue conductivity and permittivity distributions at the cellular or subcellular microscale level. Although it is commonly accepted now that, at the physiological frequency range (0–5 kHz), the inductive, magnetic, and propagative (wave) effects of the bioelectrical signals in the extracellular space can be neglected (Logothetis, Kayser, & Oeltermann, 2007; Aidley, 1998; Robinson, 1968), the presented mechanism can be used to describe the low-amplitude field potential oscillations and their propagation in the vicinity of the inhomogeneous/anisotropic cellular interfaces as well as propagation in the axon oscillations at the subthreshold levels and not only large-amplitude spikes. It is interesting that in vitro activity produced by shocking at different regions, layers, or fiber tracts, using slices with different cutting angles, does not appear to be noticeably different (Wu et al., 2001), suggesting that the activity may be sustained by nonspecific intracortical connections; thus, the wave mechanism presented here can be appropriate for its description. We have also demonstrated that the peaks these wave loops would induce in EEG spectra are consistent with those typically observed EEG data.

Although we found that the amplitude estimates for an averaged scalar potential provide physiologically important maximum amplitude variations of up to 5–10 mV, it is important to recognize the inherent difficulties in the estimation of signal amplitudes of individual wave packets (contrary to predicting an averaged local potential created by an ensemble of waves). The spherical shell wave propagation example (as well as cortical wave simulations) is not able to predict the exact values for the wave amplitude in a single packet as it requires calculation of the balance between excitation and dissipation of these waves at various frequency ranges and the simple estimates using amplitudes of spiking are not that particularly informative as they simply predict the values that are in the range of typically observed field potentials. Nevertheless, even without the amplitude estimation for an individual loop wave packet, our analysis is able to predict the preferred directions of wave propagation and loop pattern formation based on geometrical and tissue properties. Therefore, our model's ability to predict regions of cortex where external wave activity can emerge and form a sustained loop pattern has the potential to be important for understanding where the neuron spiking synchronization will have better chances to be achieved, hence it can provide substantial input in understanding effects on neural information processing and plasticity.

It is also important to recognize that direct experimental results show that, although a large-amplitude single-spike event can be visible only at a very short distance ≈ 10 μm (Bédard, Kröger, & Destexhe, 2006b), nevertheless, even subthreshold-coherent external field oscillations with an amplitude under 0.2 mV can produce phase locking and synchronization of spikes (Anastassiou, Perin, Markram, & Koch, 2011). The difference in amplitude for these events is more than 500-fold; nevertheless, the easily detectable (relatively) high-amplitude spike only influences the very limited vicinity, whereas the very low amplitude but coherent and persistent wave loop signal results in very profound consequences, although its amplitude may be below the threshold of current detecting devices or may be buried in the surrounding noise. We would like to mention that different mechanisms can be employed for the prediction of frequency scaling in the EEG spectrum (Baranauskas et al., 2012; Milstein, Mormann, Fried, & Koch, 2009; Bédard, Kröger, & Destexhe, 2006a), and the purpose of Figure 6 is not to predict this scaling but rather to show that simple estimates using typical tissue conductivity/permittivity parameters and typical values for both small scales of inhomogeneity/anisotropy and large scales of brain geometry result in noncontradictory estimates for both frequency ranges of persistent loop-like patterns that can possibly emerge from these typical tissue/geometry values, as well as noncontradictory estimates for the power spectrum slope outside (not inside) these EEG frequency ranges.

However, we would like to emphasize an important point that, although simple single-exponent power laws are ubiquitous in the EEG literature, this must be recognized as primarily a matter of simplicity and convenience, for it is highly unlikely that such simplistic descriptions adequately characterize the EEG spectrum for either task-based or resting state brain activity experiments. As such, simple power law estimates, although perhaps useful in a broad sense, can be misleading in that the very low dimensional (i.e., one parameter) description of such a complex system makes it impossible to adequately characterize the different sources and their amplitudes. The modus operandi of our paper is to develop such a framework. However, having done so, we still need to demonstrate that our model produces results consistent with these simplistic models, that is, that our results using our model and typical tissue parameters produced “noncontradictory” estimates. However, of course, it would be nonsensical to suggest that the simple single power law is a confirmation of our model, because the much larger parameter space of our model has a huge number of possible configurations that can give essentially the same exponent for the simple power law description. Our model provides characterization of the more complex problems directly related to tissue architecture and physiology, such as the distribution of power among the different bands. As an example, we demonstrate that an excess of power in different bands (from alpha to gamma) at different conditions is what distinguishes the brain activity from randomly statistically forced power decay. Our model provides a clear prediction of power spectrum for this randomly forced waves as well as predicts frequency ranges where this simple single-exponent power decay cannot be used, and these ranges are also in agreement with the observational EEG data.

### Conclusions: Cognitive and Clinical Applications

In conclusion, in this paper, we have presented an inhomogeneous anisotropic physical model of wave propagation in the brain cortex. The model predicts that, in addition to the well-known damped oscillator-like wave activity in brain fibers, there is another class of brain waves that is not directly related to major fibers but instead propagates in directions mediated by the local variations in the inhomogeneity and anisotropy of tissue conductivity and thus can potentially propagate in any direction, including the direction along the fibers or in the direction of the inhomogeneity gradient. However, the dissipation of the waves is smallest when they propagate cross fibers so that, on average, the cross-fiber direction of propagation predominates. Moreover, this dissipation can be significantly less than that of the well-known waves along brain fibers. The end result is the existence of waves perpendicular to the fibers along the highly folded cortical regions in a weakly evanescent manner that results in their persistence on time scales long compared to waves along brain fibers. For the first time, we have obtained the dispersion relation for these surface cortex waves and have shown, both analytically and numerically, a plausible physical argument for their existence. Through numerical analysis, we have developed a procedure to generate an inhomogeneous and anisotropic distribution of conductivity tensors using anatomical and diffusion brain MRI data. Although the detailed numerical studies of effects and importance of these waves and their possible biological role are out of the scope of this article, we have presented preliminary results that suggest that the time of life for these wave-like cortex activity events may significantly exceed the decay time of the typical axon action potential spikes. Thus, they can provide a persistent neuronal response in the cortical areas generated as a result of “along-the-axon” spiky activations, the mechanism important for understanding various aspects of brain dynamics, from propagation of epileptiform activity to memory consolidation.

The presented model is directly related to a very interesting research topic of including effects of realistic human connectivity in studies of brain dynamics, establishing brain spatio-temporal activity patterns using realistic brain connectivity (Robinson et al., 2016). Our approach allows straightforward incorporation and integration of the realistic human brain morphological features and structural connectivity with brain functional dynamics using, presented here, a realistic inhomogeneous/anisotropic tissue framework derived from Maxwell equations. A study based on this mechanism (Galinsky et al., 2018) applied to multimodal MRI/EEG/MEG resting and task-based datasets showed good repeatability and similarity between subjects, with the intraclass correlation coefficient ranging from .9 to higher than .99, thus demonstrating potentially broad implications for both basic neuroscience and clinical studies by enhancing the spatial–temporal resolution of the estimates derived from current neuroimaging modalities.

We emphasize that the nonlinear extension is not only an abstract model of the Maxwell equations in an inhomogeneous medium. It is directly based on conductivity/permittivity values obtained from local tissues and the level of their spatial and angular variations appropriately selected from the relevant neurophysiological brain properties. The complex interaction of the brain geometry with the fields is critical to the wave loops and their coupling with the spiking behavior. The computational environment using actual structural/dMRI data provides a very novel framework for studying this complex problem. As such, it is a practical model that offers the possibility of linking observed brain activity (from EEG, for example) with structural data (from HRA MRI or dMRI) and functional data such as fMRI. This capability would be of great importance in both research and clinical neuroscience applications.

The theory and computational framework presented in this article can be used in a joint estimation diffusion imaging approach (Frank et al., 2020) to provide further diffusion sequence improvements through robust characterization of voxel diffusion microscopic properties inferred from wave processes.

This weakly evanescent cortical wave theory and the associated computational framework thus provide a basis for relating quantitative tissue microstructural properties (such as anisotropy and inhomogeneity) and measurable larger-scale architectural features (e.g., cortical thickness) directly to electrophysiological measurements being performed with increasingly sensitive techniques (such as EEG) within a wide range of important basic and clinical research programs. The ability of this new physical model for the generation, propagation, and maintenance of brain waves may have significant implications for the analysis of electrophysiological brain recording and for current theories about human brain function. Furthermore, the dependence of these waves on brain geometry, such as cortical thickness, has potentially significant implications for understanding brain function in abnormal states, such as Alzheimer's disease, where cortical thickness changes are evident, and the dependence on tissue status may be important in conditions such as traumatic brain injury, where tissue damage may alter its anisotropic and inhomogeneous properties.

## Acknowledgments

L. R. F. and V. L. G. were supported by National Science Foundation grants DBI-1143389, DBI-1147260, EF-0850369, PHY-1201238, ACI-1440412, and ACI-1550405 and National Institutes of Health (NIH) grants R01 MH096100 and R01 AG054049. Data were provided (in part) by the Human Connectome Project (HCP), WU-Minn Consortium (principal investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH institutes and centers that support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University. Data collection and sharing for this project were provided by the HCP (principal investigators: Bruce Rosen, MD, PhD; Arthur W. Toga, PhD; and Van J. Weeden, MD). HCP funding was provided by the National Institute of Dental and Craniofacial Research, the National Institute of Mental Health, and the National Institute of Neurological Disorders and Stroke. HCP data are disseminated by the Laboratory of Neuro Imaging at the University of Southern California. The EEG data are courtesy of Antigona Martinez of the Nathan S. Kline Institute for Psychiatric Research.

Reprint requests should be sent to Vitaly L. Galinsky, Center for Scientific Computation in Imaging, University of California at San Diego, 9500 Gilman Dr., MC 854, La Jolla, CA 92093-0854, or via e-mail: vit@ucsd.edu.

## REFERENCES

REFERENCES
Abrams
,
D. M.
, &
Strogatz
,
S. H.
(
2004
).
Chimera states for coupled oscillators
.
Physical Review Letters
,
93
,
174102
.
Aidley
,
D. J.
(
1998
).
The physiology of excitable cells
(4th ed.).
Cambridge
:
Cambridge University Press
.
Anastassiou
,
C. A.
,
Perin
,
R.
,
Markram
,
H.
, &
Koch
,
C.
(
2011
).
Ephaptic coupling of cortical neurons
.
Nature Neuroscience
,
14
,
217
223
.
Bao
,
W.
, &
Wu
,
J.-Y.
(
2003
).
Propagating wave and irregular dynamics: Spatiotemporal patterns of cholinergic theta oscillations in neocortex in vitro
.
Journal of Neurophysiology
,
90
,
333
341
.
Baranauskas
,
G.
,
Maggiolini
,
E.
,
Vato
,
A.
,
Angotzi
,
G.
,
Bonfanti
,
A.
,
Zambra
,
G.
, et al
(
2012
).
Origins of 1/f2 scaling in the power spectrum of intracortical local field potential
.
Journal of Neurophysiology
,
107
,
984
994
.
Barbour
,
B.
(
2017
).
Analysis of claims that the brain extracellular impedance is high and non-resistive
.
Biophysical Journal
,
113
,
1636
1638
.
Bassett
,
D. S.
, &
Bullmore
,
E. T.
(
2017
).
Small-world brain networks revisited
.
Neuroscientist
,
23
,
499
516
.
Baxter
,
G. T.
, &
Frank
,
L. R.
(
2013
).
A computational model for diffusion weighted imaging of myelinated white matter
.
Neuroimage
,
75
,
204
212
.
Bédard
,
C.
,
Kröger
,
H.
, &
Destexhe
,
A.
(
2006a
).
Does the 1/f frequency scaling of brain signals reflect self-organized critical states?
Physical Review Letters
,
97
,
118102
.
Bédard
,
C.
,
Kröger
,
H.
, &
Destexhe
,
A.
(
2006b
).
Model of low-pass filtering of local field potentials in brain tissue
.
Physical Review E
,
73
,
051911
.
Berry
,
D. B.
,
Regner
,
B.
,
Galinsky
,
V. L.
,
Ward
,
S. R.
, &
Frank
,
L. R.
(
2018
).
Relationships between tissue microstructure and the diffusion tensor in simulated skeletal muscle
.
Magnetic Resonance in Medicine
,
80
,
317
329
.
Betzel
,
R. F.
, &
Bassett
,
D. S.
(
2017
).
Multi-scale brain networks
.
Neuroimage
,
160
,
73
83
.
Buzsáki
,
G.
(
2002
).
Theta oscillations in the hippocampus
.
Neuron
,
33
,
325
340
.
Buzsáki
,
G.
(
2006
).
Rhythms of the brain
.
New York
:
Oxford University Press
.
Buzsáki
,
G.
, &
Draguhn
,
A.
(
2004
).
Neuronal oscillations in cortical networks
.
Science
,
304
,
1926
1929
.
Chiang
,
C.-C.
,
Shivacharan
,
R. S.
,
Wei
,
X.
,
Gonzalez-Reyes
,
L. E.
, &
Durand
,
D. M.
(
2019
).
Slow periodic activity in the longitudinal hippocampal slice can self-propagate non-synaptically by a mechanism consistent with ephaptic coupling
.
Journal of Physiology
,
597
,
249
269
.
Czurkó
,
A.
,
Huxter
,
J.
,
Li
,
Y.
,
Hangya
,
B.
, &
Muller
,
R. U.
(
2011
).
Theta phase classification of interneurons in the hippocampal formation of freely moving rats
.
Journal of Neuroscience
,
31
,
2938
2947
.
Destexhe
,
A.
,
Contreras
,
D.
, &
,
M.
(
1999
).
Spatiotemporal analysis of local field potentials and unit discharges in cat cerebral cortex during natural wake and sleep states
.
Journal of Neuroscience
,
19
,
4595
4608
.
Fiser
,
J.
,
Chiu
,
C.
, &
Weliky
,
M.
(
2004
).
Small modulation of ongoing cortical dynamics by sensory input during natural vision
.
Nature
,
431
,
573
578
.
FitzHugh
,
R.
(
1961
).
Impulses and physiological states in theoretical models of nerve membrane
.
Biophysical Journal
,
1
,
445
466
.
Fonov
,
V.
,
Evans
,
A. C.
,
Botteron
,
K.
,
Almli
,
C. R.
,
McKinstry
,
R. C.
,
Collins
,
D. L.
, et al
(
2011
).
Unbiased average age-appropriate atlases for pediatric studies
.
Neuroimage
,
54
,
313
327
.
Fonov
,
V.
,
Evans
,
A. C.
,
McKinstry
,
R. C.
,
Almli
,
C. R.
, &
Collins
,
D. L.
(
2009
).
Unbiased nonlinear average age-appropriate brain templates from birth to adulthood
.
Neuroimage
,
47(Suppl. 1)
,
S102
.
Fox
,
S. E.
,
Wolfson
,
S.
, &
Ranck
,
J. B.
, Jr.
(
1986
).
Hippocampal theta rhythm and the firing of neurons in walking and urethane anesthetized rats
.
Experimental Brain Research
,
62
,
495
508
.
Frank
,
L. R.
,
Zahneisen
,
B.
, &
Galinsky
,
V. L.
(
2020
).
JEDI: Joint estimation diffusion imaging of macroscopic and microscopic tissue properties
.
Magnetic Resonance in Medicine
,
84
,
966
990
.
Frank
,
T. D.
,
Daffertshofer
,
A.
,
Peper
,
C. E.
,
Beek
,
P. J.
, &
Haken
,
H.
(
2000
).
Towards a comprehensive theory of brain activity: Coupled oscillator systems under external forces
.
Physica D: Nonlinear Phenomena
,
144
,
62
86
.
Fries
,
P.
(
2015
).
Rhythms for cognition: Communication through coherence
.
Neuron
,
88
,
220
235
.
Fukuda
,
M.
,
Hata
,
Y.
,
Ohshima
,
M.
, &
Tsumoto
,
T.
(
1998
).
Role of NMDA receptors in the propagation of excitation in rat visual cortex as studied by optical imaging
.
Neuroscience Research
,
31
,
9
21
.
Gabriel
,
S.
,
Lau
,
R. W.
, &
Gabriel
,
C.
(
1996a
).
The dielectric properties of biological tissues: II. Measurements in the frequency range 10 Hz to 20 GHz
.
Physics in Medicine and Biology
,
41
,
2251
2269
.
Gabriel
,
S.
,
Lau
,
R. W.
, &
Gabriel
,
C.
(
1996b
).
The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues
.
Physics in Medicine and Biology
,
41
,
2271
2293
.
Galinsky
,
V. L.
, &
Frank
,
L. R.
(
2014
).
Automated segmentation and shape characterization of volumetric data
.
Neuroimage
,
92
,
156
168
.
Galinsky
,
V. L.
, &
Frank
,
L. R.
(
2015
).
Simultaneous multi-scale diffusion estimation and tractography guided by entropy spectrum pathways
.
IEEE Transactions on Medical Imaging
,
34
,
1177
1193
.
Galinsky
,
V. L.
, &
Frank
,
L. R.
(
2016
).
The lamellar structure of the brain fiber pathways
.
Neural Computation
,
28
,
2533
2556
.
Galinsky
,
V. L.
, &
Frank
,
L. R.
(
2019
).
Symplectomorphic registration with phase space regularization by entropy spectrum pathways
.
Magnetic Resonance in Medicine
,
81
,
1335
1352
.
Galinsky
,
V. L.
,
Martinez
,
A.
,
Paulus
,
M. P.
, &
Frank
,
L. R.
(
2018
).
Joint estimation of effective brain wave activation modes using EEG/MEG sensor arrays and multimodal MRI volumes
.
Neural Computation
,
30
,
1725
1749
.
Gerstner
,
W.
,
Kistler
,
W. M.
,
Naud
,
R.
, &
Paninski
,
L.
(
2014
).
Neuronal dynamics: From single neurons to networks and models of cognition
.
New York
:
Cambridge University Press
.
Goel
,
P.
, &
Ermentrout
,
B.
(
2002
).
Synchrony, stability, and firing patterns in pulse-coupled oscillators
.
Physica D: Nonlinear Phenomena
,
163
,
191
216
.
Gomes
,
J.-M.
,
Bédard
,
C.
,
Valtcheva
,
S.
,
Nelson
,
M.
,
Khokhlova
,
V.
,
Pouget
,
P.
, et al
(
2016
).
Intracellular impedance measurements reveal non-ohmic properties of the extracellular medium around neurons
.
Biophysical Journal
,
110
,
234
246
.
Hallez
,
H.
,
Vanrumste
,
B.
,
Van Hese
,
P.
,
D'Asseler
,
Y.
,
Lemahieu
,
I.
, &
Van de Walle
,
R.
(
2005
).
A finite difference method with reciprocity used to incorporate anisotropy in electroencephalogram dipole source localization
.
Physics in Medicine & Biology
,
50
,
3787
3806
.
Haueisen
,
J.
,
Tuch
,
D. S.
,
Ramon
,
C.
,
Schimpf
,
P. H.
,
Wedeen
,
V. J.
,
George
,
J. S.
, et al
(
2002
).
The influence of brain tissue anisotropy on human EEG and MEG
.
Neuroimage
,
15
,
159
166
.
Hodgkin
,
A. L.
, &
Huxley
,
A. F.
(
1952
).
A quantitative description of membrane current and its application to conduction and excitation in nerve
.
Journal of Physiology
,
117
,
500
544
.
Ingber
,
L.
, &
Nunez
,
P. L.
(
2011
).
Neocortical dynamics at multiple scales: EEG standing waves, statistical mechanics, and physical analogs
.
Mathematical Biosciences
,
229
,
160
173
.
Javitt
,
D. C.
(
2009
).
When doors of perception close: Bottom–up models of disrupted cognition in schizophrenia
.
Annual Review of Clinical Psychology
,
5
,
249
275
.
Javitt
,
D. C.
, &
Sweet
,
R. A.
(
2015
).
Auditory dysfunction in schizophrenia: Integrating clinical and basic features
.
Nature Reviews Neuroscience
,
16
,
535
550
.
Kartashov
,
I. N.
, &
Kuzelev
,
M. V.
(
2014
).
Dissipative surface waves in plasma
.
Plasma Physics Reports
,
40
,
650
664
.
Kuramoto
,
Y.
(
2002
).
Reduction methods applied to nonlocally coupled oscillator systems
. In
J.
Hogan
,
A.
Champneys
,
B.
Krauskopf
,
M.
,
E.
Wilson
,
H.
Osinga
, &
M.
Homer
(Eds.),
Nonlinear dynamics and chaos: Where do we go from here?
(pp.
209
227
).
Boca Raton, FL
:
CRC Press
.
Kuramoto
,
Y.
, &
Battogtokh
,
D.
(
2002
).
Coexistence of coherence and incoherence in nonlocally coupled phase oscillators
.
Nonlinear Phenomena in Complex Systems
,
5
,
380
385
.
Li
,
X. P.
,
Xia
,
Q.
,
Qu
,
D.
,
Wu
,
T. C.
,
Yang
,
D. G.
,
Hao
,
W. D.
, et al
(
2015
).
The dynamic dielectric at a brain functional site and an EM wave approach to functional brain imaging
.
Scientific Reports
,
4
,
6893
.
Logothetis
,
N. K.
,
Kayser
,
C.
, &
Oeltermann
,
A.
(
2007
).
In vivo measurement of cortical impedance spectrum in monkeys: Implications for signal propagation
.
Neuron
,
55
,
809
823
.
Lubenov
,
E. V.
, &
Siapas
,
A. G.
(
2009
).
Hippocampal theta oscillations are travelling waves
.
Nature
,
459
,
534
539
.
Lynn
,
C. W.
, &
Bassett
,
D. S.
(
2019
).
The physics of brain network structure, function and control
.
Nature Reviews Physics
,
1
,
318
332
.
Martínez
,
A.
,
Gaspar
,
P. A.
,
Hillyard
,
S. A.
,
Bickel
,
S.
,
Lakatos
,
P.
,
Dias
,
E. C.
, et al
(
2015
).
Neural oscillatory deficits in schizophrenia predict behavioral and neurocognitive impairments
.
Frontiers in Human Neuroscience
,
9
,
371
.
Milstein
,
J.
,
Mormann
,
F.
,
Fried
,
I.
, &
Koch
,
C.
(
2009
).
Neuronal shot noise and Brownian 1/f2 behavior in the local field potential
.
PLoS One
,
4
,
e4338
.
Monai
,
H.
,
Inoue
,
M.
,
Miyakawa
,
H.
, &
Aonishi
,
T.
(
2012
).
Low-frequency dielectric dispersion of brain tissue due to electrically long neurites
.
Physical Review E
,
86
,
061911
.
Morris
,
C.
, &
Lecar
,
H.
(
1981
).
Voltage oscillations in the barnacle giant muscle fiber
.
Biophysical Journal
,
35
,
193
213
.
Muller
,
L.
,
Chavane
,
F.
,
Reynolds
,
J.
, &
Sejnowski
,
T. J.
(
2018
).
Cortical travelling waves: Mechanisms and computational principles
.
Nature Reviews Neuroscience
,
19
,
255
268
.
Muller
,
L.
,
Piantoni
,
G.
,
Koller
,
D.
,
Cash
,
S. S.
,
Halgren
,
E.
, &
Sejnowski
,
T. J.
(
2016
).
Rotating waves during human sleep spindles organize global patterns of activity that repeat precisely through the night
.
eLife
,
5
,
e17267
.
Nagumo
,
J.
,
Arimoto
,
S.
, &
Yoshizawa
,
S.
(
1962
).
An active pulse transmission line simulating nerve axon
.
Proceedings of the IRE
,
50
,
2061
2070
.
Nazarenko
,
S.
(
2011
).
Wave turbulence
(
Lecture Notes in Physics 825
.
Berlin, Germany
:
Springer
.
Niethammer
,
M.
,
Estepar
,
R. S. J.
,
Bouix
,
S.
,
Shenton
,
M.
, &
Westin
,
C.-F.
(
2006
).
On diffusion tensor estimation
. In
2006 International Conference of the IEEE Engineering in Medicine and Biology Society
(pp.
2622
2625
).
New York
:
IEEE
.
Nunez
,
P. L.
, &
Srinivasan
,
R.
(
2014
).
Neocortical dynamics due to axon propagation delays in cortico-cortical fibers: EEG traveling and standing waves with implications for top–down influences on local networks and white matter disease
.
Brain Research
,
1542
,
138
166
.
Qiu
,
C.
,
Shivacharan
,
R. S.
,
Zhang
,
M.
, &
Durand
,
D. M.
(
2015
).
Can neural activity propagate by endogenous electrical field?
Journal of Neuroscience
,
35
,
15800
15811
.
Rayleigh
,
L.
(
1885
).
On waves propagated along the plane surface of an elastic solid
.
Proceedings of the London Mathematical Society
,
17
,
4
11
.
Rinzel
,
J.
(
1987
).
A formal classification of bursting mechanisms in excitable systems
. In
E.
Teramoto
&
M.
Yamaguti
(Eds.),
Population biology, morphogenesis and neurosciences
(pp.
267
281
).
Berlin, Germany
:
Springer
.
Robinson
,
D. A.
(
1968
).
The electrical properties of metal microelectrodes
.
Proceedings of the IEEE
,
56
,
1065
1071
.
Robinson
,
P. A.
,
Zhao
,
X.
,
Aquino
,
K. M.
,
Griffiths
,
J. D.
,
Sarkar
,
S.
, &
Mehta-Pandejee
,
G.
(
2016
).
Eigenmodes of brain activity: Neural field theory predictions and comparison with experiment
.
Neuroimage
,
142
,
79
98
.
Sheremet
,
A.
,
Qin
,
Y.
,
Kennedy
,
J. P.
,
Zhou
,
Y.
, &
Maurer
,
A. P.
(
2019
).
Wave turbulence and energy cascade in the hippocampus
.
Frontiers in Systems Neuroscience
,
12
,
62
.
Shivacharan
,
R. S.
,
Chiang
,
C.-C.
,
Zhang
,
M.
,
Gonzalez-Reyes
,
L. E.
, &
Durand
,
D. M.
(
2019
).
Self-propagating, non-synaptic epileptiform activity recruits neurons by endogenous electric fields
.
Experimental Neurology
,
317
,
119
128
.
Smith
,
M. A.
, &
Kohn
,
A.
(
2008
).
Spatial and temporal scales of neuronal correlation in primary visual cortex
.
Journal of Neuroscience
,
28
,
12591
12603
.
Squire
,
L. R.
,
Genzel
,
L.
,
Wixted
,
J. T.
, &
Morris
,
R. G.
(
2015
).
Memory consolidation
.
Cold Spring Harbor Perspectives in Biology
,
7
,
a021766
.
Stewart
,
M.
,
Quirk
,
G. J.
,
Barry
,
M.
, &
Fox
,
S. E.
(
1992
).
Firing relations of medial entorhinal neurons to the hippocampal theta rhythm in urethane anesthetized and walking rats
.
Experimental Brain Research
,
90
,
21
28
.
Tuch
,
D. S.
,
Wedeen
,
V. J.
,
Dale
,
A. M.
,
George
,
J. S.
, &
Belliveau
,
J. W.
(
1999
).
Conductivity mapping of biological tissue using diffusion MRI
.
Annals of the New York Academy of Sciences
,
888
,
314
316
.
Van Essen
,
D. C.
,
Smith
,
S. M.
,
Barch
,
D. M.
,
Behrens
,
T. E. J.
,
Yacoub
,
E.
,
Ugurbil
,
K.
, et al
(
2013
).
The WU-Minn Human Connectome Project: An overview
.
Neuroimage
,
80
,
62
79
.
Weiss
,
S. A.
, &
Faber
,
D. S.
(
2010
).
Field effects in the CNS play functional roles
.
Frontiers in Neural Circuits
,
4
,
15
.
Wilting
,
J.
, &
Priesemann
,
V.
(
2018
).
Inferring collective dynamical states from widely unobserved systems
.
Nature Communications
,
9
,
2325
.
Woldorff
,
M. G.
,
Liotti
,
M.
,
Seabolt
,
M.
,
Busse
,
L.
,
Lancaster
,
J. L.
, &
Fox
,
P. T.
(
2002
).
The temporal dynamics of the effects in occipital cortex of visual–spatial selective attention
.
Cognitive Brain Research
,
15
,
1
15
.
Wu
,
J.-Y.
,
Guan
,
L.
,
Bai
,
L.
, &
Yang
,
Q.
(
2001
).
Spatiotemporal properties of an evoked population activity in rat sensory cortical slices
.
Journal of Neurophysiology
,
86
,
2461
2474
.
Wu
,
J.-Y.
,
Guan
,
L.
, &
Tsau
,
Y.
(
1999
).
Propagating activation during oscillations and evoked responses in neocortical slices
.
Journal of Neuroscience
,
19
,
5005
5015
.
Zakharov
,
V. E.
,
L'vov
,
V. S.
, &
Falkovich
,
G.
(
1992
).
Kolmogorov spectra of turbulence I: Wave turbulence
.
Berlin, Germany
:
Springer
.
Zalesky
,
A.
,
Fornito
,
A.
,
Cocchi
,
L.
,
Gollo
,
L. L.
, &
Breakspear
,
M.
(
2014
).
Time-resolved resting-state brain networks
.
Proceedings of the National Academy of Sciences, U.S.A.
,
111
,
10341
10346
.
Zhang
,
H.
,
Watrous
,
A. J.
,
Patel
,
A.
, &
Jacobs
,
J.
(
2018
).
Theta and alpha oscillations are traveling waves in the human neocortex
.
Neuron
,
98
,
1269
1281
.
Zhang
,
M.
,
,
T. P.
,
Qiu
,
C.
,
Shivacharan
,
R. S.
,
Gonzalez-Reyes
,
L. E.
, &
Durand
,
D. M.
(
2014
).
Propagation of epileptiform activity can be independent of synaptic transmission, gap junctions, or diffusion and is consistent with electrical field transmission
.
Journal of Neuroscience
,
34
,
1409
1419
.