## 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

**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

*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:

_{GM}= 4.07 · 10

^{7}ε

_{0}, ε

_{WM}= 2.76 · 10

^{7}ε

_{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.

*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}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).

_{ij}(

**r**) scaled by an inhomogeneous density ρ(

**r**), such that the total conductivity tensor Σ

_{ij}(

**r**) is defined by

_{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-mm

^{3}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

^{(1)}, σ

^{(2)}, and σ

^{(3)}, fixed (location-independent) anisotropic tensors were assigned to every location in the brain.

^{(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

_{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 $\lambda d1$. The anisotropy tensor σ

_{ij}(

**r**) was defined as

**R**is a rotation matrix between directions (0, 0, 1) and

**d**

^{(1)}≡ (

*x*,

*y*,

*z*) that can be expressed as

#### 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 $\lambda d1$ was used to define the anisotropy tensor major axis. The position-dependent anisotropy value was defined by introducing parameter ε as either $\lambda d2$/$\lambda d1$ or ($\lambda d2$ + $\lambda d3$)/2/$\lambda d1$ (in both cases resulting in axisymmetric form of the conductivity tensor) and also using two different values $\lambda d2$/$\lambda d1$ and $\lambda d3$/$\lambda 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.

*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

**R**

_{k}orientation matrix for every fiber

*k*, that is,

### Wave Trajectory Integration

**r**

_{0}with an initial wave vector

**k**

_{0}, we substitute the inhomogeneous anisotropic conductivity tensor Σ in the wave dispersion relation (imaginary part of Equation 3)

*t*= 0,

**r**=

**r**

_{0}, and

**k**=

**k**

_{0}to trace the wave trajectory

*t*(τ),

**r**(τ), and

**k**(τ), respectively.

**k**into magnitude

*k*and direction $k\u02dc$ parts (

**k**=

*k*$k\u02dc$, and |$k\u02dc$|

^{2}= 1) and rewrite the last two equations as

^{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

*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

_{⊥}< Σ

_{∥}) will reduce this estimate even further,

_{zz}≡

*S*(

*r*), the wave frequency ω can be expressed as

*S*$kz2$ because

*kh*≪ 1. Then, from the geometrical optics ray equations

*d*ω/

*dr*< 0. For

*k*

_{z}=

*k*/$2$ and

*z*= $\u2212\omega r/d\omega /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

*W*is as a path integral along its trajectory, parameterized by

*t*(τ),

**r**(τ), and

**k**(τ):

_{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

_{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

_{d}= Σ

_{xx}and Ω = ∂

_{y}Σ

_{xy}.

_{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),

_{k}and wave numbers

*k*(where “c.c.” denotes complex conjugate), which results in a set of coupled equations for time-dependent complex amplitudes

*a*

_{k}(

*t*) ≡

*a*(

*k*,

*t*)

#### Resonant Coupling of Wave Modes

The nonlinear terms 𝒩_{k} will include a sum of inputs from multiple waves, that is, *k* = $\Sigma in$*k*_{i} where *n* is the order of the nonlinearity. Those resonant conditions will give rise to coupling terms that include various combinations of exp(*i*(ω_{k} − $\Sigma in$ ± ω_{ki})*t*). Additional requirements for frequency resonances (ω_{k} = $\Sigma 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 (ω_{k} ∼ *k*^{m}, 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.

*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

*a*

_{−k}= $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}− ω

_{k′}.

_{k}= ±ω

_{k′}± ω

_{k″}allow the expression of the nonlinear resonant coupling $\mathcal{N}kR$ by extraction of only the relevant terms as

*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,

*k*

_{1}=

*k*(−1 − $5$)/2, and

*k*

_{2}=

*k*(3 + $5$)/2 are the real solutions of quadratic equations 1/

*k*± 1/

*k*′ ± 1/|

*k*−

*k*′| = 0.

#### Nonresonant Coupling of Wave Modes

*k*→ ∞, an effect that is absent for coupling of waves with directly proportional dispersion. To illustrate this, we will estimate the nonresonant nonlinear input $\mathcal{N}k0nR$ to the

*k*

_{0}wave mode.

_{0}(

*k*) = ω

_{k0}− ω

_{k}.

*k*=

*k*

_{0}(the forcing term γ

_{e}/

*k*

^{2}in Equation 27 is largest when

*k*=

*k*

_{0}). Therefore, we will derive the nonresonant input term $\mathcal{N}k0nR$ only for

*k*=

*k*

_{0}, thus neglecting nonlinear and damping terms for any

*k*>

*k*

_{0}(more correctly, for any

*k*that is not in resonance with

*k*

_{0}or for

*k*>

*k*

_{2}=

*k*

_{0}(3 + $5$/2) ≈ 2.618

*k*

_{0}). At the limit

*k*→ ∞, all frequency deltas δω

_{1–4}(

*k*) → ω

_{k0}≡ ω

_{0}and

*k*–

*k*

_{0}≈

*k*, hence approximately, we can estimate the nonresonant term $\mathcal{N}k0nR$ as

*k*solution (

*a*

_{k−k0}≈

*a*

_{k}for

*k*≫

*k*

_{0}) and keeping only terms that include

*a*

_{k0}(assuming that the amplitude

*a*

_{k0}is small and can be considered constant relative to any of the δω(

*k*) terms), we can approximately write that

*C*

_{1…4}(

*k*) and

*C*(

*k*) = Σ

*C*

_{i}(

*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

*k*

_{0}are left. Therefore, the nonresonant input term $\mathcal{N}k0nR$ in Equation 34 depends on

*a*

_{k0}and

*t*as

*C*(

*k*)

*C**(

*k*) = $C\u02dc$.

*e*

^{−iωk0t}asymptotic behavior for

*t*→ ∞.

*a*

_{k0}that integrates the nonlinear nonresonant input from smaller spatial scales can be written as

*C*

_{0}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

*A*and phase

*B*(

*a*

_{k0}=

*Ae*

^{iB}) as

_{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

*a*

_{k}for

*k*=

*k*

_{0}…

*k*

_{N}

*a*

_{kn}= 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 Σ_{ij}*k*_{i}, 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 |Σ_{ij}*k*_{i}| 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.

_{ij}will simply be

*S*δ

_{ij}, where

*S*is a scalar inhomogeneous conductivity. Therefore, both the phase velocity

*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 γ_{eff}*L*_{loop}/υ_{gr} ≤ 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 *L*_{loop} 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.

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.

### 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, 𝓕 = ∑_{i}*A*_{i}δ(*t* − *t*_{i})δ(**r** − **r**_{i}), corresponding to a flat forcing frequency spectrum in the Fourier domain. Then, from the last term in Equation 2 |(∂_{i}Σ_{ij})(∂_{j}φ)| ∼ ω$\phi \omega 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* = $\u2212\omega r/d\omega /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).

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.

### 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.

_{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*→ ∞,

*B*

_{0}is some arbitrary constant phase. Therefore, the nonoscillatory state requires that γ satisfies to

_{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 *I*^{nR} integrals in Equation 39 instead of a single *e*^{−iωk0t} exponent input. Figure 8 shows simulation results for several parameter sets with the same values as were used for plots of Figure 7 (ω_{k0} = β = *k*_{0} = 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).

#### 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.

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

*x*,

*t*) of the waves with the dispersion relation in Equation 4, one can write the two-point spatial

*C*

_{x}(ζ,

*t*

_{1},

*t*

_{2}) and temporal

*C*

_{t}(

*x*

_{1},

*x*

_{2},

*τ*) correlation functions as

*k*

_{1}and

*k*

_{2}select an appropriate spectral range and |

*a*

_{k}|

^{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.

## 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

*f*

^{2}scaling in the power spectrum of intracortical local field potential

*f*frequency scaling of brain signals reflect self-organized critical states?

*f*

^{2}behavior in the local field potential

DOI:https://doi.org/10.1103/PhysRevLett.93.174102,PMID:15525081