## Abstract

The nonlinear, metastable dynamics of the brain are essential for large-scale integration of smaller components and for the rapid organization of neurons in support of behavior. Therefore, understanding the nonlinearity of the brain is paramount for understanding the relationship between brain dynamics and behavior. Explicit quantitative descriptions of the properties and consequences of nonlinear neural networks, however, are rare. Because the local field potential (LFP) reflects the total activity across a population of neurons, nonlinearites of the nervous system should be quantifiable by examining oscillatory structure. We used high-order spectral analysis of LFP recorded from the dorsal and intermediate regions of the rat hippocampus to show that the nonlinear character of the hippocampal theta rhythm is directly related to movement speed of the animal. In the time domain, nonlinearity is expressed as the development of skewness and asymmetry in the theta shape. In the spectral domain, nonlinear dynamics manifest as the development of a chain of harmonics statistically phase coupled to the theta oscillation. This evolution was modulated across hippocampal regions, being stronger in the dorsal CA1 relative to more intermediate areas. The intensity and timing of the spiking activity of pyramidal cells and interneurons was strongly correlated to theta nonlinearity. Because theta is known to propagate from dorsal to ventral regions of the hippocampus, these data suggest that the nonlinear character of theta decreases as it travels and supports a hypothesis that activity dissipates along the longitudinal axis of the hippocampus.

**SIGNIFICANCE STATEMENT** We describe the first explicit quantification regarding how behavior enhances the nonlinearity of the nervous system. Our findings demonstrate uniquely how theta changes with increasing speed due to the altered underlying neuronal dynamics and open new directions of research on the relationship between single-neuron activity and propagation of theta through the hippocampus. This work is significant because it will encourage others to consider the nonlinear nature of the nervous system and higher-order spectral analyses when examining oscillatory interactions.

## Introduction

The nonlinear character of the brain (Buzsáki, 2006) has been recognized for >50 years (Ashby, 1947; Wiener, 1966) as the foundation of large-scale integration across local neuron structures (Steriade, 2001; Buzsáki and Draguhn, 2004) and for the metastable dynamics critical for rapid neuron organization in support of behavior (Engström et al., 1996; Friston, 1997; Tognoli and Kelso, 2014). Understanding brain nonlinearity is therefore fundamental for determining how brain dynamics translate to behavior (Hasselmo, 2015; Marder, 2015), yet quantitative descriptions of brain nonlinearity are scarce (Buzsáki, 2006). Using a thermodynamics analogy, “understanding the brain” can be achieved by describing either its microscopic states (complete description of the state of each neuron) or its macroscopic states (sets of “statistically equivalent” microscopic states). Although the former does not seem achievable (Marder, 2015), even with new high-density measures of real-time neural activity (Ziv et al., 2013), a macroscopic model of brain dynamics could be built on existing measurement techniques, such as the local-field potential (LFP). The LFP reflects the activity of a larger number of neurons (Buzsáki, 2002; Buzsáki et al., 2012; Schomburg et al., 2012), providing a macroscopic signature of microscopic activity.

Although the characteristics of the hippocampal LFP in relation to behavior are well documented (Buzsáki, 2005), its nonlinearity has not been examined explicitly. The goal of the current investigation was to quantify the nonlinearity of the hippocampal theta rhythm in relation to movement. Theta is a dominant feature of hippocampal LFP commonly reported as a 4–12 Hz oscillation associated with active exploration and REM sleep (Jung and Kornmüller, 1938; Green and Arduini, 1954; Vanderwolf, 1969) that may serve to organize hippocampal neuron firing. During movement, principal cells of the hippocampus express spatial receptive fields (O'Keefe and Dostrovsky, 1971). As an animal passes through a neuron's place field, the spike timing systematically shifts to earlier phases of theta between adjacent cycles (O'Keefe and Recce, 1993). The rate at which this precession occurs varies along the longitudinal axis of the hippocampus (Maurer et al., 2005). Interestingly, differences between dorsal and more ventral regions of the hippocampus are also observed in theta amplitude and the firing rate of neurons, such that firing rates and theta power in intermediate CA1 show less velocity modulation than in dorsal CA1 (Maurer et al., 2005). This suggests that dynamic interactions between hippocampal neurons and LFP may reflect the behavioral state of an animal, but this varies according to longitudinal position.

In addition to amplitude, the shape of the theta wave is sensitive to movement, with a transition from a sinusoid to sawtooth shape at faster running speeds (Buzsáki et al., 1983; Terrazas et al., 2005). Moreover, a 16 Hz oscillation has been observed to develop as a function of velocity (Czurko et al., 1999; Terrazas et al., 2005), described as the second harmonic of theta (Harper, 1971; Coenen, 1975; Leung et al., 1982; Leung, 1982). The standard technique used to investigate theta, however, has typically been based on Fourier power-spectrum analysis. Because it ignores phases, this approach is fundamentally linear. Therefore, the nonlinear mechanisms underlying the relationship between the 16 Hz oscillation and the sawtooth shape of theta in relation to running speed (Skaggs, 1995) and dorsoventral position is not well understood.

The theta rhythm is also known to propagate along the hippocampal longitudinal axis (Lubenov and Siapas, 2009; Patel et al., 2012), losing amplitude as it moves into more ventral regions (Maurer et al., 2005; Royer et al., 2010). It is not known, however, how this relates to longitudinal differences in the modulation of neuron activity by velocity. To understand the nature of the sawtooth shape of theta and the 16 Hz frequency component in the hippocampus, we investigated the nonlinear character of the theta rhythm in the dorsal and intermediate regions of the hippocampal CA1 subregion. The nonlinearity was quantified in terms of phase coupling between different Fourier frequency bands that were estimated using third-order Fourier statistics.

## Materials and Methods

##### Subjects and behavior.

Neurophysiology data were obtained from three Brown Norway/Fisher-344 hybrid male rats between 8 and 12 months of age. The rats were housed individually and maintained on a 12 h light/dark cycle. Recordings took place during the dark phase of the cycle. Surgery was conducted according to the National Institutes of Health's guidelines for rodents and approved Institutional Animal Care and Use Committee protocols.

Before surgery, the animals were food deprived to 85% of their *ad libitum* weight. During this time period, the rats were trained to run on circular tracks for food reinforcement. Food was given on either side of the barrier and at the 180° opposite point. Rats ran for ∼20 min, resulting in a variable number of laps per session depending on motivation and other factors. Each track running session was flanked by a rest period in which the rat rested in a towel lined pot located near the track. For the present analyses, only data obtained during running conditions are presented.

##### Surgical procedures.

Neuronal recordings were made in the dorsal and intermediate CA1 subregion of the hippocampus. Before surgery, rats were administered ampicillin (Bicillin; Wyeth Laboratories; 30,000 U, i.m., in each hindlimb). The rats were implanted, under isofluorane anesthesia, with an array of 14 separately movable microdrives (“hyperdrive”). This device, implantation methods, and the parallel recording technique have been described in detail previously (Gothard et al., 1996). Briefly, each microdrive consisted of a drive screw coupled by a nut to a guide cannula. Twelve guide cannulae contained tetrodes (McNaughton et al., 1983; Recce and O'Keefe, 1989), four-channel electrodes constructed by twisting together four strands of insulated 13 μm nichrome wire (H.P. Reid). Two additional tetrodes with their individual wires shorted together served as an indifferent reference and an EEG recording probe. A full turn of the screw advanced the tetrode 318 μm.

For 1 rat, the tetrodes were divided into 2 groups of 7 (“split bundle drive”), permitting recording simultaneously from the septal (3.0 mm posterior, 1.8 mm lateral to bregma) and middle (6 mm posterior, 5.0 mm lateral to bregma) regions of the hippocampus. For the other 2 rats, recordings were made sequentially, first from the intermediate (5.7 mm posterior, 5.0 mm lateral to bregma). After intermediate recordings were completed, the implant was removed. Rats were then implanted with a new array in dorsal/septal (3.0 mm posterior, 1.4 mm lateral to bregma) regions. In all cases, the implant was cemented in place with dental acrylic anchored by dental screws. A ground lead was connected to one of the jeweler's screws placed in the skull. After surgery, rats were orally administered 26 mg of acetaminophen (Children's Tylenol *Elixir; McNeil). They also received 2.7 mg/ml acetaminophen in the drinking water for 1–3 d after surgery and oral ampicillin (Bicillin; Wyeth Laboratories) on a 10 d on/10 d off regimen for the duration of the experiment. Data from these animals have been used in other unrelated analyses that have been published previously (Maurer et al., 2005; Maurer et al., 2006a, 2006b; Maurer et al., 2012), drawing from a database of ∼900 well isolated pyramidal neurons.

##### Electrophysiological recording.

Twelve tetrodes were lowered after surgery into the hippocampus, allowed to stabilize for several days just above the CA1 hippocampal subregion, and then gradually advanced into the CA1 stratum pyramidale. Another probe was used as a neutral reference electrode and was located in or near the corpus callosum. The final probe was used to record theta field activity from the vicinity of the hippocampal fissure. Each tetrode was attached to four separate channels of a 50-channel unity-gain head stage (Neuralynx). A multiwire cable connected the head stage to digitally programmable amplifiers (Neuralynx). The spike signals were amplified by a factor of 1000–5000, band-pass filtered between 600 and 6 kHz, and transmitted to the Cheetah Data Acquisition system (Neuralynx). Signals were digitized at 32 kHz and events that reached a predetermined threshold were recorded for a duration of 1 ms. Spikes were sorted offline on the basis of the amplitude and principal components from the four tetrode channels by means of a semiautomatic clustering algorithm (BBClust; P. Lipa, University of Arizona, Tucson, AZ, and KlustaKwik; K.D. Harris, Rutgers University, Newark, NJ). The resulting classification was corrected and refined manually with custom-written software (MClust; A.D. Redish, University of Minnesota, Minneapolis, MN), resulting in a spike-train time series for each of the well isolated cells. No attempt was made to match cells from one daily session to the next, so the numbers of recorded cells reported does not take into account possible recordings from the same cells on consecutive days. However, because the electrode positions were adjusted from one day to the next, recordings from the same cell over days were probably relatively infrequent. Putative pyramidal neurons were identified by means of the standard parameters of firing rate, burstiness, spike waveform shape characteristics (Ranck, 1973), and the first moment of the autocorrelation (Csicsvari et al., 1998).

Theta activity in the CA1 layer was taken from the tetrode that collected the most pyramidal neurons while the fissure LFP was recorded from a separate probe that was positioned 0.5 mm below the CA1 pyramidal layer. LFP signals were band-pass filtered between 1 and 300 Hz and sampled at 2.4 kHz, amplified on the head stage with unity gain, and then amplified again with variable gain amplifiers (up to 5000). Several light-emitting diodes were mounted on the head stage to allow position tracking. The position of the diode array was detected by a television camera placed directly above the experimental apparatus and recorded with a sampling frequency of 60 Hz. The sampling resolution was such that a pixel was 0.3 cm.

##### Time-series analysis.

Throughout the current study, the phrase “time series” will be used in its traditional sense of a sequence of unprocessed measurements (numbers) indexed in time, synonymous to raw, unfiltered voltage traces of the EEG. This study assumes that EEG time series express essential properties of the underlying physical processes of the brain; that is, the physiological interactions contributing to the extracellular transmembrane currents. In principle, these processes could be represented using a mathematical formalism that uses a set of differential equations to describe their time evolution. The term “system” will therefore be used as it is in “dynamical systems” theory, to denote the differential equations that describe the brain processes and the time evolution of their solutions. Therefore, “system” can be thought of as a mathematical abstraction completely describing the brain processes.

In the mathematical description introduced above, EEGs are functions of, and carry information about, the state of the underlying system. The central focus of this study is the nonlinear character of the system and the central observation that the nonlinearity of the brain is reflected in the EEG time series. If the system is linear, then known solutions can be added to construct new solutions; that is, the solution space is a linear space. Under quite general uniformity conditions, sinusoids form a basis of the solution space, the Fourier basis. The general solution is a superposition of sinusoids with different amplitudes (which include the initial phase) and thus completely defined by the distribution of amplitudes. The stochastic general solution is a superposition of sinusoids with random amplitudes. The joint probability density of the amplitude set defined completely the statistics of the solution. Therefore, the geometry of the solution space of a linear system is trivial: any point in the space spanned by the Fourier basis is a solution to the system and the general solution is completely characterized by the projections on the elements of the basis (amplitudes).

In contrast, in the case of nonlinear systems, two solutions cannot be added to construct a new solution, which makes them rather intractable unless some simplifying assumptions can be made. Our central assumption is that the brain is a weakly nonlinear system and, for weakly nonlinear systems, the general solution can still be represented, in the leading order, as a superposition of sinusoids. However, the amplitudes cannot be constant (otherwise, the system would be linear) and therefore have to evolve. Indeed, decomposition on a linear basis (e.g., the Fourier one) yields a system of equations that describes the evolution of amplitudes through mutual interactions.

The geometry of the nonlinear solution space becomes nontrivial: a general solution is a trajectory in the space spanned by the Fourier basis. The stochastic solution will therefore be characterized by statistically correlated amplitudes and phases. Phase correlations are expressed in the peculiar appearance of the shape of the solution, for example, in the development of time series asymmetries. It follows that such a solution is not completely defined by its power spectrum and knowledge about the phase correlations is essential.

The analysis of the theta spectral band EEG time series used in the current study was based on standard techniques used for stationary signals (Priestley, 1981; Papoulis and Pillai, 2002). We assume that the EEG time series *g*(*t*) is a stochastic process, stationary in the relevant statistics, and decompose it using the discrete Fourier transform (DFT) as follows:
Where *g*_{j} = *g*(*t*_{j}) is a sequence of *N* points collected at times *t*_{j} = *j*Δ*t*, with *j* = 1, 2, …, *N*, Δ*t* is the sampling interval, and *G*_{n} = *G*(*f*_{n}) is the sequence of complex Fourier coefficients corresponding to the following frequencies:
The family of sequences *s*_{n}(*t*_{j}) = *exp* (−2π*f*_{n}*t*_{j}) form the Fourier basis, also called modes. The second-order (linear) statistics of the Fourier spectrum of time series *g* are characterized by the spectral density as follows:
Where *G*_{n} is the series of complex Fourier coefficients of process *g*, *f*_{n} is a frequency band in the Fourier representation (Eq. 2), and **E**[…] is the expected-value operator. Spectral densities describe the frequency distribution of variance. If the process *g* is linear, it is completely characterized by its variance and consequently its variance density *S _{n}*. A realization of the process can be derived by assigning a set of random phases (uniformly distributed in [−π, π]) to modal amplitudes defined e.g., as √(

*S*)/2. Note that the spectral density contains no information about phases, their correlation, and, therefore, about the nonlinearity of the process.

_{n}The nonlinear character of the system, expressed in phase correlations across spectral components, is described in the lowest order by the bispectrum, first proposed for ocean waves by Hasselmann et al. (1963) and further developed by Rosenblatt and Van Ness (1965) and others (Kim and Powers, 1979; Masuda and Kuo, 1981; Elgar, 1987), as follows:
Where *f*_{n} and *f*_{m} are two frequencies in the Fourier sequence in Equation 2. In other words, in this notation, *f*_{n} and *f*_{m} are waves of different frequencies. When the bispectrum is derived from a single time series, the result is described as the auto-bicoherence (simply referred to as bicoherence in the present study). The bispectrum is statistically zero if the Fourier coefficients are mutually independent; that is, for a linear system. For nonlinear systems, the bispectrum will exhibit peaks at triads (*f*_{n}, *f*_{m}, *f*_{n+m}) that are phase correlated, measuring the degree of three-wave coupling. The real and imaginary part of the bispectrum *B* provide measures of third-order statistics (Masuda and Kuo, 1981; Elgar, 1987), as follows:
Where ζ is the skewness and *A* is the asymmetry of the process (Haubrich and MacKenzie, 1965; Masuda and Kuo, 1981) and σ is its SD, with the sign of ζ and *A* (±) indicating the direction of the skew or asymmetry. Both ζ and *A* are “global” measures of the nonlinearity of the process (Fourier modes obviously have zero skewness and asymmetry). Meaningful skewness/asymmetry estimates can be defined for a frequency band only if it is wide enough (see below) to contain all of the relevant phase-coupled modes.

To eliminate the distortion induced by the variance distribution, the bispectrum can be normalized (Haubrich and MacKenzie, 1965) as follows:
The squared modulus and the phase of the normalized bispectrum are called bicoherence and biphase. Although definition 5 does not insure that |*b*_{n,m}| < 1 (see, e.g., Kim and Powers, 1979; Elgar and Guza, 1985), it was adopted here for its simplicity and ease of implementation. Equation 4 implies that the real and imaginary part of the normalized bispectrum can be interpreted as the “frequency distribution” of skewness and asymmetry. In the sequel, we will refer to the real and imaginary parts of the normalized bispectrum as skewness and asymmetry distributions as follows:
Where **R**{*z*} and **I**{*z*} are the real and imaginary parts of the complex number *z*, respectively.

The bispectrum (and thus the normalized bispectrum) has well known symmetries with respect to its arguments (Rosenblatt and Van Ness, 1965; Fig. 1), which is equivalent to saying that different regions in the plane spanned by (*f*_{n}, *f*_{m}) contain equivalent (redundant) bispectral information, as follows:
From rules 8–10, it follows that the first octant of the plane (*f*_{n}, *f*_{m}) contains all nonredundant bispectral information. For the DFT, the nonredundant domain reduces to a segment of the first octant bounded by a line parallel to the second diagonal passing through the maximum (Nyquist) frequency. The nonredundant domain in the (*f*_{n}, *f*_{m}) plane is illustrated in Figure 1. The figure contains a basic, minimal receipt for understanding bispectral distributions. At any point in the (*f*_{n}, *f*_{m}) plane, the value of the bispectrum *B*(*f*_{n}, *f*_{m}) represents the phase correlation between the Fourier modes with frequencies *f*_{n}, *f*_{m}, and *f*_{n+m}= *f*_{n} + *f*_{m},. In other words, the triad contains the two coordinates of the point and the frequency defined by the intersection of the second diagonal with the horizontal axis.

##### Implementation of bispectral analysis.

To estimate second-order Fourier statistics, the time series were de-meaned, linearly de-trended, and divided into 50% overlapping segments of ∼4 s windows (2^{13}-point windows for the sampling rate *f*_{s} = 1988 Hz for rats 7951 and 8042; 214 for *f*_{s} = 936 Hz for rat 7805), with a frequency resolution of ∼Δ*f* = 0.25 Hz. The total time intervals for analysis were chosen so as to yield a number of degrees of freedom (DOF) >170 for all data analyzed and DOF = 300 for the data presented. The statistics of the bispectrum of stationary processes are well understood (Haubrich and MacKenzie, 1965; Elgar and Guza, 1985; Elgar and Sebert, 1989; and many others). Briefly, for the normalization used here (Eq. 6), the probability density function (PDF) is approximated by the noncentral χ^{2} distribution, where |*b*|^{2} is the mean value of the bicoherence, and the parameters *n* and α are given by the following (equations 6–8 in Elgar and Sebert, 1989).
Knowledge of the theoretical PDF allows for estimating confidence limits for any mean value of the bicoherence. Because zero mean that bicoherence is meaningful for distinguishing between linear and nonlinear stochastic processes, an important consequence of the above discussion is that the zero-mean bicoherence is χ^{2} distributed with *n* = 2. From this, a confidence level can be derived. With DOF = 300, the zero-mean bicoherence |*b*| < 0.05 at 90% confidence level and bicoherence |*b*| < 0.1 at 95%. This is consistent with Elgar and Guza's (1985) estimate of |*b*| <

The time–frequency representation (windowed Fourier transform; Mallat, 1999) was used to represent the time evolution of the EEG frequency content (2^{12}-point window at 90% overlap). Bispectral analysis was used to examine the nonlinearity of the EEG time series and was implemented using a lower-frequency resolution and increased DOF by using 210-point windows (0.5 s) segments.

To compare the nonlinear character of theta across hippocampal regions and velocity conditions, we defined a global nonlinearity measure as the square root of the bicoherence integrated over the frequency domain (Sheremet et al., 2002), as follows:
where
The value 0.1 is the 95% confidence level of |*b̄*_{n,m}| = 0, where |*b*_{n,m}| is the normalized bispectrum matrix (Eq. 6); Δ*f* = *b*_{n,m}|, where the coefficient *b*_{n,m}, that is, for the diagonal appearing only four times in the full matrix.

All calculations were coded in MATLAB using its implementation of the DFT. The bispectral estimates were computed using code based on modified functions of the HOSA toolbox (Swami et al., 2000).

##### Spike spectrogram.

To determine the frequency in which neuronal spiking occurred, we implemented spectral analyses on spike trains (Leung and Buzsáki, 1983). First, action potentials were sorted based upon the velocity of the rat. Spikes between either 10–20 cm/s or 60–70 cm/s were analyzed separately, placing the spikes for each velocity segment into 1 ms bins for the duration of the recording. This often resulted in a sparse, but periodic vector of spike counts. The power spectra of the spike trains were then calculated via the Thomson multitaper method.

## Results

### Simplified nonlinear EEG model

Here, we constructed a simplified model that illustrates the bispectral signature of stochastic processes with measurable global skewness and asymmetry. The model consists of two nonlinear time series (Figs. 2, 3) constructed using periodic analytical functions (period 1 s) based on elliptic Jacobi functions (Whitham, 2011) deformed to have controllable asymmetry (Figs. 2*a*, 3*a*). The process analyzed in Figure 2 has negative skewness and zero asymmetry; the process analyzed in Figure 3 has zero skewness and positive asymmetry (Figs. 2*a*, 3*a*). For a realistic aspect, a Gaussian noise process with a spectral law of *f*^{−2} (Figs. 2*b*,*c*, 3*b*,*c*) is added to the analytic functions. The spectra of the time series (both analytic and “noisified”) exhibit peaks at the harmonics of the fundamental frequency of 10 Hz. Note that the DOF values for the model are arbitrary, so the statistics of the processes represented in Figures 2, 3, and 4 are arbitrarily close to the theoretical ones.

Different components of the bispectrum of the noisified time series are shown in Figures 2*d–f*, 3*d–f*, and 4*d–f*. The modulus of the normalized bispectrum (Eq. 6, Fig. 2*d*, 3*d*) exhibits seven, but potentially more, distinct peaks. Based on the diagram in Figure 1, these indicate the frequency triplets that are strongly phase coupled. One can distinguish coupling up to the sixth harmonic of the peak: (*f*_{p}, *f*_{p}, 2*f*_{p}), (*f*_{p}, 2*f*_{p}, 3*f*_{p}), (*f*_{p}, 3*f*_{p}, 4*f*_{p}), (*f*_{p}, 4*f*_{p}, 5*f*_{p}), (2*f*_{p}, 2*f*_{p}, 4*f*_{p}), (2*f*_{p}, 3*f*_{p}, 5*f*_{p}), (2*f*_{p}, 4*f*_{p}, 6*f*_{p}), as well as (3*f*_{p}, 3*f*_{p}, 6*f*_{p}). For the first nonlinear time series (Fig. 2), the skewness distribution γ (real part of the normalized bispectrum, Eq. 6, Figs 2*e*, 3*e*) shows the same peaks, but negative, in accordance with the global skewness of the time series. Because the time series has no asymmetry, the asymmetry distribution α is statistically zero (Fig. 2*f*). Alternatively, for the nonskewed but asymmetric time series (Fig. 3), the skewness distribution is statistically zero (Fig. 3*e*), but the asymmetry distribution is positive (Fig. 3*f*).

The contrast spectral and bispectral analysis and the effectiveness of the bispectral components in detecting nonlinearity is illustrated in Figure 4. Here, a linear time series is constructed by assigning Fourier modes random phases uniformly distributed in [−π,π] to amplitudes derived from the spectral density of the negatively skewed, symmetric time series (Fig. 2*b*), as described above. Importantly, a simple spectral density analysis does not distinguish between the linear and the nonlinear time series because these two oscillations have the same spectral density characteristics (Fig. 4*c*). This illustrates the point that spectral density estimators contain no information about phase correlation across the spectrum, and thus no information about the nonlinear character of the process. Bispectral analysis (cf. Figs. 2*d–f*, 4*d–f*), however, reveals clearly the linear character of the linear time series in Figure 4*b*. The normalized bispectrum and all of its components are statistically zero.

This simple example makes several points: (1) the bispectrum can identify phase coupling, (2) the different components of the bispectrum are useful to identify the effect of the phase coupling on the skewness and asymmetry of the process, and (3) the presence of phase coupling can be detected even in the presence of Gaussian noise and below the noise level. Finally, this analysis suggests a more general meaning for the concept of harmonics. Harmonics can be defined in several ways. A mode *G*_{kn} with frequency *f*_{kn} = *kf*_{n}, with *k* integer, is called the *k*th harmonic of mode *Gn*. This definition is trivial and not useful: because spectral estimates are never exactly zero at any frequency, such harmonics always exist but have no statistical relationship with the mode *G*_{n}. It is tempting to describe spectral distributions such as that in Figure 3 as fundamental frequency and its harmonics.

The inability of the spectral density estimators to distinguish between linear and nonlinear systems makes another important point: the presence of peaks in spectral densities at multiples of a given frequency (harmonics) might have different meanings. For linear systems, harmonics are statistically independent and their presence does not modify the statistics of the process. For nonlinear systems, the presence of harmonics might be accompanied by phase correlations that have significant and fundamental effects on the shape and statistics of the process.

Based on this discussion, hereafter, we reserve the term “harmonic” only for cases in which there is a phase correlation between the harmonics and the fundamental frequency (mode). The determination of whether certain modes are or not harmonics of a spectral peak has to involve an examination of the bispectral characteristics.

### Relation between theta and velocity in EEG measured in the dorsal and intermediate CA1

An examination of LFP time series showed obvious changes in the shape of the hippocampal theta (Fig. 5), in agreement with previous studies (Harper, 1971; Coenen, 1975; Leung et al., 1982; Leung, 1982; Terrazas et al., 2005). The unfiltered LFP trace was well approximated by a near-sinusoidal 6–8 Hz band-pass signal for running speeds <20 cm/s (Fig. 5*a*,*c*). At running speeds >20 cm/s (Fig. 5*b*), the dorsal LFP trace became asymmetric and skewed, departing significantly from the 6–8 Hz signal. Note that the departure of the raw LFP from the 6–8 Hz filter was attenuated in the intermediate hippocampus, even at fast running speeds (Fig. 5*d*). These data suggest the development of cycle skew and asymmetry may be associated with the 16 Hz frequency noted in prior investigations.

### Velocity and the CA1 spectrogram

To examine the effect of running speed on the LFP, the average shape of the power spectrogram as a function of running speed was calculated (Fig. 6). At low speeds, the 4–60 Hz frequency band was dominated by a single spectral peak located approximately between 6 and 7 Hz (Fig. 6*a–d*). As the speed increased, the peak quickly narrowed, became more prominent, and shifted toward 8 Hz while additional peaks developed at frequencies *n f*_{θ} (*n* > 1 integer) where *f*_{θ} = 8 Hz. For simplicity, we will call the *n f*_{θ} peaks putative theta “harmonics” (see simulation above). The dominance of the putative harmonics increased with movement speed at the expense of the bands separating them. This effect was more pronounced in the dorsal region of the hippocampus, where harmonics could be identified up to *n* = 6 (Fig. 6*a*,*b*), whereas the intermediate region developed no more than *n* = 2 (Fig. 6*c*,*d*).

### Phase dependence of spectral peaks

The development of spectral peaks in the 16 Hz range and higher have considerable overlap with beta (10–30 Hz) (Penttonen and Buzsáki, 2003) and low-gamma (25–50 Hz) (Colgin et al., 2009; Belluscio et al., 2012). Therefore, the development of higher frequencies at faster velocities may be due to the appearance of rhythms other than the evolution of theta harmonics. Beta and low-gamma, however, are meaningful as distinct rhythms only if they are statistically independent of the theta rhythm. Therefore, the identity of harmonics can be tested simply by checking their statistical relation to theta. In the lowest order, this can be achieved by using bispectral analysis. Compared with other methods for examining frequency coupling between two different frequencies (Lachaux et al., 1999), this approach has the advantage that it provides direct measures of phase coupling between Fourier frequency bands with no *a priori* (arbitrary) band identification (filtering). Bispectrum components provide complementary information about the stochastic process analyzed. The absolute value of the normalized bispectrum (bicoherence; Fig. 2) is a measure of the nonlinearity of the signal and is statistically zero for independent frequency bands (stochastic process with uncorrelated phases). The real and imaginary parts provide measures of the contribution of different bands to the skewness (Fig. 3) and asymmetry (Fig. 4) of the total signal. For low-velocity conditions (Fig. 7*a*), dorsal CA1 bicoherence estimates (Fig. 7*a*,*d*) were statistically zero over the entire Fourier space (no phase correlation). This implies that the asymmetry and skewness measures were also statistically zero (Fig. 7*b*,*c*). That is, the LFP process was nearly sinusoidal. This agreed with the general aspect of the trace (Fig. 5*a*,*c*). In contrast, bispectral measures estimated for high-velocity conditions exhibited a distinct pattern of peaks, coinciding with the location and development of the harmonics (see Fig. 1 for instructions on how to read the plots). The bicoherence showed significant phase coupling between theta and the harmonics, with correlations distinguishable for up to the sixth harmonic in dorsal CA1. The bicoherence peaks were reproduced in the skewness and asymmetry distributions shown in Fig. 7, *e* and *f*, where the negative sign indicates the type of skewness and asymmetry produced by the peaks. Remarkably, the nonlinearity of theta was modulated across the hippocampal regions. The same analysis performed for intermediate CA1 (Fig. 7*j–r*) exhibited weaker nonlinearity, with coupling only between theta and the second harmonic clearly distinguishable.

To compare the magnitude of the EEG nonlinearity across hippocampal regions and velocity conditions directly, the mean global nonlinearity measure (Eq. 12–14) for each rat in the dorsal and intermediate hippocampus for the low and high velocities was calculated. Figure 8 illustrates the effects of speed on theta nonlinearity for the three rats tested. The nonlinearity measure ϕ is consistently higher for higher speeds, with the nonlinearity of the dorsal CA1 region dominating the nonlinearity of the intermediate CA1 region. This comparison revealed that there was a significant main effect of hippocampal region (*F*_{(1,4)} = 30.6, *p* < 0.01; repeated measures) on the magnitude of the nonlinearity such that the EEG in the dorsal hippocampus was more nonlinear than the intermediate. Although there was no overall main effect of velocity (low vs high) on nonlinearity (*F*_{(1,4)} = 5.9, *p* = 0.08; repeated measures), the interaction between hippocampal region and velocity was significant (*F*_{(1,4)} = 9.1, *p* < 0.05; repeated measures). *Post hoc* analysis indicated that this interaction effect was due to a significant difference in EEG nonlinearity across velocity conditions in the dorsal hippocampus (*p* < 0.05), but not in the intermediate hippocampus (*p* = 0.2). Together, these data indicate that EEG in the dorsal hippocampus is more nonlinear and more sensitive to changes in speed of movement relative to the intermediate hippocampus.

### Neuronal activity and theta nonlinearity

Because theta has been hypothesized to organize spike timing of neurons (O'Keefe and Recce, 1993; Lisman and Idiart, 1995; Skaggs et al., 1996; Harris et al., 2003), we examined action potential frequency as a function of velocity for single cells in dorsal and intermediate CA1 using spectral density estimates (Fig. 9). The spectral power of spiking in the theta band demonstrated that dorsal neurons were sensitive to changes in velocity. There was a significant increase in both frequency (*t*_{(5)} = 7.75, *p* < 0.001) and amplitude (*t*_{(5)} = 3.58, *p* < 0.02) between slow and fast velocities, which is consistent with a previous report (Geisler et al., 2010). This was not the case in the intermediate hippocampus, in which neither the frequency (*t*_{(5)} = 0.48, *p* = 0.66) nor the amplitude (*t*_{(5)} = 1.12, *p* < 0.33) was significantly affected by velocity. Second harmonic peaks were seen to develop for the pyramidal cell action potential spectrograms in the dorsal and intermediate region of the hippocampus, but not for the interneurons. The arrows in Figure 9 indicate a second harmonic peak that is significantly different from 0 (*p* < 0.05 for both comparisons) for the pyramidal cell spike spectrogram in high-velocity conditions only. Given the appearance of harmonics exclusively in principal neurons, it suggests that harmonic activity may be related to rebound activity after inhibition (Cobb et al., 1995; Diba et al., 2014). Specifically, velocity-dependent increases of interneuron activity may invoke rebound excitation, leading to pyramidal cells firing two distinct bursts within a single theta cycle (see Discussion).

## Discussion

Using bispectral analysis, the current study documents the novel findings that (1) the CA1 subregion of the hippocampus has harmonics that are statistically dependent on theta; (2) higher-order harmonics emerge as a function of running speed in dorsal CA1 and, to a lesser extent, in intermediate CA1; and (3) harmonics were the source of the overall deformation of the LFP shape (negative skewness and asymmetry). Therefore, theta harmonics are in fact the spectral signature of the shape change of the LFP and express the increased nonlinearity of the theta rhythm in response to faster running speeds and increased neural activity. The harmonics cannot be regarded as a distinct (statistically independent) set of oscillations and are therefore not related to the beta and low-gamma rhythms.

These results show that the theta rhythm changed in response to running behavior, becoming more nonlinear as the animal's movement velocity increased. In the time domain, this is expressed as the development of skewness and asymmetry in the theta shape; in the spectral domain, it is expressed as the development of a chain of harmonics statistically phase coupled to the theta oscillation. This evolution is reflected in the spiking activity of principal cells (Maurer et al., 2005) and modulated across regions of the hippocampus: it is stronger in the dorsal CA1 than in the intermediate CA1. The current findings add to our understanding of how dynamic activity patterns vary as a function of anatomical position along the longitudinal axis of the hippocampus. Because theta is known to propagate from dorsal to ventral regions (Petsche and Stumpf, 1960; Lubenov and Siapas, 2009; Patel et al., 2012), these data suggest that the nonlinear character of theta decreases as it travels and suggests a hypothesis that activity during movement dissipates along the longitudinal axis of the hippocampus. Future studies will need to test this idea directly.

Although the potential mechanisms that attenuate theta nonlinearity as it travels are not known, the dorsal and ventral regions of the hippocampus are dissociated in terms of patterns of gene expression and the behavioral impact of lesions (Nadel, 1968; Stevens and Cowey, 1973; Sinnamon et al., 1978; Moser et al., 1993; Moser et al., 1995; Long and Kesner, 1996; Moser and Moser, 1998; Kjelstrup et al., 2002; Burton et al., 2009; Wang et al., 2015). In fact, it has been proposed that the hippocampal longitudinal axis is functionally organized along a gradient (Strange et al., 2014; Long et al., 2015) and others have even proposed that the dorsal and ventral hippocampus should be considered distinct structures (Fanselow and Dong, 2010). The observation that the theta rhythm is more nonlinear in the dorsal relative to intermediate hippocampus supports the idea that information may be processed differently along the dorsoventral axis. It also leaves open the question of whether different behavioral parameters could more strongly modulate theta nonlinearity in the intermediate and ventral areas. Specifically, whereas spatial cognitive behaviors are attributed to the dorsal hippocampus, the ventral hippocampus is often associated with emotional processing, stress and affect (for review, see Fanselow and Dong, 2010). This presents the possibility that LFP nonlinearity in more ventral regions of the hippocampus might be influenced by the emotional valence of an experience rather than by the animal's movement.

The development of the theta harmonic with faster running speeds was reflected in the spectrogram of principal cell firing. Although the correlation between theta nonlinearity and neuron spiking is obvious, establishing causality is more difficult. One can make the argument that spiking activity contributes to the background field (Geisler et al., 2010) and, therefore, that altered spiking dynamics with velocity would change the LFP shape. This argument, however, may be circular in that the evolution of theta with velocity alters the spike timing by mechanisms such as ephatic coupling (Buzsáki et al., 1991; Anastassiou et al., 2011). Perhaps more importantly, theta is not carried by the spiking of neurons, but rather is the summed activity of excitatory and inhibitory postsynaptic potentials (for review, see Buzsáki, 2005). Because hippocampal pyramidal neurons project to local interneurons (Csicsvari et al., 1998), which in turn provide dense feedback inhibition to several hundred principal neurons (Sik et al., 1995), local recurrent dynamics have been proposed to support the local generation of theta (Leung, 1998). In the model of Leung (1998), two inputs are responsible for the theta dipole across the CA1 layer: one onto the soma and proximal dendrites of pyramidal neurons carried by rhythmic basket cells and the other from input onto distal apical dendrites. Moreover, the interneurons across these regions exhibit differences in their phase preference (Klausberger et al., 2003; Somogyi and Klausberger, 2005), providing structured IPSPs along different pyramidal neuron domains (Allen and Monyer, 2015). Therefore, the consequence of two inhibitory rhythmic inputs at ∼8 Hz across distinct subcompartments of CA1 pyramidal neurons will interact, providing the basis for the 16 Hz modulation seen in Figure 9. Because there are fewer neurons active in the intermediate region of the hippocampus at any given moment in time (Maurer et al., 2006a) and these cells are less sensitive to changes in velocity (Maurer et al., 2005), fewer excitatory and inhibitory events contribute to the theta rhythm and the harmonic. Nonetheless, the intermediate pyramidal cells are still influenced by the structured inhibition, yielding a harmonic in the firing patterns.

This brings to the forefront the question of the role of the theta oscillation in the architecture of the brain. Although axons tend to propagate in a hierarchical manner (Felleman and Van Essen, 1991), the brain's architecture is more akin to a recurrent network in which neurons generate, as well as add, their own information to the activity propagating through the networks (Buzsáki, 2006). Because the theta rhythm is theorized to play a role in the coordination of brain activity, the nonlinear nature of wave propagation may govern the organization of spike timing. This idea is supported by the brain's capacity for large-scale integration of small components (Steriade, 2001; Buzsáki and Draguhn, 2004) and self-organized dynamics (Ashby, 1947; Kelso, 1997), as well as the simple observations that neurons fire during sleep and quiescent periods, when integration is largely disabled. The self-organized emergence of these oscillatory patterns, which reflects the excitatory (or inhibitory) state of the local groups of neurons (Haken, 1984), allows the firing patterns among assemblies of neurons to be “constrained” or timed via fluxes in the extracellular ionic concentration that modifies ionic driving force (Anastassiou et al., 2011). Although there have been significant theoretical advancements in describing self-organized, nonlinear dynamics in the nervous system (Amari, 1982; Haken, 1984; Freeman, 1992; McKenna et al., 1994; Kelso, 1997), the majority of our knowledge has been gained via linear methods, which ignore important elements of the complexity of the brain (Buzsáki, 2006). Future research attempting to untangle the dynamics of the nervous system should integrate nonlinear methods to understand the interaction between anatomy and activity.

## Footnotes

This work was supported by the McKnight Brain Research Foundation, the University of Florida Research Seed Opportunity Fund, the Claude D. Pepper Older Americans Independence Center, and The National Institutes of Health (Grant R03AG049411). We thank the reviewers for their thoughtful comments in revision.

The authors declare no competing financial interests.

- Correspondence should be addressed to Andrew P. Maurer, Department of Neuroscience, University of Florida College of Medicine, P.O. Box 100244, 1149 Newell Drive, McKnight Brain Institute, L1-100 Gainesville, FL 32610. drewmaurer{at}ufl.edu