Universal theory of brain waves: From linear loops to nonlinear synchronized spiking and collective brain rhythms

An inhomogeneous anisotropic physical model of the brain cortex is presented that predicts the emergence of nonevanescent (weakly damped) wavelike 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 standard paradigm that the neural fibers determine the pathways for signal propagation in the brain, the signal propagation due to the cortex wave modes in highly folded areas will exhibit no apparent correlation with the fiber directions. The results are consistent with numerous recent experimental animal and human brain studies demonstrating the existence of electrostatic field activity in the form of traveling waves (including studies where neuronal connections were severed) and with wave-loop-induced peaks observed in EEG spectra. In addition, we demonstrate that the resonant and nonresonant terms of the nonlinear coupling between multiple modes produce both synchronous spiking-like high-frequency wave activity as well as low-frequency wave rhythms as a result of their unique dispersion properties. Numerical simulation of forced multiple mode dynamics shows that as forcing increases there is a transition from a damped to an oscillatory regime that subsequently decays away as over-excitation 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 cortical wave modes, and this mechanism of their collective input to the generation of synchronized spiking and low-frequency rhythms, have significant implications in particular for neuroimaging methods that detect electromagnetic physiological activity, such as EEG and MEG, and in general for the understanding of brain activity, including mechanisms of memory.


I. INTRODUCTION
The majority of 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 spatiotemporal scales. Experimentally, data on which this assumption is employed range from high temporal resolu-tion neural oscillations detected at low spatial resolution by EEG/MEG [1] to high spatial resolution functional MRI resting state modes oscillation detected at low temporal resolution [2]. 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 [3]) 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 [4], i.e., by allowing input from different scales or introducing axon propagation delays [5].
However, recent detection [6] of cortical wave activity spatiotemporally organized into circular wavelike 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. Additionally, many studies show evidence that electrostatic field activity in animal or human hippocampus (as well as cortex) are traveling waves [7] that can affect neuronal activity by modulating the firing rates and may possibly play an important functional role in diverse brain structures [8]. And perhaps more importantly it has been experimentally shown [9] that periodic activity can self-propagate by endogenous electric fields even 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 due to the interplay between tissue inhomogeneity and anisotropy in the thin surface cortex layer. This mechanism has been overlooked by previous models of brain wave characterization and thus is absent from current network pathway reconstruction and analysis approaches.
The abundance of oscillatory patterns across a wide range of spatial and temporal scales of brain electromagnetic activity has generated much discussion and debate in the literature on their interaction [10]. The standard approach involves representing the brain as a large network of coupled oscillators [11] and using this as a test bed for the study of network wave propagation, mechanisms of synchrony, possibly deriving some mean field equations and properties, etc. However, such models are necessarily descriptive, and their relationship to actual physical properties of actual brain tissue or the electromagnetic waves they support is tenuous.
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 wavelike 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 "alongthe-axon" spiky activations. This mechanism may provide an alternative approach for the integration of microscopic brain properties and for the development of a "physical" model for memory.
The 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, timescales, and directions and supports extended and highly structured regions of existence in dissipative media using simple one-and two-dimensional anisotropic models, as well as more realistic full brain model 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., a thin inhomogeneous layer around a sphere with a 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 scales similar to the scales of global rotation recently detected experimentally [6].
An important aspect of the cortical waves model we present 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. While 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.
In addition, we demonstrate that 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 that have been observed experimentally throughout a wide range of scales of brain activity. 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. This provides a physical mechanism to explain the flexibility and adaptability of the human brain across a wide range of tasks.
In order to describe this complex behavior we show for the first time that the inverse proportionality of frequency and wave number in 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 a nonlinear equation that integrates and includes the collective nonresonant input from multiple wave modes. This general mechanism explains the emergence of synchronized spiking from an elegant perspective of nonlinear wave physics, rather than relying on empirical fitting of a single measured quantity to a set of ad hoc multiparametric differential equations using 20 or more dimensional and dimensionless fitting parameters, as is typically employed by a multitude of single-neuron spiking models [10,12]. Following the ideas of wave turbulence [13] 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 Sec. II we present the linear theory that explains the existence of evanescent wave loops. In Sec. III we extend this to the nonlinear regime and demonstrate the emergence of spiking behavior and its interaction with the cortical loops of the linear theory.

A. Wave dispersion
We start with the most general form of brain electromagnetic activity using Maxwell equations in a general (i.e., inhomogeneous and anisotropic) medium: Using the electrostatic potential E = −∇φ, Ohm's law J = σ · E (where σ ≡ {σ i j } is an anisotropic conductivity tensor), a linear electrostatic property for brain tissue D = εE, assuming that the scalar permittivity ε is a "good" function (i.e., it does not go to zero or infinity everywhere), and taking the change of variables ∂x → ε∂x , the charge continuity equation for the spatial-temporal evolution of the potential φ can be written in terms of a permittivity scaled conductivity where we have included a possible external source (or forcing) term F. For brain fiber tissues the conductivity tensor might have significantly larger values along the fiber direction than across them. The charge continuity without forcing, i.e., (F = 0) can be written in tensor notation as where repeating indices denote summation. Simple linear wave analysis, i.e., substitution where k is the wave number, r is the coordinate, is the frequency, and t is the time, gives the following complex dispersion relation: which is composed of the real and imaginary components: The condition for non-(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 γ : i.e., that 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 Ref. [14]) ε GM = 4.07 × 10 7 ε 0 , ε WM = 2.76 × 10 7 ε 0 , σ GM = 2.75 × 10 −2 S/m, σ 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.

B. Effects of anisotropy
In order to better understand the effects of brain tissue micro-and macrostructure on the manifestation of propagating brain waves, it is instructive to consider two idealized tissue models. In the first model [ Fig. 1(a)] all brain fibers The uniform area of fibers oriented along z direction (shown in green) is bounded by a thin transitional area (magenta) where the conductivity gradient may be important (a sketch of one possible conductivity profile is shown at the bottom of the panel). (b) Schematic picture that can be used as a crude two-dimensional approximation of a fold. The direction of fiber conductivity has only x and z components, and all quantities are assumed to be uniform in the y direction.
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. We assume that small cross-fiber currents can be characterized by a small parameter and introduce the conductivity tensor as 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 Fig. 1(a)] and the conductivity gradient is directed along x axis. Then we will look for a solution for the potential φ located in the thin layer of inhomogeneity (that is, we substitute x → x) that depends 023061-3 on z and y only as φ = φ (z) + φ ⊥ (y): where 0 and 1 denote the zeroth and the first orders of power. Equation (6) describes a potential along the fiber direction and is a damped oscillator equation that has a decaying solution, but Eq. (7) describes a potential perpendicular to the fiber direction and does not include a damping term, hence it describes a pure wavelike solution that propagates in the thin layer transverse to the main fiber direction. Thus although this wavelike solution φ ⊥ has a smaller amplitude than along the fiber action potential φ , but it can nevertheless have a much longer lifetime.
We would like to stress that Eqs. (6) and (7) are here only for illustrative purposes to emphasize a relatively obvious but often overlooked consideration that under anisotropic inhomogeneous conditions some directions may happen to be better suited for wave propagation than the others, and, in particular, those directions are not necessarily along the direction of the fibers. Below we will develop the more complete theory of brain wave propagation using the more general dispersion relation (3).
In order to account for geometric variations, we construct a slightly more complex two-dimensional model that can be viewed as a very crude approximation of a cortical fold [ Fig. 1(b)]. It is easy to derive an analogue of Eq. (7) for the φ ⊥ (y) in the case when υ ≡ υ(x, z). This equation for the φ ⊥ (y) will again be of a wave type, similar to Eq. (7), with the addition of a z component of the conductivity gradient and with a similar wavelike 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 wavelike 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 [15]. The existence of surface waves at the dissipative medium boundary is also known [16].

C. Effects of brain composition and architecture
In order to extend the above analysis to a more realistic model of brain tissue architecture we extracted volumetric structural brain parameters from high-resolution anatomical MRI data sets as well as brain fiber anisotropy from diffusion weighted MRI data sets. All anatomical and diffusion MRI data sets are from the Human Connectome Project [3]. The details for all of the processing steps can be found in Refs. [17][18][19]. 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.
To provide an illustration of where the conductivity anisotropy and inhomogeneity can form appropriate The isosurface that is located in the outer cortex area of gray matter. The color scheme uses shades of green to mark regions where dissipative term dominates, i.e., | i j k i | |∂ i i j |, shades of red where | i j k i | < |∂ i i j | 2| i j k i |, and shades of blue where the wavelike term is more than two times dominant, i.e., 2| i j k i | < |∂ i i j |. The inner cortex shown in (a) clearly display prevalence of dissipation, whereas the outer cortex shown in (b) allows for wavelike cortex activity in a majority of the locations.
conditions for cortex surface waves generation we created plots that can be used to characterize the ratio |γ /ω|. To construct the ratio we calculated two vectors, ∂ i i j and i j k i , and compared their norms. For wave vector k we used a vector with the same direction as in the ∂ i i j vector and magnitude |k| = |∇ |/| |; i.e., our intention is to compare norms of dissipative and wavelike terms at kh ≈ 1 where h is the cortical thickness. Figure 2 shows plot of the ratio of | i j k i | and |∂ i i j | at two different depths inside the brain with a dissipative regime in the inner cortex [ Fig. 2(a)] and wavelike conditions in the outer cortex [ Fig. 2(b)].
Both anisotropy and inhomogeneity are important for the existence of the cortex surface waves. For example, for inhomogeneous but isotropic tissue the conductivity tensor i j will simply be Sδ i j , where S is a scalar inhomogeneous conductivity. Therefore both the phase velocity 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 (due to tensor products in ∇ · and ∇ · · k).

D. Spherical cortex shell model
To provide some (possibly overly optimistic) estimates based on a typical human brain dimensions we can consider a simple spherical cortex shell model with a cortical layer of fixed thickness h ≈ 1.5-3 mm 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 provides from Eq. (4) for the longest waves (with the smallest amount of damping) with Anisotropy ( ⊥ < ) will reduce this estimate even further, 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 (3) with only a single component of the conductivity tensor zz ≡ S(r), the wave frequency ω can be expressed as where we have neglected the term Sk 2 z because kh 1. Then from the geometrical optics ray equations where r l and k l are l component of the coordinate r and wave number k vectors, respectively, we can get 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 k z = 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 (more details of spherical shell model are available in Appendix A).

E. Brain waves power spectra
In order to determine a possible energy distribution across different frequencies for these cortical waves some knowledge about the forcing term F in Eq. (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, corresponding to a flat forcing frequency spectrum in the Fourier domain. Then, from the last term in Eq.
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 Eq. (8) 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 1.2-92 Hz (v ph ≈ 0.002-7 m/s). The large-scale circular cortical waves of Ref. [6] (9-18 Hz, 2-5 m/s) are clearly inside this range. The above spherical shell wave propagation example as well as wave simulations presented below are not able to predict the exact values for wave amplitude 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 our analysis is able to predict the preferred directions of wave propagation and loop pattern formation based on geometrical and tissue properties. Moreover, this framework provides the mechanisms to incorporate a more complete quantitative description of axonal spiking, which we address in the next section.
Evidence for the existence of these wave loop induced spectral peaks is shown in Fig. 3, which shows the spectral power of the EEG signal for six subjects [20], 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.

F. Linear brain wave propagation
The effect described theoretically above can be demonstrated through numerical simulation of wave propagation in a thin dissipative inhomogeneous and anisotropic cortical layer. As a starting point for numerical study we include Fig. 4, which shows spatial snapshots of the dynamical behavior of randomly generated wave trajectories (the movies are in [21] with additional technical details in Appendix G) in a regime with γ eff L loop /v gr 1 (where γ eff is an effective wave dissipation rate, i.e., a difference between average dissipation and spiky activations rates for a wave packet propagating with group velocity v gr along some characteristic loop of L loop length). Wave packets are simulated using the ray equations (9) where the general form of anisotropic dispersion relation (3) and (4) is used.
In the idealized spherical model some of the emergent persistent localized cortical wave patterns are precisely the simple loop pattern predicted by Eq. (10), despite the complex initial spatiotemporal 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 Fig. 2 (middle and bottom panels) also produce stable loop structures, now embedded within the complex geometry of the cortical folds.
All wave simulations were initialized with wave packets of random parameters (frequency, wave number, location, etc.), but clearly show an emergence of localized persistent closed loop patterns at different spatial and temporal scales, from scales as large as the whole brain where rotational wave activity has been experimentally detected in Ref. [6] to as small as the resolution used for the cortical layer thickness detection (more details of cortical fold model and wave integration are available in Appendixes B-D). . All wave packets were initialized with random parameters assuming the presence of spiky activation sources (not shown). Movie files for these and additional loop examples can be found in [21].

A. Reduced equation
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 two-dimensional symmetric form of the conductivity tensor with constant diagonal terms xx and yy (where yy is along the fibers 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 be interested in only a one-dimensional solution for the potential φ(x) located in this thin layer of inhomogeneity that can be described by the reduced equation where γ d = xx and = ∂ y xy .

023061-6
The source term F can be assumed to have a frequencyindependent 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 [7,22], The solution φ can be sought as a Fourier integral expansion, 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 a k (t ) ≡ a(k, t ), where from Eq. (3) the wave mode frequencies are inversely proportional to the wave number and where all spectral and spatial integrals are assumed to be taken on the infinite domain [−∞, ∞], and the spectral integrals are actually converted to the half infinite domain [0, ∞] as a consequence of the symmetry of a k due to φ(x, t ) being real (from this point we will omit repeating explicit limits of integration where it does not create ambiguity). The condition of wave propagation (ω/γ > 1), corresponding to small wave growth or damping during a single wave period, provides natural limits for the spectral (k 0 , k ∞ ) or spatial (l, L) domains to be taken as cutoffs in a case of presence of large-or small-scale divergences, i.e., The existence of spatial scales L ≡ 2π /γ e and l ≡ 2πγ d / does not mean that the waves propagate only inside some fixed [l, L] domain that requires additional boundary conditions. Both the spectral (k 0 , k ∞ ) and the spatial (l, L) scales simply emphasize the fact that for any frequency it is the combination of the excitation γ e and the damping γ d that defines an extent where any disturbance of the scalar potential φ may be assumed to have wave properties.

B. Resonant coupling of wave modes
The nonlinear terms N k will include a sum of inputs from multiple waves, i.e., k = n i 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 − n i ±ω k i )t ). Additional requirements for frequency resonances (ω k = n i ±ω k i ) produce wave turbulence-like [13] selection rules for the coupling terms that are similar to phase coupling terms in a ring of connected oscillators [23].
For waves having typical dispersion properties, that is, with the frequencies directly proportional to the wave numbers (ω k ∼ k α α > 0), the maximum oscillatory frequency is increasing and going to infinity with the increase of wave numbers. In this case the nonlinear terms produce a direct cascade of wave energy [13] constantly generating larger and larger frequencies. 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 the field potential and firing rate [7,22]. The feedback can be represented through nonlinear coupling. This will be demonstrated using the simplest quadratic form N (φ) = φ(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 Using symmetry conditions a −k = a * k and ω −k = −ω k [a consequence of φ * (x) = φ(x)] this can be rewritten as The corresponding conditions for the frequency resonances ω k = ±ω k ± ω k allow the expression of the nonlinear resonant coupling N R k by extraction of only the relevant terms as where only three out of four wave-number resonances appear, as the resonance k + k + k = 0 is not possible [13],

C. 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 (15). 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 023061-7 dispersion. To illustrate this, we will estimate the nonresonant nonlinear input N nR k 0 to the k 0 wave mode: where The approximate expression for the forced oscillations solution can be obtained assuming that all the forcing input originates from the scales of k = k 0 [the forcing term γ e /k 2 in Eq. (14) is largest when k = k 0 ]. Therefore we will derive the nonresonant input term N nR k 0 only for k = k 0 , thus neglecting nonlinear and damping terms for any k > k 0 [more correctly, for any k that are not in resonance with k 0 or for k > k 2 = k 0 (3 + √ 5/2) ≈ 2.618k 0 ]. At the limit k → ∞ all frequency deltas δω 1−4 (k) → ω k 0 ≡ ω 0 and k − k 0 ≈ k, hence we can estimate approximately the nonresonant term N nR k 0 as To estimate forced oscillations terms required in evaluation of the integral Eq. (23), one can write from Eqs. (14), (20), and (22) that Looking again for an approximate large k solution (a k−k 0 ≈ a k for k k 0 ) and keeping only terms that include a k 0 (assuming that the amplitude a k 0 is small and can be considered constant relative to any of the δω(k) terms), we can approximately write that where C 1,...,4 (k) and C(k) = C i (k) are some complex integration constants that we assume to have random phases with the amplitude independent of k, and, hence, we can use that C(k)C * (k) =C. Therefore, the nonresonant input term N nR k 0 Eq. (23) depends on a k 0 and t as where terms with C(k)C(k) and its complex conjugate vanish because of the randomness of the phases. More accurate estimation of N nR k 0 will require evaluation of integrals similar to that, although resulting in more complex expressions, nevertheless have the same e −iω k 0 t asymptotic behavior for t → ∞ (the details of the derivation of I nR in general form without using the k k 0 condition and an evaluation of the asymptotic limit t → ∞ are presented in Appendixes E and F).
Because the exact numeric coefficients that appear in Eq. (27) [or the exact form of algebraic time dependence in Eqs. (F8) to (F11) for the leading asymptotic terms] may depend on some arbitrary initial background wave spectra, in the following analysis we will introduce some free parameters that can be used to study and characterize possible effects of various linear and nonlinear (e.g., resonant and nonresonant) terms.
Therefore, an equation for the longest wave length brain mode a k 0 that integrates the nonlinear nonresonant input from smaller spatial scales can be written as 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 Eq. (29) (where C 0 is a constant) has the same 1/k 2 0 asymptotic behavior for t → ∞ as the solution of Eq. (11) obtained with timeand space-scale-independent forcing.
Equation (29) can be converted to a system of equations for the amplitude A and phase B (a k 0 = Ae iB ) as where δ A and δ B were added to introduce tunable phase delays [24].

D. Solution for nonresonant coupling
The system of Eqs. (31) and (32) includes nonresonant input from all wave modes, but the resonant term Eq. (21) should also be added together with additional equations for resonant wave amplitudes that participate in resonant coupling. We will consider this complete system later, but first we will investigate the behavior of the nonlinear nonresonant part only. Figure 5 shows the results of numerical solution of the system Eqs. (31) and (32) for several different sets of parameters. The time evolution of the highest frequency, longest wave length 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 [10]. Increase in the level of activation γ produces a nonlinear signal with the spikelike shape of a single-neuron firing.

E. Critical point and synchronized spiking
It is interesting that this spiking-like solution of system Eqs. (31) and (32) appears near the critical point, and the oscillatory state undergoes bifurcation and transitions to a nonoscillatory regime as γ reaches the value above some critical point. To illustrate the reason for this transition we will consider the simplest case of = 1, k 0 = 1, δ A = 0 and δ B = π/2 (although different parameter values can be used for a similar analysis as well). The nonoscillatory regime can be reached if dA/dt → 0 and dB/dt → −ω k 0 as t → ∞. Then from Eqs. (31) and (32) one can write that at t → ∞, where B 0 is some arbitrary constant phase. Therefore, the nonoscillatory state requires that γ satisfies to Hence for the nonoscillatory solution is not possible. The simulations shown in the bottom panels of Fig. 5 confirm that the critical γ value is indeed 3 when ω k 0 = α = β = 1. Similar analysis when δ A = δ B gives the critical γ value equal to [2 − cos (π/3)]/ sin (π/3) ≈ 1.732. 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 [12], and finishing with many dynamical integrate-and-fire models of neuron [10], are based on an approximation of several local neuron variables, e.g., membrane currents, gate voltages, etc., and defining the relations between these local properties. Contrary to this and rather unexpectedly, Eq. (29) is obtained through an integration of a large number of oscillatory brain wave modes nonresonantly interacting in an 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 that Eq. (29) cannot be separated into "fast" and "slow" parts as typically required for 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 singleneuron models [25]. It describes a mechanism for 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 [26] of brain activity.

F. Bursting activity
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 Eq. (28) instead of a single e −iω k 0 t exponent input. Figure 6 shows simulation results for several parameter sets with the same values as were used for plots of Fig. 5 (ω k 0 = β = k 0 = 1, δ A = 0). The numerical solution shows more complex behavior that now includes modulation of the spiking rate with lower frequency and emergence of a burstlike train of spikes, effects often observed in different types of neuronal activity [10].  31) and (32) when the exponential term was replaced by I nR integrals Eq. (28) with the region of integration set to 50k 0 < k < 1000k 0 . For all plots the values of ω k 0 , k 0 , α, and β were set to be equal to 1, δ A = 0, and γ and δ B were varied. The top and bottom rows show plots for phase delay δ B equal to 3π/4 and π/2, respectively. The left column displays modulation of spiking rate for γ = 4.5. The right column shows the nonlinear bursting of spikes for γ = 5.1 (time and amplitude units are arbitrary).

G. Emergence of low-frequency brain rhythms
Finally, we considered a model that combines a chain of inputs from nonlinear modes generated due to resonant terms Eq. (21) into a set of nonresonant mode equations (29), which results in a system of equations for mode amplitudes a k for k = k 0 · · · k N : · · · da k n dt = a k n k 2 n γ + β 2 |k n | e −iω kn t+δ a * k n − 2α a k n a k n + λ k 2 n a k n−2 a * k n−1 + a k n−1 a k n+1 + a * k n+1 a k n+2 , where the parameter λ describes the strength of resonant coupling between modes.
We would like to mention two new, rather important, and not entirely obvious features that appear in the nonlinear system (35) but are absent for phase-coupled oscillator models (e.g., Ref. [23]). First, the system (35) may show multiple critical point transitions corresponding to multiple linear resonant frequencies ω k i 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 emergence of unexpected oscillations difficult to explain by more simplistic models. Figures 7 and 8 show results of numerical simulation of the system (35), clearly indicating that weak nonlinear resonant coupling between just three modes with frequencies of ω k 0 , The results of numerical integration of the system (35). For all plots the values of ω k 0 , k 0 , α, and β were set to be equal to 1, and the resonant coupling λ was 0.05. Different values of δ were again used in real and imaginary parts [as in Eqs. (31) and (32)] with δ A = 3π/4, 0, 0, δ B = 3π/4, 3π/4, π/2 and close to the critical values of γ = 1.731, 2.575, 2, 9969 for the left, middle, and right columns, respectively. The total potential φ is plotted with the black and different colors showing the oscillations of the individual modes. All plots show that when γ is sufficiently close to criticality a weak coupling produces jumps from subcritical to supercritical regimes with amazingly regular low-frequency quasiperiodicity (time and amplitude units are arbitrary).
an emergence of periodic activity with frequencies up to 100-1000 times lower then the linear frequencies of participating modes. We would like to emphasize again that the system (35) cannot be separated into traditional "slow" and "fast" subsystems, hence the low-frequency component cannot be explained by a modulation [27] of the "fast" subsystem with oscillations of the "slow" part.
In Fig. 7 the high-frequency spiking is generated with the level of activation γ = 1.535. This activation level is still 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 rows clearly show that small increase of the resonant coupling (0.001, 0.01, and 0.05, respectively) results in the appearance of component with significantly lower frequency. Figure 8 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 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, e.g., distances from the critical points and other resonances, phase delays, etc., is an interesting open question that will be considered in future work.

IV. DISCUSSION AND CONCLUSIONS
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 are not directly related to major fibers, but instead propagate perpendicular to the fibers along the highly folded cortical regions in a weakly evanescent manner that results in their persistence on timescales long compared to waves along brain fibers. The waves 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 perpendicular to both the inhomogeneity gradient and the direction of anisotropy, therefore, on average the cross-fiber direction of propagation should be seen more often. We have obtained the dispersion relation for those surface cortex waves and have shown, both analytically and numerically, a plausible 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. While the detailed numerical studies of effects and importance of these waves and their possible biological role are beyond the scope of this paper, we have presented preliminary results that suggest that the time of life for these wavelike cortex activity events may significantly exceed the decay time of the typical axon action potential spikes. Thus, they can provide an persistent neuronal response in the cortical areas generated as a result of "along-the-axon" spiky activations.
The ranges of parameters for the waves produced by our model are in agreement with those presented by several studies [7] that show 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 [7] and may possibly play an important functional role in diverse brain structures [8]. The natural selforganization of these traveling waves into looplike structures that our model produces agrees well with recently detected [6] cortical wave activity spatiotemporally organized into circular wavelike 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 have been destroyed [9] can potentially be attributed to these waves as well. We have also demonstrated that the peaks these wave loops would induce in EEG spectra are consistent with those typically observed EEG data.
Direct experimental results have shown that even despite the small amplitudes of the external field potentials relative to the threshold of a spiking neuron, external fields can play a substantial role in the spiking activity, and "even very small and slowly changing fields that triggered V e changes under 0.2 mV led to phase locking of spikes to the external field and to a greatly enhanced spike-field synchrony" [28]. 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.
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 due to the 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. This theory thus provides 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 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.

APPENDIX A: SPHERICAL SHELL CORTEX MODEL
We first provide details about the procedures used for generating inhomogeneous and anisotropic components of the permittivity-scaled conductivity tensor .
A spherical shell cortex model is represented by fixed anisotropy tensor σ i j (r) ≡ σ i j scaled by a radial inhomogeneous density ρ(r) ≡ ρ(r) (see Fig. 9), such that the total conductivity tensor i j (r) is defined by where σ i j is a constant, anisotropic, positive semidefinite symmetric tensor and ρ(r) is piecewise continuous function of radius r (0 r 1): The spherical shell cortex wave simulations include three different anisotropic tensor σ i j choices, where (for ε = 0) σ (1) represents the conductivity tensor used in analytical solution of Eqs. (10), i.e., currents only in the z direction, σ (2) represents different orientation, where currents are allowed in the direction of 45 • relative to all axis, and σ (3) represents more complicated current anisotropy, with crossing currents flowing in the x and z directions, but with no currents in the y direction. Parameters n, r 0 , and r 1 used to control the thickness of the inhomogeneous "cortical" layer (r 0 = 0.5, r 1 = 0.9, α ∞ = 500, and n = 0, 2, 14, and 25 used for Figs. 10 to 13).

APPENDIX B: CORTICAL FOLD MODEL
For a cortical fold model the conductivity tensor i j (r) is again defined as a product of inhomogeneous ρ(r) and anisotropic σ i j (r) parts Eq. (A1), but both parts are now functions of brain locations.

Inhomogeneity estimation
The inhomogeneous density function ρ(r) is estimated from high-resolution anatomical MRI data by processing the data with SWD [17] (skull stripping, field of view normalization, noise filtering) and then registering to MNI152 space [29,30] with SYM-REG [19]. The final 1 mm 3 182 × 218 × 182 volume is used as the inhomogeneous density ρ(r).

Anisotropy estimation
For the anisotropy tensor σ i j (r) several different test cases were employed.

a. Fixed anisotropy orientation and value
The same σ (1) , σ (2) , and σ (3) fixed (location-independent) anisotropic tensors as in the spherical shell cortex model Eq. (A2) were assigned to every location in the brain. This is the simplest case that may allow to separate the effects of inhomogeneity and anisotropy on loop pattern formation.

b. Varying anisotropy orientation and fixed anisotropy value
Multiple diffusion direction and diffusion gradient strength MRI data were used to estimate diffusion tensor D i j [31].
The anisotropy orientation was estimated using eigenvector d (1) of the diffusion tensor D i j with the largest eigenvalue λ (1) d . The anisotropy tensor σ i j (r) was defined as 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

c. Varying anisotropy orientation and value
Similarly to the previous case the diffusion tensor D i j was estimated for each voxel, and eigenvector d (1) with the largest eigenvalue λ (1) d was used to define the anisotropy tensor major axis. The position-dependent anisotropy value was defined by introducing parameter ε as either λ (2) d /λ (1) d or (λ (2) d + λ (3) d )/2/λ (1) d (in both cases resulting in an axisymmetric form of the conductivity tensor) and using two different values λ (2) d /λ (1) d and λ (3) d /λ (1) d for σ 22 and σ 11 .

d. 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 diffusion MRI estimates. Recent measurements of effective impedance that involve both intercellular and extracellular electrodes [32] 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 the extracellular medium may be questionable [33], this clearly confirms high anisotropy for the effective conductivity that should be used for analysis of effects of intercellular and membrane sources (and these are the most important type of sources for 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 anisotropic form of σ (1) i j 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 the z direction corresponds to the direction of the fiber). Using full brain tractography results [34] 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, i.e., with the complete form of the conductivity tensor again given by Eq. (A1).

APPENDIX C: WAVE TRAJECTORY INTEGRATION
To obtain the trajectory of a brain wave with frequency ω emitted at a point r 0 with an initial wave vector k 0 we substitute the inhomogeneous anisotropic conductivity tensor in the wave dispersion relation [imaginary part of Eq. (3)] and obtain wave ray equations as We integrate these equations starting at τ = 0, t = 0, r = r 0 , and k = k 0 to trace the wave trajectory t (τ ), r(τ ), and k(τ ). For numerical integration it is beneficial to split k into magnitude k and directionk parts (k = kk, and |k| 2 = 1), and 023061-13 rewrite the last two equations as where in the last equation the right-hand side is orthogonal tok to guarantee that |k| 2 = 1, and the algebraic expression Eq. (C1) was substituted for the differential equation for wave vector magnitude k.

APPENDIX D: WAVE DISSIPATION AND EXCITATION
An integral along the wave trajectory t (τ ), r(τ ), and k(τ ), with any appropriate model form of wave excitation γ exc describes a change of wave energy W along its path. This may allow the study of many interesting questions of brain wave dynamics, i.e., 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, to find when and/or how these waves may potentially become important in triggering neuronal firing, that is, to study possible mechanisms of synchronization and feedback, etc.

APPENDIX E: DERIVATION OF I nR INTEGRAL
We assume that the initial amplitude values for all wave modes are set at some random level (k), and we can write for a k (0) where for the initial phase ψ (k) we assume that it is a continuously differentiable function up to at least the second derivative with dψ (k)/dk → 0 at both k → 0 and k → ∞ and with the derivatives bounded for any k (including at k → ∞ and k → 0), i.e., |dψ (k)/dk| < M and |d 2 ψ (k)/dk 2 | < M, where M is a constant. For the initial amplitude (k) will assume a power-law-type dependence (k) ∼ k ν . The linear solution a l k (t ) of Eq. (14) then will be a l As the fastest growing linear solution for any k ∞ k k 0 corresponds to k = k 0 , we will derive a nonlinear equation for a k 0 evaluating N nR k 0 Eq. (23) and including only input from this a k 0 mode in the evaluation of nonlinear terms for the rest of a k modes, hence for k = k 0 , and for k > k 0 where in Eq. (E4) we don't need to keep linear terms as either they result in a decaying solution or their linear growth is exponentially slower than the linear growth due to forcing by the a k 0 term (a l k (t )/a l k 0 (t ) ∼ (k)/ (k 0 ) exp(−γ e /k 2 0 t + γ e /k 2 t ) < 1 when k > k 0 and t > ln [ (k 0 )/ (k)]k 2 0 k 2 /γ e (k 2 − k 2 0 )), therefore we ignored the linear terms obtaining Eq. (24). We also don't need to keep them in nonlinear parts as, due to our propagating waves requirements Eqs. (17) and (18), they result only in the appearance of small imaginary components of frequencies ω k .
Therefore we can write an approximate forced solution for Eq. (E4) (without using the k k 0 condition this time) as Now we can substitute Eq. (E5) into a k a k±k 0 and a k a * k±k 0 products of Eq. (E3) to obtain the nonlinear equation for a k 0 . Substituting it only into a single part of the products and using Eq. (E1) for the second part will result in a set of integrals that may be divergent at k → 0 or k → ∞ (therefore the k = k 0 and/or k = k ∞ cutoffs should be used) and with any possible resonances in δω j (k) ignored (i.e., taken as a principal value integrals). The obvious results of those integrals is just a small nonlinear modification of ω k 0 frequency (a nonlinear frequency shift). If we substitute Eq. (E5) into both parts of a k a , * k±k 0 products we will obtain expressions similar to Eq. (28) integrals where we now also included initial amplitudes (k) and phases ψ (k) where a Taylor expansion was used to replace ±[ψ (k) − ψ (k ± k 0 )] with (k) ≡ k 0 dψ/dk and the variable of integration for the second expression has been changed to k/k 0 . The integral limits are not shown (assumed to be from 0 to ∞), but as the integrals may be divergent at k = 0 and/or k = ∞, the cutoffs at k = k 0 and/or k = k ∞ may be applied if needed.
All the above expressions show that the oscillatory mode e −iω k 0 t emerges as result of integration of multimode nonresonant terms, but the conditions on the initial amplitudes and phases of the modes determine if this oscillation mode slowly (algebraically not exponentially) grows or decays at large times.
The wave-number spectral power |a k | 2 ∼ k μ can be related to the frequency spectral power |a ω | 2 ∼ ω τ as therefore, τ = −μ − 2. Of course, this equality gives the exact relation between frequency and wave-number spectra only if just the wave modes with ω k = /k dispersion contribute to all temporal and spatial oscillation modes, and this is not necessarily true as other propagating and nonpropagating modes can be present. Nevertheless, it outlines the general trend that the spectra with a steeper frequency dependence (with smaller τ or larger μ ≡ 2ν, that is, with a greater "energy" of lower frequency modes relatively to higher frequency modes) will generate asymptotically growing high-frequency e −iω k 0 t oscillations, whereas more shallow frequency spectra (with large τ or smaller μ ≡ 2ν, that is with increased presence of high-frequency modes) will show asymptotic decay of the high-frequency e −iω k 0 t oscillations. Thus the asymptotic expressions show the existence of nonlinear processes that provide a mechanism for spectral regulation.

APPENDIX G: PERSISTENT WAVE LOOP PATTERNS
One of the interesting questions is the possibility of persistent pattern formation/self-organization from a randomly emitted distribution of brain waves influenced by an inhomogeneous anisotropic structure of their dispersion. Considering a simplest case of fixed difference between wave excitation and dissipation, i.e., assuming propagation of wave packets of fixed width (exponentially decaying) and searching for any closed parts of wave trajectories (loops) that are shorter than the packet width can provide a clue about possible answer. Figure 4 as well as Figs. 10 to 25 (and the movies found in [21]) show formation of persistent loops for a variety of wave packet initial conditions as well as for different models of conductivity tensor inhomogeneity and anisotropy constructed from dMRI or whole-brain tractography.
Examples of wave trajectories and emergent persistent loop patterns for the spherical shell cortex model with varying amounts of tensor anisotropy and inhomogeneous shell layer thickness (Figs. 10 to 13).
Examples of wave trajectories and emergent persistent loop patterns for cortical fold geometry with different approaches used for estimation of inhomogeneity and anisotropy are shown in Figs. 14 to 25. Among those examples are several simple cases with variable inhomogeneity and fixed anisotropy (similar to the above spherical shell cortex model) as well as with more complex estimates of anisotropy based on multiple diffusion gradients MRI (dMRI) acquisitions. Assuming a linear dependence between diffusion and conductivity tensors [35] several combinations of the diffusion tensor eigenvectors were used to describe the conductivity tensor anisotropy. In a more complex approach, the full brain tractography [18] generated fiber distributions [34] were used to infer anisotropy through direct integration of single-fiber anisotropy parameters in each voxel.  = (142, 142, 142). High-resolution and low-resolution movie links can be found online in [21].
FIG. 14. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with diffusion MRI-derived position-dependent anisotropy tensor R T σ (1) The trajectory was initialized with wave vector k = (−0.75, 0.23, −0.62) inside the inhomogeneous layer with voxel coordinates r=(55, 130, 130). High-resolution and low-resolution movie links can be found online in [21].

023061-18
FIG. 17. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with fixed anisotropy tensor σ (1) (ε = 0.1). The trajectory was initialized with wave vector k = (0.56, 0.58, −0.59) inside the inhomogeneous layer with voxel coordinates r = (72, 156, 120). Highresolution and low-resolution movie links can be found online in [21].
FIG. 18. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with diffusion MRI-derived position-independent anisotropy tensor R T σ (1) R (ε = 0.1). The trajectory was initialized with wave vector k = (−0.13, −0.97, −0.21) inside the inhomogeneous layer with voxel coordinates r = (111, 103, 110). High-resolution and low-resolution movie links can be found online in [21].
FIG. 19. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with fixed anisotropy tensor σ (1) (ε = 0.01). The trajectory was initialized with wave vector k = (0.42, −0.41, 0.81) inside the inhomogeneous layer with voxel coordinates r = (114, 85, 22). High-resolution and low-resolution movie links can be found online in [21].
FIG. 20. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with fixed anisotropy tensor σ (2) (ε = 0.1). The trajectory was initialized with wave vector k = (−0.90, 0.24, 0.36) inside the inhomogeneous layer with voxel coordinates r = (66, 86, 126). High-resolution and low-resolution movie links can be found online in [21].
FIG. 21. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with diffusion MRI-derived position-independent anisotropy tensor R T σ (1) R (ε = 0.01). The trajectory was initialized with wave vector k = (−0.90, 0.24, 0.36) inside the inhomogeneous layer with voxel coordinates r = (66, 87, 126). High-resolution and low-resolution movie links can be found online in [21].
FIG. 22. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with tractography-derived anisotropy tensor 1/N k R T k σ (1) R k (ε = 0.001). The trajectory was initialized with wave vector k = (0.4, −0.83, 0.4) inside the inhomogeneous layer with voxel coordinates r = (75, 61, 90). High-resolution and low-resolution movie links can be found online in [21].

023061-20
FIG. 23. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with tractography-derived anisotropy tensor 1/N k R T k σ (1) R k (ε = 0.01). The trajectory was initialized with wave vector k = (0.66, 0.5, 0.56) inside the inhomogeneous layer with voxel coordinates r = (33, 77, 56). High-resolution and low-resolution movie links can be found online in [21].
FIG. 24. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with tractography-derived anisotropy tensor 1/N k R T k σ (1) R k (ε = 0.1). The trajectory was initialized with wave vector k = (0.36, 0.67, −0.65) inside the inhomogeneous layer with voxel coordinates r = (77, 123, 64). High-resolution and low-resolution movie links can be found online in [21].
FIG. 25. An example of randomly initialized complete wave trajectory (a) and emergent loop pattern (b) for the cortical fold model with inhomogeneity extracted from HRA volume registered to MNI152 space and with tractography-derived anisotropy tensor 1/N k R T k σ (1) R k (ε = 0). The trajectory was initialized with wave vector k = (−0.06, 0.38, 0.92) inside the inhomogeneous layer with voxel coordinates r = (132, 113, 96). High-resolution and low-resolution movie links can be found online in [21].