High-order theta harmonics account for the detection slow gamma

Local field potential (LFP) oscillations are the superposition of excitatory/inhibitory postsynaptic potentials. In the hippocampus, the 20-55 Hz range (‘slow gamma’) is proposed to support cognition independent of other frequencies. However, this band overlaps with theta harmonics. We aimed to dissociate the generators of slow gamma versus theta harmonics with current source density and different LFP decompositions. Hippocampal theta harmonic and low gamma generators were not dissociable. Moreover, comparison of wavelet, ensemble empirical-mode (EEMD), and Fourier decompositions produced distinct outcomes with wavelet and EEMD failing to resolve high-order theta harmonics well defined by Fourier analysis. The varying sizes of the time-frequency atoms used by wavelet distributed the higher-order harmonics over a broader range giving the impression of a wide band, low frequency burst. The absence of detectable slow gamma refutes a multiplexed model of cognition in favor of the energy cascade hypothesis in which dependency across oscillatory frequencies exists.


14
The effort of understanding the connection between brain activity and behavior (started as early 15 as [Jung and Kornmüller, 1938, Green and Arduini, 1954, Vanderwolf, 1969, O'keefe and Nadel, 16 1978, Buzsáki, 2005]) may be described, in essence, as collecting a variety of measurements of 17 brain activity and associating (often statistically) data "features" to classes of behavior. For ex-18 ample, local field potential (LFP) measurements in the hippocampus exhibit well-known, unique, 19 and easily identifiable oscillatory patterns that associate with obvious behavior. The theta rhythm 20 (4-12 Hz) and associated gamma oscillations (approx. 40-100 Hz) [Stumpf, 1965, Buzsáki et [Buzsaki,28 1986], and is believed to be responsible for the reactivation of memories during off-line periods 29 [Wilson and McNaughton, 1994]. 30 While this associative approach seems simple enough, it is in fact very difficult to implement in 31 practice, due to both the complexity of the system observed (the brain) and difficulty of defining 32 behavior. The brain controls a staggering range of biological functions, in a manner that seems to 33 be at once distributed and tightly integrated. This makes it hard to identify the precise brain region 34 and type of activity "in charge" of a given function (assuming the function is well defined), and, 35 if identified, to select a clear measurement "feature" to associate to behavior. Moreover, defining 36 behavior beyond broad characterizations such as "awake", "active", or "sleep", is itself a challenge, 37 as one can seldom assign clear causal relations to sequences of behavior actions. 38 As data collection and analysis techniques refine, the need of robust methods becomes more 39 stringent, both for identifying data "features", and behavior markers. A number of recent studies 40 report that rat running-speed is monotonically correlated to the power and frequency of both theta probability distribution of rat speed and LFP variance (figure 2). Importantly, variance was chosen 87 over standard post-filtering amplitude method as it is insensitive to a specific frequency range of 88 interest, but took advantage of the power-law shape of the spectra, where lower frequencies (e.g., 89 theta) would carry the majority of the variance. Two speed domains are obvious in figure 2. At 90 low rat speeds (~0-3 cm/s, region A) variance and speed are statistically independent: speed is 91 not a marker for behavior. This region corresponds to low activity, quiescent behavior, or sleep 92 (e.g., REM or non-REM). At speeds larger than, say, 5 cm/s (region B), variance shows a clear 93 correlation to speed. Remarkably, the relation between variance and speed appears to follow the

99
FIGURE 2. Example of joint probability distribution of the variance of hippocampal LFP and velocity (logarithmic scale) for rat r695 ♀ , LM layer. Each gray dot represents one 1-s segment of LFP recording (one realization). Two region may be identified: in region A velocity and variance are approximately statistically independent; in region B, they are correlated. The red line is a linear regression V = a log 10 v +b, where v is the velocity, V is the variance, and a = 0.2, and b = −0.15. The transition region (blue) is excluded from this classification.
Below, we use these two regions to classify behavior and discuss results of different analysis 100 techniques. While this classification is obviously crude and broad, with possibly multiple behavior FIGURE 3. Distribution of cross-spectrum over hippocampal layers during active behavior (region B in figure 2), at low speed. Cross-spectra shown are estimated with reference to the DG layer. For each rat, the observed regions of the hippocampus Top: spectral density of LFP variance. Middle: coherence (cross-spectrum modulus). Bottom: phase lag (cross-spectrum argument). On each sub-panel, the vertical axis is channel number (not shown); the panels are aligned according to the layer identification procedure (figure 1).  . Cross-spectra shown are estimated with reference to the DG layer. For each rat, the observed regions of the hippocampus Top: spectral density of LFP variance. Middle: coherence (cross-spectrum modulus). Bottom: phase lag (crossspectrum argument). On each sub-panel, the vertical axis is channel number (not shown); the panels are aligned according to the layer identification procedure (figure 1).

RESULTS
cross-spectrum matrix (appendix ) over spectral densities, coherence, and phase lag (details of the 128 methodology for estimating cross-spectra are given in appendix, section B). The lamina identifica-129 tion map discussed in section 2.1 (figure 1) provides a convenient framework for the examination 130 of LFP power and phase distribution across hippocampal layers and across rat population.

131
At low speed, the distribution of power across layers is dominated by the LM and DG layers 132 (figure 3). In general, the spectra are rather featureless, exhibiting a broad peak at approximately 133 8 Hz and,in some cases, extending into lower frequencies (figure 3, rat 539). A weak second 134 harmonic of theta is visible in the LM layer at 16 Hz. The coherence of LM with DG is high 135 throughout the entire frequency range. Across layers (with respect to DG) the coherence is high 136 only in theta band (around 8 Hz; the phase lag is not meaningful unless the coherence is high).  [Harper, 1971, Coenen, 1975, Leung et al., 1982, Leung, 1982. Overall, spectral evolution is re-141 markably consistent across the examined rat population revealing that the lower frequencies of the 142 hippocampal spectrum are dominated by theta harmonics that develop with velocity. By definition, 143 harmonics have an integer relationship to the fundamental and are phase-coupled. Therefore, un-144 surprisingly, the phase shift in the harmonics with depth mirror resemble that of the fundamental 145 with the caveat that a substantial amount of the effect is carried by the amplitude/detectability of the 146 harmonic in the decomposition. For example, the pyramidal layer exhibits a low coherence with 147 the DG layer in higher-order harmonics at 24 Hz. It is also worth noting that there is little to no co-148 herence higher frequencies in the dentate gyrus with other lamina of the hippocampus, suggesting 149 that traditional gamma (40-100Hz) can be described as largely "incoherent" across lamina.

150
A closer examination of the spectra estimates for the Pyr, Rad, LM, and DG layers (figure 5, 151 152 panel B) allows for some quantification of the evolution of spectral characteristics with speed. All 153 layers exhibit a significant growth of theta amplitude as speed increases, also associated with a 154 narrowing of the theta peak and a slight shift in peak frequency (roughly from 7 Hz to 9 Hz; the 155 resolution of the spectra is 1 Hz; see appendix B). Theta growth also coincides with the growth 156 of additional spectral peaks at harmonic frequencies of theta, most evident in LM and DG spectra, 157 but also, albeit weaker, in the Pyr and Rad layers. In the frequency range occupied by theta and harmonics (up to ≈ 40 Hz), all spectra shown (figure 5) exhibit a clear power tendency to align to a 158 power law shape f −α , but slopes varying with layer. Approximate slopes for the Pyr, Rad and LM 159 are similar (α Pyr ≈ 1.9, α Rad ≈ 2.1, α LM ≈ 2.1); the DG slope stands out as significantly smaller, 160 α DG ≈ 1.1. The LFP variance in the gamma range (say, between 40 and 100 Hz) also increases 161 with speed, especially in the LM and DG layers. While the gamma range appears to deviate only 162 slightly from the theta-range slope in the Pyr and Rad layers, and only at large speeds, the LM and 163 DG spectra exhibit significantly steeper gamma slopes.  with the Fourier spectral estimates (figure 5), but are much smoother (featureless). In particular, 181 the deteriorating frequency resolution at small scales completely removes higher theta harmonics 182 evident in the Fourier spectra (e.g., third harmonic and higher). The situation is more severe in the 183 DG layer, where only the theta peak is discernible.

184
The examples of scalograms shown in figure 6, middle row, exhibit peaks at theta and its second    The two sharp-wave-ripple complexes shown in figure 7 cover a frequency range 4 to 150 Hz.

199
In DG, the event synchronous to the sharp-wave/ripples is a 10-ms dentate spike [Bragin et al.,  The concept of a stochastic process assumes that such time series are in fact random "realiza-213 tions" of a "virtual" process. When this assumption is made, we expect that each experiment will 214 return a different time series, but as a realization of a unique virtual process, all time series should 215 share some average features that define the individuality of the stochastic process.   ing that is often encountered: while they are both equally applicable, the results are not inter-246 changeable. We stress that the wavelet ordering parameter is scale (the value of the dilation param-247 eter) and not frequency. While scales are often translated to frequencies for ease of interpretation, 248 one should not be confused: Wavelet "frequencies" are, strictly speaking, just labels. Wavelet 249 have a non-trivial frequency distribution of power, i.e., "comprise" an entire range of (true) Fourier 250 frequencies (appendix, subsection B.1.2). Referring to them as if they had the Fourier meaning is 251 incorrect, and may lead to erroneous conclusions. 252 3.3.2. Stochastic processes. The dissimilarity of the two transforms comes into focus when con-253 sidering observations as realizations of a stochastic process. In estimating any quantity that is 254 derived by applying an averaging operator, one implicitly assumes that 1) the data are realizations 255 of a stochastic process; 2) that the quantities derived by averaging exist for that stochastic process. 256 In practice, as the stochastic process is a virtual object, not entirely accessible, it is not possible 257 to guarantee the existence of a particular average quantity. Therefore one has to assume a certain  (1) The existence of the spectral density of the process is guarantied by the celebrated Wiener-Khinchin 273 theorem [Priestley, 1981], which states that the power spectral density of process x is the Fourier  therefore, it appears that it makes little sense to look for a stochastic scope of the wavelet transform,

326
The present study aims to clarify the usefulness and limitations of Fourier and wavelet transform 327 in the analysis of hippocampal LFP recordings. The Fourier transform has been widely adopted as . Left: LM layer, active exploration at high speed motion; arrows mark the areas in the bicoherence map corresponding to theta/harmonics, and theta/gamma phase-coupling. Right: Rad layer, non-REM sleep; arrows mark the areas in the bicoherence map corresponding to the phase-coupling between Fourier components that assemble together into the solitary wave, and those that correspond to sharp wave and sharp-wave/ripples phase coupling (see also discussion in text).
there exists a stochastic process model that guaranties that averaged quantities exist, and provides  This study suggests perhaps, a more general conclusion. Referring back to what we call the 381 "associative approach" for studying the relationship between brain activity and cognitive behavior, 382 it should be quite clear that, past a certain boundary, associations become a rather sterile game.

383
Powerful data analysis techniques, used or misused, will generate ever finer, multi-dimensional 384 features and details. As we increase our details in the the levels of analysis, steps need to be taken 385 to carefully interpret and describe the results in the context of a quantitative brain models. We hope 386 to take up this issue in a subsequent paper. (2) The "coefficients" G(f ) of the decomposition (also referred to as the transform of g) are obtained 464 by taking the inner product of the function g with the basis elements (i.e., projecting g(t) onto the 465 basis), The pair of equations 2 are best known as direct and inverse transforms, sometimes also 466 called analysis and synthesis.

467
The Fourier transform pair is obtained letting Unless the set of functions is carefully defined, equations 3 are only formal, and in most cases  The wavelet transform and its inverse are given by where the functions ψ s,τ (t) = αφ ( t−τ s ) are "copies" of the "mother wavelet" φ(t),  [Strang, 1986, Briggs 485 and Henson, 1995, Mallat, 1998].

494
The Parseval relation insures that equations 5 conserve the variance of the time sequence and its 495 transform [e.g., Briggs and Henson, 1995], i.e., σ g = σ G , where σ g = ∑ n |g n | 2 is the variance of 496 g. The discrete wavelet transform does conserve variance (σ g ̸ = σ G ) and in general, the variance 497 ratio depends on g, which means that a universal correction factor does not exist.
or, alternatively by projecting g onto functions φ(t) = w(t)ψ(t). The function φ could be de-509 scribed as a localized oscillation. If one constructs in the time-frequency plane a rectangle of sides 510 T and ∆f centered, say, at t 0 = T 2 and f 0 = T ∆f 2∆t , equation 7 states that the area of this rectan-511 gle is constant, regardless of the value of T . This rectangle is sometimes called Heisenberg box 512 [Mallat, 1998]. Equation 7 is in fact an example of the application of the general Heisenberg un-513 certainty principle, which states that the area of a Heisenberg box cannot be made arbitrarily small.

514
For an arbitrarily-shaped localized oscillation φ(t) with Fourier transform ϕ(f ) and unit variance , defining the time and frequency widths as one can show [Gabor, 1946, Mallat, 1998, Percival and Walden, 2009] that 517 (10) In other words, it is impossible to achieve simultaneous arbitrary resolutions both in time and 518 frequency.

519
While the time-frequency resolution (area of Heisenberg boxes) cannot be made arbitrarily small, 520 it can be minimized. [Gabor, 1946] showed that the minimal area for a Heisenberg box is achieved 521 by localized oscillation φ(t) = w(t)e 2πif t where w is a a Gaussian-shaped window. Following using wavelets as a "frequency" representation results in a non-uniform frequency resolution, with 535 resolution degrading at higher frequencies.

536
Treating the wavelet transform as a frequency decomposition amounts to assigning a unique 537 frequency to each wavelet scale, in effect transforming the dual space from "scale space" to a 538 "frequency space". In the case of the Morlet wavelet, because the Gaussian is self-similar under 539 the Fourier transform, the frequency distribution of the wavelet is also a Gaussian, is symmetric and 540 has a maximum, therefore it seems logical to choose the peak frequency as the nominal frequency.

541
It should be clear that this scale-to-frequency transformation is arbitrary in the general case and 542 most likely meaningless.
FIGURE 11. A schematic of the Heisenberg boxes of the wavelet transform. The mother wavelet is obtained by the multiplication of a sine function with a window (panels a-c). d) Wavelets at three scales resulting by the scaling of the mother wavelet. e) Frequency distribution of the power of wavelets shown in e). f) Heisenberg boxes of the wavelet transform.
B.2. Spectral analysis. The spectral analysis of the LFP in the current study was based on standard techniques used for stationary signals [Priestley, 1981, Papoulis andPillai, 2002] as previously described in [Sheremet et al., 2016a]. Assume the LFP recordings g(t) and h(t) are realizations of zero-mean stochastic processes, stationary in the relevant statistics, with Fourier transforms G(f n ) and H(f n ), n = 1, · · · , N . The second and third order spectral statistics are estimated using cross-spectrum and bispectrum, defined as where the angular brackets denote the ensemble average, the asterisk denotes complex conjugation, 544 and we omit the superscript when a single time series is involved. The diagonal S g,g n of the cross- where the discrete two-dimensional interval I is typically defined as all the indices in the bis-560 pectrum triangle, i.e., pairs of integer indices (m, n), such that 1 ≤ m ≤ N and 1 < n ≤ 561 min(m, N − m + 1) ; due to the symmetries of the bispectrum, diagonal terms should be halved.

562
The proposed measure 16 is not independent of frequency resolution and and number of realiza-563 tions used, therefore the interpretation is somewhat loose: one should exercise care and discern-564 ment using it. Figure 12 illustrates this idea. The authors declare that the research was conducted in the absence of any commercial or finan-577 cial relationships that could be construed as a potential conflict of interest.