Brain waves are manifestations of synchronized neuronal currents widely used for describing neurophysiological activity (Fries, 2005; Buzsáki, 2011; Thut et al., 2012; Cannon et al., 2014). However, our understanding of these phenomena depends fundamentally on mathematical and computational tools used for analyzing the recorded Local Field Potentials (LFPs). Most computational methods are based on breaking the signal into a combination of basic components suggested by the study specifics, e.g., wavelet analysis is most appropriate for studying time-localized events, such as ripples or spindles (Battaglia et al., 2004; Bosnyakova et al., 2006; Sitnikova et al., 2009; Luijtelaar et al., 2011), whereas Fourier decomposition is used for describing the oscillatory patterns of LFPs (Roopun et al., 2008; Aru et al., 2015; Lozano-Soldevilla et al., 2016; Cole and Voytek, 2017). Since most techniques are backed up by a completeness theorem, it may appear that selecting a specific decomposition is only a matter of convenience. This, however, is not the case: given that physiological mechanisms of the LFP oscillations and their functions are not yet fully understood, the task of establishing a physically adequate description of the signal's structure is not idle (Kopell et al., 2010; Buzsáki at al., 2012; Sreenivasan and D'Esposito, 2019). One may draw here a historical parallel with the use of the Ptolemaic system, in which every movement of a celestial object could be decomposed into a sufficient system of epicycles (Hanson, 1960; Van der Waerden, 1974, 1982; Babb, 1977). However, it was the discovery of the heliocentric system that eventually revealed the physical laws governing planetary motion (Gallavotti, 2001).

Discrete Fourier Transform (DFT) converts data series into a superposition of discrete harmonics with fixed frequencies, proportional to a certain base frequency ω₀ (Brigham, 1988). This built-in rigidity of the Fourier spectra leads to a well-known conflict between the temporal and the frequency resolutions, manifested in many fields, from biology to Quantum Mechanics, which limits the method's resolution (Folland and Sitaram, 1997; Grünbaum, 2003). In the following, we use an alternative technique—Discrete Padé Transform (DPT)¹, which also converts data points into a superposition of harmonics. However, the DPT harmonics are free to change frequencies independently, adapting their values on a moment to moment basis through the Padé Approximation Theory algorithms (Baker and Graves-Morris, 1996; Bessis, 1996; Bessis and Perotti, 2009; Perotti et al., 2013, 2019; DeVito and Dabaghian, 2014; Perotti and Wojtylak, 2018).

The spectrograms of the LFPs recorded in the CA1 area of the rat's hippocampus, built using a “sliding window” version of DPT (refer to section 4), reveal patterns that open a novel perspective on the analyses of extracellular field dynamics. First, there appear to be two types of reconstructed frequencies. The first kind changes regularly across time, leaving distinct traces in the spectrogram—the spectral waves (Figure 1A). The frequencies of the second kind assume sporadic values from moment to moment and correspond to instantaneous “irregular” harmonics with much lower amplitudes. The nature of these two classes of harmonics can be explained based on several subtle theorems of Complex Analysis (Steinhaus, 1929; Froissart, 1973; Gilewicz and Pindor, 1997; Gilewicz and Kryakin, 2003). In essence, it turns out that the irregular harmonics represent the signal's noise component, ξ(t), whereas the regular, stable frequencies define its genuine oscillatory part, r(t) (Bessis, 1996; Bessis and Perotti, 2009; Perotti et al., 2013, 2019; DeVito and Dabaghian, 2014; Perotti and Wojtylak, 2018). Interestingly, the superposition of the regular harmonics, which typically constitute only 1 − 5% of the full set, captures the shape of the original signal with over 90 − 95% precision. Correspondingly, the contribution of the remaining 95 − 99% harmonics is small, typically less than 5 − 10% of the signal's amplitude (Figure 1B). Thus, according to the DPT, the brain waves consist of a few phase-modulated waves embedded into a weak noise background.

The individual oscillatory terms in (1), ϑq(t)=AqeiΦq(t), are referred below as brain wave oscillons (DeVito and Dabaghian, 2014; Perotti et al., 2019).

Since the decomposition (1) emerges through empirical analyses, with no a priori assumptions or ansatzs, the oscillons may capture the physical organization of synchronized neuronal activity and help link empirical observations to theoretical models (Berger, 1933; Hoppensteadt and Izhikevich, 1997; Boashash, 2003; Vugt et al., 2007; Colgin, 2016).

Second, higher temporal resolutions reveal a quasiperiodic pattern of the reconstructed frequencies,

where ω_q,0 is the mean frequency, and ω_q,i are the magnitudes of the embedded undulations with frequencies Ω_q,i and phases φ_q,i (Figure 1C). Importantly, the mean frequencies of the spectral waves dovetail with the mean frequencies of the traditional (i.e., Fourier-defined) rhythms. For example, the mean frequency of the lowest spectral wave (about 8 Hz) matches the mean θ-frequency and the mean frequency of the next spectral wave (about 32 Hz) aligns with the characteristic slow-γ frequency. Furthermore, the spectral undulation magnitudes are consistent with the widths of the corresponding Fourier bands (Senior et al., 2008; Carr et al., 2012; Colgin, 2015), which allows using the standard nomenclature, e.g., ω_θ(t) for the spectral θ-wave, ω_γ1(t) for the spectral slow-γ wave, and to write the oscillon decomposition (1) in the form

The analysis of the spectral waves carried in Perotti et al. (2019) was motivated by the assumption that, in a given physiological state, the magnitudes ω_q,i and the embedded frequencies Ω_q,i in the expansion (2) are relatively stable and extractable through Fourier-based analyses, such as Welch's transform (Welch, 1967; Proakis and Manolakis, 1996). Indeed, the power profiles of approximately 1 s long segments of spectral waves exhibit consistent series of isolated peaks, suggesting that the hippocampal oscillons are driven by a discrete and comparatively scarce set of spectral harmonics (Figure 1D). However, further analyses revealed that the spectral dynamics are substantially more complex, as discussed below.

To validate the qualitative conclusion drawn from Figure 2C, we simulated a superposition of two oscillons with the spectral wave parameters derived from the recorded data. For example, the waves illustrated in Figures 2A,B were generated for the mean θ and γ frequencies, ω_θ,0 ≈ 8 Hz and ω_γ,0 ≈ 32 Hz, along with various specific sets of the reconstructed embedded frequencies, e.g., Ω_θ,* ≈ {1.9, 4.4, 6.2, 8.2, …} Hz and Ω_γ,* = {1.4, 3.9, 6.9, 9.2, …} Hz. These values were used to build “synthetic” θ and γ oscillons with spectral frequencies

which were then processed using DPT algorithms.

As expected, the stable frequencies reconstructed through DPT procedures produce clear point traces across the spectrogram, and the reconstructed contiguous spectral waves match the input data (Figures 3A,B). Correspondingly, the peaks representing the embedded frequencies appear in the Welch's spectrogram in correct positions and remain steady, nearly unchanged over the entire duration of the signal (Figure 3C). Furthermore, numerous computational experiments with synthetic oscillons produced no spurious peaks or other artifacts suggestive of the patterns visible in Figure 2C.

The ostensible difference between the spectrogram produced by the simulated oscillons with constant spectral waves (Figure 3C) and the ones reconstructed from the recorded LFP data (Figure 2C) suggests that the hippocampal extracellular field dynamics may not be described by quasiperiodic series (4a) with steady coefficients. The time-localized peaks visible in Figure 2C suggest that the hippocampal frequency spectra are disturbed by rapid, transient processes that appear for a short time and rapidly disappear. To verify this possibility, we applied the DPT analyses to a numerically generated signal in which the spectral waves with constant coefficients (4) were replaced by a superposition of harmonics with time-localized spectral magnitudes,

where ω^θ,i(t) are narrow (σt2≈70-80 ms) Gaussian pulses localized at a few discrete moments (Figure 4A). These “spectral kicks” are clearly manifested on Welch's spectrogram computed directly for the simulated spectral θ- and slow-γ waves (Figure 4B), but they do not significantly alter the reconstructed stable frequency traces (Figure 4C). Applying DPT analyses to the corresponding oscillons produces Welch's spectrograms that bear an uncanny resemblance to the spectrograms obtained for hippocampal LFP spectral dynamics (Figure 4D). The latter result suggests that the hippocampal oscillons may exhibit elaborate behaviors that include rapid, nonstationary spectral modulations that may be due to the extracellular field's endogenous dynamics or to inputs from parahippocampal or cortical networks.

Discrete Fourier Transform is produced by convolving the data values, s₁, s₂, …, s_N, with a discrete set of harmonics with fixed frequencies,

arranged uniformly over the unit circle in the complex plane. The closer is the discrete frequency ω_l = 2πl/N to the frequency of the signal's constituent waves, r(t)=∑pApeiωpt, the bigger is the contribution of the corresponding harmonic into the decomposition (Brigham, 1988). If the data are sampled from a combination of harmonic oscillations and a noise background,

then each frequency ω_p produces a Fourier-peak, broadened and lowered by the noise ξ(t) (Newland, 2005; Perotti et al., 2014).

The DPT extends the expansion (8) from the unit circle into the complex plane, ei2πlN→z,

where z is a generic complex number. For the oscillatory component of s(t), the sum (9), extended to infinity, yields a meromorphic function,

whose poles, zp=e-iωpτ, and residues, apeiφp, define the frequencies, the amplitudes, and the phases of the contributing harmonics (Bessis, 1996; Bessis and Perotti, 2009; Perotti et al., 2013, 2019; DeVito and Dabaghian, 2014; Perotti and Wojtylak, 2018).

The sub-diagonal Padé approximant to (10),

rapidly approaches R(z) as the degree N of the polynomials P_N(z) and Q_N+1(z) grows, R(z)=RN(z)+O(z2N) (Baker and Graves-Morris, 1996). In particular, the poles z_p of R(z) are approximated by the roots ζ_q of the denominator in (11), Q_N+1(ζ_q) = 0 (Bessis, 1996; Bessis and Perotti, 2009; Perotti et al., 2013).

As for the z-transform of the noise component,

the Steinhaus theorem establishes that its poles appear at the unit circle with unit probability (Steinhaus, 1929). The manifestation of this effect in the Padé approximations to Ξ(z) is subtle: the “noisy” poles are the ones that not only cluster around the unit circle, but also pair with the zeroes of Ξ_N(z), thus forming the so-called “Froissart doublets” (Froissart, 1973; Gilewicz and Pindor, 1997; Gilewicz and Kryakin, 2003). A typical pole-zero distance in these pairs is smaller than 10⁻⁶ − 10⁻⁷ in the standard Euclidean metric in C¹. Furthermore, the Froissart doublets are unstable with respect to variations of the algorithm's parameters, in contrast with the unpaired, stable poles produced by the regular part of the signal (Froissart, 1973; Bessis, 1996; Bessis and Perotti, 2009). These qualitative differences allow the separation of the regular component of the signal from its noise background, as expressed by the decomposition (1). The original study of Steinhaus (1929) presumed uniformly distributed noise series; subsequent works cited above allow generic, continuous noise distributions.

Sliding window or the Short Time Padé Transform (STPT) uses a segment of the signal of length T_W, centered at time t_i, to extract the time-localized spectra—in full analogy with the Short Time Fourier Transform, STFT (Howell, 2001; Jacobsen and Lyons, 2003). Plotting the reconstructed frequencies along the vertical axis and arranging the times t_i horizontally yields the Padé spectrogram, which we use to illustrate spectral dynamics, in direct analogy with the standard Fourier spectrograms.

The mean amplitude of the input data was normalized to s-(t)=2. The LFPs were originally recorded at the rate S_r = 8 kHz. To increase time resolution in the biologically relevant range of frequencies (f <300 Hz), we interpolated the signal to higher rates (S~r=30 kHz, S~r=36 kHz or S~r=44 kHz), which did not alter the shape of the studied spectral patterns but significantly improved stability and sharpness of the results. The oversampled time series were then downsampled 2 ≤ m ≤ 4 times, which produced m interlaced subseries that were independently studied with DPT. As one would anticipate, the stable frequencies generated by each subsequence form tight clusters of m points, grouping around the frequency produced by the original sequence, while the Froissart doublets exhibit erratic behavior (Bessis, 1996; Bessis and Perotti, 2009; Perotti et al., 2013, 2019; DeVito and Dabaghian, 2014). These procedures allow using time windows as short as T_W = 10 − −20 ms while keeping the order of the Padé approximants high, N = 100 or more. Shifting the time windows by a single data point ensures maximal contiguity of the reconstructed spectral waves and the oscillons' amplitudes. The Froissart distance used to identify close pole-zero pairs (Froissart doublet) is dF=10-6. To increase stability, the signals were filtered between 1 and 40 Hz