High-order theta harmonics account for the detection of slow gamma

ABSTRACTLocal 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 slow 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 low frequency burst (“slow gamma”). 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.


Introduction 56
In vivo neurophysiology has been a major tool in attempting to understand how the brain 57 organizes behavior, with local-field potential (LFP) oscillations taking center stage. The large 58 rhythmic oscillation in the hippocampus, theta (4-12Hz), was among the first to be characterized 59 in the freely behaving animal (Jung and Kornmüller, 1938;Green and Arduini, 1954; Green and 60 Machne, 1955; Vanderwolf and Heron, 1964;Vanderwolf, 1969; for review, see Buzsáki, 2005). 61 The faster gamma oscillation (Stumpf, 1965;Leung, 1992;Bragin et al., 1995;Chrobak and 62 Buzsáki, 1998;Penttonen et al., 1998), initially defined as a single, broad frequency range (e.g., 63 40-100 Hz, Bragin et al., 1995), has recently been subdivided into two, sometimes three 64 independent ranges referred to as slow-, medium-and fast gamma (Colgin et  identified as 20-50Hz (slow), 50-90Hz (medium), and fast (90-120Hz). The identification of more 69 than one gamma band has led to the multiple-gammas hypothesis of cognition, which proposes 70 that different gamma frequencies reflect different synaptic inputs to the hippocampus and 71 possibly distinct cognitive processing states (for review, see Colgin, 2015). Notably, the data 72 have been equivocal on what those states might be, with both slow and fast gamma being 73 Therefore, the exact relationship between the different gamma ranges and behavior is 80 not straight forward. The problem is further complicated by the existence of theta harmonics 81 (Harper, 1971;Coenen, 1975;Leung et al., 1982;Buzsaki et al., 1983;Leung and Buzsáki, 82 1983; Buzsáki et al., 1985;Ning and Bronzino, 1993;Czurkó et al., 1999;Terrazas et al., 2005), 83 which are frequency integer, phase coupled oscillations that precipitate in spectral 84 decomposition as a consequence of the fundamental 8Hz rhythm becoming both skewed and 85 asymmetric. As the harmonics of theta have been observed to go as high as 48 Hz (Sheremet 86  In order to address the overlap between low gamma and high-order theta harmonics in 91 the power spectra, we sought to disambiguate the theta harmonics from slow gamma using 92 current source density analysis across hippocampal lamina and spectral decomposition. 93 Contemporary theories of slow gamma suggest that the oscillation is generated via CA3 to CA1 94 projections whereas the fast gamma is generated by entorhinal to CA1 projections (Colgin et al.,95 written code as well as code imported from the HOSAtoolbox. Raw LFP records sampled at 24 205 kHz (Tucker-Davis system) were low-pass filtered down to 2 kHz and divided into fragments of 206 2048 time samples (approximately 1 second). To eliminate the effects of anatomic variations 207 between individuals and because the position and orientation of the hippocampal LFP probes 208 cannot be controlled precisely, the position of the hippocampal layers with respect to the 209 recording channels was determined by estimating the distribution of current-source density 210 (Rappelsberger et al., 1981;Mitzdorf, 1985;Buzsaki et al., 1986;Bragin et al., 1995). Current-211 source density exhibits a unique distribution that is largely independent of anatomic details and 212 that may be used to accurately identify hippocampal layers. 213

Spike autocorrelation and cross-correlation spectrogram analyses. 214
Often the presence or absence of an oscillation is buttressed by whether or not single units are 215 modulated at a specific frequency, accomplished for examining if the spike times of the 216 individual cells prefer a specific phase of an oscillation. However, nonlinear issues can become 217 a factor in these analyses (Belluscio et al., 2012), necessitating a different approach. Therefore, 218 to determine the frequency in which neuronal spiking occurred, we implemented spectral 219 analyses on spike trains (Leung and Buzsáki, 1983;Sheremet et al., 2016). Single unit data was 220 generously provided by the Buzsaki laboratory and curated by the Collaborative Research in 221 Computational Neuroscience Pastalkova et al., 2015;datasets: 222 i01_maze06.005, i01_maze08.001, and i01_maze08.004). For these dataset, the rat performed 223 a delayed alternation task on a figure-8 maze, running in a wheel during the delay. Only 224 bilaterally recorded CA1 neurons were used for the analysis. Action potentials for pyramidal 225 cells and interneurons were initially sorted in to velocity bins, analyzing the 5-15 or 35+ cm/sec 226 conditions separately. These spike times were converted into a binary time series with a 227 sampling frequency of the 1250 Hz. As this binsize is slightly smaller than the traditional window 228 of 1 msec used for waveforms, a box car convolution was performed such that a single spike 229 registers as 1 across 3 adjacent bins. Each spike train was then passed through Matlab®'s 230 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint xcorr.m function which finds the correlation at specified lag times. Lags we only calculated up to 231 3 seconds. Each output was then passed through the pmtmPH 232 (https://www.mathworks.com/matlabcentral/fileexchange/2927-pmtmph-m) function to calculate 233 the power spectral density (PSD) based on multitaper analysis. Each unit's power was 234 normalized by the total power and then sorted by peak theta frequencies to align the harmonics 235 between units. The average and standard error power spectra is presented.  The "coefficients" ( ) of the decomposition (also referred to as the transform of ) are obtained 244 by taking the inner product of the function with the basis elements (i.e., projecting ( ) onto 245 the basis), The pair of equations are best known as direct and inverse transforms, sometimes 246 also called analysis and synthesis. 247 The Fourier transform pair is obtained letting ( ) = 2 248 249 Unless the set of functions is carefully defined, these equations are only formal, and in most 252 cases they do not have any elementary mathematical meaning; when they do, their usefulness 253 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The wavelet transform and its inverse are given by 256 Where the constant depends only on 259 The functions , ( ) = � − � are "copies" of the "mother wavelet" ( ), shifted in time by 261 and scaled by , with a normalization constant. And � ( ) is the frequency space 262 representation of the mother wavelet. The full set of wavelets � , ( )� , ∈ is in general 263 complete but not independent (i.e., larger than a basis). 264

II. Discrete Transforms 266
Evolving the generic Fourier and Wavelet equations into useful analysis tools has complexities 267 that required the development of full theories (e.g., the theory of distributions, or generalized 268 functions, e.g., Schwartz, 1950;Lighthill, 1958). For example, in Fourier equations ( ) = 269 2 are orthogonal in the sense that ∫ 2 ∞ −∞ = ( ), where is the Dirac delta function 270 (Vladimirov, 2002;Strichartz, 2003). A complete discussion of these is far beyond the scope of 271 this study, and is also unnecessary, because time series measured in practical applications are represented by finite-dimensional linear operators, i.e., × matrices. These discretized 276 versions of the Fourier and wavelet representations are called the discrete transforms (Briggs, 277 1995;Mallat, 1999;Strang, 2006). 278

279
The discrete Fourier transform pair is (in Matlab® convention) 280 where Δ is a time-shift increment. The wavelets form orthogonal only for compact-support 292 wavelet shapes (e.g., Haar, and Daubechies wavelets; Daubechies, 1988Daubechies, , 1992Mallat, 293 1999). 294 The Parseval relation insures that the discrete Fourier equations conserve the variance of the 295 time sequence and its transform (e.g., Briggs, 1995), i.e., = , where = ∑ | | 2 is the 296 variance of . The discrete wavelet transform does conserve variance ( ≠ ) and in general, 297 the variance ratio depends on , which means that a universal correction factor does not exist. 298 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. To investigate the evolution of spectra over time, windowed Fourier transform (WFT) is 301 introduced where the original signal is multiplied by a window function which is nonzero 302 for only a short period of time where w is the window function (Priestley, 1981;Roads, 2004). The window function 305 slides along time axis, which gives rise to a time-frequency representation of the original 306 time series. This time-frequency representation can be plotted as the spectrogram.

307
Typically, window functions are smooth, "bell shaped", time-localized curves, the Fourier 308 transform of the window function is In frequency space, the window function W is usually composed of a "bell shaped" main 311 lobe and symmetric side lobes. The window bandwidth of a window function is defined . 314 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint The spectrum estimation depends critically on the window bandwidth as it will smooth 315 the spectrum and influence the frequency resolution (Fig. 1). The time-frequency 316 resolution of transforms will be further discussed in the following section.

317
In order to approach LFP interactions across regions, we implemented the method of 318 calculating the correlation coefficients of the spectrogram as outlined by Masimore and  Mode Functions (IMFs). White noise was added to the original signal before decomposition to 331 alleviate mode mixing, and was canceled out after the decomposition by averaging over 332 ensemble (Wu and Huang, 2009). In this study, the ratio between the variance of added white 333 noise and the original LFP was 0.5, and the ensemble number was 200. After applying the 334 EEMD, the power spectrum for each IMF was used to identify the peak frequency. The supra-335 theta signal is defined as the sum of all the IMFs with peak frequency above 12 Hz. The time- Heisenberg box (Mallat, 1999). This is in fact an example of the application of the general 355 Heisenberg uncertainty principle, which states that the area of a Heisenberg box cannot be one can show (Gabor, 1946;Percival and Walden, 1993;Mallat, 1999 In other words, it is impossible to achieve simultaneous arbitrary resolutions both in time and 361 frequency. While the time-frequency resolution (area of Heisenberg boxes) cannot be made 362 arbitrarily small, it can be minimized. Gabor (1946) showed that the minimal area for a 363 Heisenberg box is achieved by localized oscillation ( ) = ( ) 2 where is a Gaussian-364 shaped window. Following his work, Goupillaud (1984) used this shape, later called the Morlet 365 wavelet, to introduce the continuous wavelet transform Mallat, 366 1999). The limitations imposed by the Heisenberg uncertainty principle are illustrated in figure  367 2. A Morlet wavelet results as a product of a sine function with a Gaussian window. The 368 resulting function is the "mother wavelet". The wavelet transform then uses scaled versions of 369 the mother wavelet as "elementary" functions. 370 It is important to note that the wavelets shown in figure 2, panel d, are not harmonic, 371 therefore they are not uniquely characterized by a single frequency value; instead, they are 372 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint characterized by a frequency interval, say, corresponding to the width of the frequency spread 373 of their power (as given by the Fourier transform). The frequency distribution of power is, as 374 expected, narrower for the larger-scale wavelets, and wider for smaller-scale wavelets. 375 Therefore, the wavelet Heisenberg boxes coverage of the frequency axis is not uniform, and 376 consequently, using wavelets as a "frequency" representation results in a non-uniform 377 frequency resolution, with resolution degrading at higher frequencies. 378 379

IV. Numerical implementation 380
Fourier and wavelet analysis procedures were implemented in the Matlab® environment 381 using available Matlab® toolboxes. The Fourier cross-spectrum was estimated by dividing the 382  The wavelet scalogram is defined as the logarithm of the squared modulus |G mn | 2 of the 391 wavelet coefficients. The power spectral density for the wavelet transform were computed as 392 the time marginal of the wavelet transform (Abry et al., 1993). The Fourier spectra were 393 estimated using the Welch method (Welch, 1967). The anatomical connectivity of the hippocampus is well conserved such that CA3 input into CA1 403 terminates into the str. radiatum layer of CA1 while entorhinal input terminates in both the oriens 404 and str. lacunosum-moleculare (Sheremet et al., 2018b). As the entorhinal cortex is has been 405 reported to convey fast gamma into the hippocampus whereas the radiatum supports low 406 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.  Figure 3 shows the current-source density plots for the theta 408 frequency band as well as the first harmonic (15-17 Hz) and slow gamma (25-55Hz). With the 409 evidence that prospective coding increases with running velocity (  Unsurprisingly, the distribution of sources and sink for the theta and the first harmonic are nearly 419 identical with both frequencies exhibiting maximum current deflection in the hippocampal fissure 420 (compare with figure 2 of Bragin et al., 1995). However, of concern is the similarity of fissure 421 sources and sinks between the traditional gamma range (25-55Hz) and both the fundamental 422 and first theta harmonic. One sure way to cleanly differentiate oscillations would be through the 423 demonstration of minimal synaptic overlap. These results, however, reveal that the synapses 424 that terminate around the fissure support both theta, the harmonic and the 25-55Hz oscillation 425 band (cf. figure 4C of Belluscio et al., 2012). Therefore, the current source density analysis is -426 at the least -inconclusive regarding the disambiguation of low gamma from theta harmonics. At 427 the worst, it suggests that precisely the same synapses contribute to the spectral frequencies 428 between 8-55 Hz (suggesting that oscillations will inevitably be interdependent and coupled; 429 evidence against the multiplexing hypothesis). 430 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. are partially at odds. Specifically, these studies identified a slow gamma rhythm in the radiatum. 433 Figure 3: Current source density (CSD) analysis of seven rats during high speed movement. LFP recordings were split into 1 second segments and only segments with mean velocity larger than 35 cm/s were kept. LFP segments were aligned based on maximum theta peak in stratum LM and then were averaged. CSD analysis were carried out for three frequency bands: 1) theta, 7-9 Hz (top row); 2) theta harmonic, 15-17 Hz (middle row); 3) slow-gamma, 25-55 Hz (bottom row). For each panel, CSD values were normalized by maximum source/sink strength. In contour plots, warm color represents current source and cold color represents current sink. Considering the differences of probe positions and hippocampus sizes among rats, the CSD data were rescaled and shifted based on hippocampus lamina. In individual rat maps the vertical axis is channel number (not shown). Phase changes across channels (vertical) of the current-source density allow for the identification of oriens (Or), pyramidal layer (Pyr), stratum radiatum (Rad), lacunosum-moleculare (LM), molecular layer (M) and upper granule layer, or dentate gyrus (DG). It may be interpreted that the dentate is capable supporting a 25-55Hz frequency oscillation such as low gamma  or Beta (15-30 Hz; Rangel et al., 2015). We offer caution here as there are multiple theta diploes in the hippocampus with theta in the dentate potentially decoupled from theta in the str. radiatum (Montgomery et al., 2009). Furthermore, the theta oscillation in the dentate, cast a high number of harmonics relative to the other layers with the consequence of higher sources and sinks in the 25-55Hz band (see Fig. 5 below). not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint In particular, Belluscio and colleagues observed a significant current source density in the 434 radiatum when the LFP was filtered in the low gamma range (see Figure 4C of Belluscio et al., 435 2012). Therefore, we considered the possibility that due to unaccounted confounds, our results 436 in Figure 3 were incorrect. In light of this, we analyzed a single shank of a 256-site silicon probe 437 generously provided by the laboratory of Gyuri Buzsaki (Fernández-Ruiz et al., 2017). Current 438 source density was calculated during ripples as well as the LFP during behavior filtered for 439 theta, the first harmonic of theta, slow gamma (25-55 Hz) and a broad gamma (defined as 50-440 120 Hz based on power spectral density cross-correlation and bicoherence analyses presented 441 in Sheremet et al., 2018b). Similar to figure 3, we found nearly identical profiles between theta, 442 the first harmonic, low gamma and the broad gamma frequency (Fig. 4). Therefore, in yet 443 another preparation, the radiatum has no more or less "slow gamma" than would be anticipated  Figs. 1 and 2), the wavelet "convolves" over the higher order harmonics (potentially accounting 459 for why studies using wavelet decomposition rarely identify a harmonic above 16Hz). 460 Therefore, as an initial benchmark comparison of Fourier decomposition versus a Wavelet 461 decomposition, we generated a synthetic time series of pink noise embedded with an oscillation 462 at 8 Hz and harmonics at 16, 24 and 32 Hz (Fig. 5). Notably the power estimate of the Fourier 463 decomposition closely matched the theoretical distribution. However, as a consequence of the 464 multi-resolution analysis of the wavelet decomposition, the power estimate significantly expands 465 . The ripple triggered CSD clearly resolves sources and sinks associated with the pyramidal layer, radiatum, lacunosum moleulcare and dentate gyrus. The source sink profile is distinct between ripple events and theta. Note well, however, that CSD plots for theta, the 16 Hz harmonic, slow gamma, and traditional gamma, the source sink profiles are highly similar. Simply, current is generated where there are synaptic contacts. As there does not appear to be separate layers that support specific frequencies, a parsimonious conclusion is that the synapses themselves are agnostic to the LFP frequency that they support. claim that, following EEMD, the remaining supra-theta signal is "harmonic free". Therefore, we 471 conducted a direct comparison of fast Fourier transform, Wavelet decomposition and EEMD 472 from the lacunosum-moleculare of an awake-behaving rat as a function of velocity (Fig. 6). This 473 comparison clearly shows that, as anticipated from figure 6, that both the wavelet and EEMD 474 analytical approaches severely deviates from the Fourier representation at movement speeds 475 greater than 15 cm/s. Notably, this results in an apparent "convolution" of the third and fourth 476 harmonics of theta raising a significant cause for concern. In particular, wavelet and EEMD 477  (Fig. 1) convolve over the 24Hz and higher order harmonics of theta. The peak of these broad bands roughly resembles what is described as "low gamma". Also, please note that in the decomposition of the str. pyramidale, str. radiatum and dentate gyrus, wavelet analyses potentially underestimates the true power of the traditional 60-100Hz gamma (Fourier decomposition power exceeds the wavelet). This most-likely is a consequence of broad range averaging as the Heisenberg box becomes larger with higher frequencies. (Fig. 1). Finally, it Is worth noting that the removal of a fundamental oscillation, as in EEMD, reduces power in the theta rhythm, but fails to demonstrate the capability of addressing higher order harmonics. Rather, it parallels the outcome of the wavelet analysis. Once again, we remind the reader that, for wavelet and EEMD, the assignment of variance to any specific frequency is ultimately arbitrary and meaningless. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint With this concern in hand, we proceeded to revisit the wavelet scalogram analysis used 484 by  to investigate the interaction between gamma frequency and the theta 485 oscillation. This was first conducted for synthetically generated data (a sawtooth wave and an 486 oscillation with integer locked harmonics) as well as data collected from the dentate gyrus and 487 str. lacunosum-moleculare (Fig. 7). For the wavelet scalogram decomposition of the sawtooth 488 wave, an infinite series of harmonics capture the instantaneous change in amplitude. The 489 wavelet decomposition achieves this through a considerable peak in high frequency mid-trace. 490 However, this is accompanied by multiple low-level harmonics within the low gamma range. 491 While this is a simulation to an extreme derivative, the time-frequency representation the 492 synthetic harmonic series (8, 16, 24, and 32 Hz) is nearly identical to the slow gamma bumps 493 that can be seen in the decomposition of either the dentate granule layer or the lacunosum-494 moleculare (compare Fig. 7 to figure 1Eb of ). This "slow gamma" effect, 495 however, can be accounted for the smoothing of harmonics observed in figure 6. 496 As the results of wavelet and EEMD are rarely portrayed as frequency versus power, but 497 rather in terms of a time frequency representation, we furthered the comparison using two 498 different windows of analysis (Fig. 8). Both in the long, three-minute fast Fourier transform and 499 short, 6 second fast Fourier transform there is the clear presence of a 3 rd order theta harmonic 500 (~24-27 Hz) and a trace of a 4 th harmonic (~32-36 Hz; Fig.8, top panels). However, as 501 anticipated from the Heisenberg boxes (Fig. 2), examination of the wavelet on both the long and 502 short time series result in a "smearing" and apparent intermittency of the ~24-27 Hz frequency 503 band. Stated differently, the 24-27 Hz component is a direct consequence of theta being 504 skewed and asymmetric (for example, more "saw tooth" shape than a sinusoid; Buzsaki et al., inherently coupled to the fundamental 8 Hz frequency, but the wavelet decomposition artificially 507 detaches the 3 rd harmonic from the fundamental by using a different decimation of time and 508 frequency (Fig. 2). This "bleed" and intermittency offers a strong resemblance to experiments 509 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . Synthetic time series have the total length of 10 seconds and were added with a background pink noise. LFP recordings were selected during high speed running with the total length of 5 seconds. The averaged time series were obtained by averaging over theta (8 Hz) periods (top row). The wavelet scalograms for averaged time series were computed with wavelet transform (middle row). Averaged time series were filtered within theta band (7.5-8.5 Hz) which represent the theta phase (bottom row). LFP data were collected from r530. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint 3) Wavelet spectrogram of supra-theta signal (3 rd row). The supra-theta signal was obtained from ensemble empirical mode decomposition (EEMD) with noise level equaling to 0.5 total variance and ensemble number equaling to 200. The supra-theta was defined as the sum of decomposed modes whose central frequencies were large than 12 Hz. The wavelet spectrogram of supra-theta was computed with the method described above. Power spectra were computed by averaging the spectrogram over time (bottom row). For both epochs, the Fourier transform identifies theta and high order harmonics while wavelet analysis tends to resolve a theta rhythm and a wide-band frequency component (16-30 Hz). By considering this band as independent from theta in this manner gives the unintentional representation that there are "bursts" of gamma. Stated differently, the 24Hz oscillation is inherently dependent on the nonlinearity of theta ("saw tooth" shape of theta) and thus is incorrectly decomposed by the wavelet. Data were collected from r782. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint second harmonic of theta as well as the third harmonic (Fig. 8). Finally, the summary power 518 spectra for each method is presented demonstrating the inability of either wavelet or EEMD to 519 clearly resolve the theta harmonics. 520 Up to this point, the current manuscript has focused on the primary methods of spectral 521 decomposition that have led to the impression of a slow gamma oscillation in the hippocampus. 522 However, this does not encompass the entirety of support for slow gamma. For example, 523 analyses of hippocampal neuron spike times have revealed modulation to low gamma (e.g., 524 Colgin et al., 2009), with the logic being that as the cells exhibit a phase preference, one can 525 sketch a theoretical mechanism by which they support the oscillation. The concern, however, 526 has been raised that modulation in the 20-30Hz band may not actually be related to the low 527 gamma, but rather coupled to the asymmetry of theta (Schomburg et al., 2014). That is, the 528 spike depth of modulation plots are not pure sinusoids, but exhibit skew and asymmetry 529 (Skaggs et al., 1996;Quilichini et al., 2010). Furthermore, spectral decomposition of spike trains 530 has revealed the presence an ~16 Hz harmonic modulation (Sheremet et al., 2016). Therefore, 531 we considered the possibility that prior publications reporting phase coupling of neurons to slow 532 gamma find significance for the reason that higher order harmonics in the LFP have a direct 533 relationship to higher order harmonics in the frequency of spikes. Therefore, the point process 534 spike trains (readily available on the CRCNS.org website; see Methods) were sorted into either 535 a low or high velocity bin, converted into a binary timeseries, and either correlated against itself 536 (autocorrelogram) or other neurons (cross-correlogram). These correlograms were then 537 spectrally decomposed to examine the burst frequency modulation (see Geisler et al., 2007 for 538 similar methods). Notably, in accordance with prior publications, the "power" (number of spikes) 539 increases from low to high velocity and the burst frequency of the population increases with 540 velocity ( Fig. 9; Maurer et al., 2005;Geisler et al., 2007). Under the assumption that slow 541 gamma modulation is highest when prospective coding is the greatest (Bieri et al., 2014;Zheng 542 et al., 2016), which also describes the farther look-ahead at higher running speed 543 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint Fig. 9: Power spectral density of neuron autocorrelograms (top two rows) or crosscorrelograms (bottom three rows). The left column are the individual power spectra for autoand cross-correlograms taken at slow velocity (5-15 cm/s) while the middle column are the same neurons or neuron pairs at high velocity (35+ cm/s). The power spectral densities are normalized by the overall power of the high velocity spectrogram. The rows are sorted by the maximum frequency between 6-10Hz. The right column is the average power spectral density of the spike auto-and cross-correlograms at low (blue) and high velocities (red); errorbars are the standard error of the mean. The red dashed vertical lines are plotted at integers of 9 Hz up to 36 Hz. Note that neither the auto-or cross-correlogram power spectral show a single, slow gamma peak. Rather, should there be a notable deflection, it coincides with harmonics of the fundamental frequency. Also note that the dispersion in frequencies of the auto-and cross-correlations that include pyramidal neurons decreases as velocity increases. Energy/input stabilizes coherent oscillatory patterns. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint (Maurer et al., 2012), it may be anticipated that there would be a single peak in the gamma 545 range of the spike correlation spectrograms. Investigations of the spectrograms, however, do 546 not support this perspective. Rather, the small observable peaks tend to coincide with harmonic 547 lines (7.5, 15, 22.5, and 30 Hz). Therefore, the most reasonable conclusion is that phase 548 coupling observed between a 25-55Hz band and spikes is carried by the harmonics evident in 549 the LFP and the spike trains rather than being a consequence of a slow gamma oscillation that 550 is distinct from theta. Note well that, when the power spectra are sorted by the frequency of the 551 maximum power, pyramidal neurons (and their pairs) exhibit a greater frequency dispersion in 552 the autocorrelograms and cross-correlograms relative to high velocity -simply, oscillatory 553 dynamics are stabilized by energy/input (Churchland et al., 2010). 554 Significant emphasis has been placed on coherence and cross-frequency analysis as a 555 justification for low gamma. Caution, however, has been emphasized with respect to low values 556 of coherence (Buzsáki and Schomburg, 2015), alluding to the idea that low values have 557 contracted meaning. Furthermore, under the hypothesis that hippocampal LFP patterns are 558 organized by a cascade of energy from low frequencies to high frequencies (Buzsaki, 2006;559 Sheremet et al., 2018b), coherence should be significantly driven by activity into the 560 hippocampus (which is approximately correlated with running speed). Therefore, we 561 investigated coherence between the dentate gyrus and other lamina of the hippocampus as a 562 function of velocity (Fig. 10). For comparison, we have also included the power spectral density 563 measure and phase offset as a function of depth. At low velocities, coherence is fairly low 564 without any notable deflections above theta. At high velocities, however, there is both an 565 increase in the power of theta harmonics along with large coherence in theta harmonics across 566 hippocampal lamina. Given the substantial increase in coherence across the theta harmonics 567 with velocity, a traditional approach of considering coherence that does not account for velocity 568 may unintentionally smear the range between 16-48 Hz, giving the impression of a single peak. 569 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint Coherence analysis, however, to a large extent is unaffected by changes in power (of 570 course, an amplitude must be present to define a phase). Therefore, we considered the 571 Figure 10: Distribution of cross-spectrum over hippocampal layers during active behavior, at low speed (top three rows) or high speed (bottom three rows). 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. Distribution of cross-spectrum over hippocampal layers during active behavior), at high speed. Note that at high velocities, the power spectra, coherence and phase measures support the observation of multiple theta harmonics. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint argument that slow gamma is not entrained through the hippocampus in terms of phase-phase 572 coupling. Rather, power is the important driver of coupling. Therefore, using the methods of 573 Masimore et al. (2004;, we investigated how changes in power within either the medial 574 entorhinal cortex or dentate gyrus correlated with power changes across the lamina of CA1. In 575 accordance with our recent publication (Sheremet et al., 2018b), we only found evidence of 576 cross-regional power interactions between theta, theta harmonics and a unitary gamma (Fig.  577   11). From these power-power analyses, there is no reason to believe a slow gamma oscillation 578

exists. 579
Finally, if the observation of slow gamma were a consequence of theta harmonics, then 580 there should be a strong dependence on velocity. Investigation of the primary literature, 581 Figure 11: Average cross correlation coefficients of Fourier transform across three rats with electrodes in the medial entorhinal cortex in order to identify fundamental frequencies of interaction. As previously described, regions of non-zero correlation potentially indicate that two oscillatory bands showed high power at similar times suggestive of an interaction (see Masimore et al., 2004;. In both cross-correlations between the DG and hippocampal lamina as well as the MEC and hippocampal lamina, there are notable "dots" of correlations below at ~40Hz and below, indicative of theta and its harmonics (Sheremet et al., 2016;. Further off-axis, there is a correlation between ~80-120 Hz and theta as well as its harmonic. Notably, there is no reason to expect the need to subdivide gamma into multiple bands. Rather, the interactions are constrained to theta, its harmonics and the broad gamma band as originally described by Bragin et al. (1995;50-120 Hz). For all plots, values below 0.03 are set to white. Or, oriens; Pyr, CA1 pyramidal layer; Rad, Radiatum; and LM, lacunosum-moleculare. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint however, is not so clear cut. For example, Chen and colleagues observed a linear change in 582 slow gamma with velocity (2011), supporting the harmonic hypothesis. Ahmed and Mehta, on 583 the other hand, found a reduction in slow gamma power with running speed (2012). Analytically, 584 these results are not incorrect, but more nuanced. Specifically, as animal running velocity 585 increases, there is an associated increase in theta power and the harmonics. These increases 586 in power can be viewed as the energy into the hippocampus increasing resonance or bringing 587 otherwise incoherent neural activity into coherent, dynamic patterns. The consequence is that 588 the incoherent activity surrounding the theta and harmonic bands: 1-6, 10-14, 19-23, and 25-40 589 Hz will effectively lose power (Fig. 12). Note well, this idea is against the idea of spectral 590 multiplexing, in which each frequency has an independent meaning. Thus, the parsimonious 591 explanation is that neurons can "drift" (exhibit a range of oscillations) when provided a low level 592 of input, but become strongly entrained with energy and velocity (Sheremet et al., 2018b). This 593 accretion of power in one or multiple bands comes at the expense of losing power in others.  contradicted by Fourier analysis. Specifically, the Fourier decomposition exhibits significant 608 peaks in the high-order theta harmonics (24, 32, and 40 Hz). Critically, these contradictory 609 outcomes are incompatible, necessitating a reevaluation of multiple concepts. 610 Figure 12: Spectral changes with velocity. It is often reported that there is a either no change or a decrease in slow gamma with increasing running velocity. While this is analytically, it is most-likely a result of the same reduction in power from 1-6, 10-14, 19-23, and 25-40Hz. That is, when a network of neurons receives a low level of input, the activity is unorganized, resulting in a smooth decrease in power from 1-100Hz. However, as velocity is associated with increasing firing rates, the network organizes as a function of increasing energy (synaptic activity). Input stabilizes oscillatory dynamics (see Fig. 9; for similar data, see Churchland et al., 2010). By acting on one another, neurons will self-organize into a pattern of activity that supports theta, its harmonics and gamma. Thus, the reduction of the 25-40 Hz band is most likely a consequence of neurons becoming entrained. This effectively challenges the argument for slow gamma (it may be considered inappropriate to suggest an oscillatory band disappears when neural entrainment occurs). Compare to Ahmed and Mehta (2012; Fig 2). not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint The fundamental issue is what defines an oscillation or rhythm? The identification of a 611 hippocampal oscillation is often based on visual examination, followed by analytical methods. 612 This has been true for theta which corresponds to a significant local maximum above the 613 background activity (Vanderwolf and Heron, 1964;Vanderwolf, 1969). Note, that in this 614 definition, the power distribution defines the frequency band. If the 6-8 Hz band deviates from 615 the noise slope, it is interpreted as the presence of an oscillation. This distinction is peculiar, 616 since in fact any frequency band with non-zero power is in the Fourier representation an 617 oscillation. 618 Moreover, the "peak above noise" definition of a rhythm does not work in general for the 619 gamma rhythm, because (in all cases but high running speed) there is no "significant local 620 maximum" in the power spectral density. Rather, gamma is an accumulation of energy that is 621 often made evident by compensating (e.g., whitening) the spectrum (n.b., finding peaks after 622 standard whitening procedures are problematic if the spectra does not follow a unitary slope 623 which appears to be true for hippocampal LFP; Sheremet et al., 2018b;Sheremet et al., 2018a). 624 This raises the question "what is so special about spectral deviations that single out those 625 frequency bands from 'noise bands' that have comparable energy?". The answer seems to be 626 the idea that there is a "background spectrum" that has a different meaning and does not really 627 matter enough to be called a rhythm. Only accumulations of energy past this background level 628 are significant. The significance of these accumulations of energy, however, is unclear in the 629 absence of a physical model of the phenomenon. Such a model should explain both the 630 processes underlying the "non-rhythm" background spectrum, and the "rhythm" spectral 631 deviations (e.g., Sheremet et al., 2018a). In failing to do so, the multiplexing model of cognition 632 seems improvised to account for the wavelet decomposition results. One must conclude that 633 the definition of a "rhythm" based on the power distribution is an arbitrary decision, because it is 634 based on the "aspect" of subjective representations and not on a fundamental, analysis-635 invariant description of the physics underlying the LFP process. 636 In contrast to the multiplexed "spectral finger print" perspective (e.g., Canolty  While these are supporting pieces of data, they neither confirm or deny the presence of an 664 oscillation. Nevertheless, we have addressed each of these in the current manuscript. Notably, 665 rather than perform a second-order analysis of investigating either phase locking to LFP or 666 spike-LFP coherence, we simply sought to determine what frequencies are evident in the action 667 potentials, using Thomson multitaper decomposition on auto-and cross-correlograms (Leung 668 and Buzsáki, 1983). We did not observe any slow gamma modulation, but did observe multiple 669 harmonics of theta ( Fig. 9; Sheremet et al., 2016). Therefore, a wide band filter of 25-55Hz will 25-55Hz will be harmonic dominated. Furthermore, we have replicated cross-regional 675 coherence, this time using velocity as a parameter, to investigate whether it is tenable to believe 676 a broad slow gamma band exists. Although prior research would suggest that the dentate 677 entrains other regions through slow gamma , our analyses here (Fig. 10) and 678 elsewhere (Sheremet et al., 2018b) reveal coherence up to the 40Hz harmonic of theta. 679 Furthermore, cross-correlation of power between the medial entorhinal cortex or dentate gyrus 680 and layers of CA1 fail to show any interactions between frequencies other than theta, its 681 harmonics, and a single 50-120 Hz gamma band (Fig. 11). Finally, we have addressed all forms 682 of cross-frequency coupling -phase-phase, phase-amplitude, and amplitude-amplitude -in 683 another manuscript (Sheremet et al., 2018b). Therefore, when considering both the primary and 684 secondary evidence supporting the existence of slow gamma, there appear to be fundamental, 685 critical shortcomings. 686 Of course, our observations do not cover the compendium of manuscripts that have 687 supported the existence of slow gamma. An attempt to refute them all may prove to be a 688 there is no reason to believe that an independent slow gamma oscillation exists in the normal 710 rat. However, our failure to observe this oscillation does not mean that the brain cannot venture 711 into a parameter space in which slow gamma is prominently detectable from noise. It is 712 therefore worth discussing what empirical data would be needed to support the assertion of 713 slow gamma. Control for theta harmonics using bicoherence analysis should be required. As the 714 argument has been made that slow gamma is coupled to theta, a strong interaction in the 715 bicoherence between slow gamma and theta would further enhance this claim. Notably, the 716 bicoherence analysis detects a strong coupling of the 50-120 Hz gamma to both theta and its 717 first harmonic (Sheremet et al., 2018b). It should be mentioned that "wavelet bicoherence" 718 analysis (Van Milligen et al., 1995) would suffer similar pitfalls to wavelet time-frequency 719 representation. Moreover, considering asymmetry alone (e.g., Amemiya and Redish, 2018) is 720 an insufficient control as oscillatory skew can also contribute to harmonics (Sheremet et al., 721 2016). Furthermore, as discussed above, oscillatory identification is partially determined by a 722 deflection that significantly deviates from the "background" activity. Therefore, a simple, 723 unwhitened, log-log power spectral density with a peak in the 25-55 Hz band using the Fourier 724 transform would be indicative that there is a uniform signal that exceeds noise. Finally, 725 controlling for animal velocity when conducting these analyses should become common 726 practice. The secondary analyses should only be implemented once due diligence has been 727 employed to ensure that a signal exists. 728 Through this presentation, we believe that fundamental, critical analytical inconsistencies 729 are responsible for "slow gamma" detection that are instead a consequence of the theta 730 harmonics. Moving forward, vigilance is necessary to ensure that the analytical treatment is 731 "true to the essence of the biological process studied" or otherwise risk generating mistakes as 732 a consequence of inappropriate methods (Marder, 2015). While quantitative data analysis 733 reflecting what can be observed qualitatively is the first step toward a model, caution should be 734 taken prior to applying analytical techniques that incorporate specific assumptions about the 735 system without a physical model. 736 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint of the properties of the stochastic process. This is, however, essential for interpreting the 789 observations and insuring that one does not elaborate on meaningless, at worst non-existent, 790 quantities. It is important to understand that just because a finite number of realizations always 791 yields an average does not necessarily imply that the stochastic average (over all realizations) 792 exists. Therefore, in practice, one is forced to assume a model for the stochastic process, which 793 should guide the effort to estimate its characteristics. 794 The dissimilarity of the two transforms becomes particularly acute when considering 795 observations as realizations of a stochastic process. For Fourier analysis, a stochastic process 796 model exists: it is the generalized harmonic process (Priestley, 1981;Percival and Walden, 797 1993), a zero-mean, variance-stationary process (variance per unit of time is conserved). The existence of the spectral density of the process is guaranteed by the celebrated 803 Wiener-Khinchin theorem (Priestley, 1981), which states that the power spectral density of 804 process is the Fourier transform of its auto-covariance function. It is important to note the 805 condition of orthogonality of increments: it means that the Fourier components of different 806 frequencies are uncorrelated. Therefore, the Fourier stochastic mode cannot have (does not 807 include) non-stationary realizations such as sharp waves. 808 As far as we are aware, there no widely-accepted stochastic-process model for transient 809 signals. therefore, it appears that it makes little sense to look for a stochastic scope of the 810 wavelet transform, as it is hard to imagine an overarching stochastic model for non-stationary 811 processes (although recent attempts, e.g., Antoniou and Gustafson, 1999, might prove this 812 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint argument wrong). While one is free to average whatever feature seems relevant, we do not in 813 fact know under what conditions a particular average exists. 814 This discussion should show that the Fourier and wavelet representations are not 815 interchangeable. For example, remapping a wavelet scalogram to frequencies and averaging 816 over time to estimate the distribution of the power over frequencies has at least two fundamental 817 problem: the frequency mapping is basically arbitrary, and the time average is not defined, if the 818 process is non-stationary. Alternatively, if the process is stationary, the Fourier transform is the 819 correct tool, because the wavelet distorts the distribution due to its peculiar shape of time-820 frequency atoms. 821

A simple stationarity test 822
Remarkably, the simple idea that a variance-stationary generalized harmonic process should 823 have at most weakly-correlated frequency components suggests a simple test for stationarity. 824 The bicoherence of a stationary stochastic process should be weak (Elgar, 1987 stationary series, although weakly asymmetric (not entirely free of cross-spectrum correlation) 827 exhibits a bicoherence that is essentially statistically zero, with the exception of the harmonic 828 phase coupling. One could classify this series as matching the generalized harmonic stochastic 829 process. In contrast, the non-stationary example in figure 13 shows large areas of strong 830 coupling, suggesting that this should indeed be classified as non-stationary. 831 not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted December 20, 2018.  The copyright holder for this preprint (which was this version posted December 20, 2018. . https://doi.org/10.1101/428490 doi: bioRxiv preprint