# Relationships between electroencephalographic spectral peaks across frequency bands

^{1}Institute of Neuroscience and Medicine (INM-6) and Institute for Advanced Simulation (IAS-6), Jülich Research Centre and Jülich-Aachen Research Alliance, Jülich, Germany^{2}School of Physics, The University of Sydney, Sydney, NSW, Australia^{3}Brain Dynamics Center, Sydney Medical School – Western, University of Sydney, Sydney, NSW, Australia^{4}Center for Integrated Research and Understanding of Sleep, Glebe, NSW, Australia

The degree to which electroencephalographic spectral peaks are independent, and the relationships between their frequencies have been debated. A novel fitting method was used to determine peak parameters in the range 2–35 Hz from a large sample of eyes-closed spectra, and their interrelationships were investigated. Findings were compared with a mean-field model of thalamocortical activity, which predicts near-harmonic relationships between peaks. The subject set consisted of 1424 healthy subjects from the Brain Resource International Database. Peaks in the theta range occurred on average near half the alpha peak frequency, while peaks in the beta range tended to occur near twice and three times the alpha peak frequency on an individual-subject basis. Moreover, for the majority of subjects, alpha peak frequencies were significantly positively correlated with frequencies of peaks in the theta and low and high beta ranges. Such a harmonic progression agrees semiquantitatively with theoretical predictions from the mean-field model. These findings indicate a common or analogous source for different rhythms, and help to define appropriate individual frequency bands for peak identification.

## 1. Introduction

Electroencephalographic (EEG) spectra are often characterized by peaks at various frequencies. Most notable is the alpha peak, which usually lies between 8 and 12 Hz in healthy adult humans. It was the first feature reliably detected in human EEG (Berger, 1933), and has often been subcategorized into variants in different regions of the cortex (Niedermeyer and Lopes da Silva, 2005). Other peaks have been widely noted, including beta peaks typically in the range 13–30 Hz in healthy adults, spatially localized gamma peaks above 30 Hz, the theta peak at 4–8 Hz, and (in sleep) spindle peaks at 11–15 Hz (Niedermeyer and Lopes da Silva, 2005). All these peaks are superposed on broadband activity that falls off with increasing frequency.

In the most common forms of quantitative EEG (qEEG), the frequency spectrum is divided into several bands, and the total absolute or relative power in each band is analyzed. While the analysis of band powers has proved to be useful, it amounts to approximating the rich structure of actual EEG spectra by just a few numbers (see Figure 1). Moreover, the bands used to calculate these powers are almost invariably based on average parameters for normal adult humans. This procedure for instance fails to capture the fact that the alpha peak rises from around 3–5 Hz in newborns (Niedermeyer, 1997; Marshall et al., 2002) to 7–13 Hz in normal adults (Van Albada et al., 2010; Chiang et al., 2011), and individual variability which can take peaks outside the normal ranges. In the present study, we perform EEG spectroscopy of a large sample of healthy individuals, characterizing spectral structure in detail, and allowing for individual variations in frequency bands.

**Figure 1. Example of an EEG spectrum (black line) with its qEEG approximation in terms of band powers, given by the areas of the gray bars**.

Various EEG rhythms have been noted to reflect different states of vigilance or independent aspects of cognitive processing (Niedermeyer and Lopes da Silva, 2005). For example, the alpha peak is most prominent in the eyes-closed condition and is associated with attentional suppression (Snyder and Foxe, 2010), while a spindle peak is associated with non-REM sleep, and theta peaks occur especially during drowsiness (Niedermeyer and Lopes da Silva, 2005). Characterizing the relationships between spectral peaks helps to refine such interpretations and sheds light on the underlying mechanisms.

Multiple suggestions have been made as to why EEG peaks have the observed frequencies:

(i) Based on spectral estimates in rats, it was suggested that successive functional frequency bands increase in center frequency by a factor *e* ≈ 2.718 (Penttonen and Buzsáki, 2003; Buzsáki and Draguhn, 2004). In rat brain slices, oscillations could be induced at relative frequencies corresponding approximately to the golden ratio, suggesting period concatenation as an underlying mechanism (Roopun et al., 2008a,b). Since the golden ratio is close to *e*^{0.5}, the second proposal is related to the first, but implies a denser packing of rhythms across frequency. Both Euler’s number and the golden ratio were proposed by offer a computational advantage by minimizing interference between rhythms (Roopun et al., 2008a,b; Pletzer et al., 2010).

(ii) A second suggestion is that rhythms are produced by groups of neurons with similar characteristic frequencies, which might synchronize and act as “pacemakers.” Despite the existence of neurons with intrinsic oscillation properties, this hypothesis suffers from a number of drawbacks (Nunez and Srinivasan, 1981); for instance, it would require a separate pacemaker to be postulated *ad hoc* for each spectral peak.

(iii) Nunez suggested that global EEG rhythms arise as spatial cortical eigen-modes, yielding a non-harmonic progression of peak frequencies (Nunez and Srinivasan, 1981; Nunez, 1995). One prediction of this hypothesis is that alpha frequency should be negatively related to head size, which was found by Nunez (1978) and Posthuma et al. (2001) but was recently challenged (Valdés-Hernández et al., 2010).

(iv) Several other models have considered purely cortical oscillations (Van Rotterdam and Lopes da Silva, 1982; Liley et al., 1999, 2002; Wright, 1999; Jirsa et al., 2002; David and Friston, 2003). For instance, networks of simulated multicompartmental cortical neurons can produce oscillations in the range 8–20 Hz (Liley et al., 1999), and in a non-linear continuum theory, peaks at various frequencies in the range 2–16 Hz were obtained depending on the parameters (Liley et al., 2002).

(v) Considerations of the importance of the thalamus in synchronized oscillations in both sleeping and waking states (Lopes da Silva et al., 1973, 1980; Steriade et al., 1993, 1996; Steriade, 2000) have motivated thalamocortical models (Lumer et al., 1997; Robinson et al., 2001b, 2002; Rennie and Robinson, 2002; Hill and Tononi, 2005; Izhikevich and Edelman, 2008). The proposed models display resonances in various ranges: Lumer et al. (1997) found mostly gamma oscillations with precise frequencies depending on the parameters, Izhikevich and Edelman (2008) found oscillations in the delta and alpha ranges, and the model of Hill and Tononi (2005) exhibited slow waves in sleep and gamma oscillations in activated states. The neural field models of Rennie and Robinson (2002) and Robinson et al. (2001b, 2002), which are further explored here, provide a unified mechanism for slow-wave and spindle oscillations in sleep, and alpha, beta, and higher-frequency oscillations in the waking state. These models predict clear relationships between peak frequencies and amplitudes, with the theta peak occurring at approximately half the alpha frequency on an individual-subject basis, and alpha and beta peaks forming part of a near-harmonic progression.

The latter prediction is consistent with a number of previous studies: Carlqvist et al. (2005) found clear frequency, power, and phase relationships between alpha and beta activity in the resting EEG. The average ratio between beta and alpha peak frequencies was 1.9–2.0, consistent with the beta peak being generated as a harmonic of alpha. Similarly, bispectral analysis of subjects with high alpha activity revealed significant phase and amplitude relationships between alpha and its second harmonic (Barnett et al., 1971). In addition, Barnett et al. (1971) observed that 10 Hz activity was significantly phase-related to third and fourth harmonics at 30 and 40 Hz in some cases, and less prominently to activity at 2 and 7 Hz. Palva et al. (2005) reported cross-frequency phase synchrony between alpha, beta, and gamma oscillations in the human MEG. Finally, some studies have revealed similarities in the scalp topographies and functional characteristics of alpha and beta activity (Chen et al., 2008; Shackman et al., 2010). The present study extends these findings using EEG spectroscopy of a large sample of healthy individuals.

Besides frequencies, we also examine the amplitudes of spectral peaks. These can provide additional evidence for the independence or interdependence of rhythms and allow the thalamocortical mean-field model to be tested further. This model has already been shown to be able to account for various aspects of evoked response potentials (Rennie and Robinson, 2002), onset and dynamics of epileptic seizures (Robinson et al., 2002), and correlation and coherence of EEG and electrocorticographic signals (Robinson, 2003). An extension of this model incorporating the basal ganglia successfully mimicked a number of electrophysiological changes in Parkinson’s disease (Van Albada and Robinson, 2009; Van Albada et al., 2009). Correspondence of amplitude relationships with model predictions would constitute additional evidence for its plausibility.

We perform the analyses partly in the light of aforementioned model of thalamocortical activity, but in a way that would allow the model to be invalidated by the data. The model is fitted to eyes-closed spectra of a large group of healthy subjects, and the model parameters are used to estimate a background spectrum without peaks or troughs. This method balances the dual goals of determining a physiologically realistic background, and not making any prior assumptions about relationships between spectral peaks. Frequencies and amplitudes are then estimated of the empirically measured peaks relative to this background, and their interrelationships are explored.

## 2. Materials and Methods

In this section we describe our data collection, peak fitting, and statistical methods. Section 2.1 describes the subject group, EEG recording procedures, and calculation of spectra. Section 2.2 gives a brief account of the model of thalamocortical activity and its predictions concerning relationships between spectral peaks. Sections 2.3 and 2.4 respectively detail the methods for peak fitting and classification.

### 2.1. Subjects and Recordings

The data were eyes-closed resting EEG spectra of 1424 healthy subjects (702 females and 722 males), a subset (95%) of those in Van Albada et al. (2010) and Chiang et al. (2011), where any subjects rejected in that study based on excessive voltage fluctuations at 14 or more electrodes were also excluded here, resulting in the removal of 39 subjects of the original 1463. Subjects’ ages ranged from 6.08 to 86.55 years (mean 26.88 years). The recordings were obtained with a NuAmps amplifier (Neuroscan) by Brain Resource, Ltd. (www.brainresource.com) and made available through the Brain Resource International Database (BRID; Gordon et al., 2005). The montage included 26 electrodes placed according to an extended International 10–20 system (Klem et al., 1999). Of these, we focus on the Cz electrode, which is relatively unaffected by muscle artifact and combines frontal and occipital influences. The sampling rate was 500 Hz and average of mastoids was used as a reference. An analog low-pass filter removed 40 dB per decade above 100 Hz. Data were corrected offline for eye movements using a method based on that of Gratton et al. (1983). The spectrum was calculated from 2 min of relatively artifact-free EEG with a resolution of 0.25 Hz by averaging the spectra of 50% overlapping 4 s epochs after multiplying each epoch’s time series by a Welch window. We compared our findings with 981 spectra that were identical except for the use of a Hann instead of a Welch window, to exclude the possibility of results depending on the particular choice of windowing function.

### 2.2. Thalamocortical Model

Background spectra and predictions of peak frequencies and amplitudes were calculated using a mean-field model of thalamocortical electrical activity (Robinson et al., 2001b, 2002, 2003a, 2005). It is beyond the scope of this paper to give a detailed mathematical account of the model, but we introduce some aspects here to clarify theoretical predictions of peak frequencies and amplitudes. Section 2.2.1 gives a brief overview of the model, Section 2.2.2 provides approximate frequencies of corticothalamic resonances, and Section 2.2.3 discusses qualitative predictions on relationships between peak amplitudes. For a more detailed treatment we refer the reader to the papers cited.

#### 2.2.1. Overview of the model

The structure of the model is illustrated in Figure 2. We here consider only the version obtained by linearizing about its fixed-point firing rates. The neural populations included are cortical excitatory (*e*), cortical inhibitory (*i*), thalamic reticular (*r*), and thalamic relay (*s*) including both primary relay and association nuclei. Each population is described by its instantaneous mean firing rate. The *e* and *i* populations connect both to themselves with gains *G _{ee}* and

*G*, and to each other with gains

_{ii}*G*and

_{ie}*G*, quantifying the change in output rate divided by the change in input rate. Similarly, the relay nuclei project to the cortical populations with gains

_{ei}*G*and

_{es}*G*. In each case, the second subscript corresponds to the sending population and the first subscript to the receiving population. Approximating connections in the cortex as random leads to

_{is}*G*=

_{ii}*G*,

_{ei}*G*=

_{ie}*G*, and

_{ee}*G*=

_{is}*G*(Braitenberg and Schüz, 1998; Robinson et al., 2001b). Besides cortical interactions, the following loops involving the thalamus are seen: a direct corticothalamic loop passing only through the relay nuclei; an indirect corticothalamic loop also passing through the reticular nucleus; and an intrathalamic loop that involves reciprocal connections between the relay and reticular nuclei. These loops are associated with gains

_{es}*G*=

_{ese}*G*

_{es}G_{se},

*G*=

_{esre}*G*, and

_{es}G_{sr}G_{re}*G*=

_{srs}*G*, respectively.

_{sr}G_{rs}**Figure 2. Schematic representation of the model, including the following populations: e, cortical excitatory neurons; i, cortical inhibitory neurons; r, thalamic reticular nucleus; s, primary and secondary thalamic relay nuclei**. In the linearized version of the model, each connection has a gain

*G*(

_{ab}*a*,

*b*=

*e*,

*i*,

*r*,

*s*). The relay nuclei receive input from the brainstem, indicated by the subscript

*n*.

Spectra can be computed from the model by approximating brainstem input as white noise, and assuming that EEG signals are proportional to the activities of the cortical excitatory neurons (Robinson et al., 1997, 2001a, 2005). Such model spectra were fitted to empirical ones using a fitting procedure that uses a Monte Carlo method with repeated random initializations to avoid finding false minimums (Robinson and Rennie, 2010; Rowe et al., 2004). The quantity minimized was a weighted sum of squared differences between log empirical and log predicted spectra at each frequency. The free parameters were a synaptodendritic time constant *α*, a cortical damping rate *γ*, the corticothalamic axonal latency *t*_{0}, an overall scale factor *p*_{0}, and the gains *G _{ee}*,

*G*,

_{ei}*G*,

_{ese}*G*, and

_{esre}*G*. For further details we refer to the papers cited.

_{srs}Model spectra consist of a background modulated by thalamocortical interactions yielding peaks and troughs. The background is calculated by retaining projections from thalamus to cortex, but setting the strengths of projections from cortex to thalamus to zero.

#### 2.2.2. Frequency estimates via approximation of the dispersion relation

The thalamocortical model uses a damped-wave equation to describe the propagation of neural activity across the cortical sheet (Robinson et al., 2001b). By Fourier transforming the spatiotemporal model equations, an expression for the activity of the cortical excitatory neurons can be obtained in terms of frequencies and wavenumbers. Equating the denominator of this expression to zero yields a dispersion relation, determining the characteristics of the damped waves making up the activity.

In this study we estimated peak frequencies for model spectra in two ways: the first is based on approximations of the dispersion relation for the linearized model, and the second refines these estimates by looking for peaks close to these approximations in background-subtracted model spectra. In the present section, we focus on the approximate frequencies, while the peaks in fitted spectra are described in Section 3.1. Results of these two methods are illustrated in Figure 3.

**Figure 3. Initial approximations (dashed vertical lines) and precise estimates (solid vertical lines) of peak frequencies in model spectra for four subjects with different values of their corticothalamic loop gains G_{ese} and G_{esre}**. For subjects with

*G*> 0 >

_{ese}*G*and |

_{esre}*G*| < |

_{ese}*G*|, we determined peaks in the theta and what we term “iota” and “xi” ranges. For all other subjects, we determined alpha and what we term “beta

_{esre}_{1}” and “beta

_{2}” peaks. For our definitions of iota, xi, beta

_{1}, and beta

_{2}, see Figure 7 and Table 1. Empirical spectra are shown in gray, model spectra in black.

In general, the dispersion relation has complex angular frequencies *ω* = *ω _{x}* +

*iω*as solutions, where

_{y}*ω*determines the oscillation frequency of the solution, and

_{x}*ω*its temporal damping rate. There are no relevant solutions with

_{y}*ω*= 0, since instabilities set in at boundaries where the dispersion relation has real solutions. Spectral peaks for real frequencies

_{y}*ω*=

*ω*occur when the dispersion relation is closest to having a zero. Since uniform modes turn out to be the least damped (Robinson et al., 1997, 1998), we consider only the dispersion relation for zero wavenumber:

_{r}where *t*_{0} is the thalamocortical axonal loop delay, γ is a damping rate for cortical activity propagation, and *L*(ω) accounts for low-pass filtering of signals in synapses and the dendritic tree,

Here, β and α are synaptodendritic rise and decay rates, respectively.

To simplify equation (1) we use the approximations (Roberts and Robinson, 2008)

valid for *ω* < *α*, *β*, *γ*. Dividing equation (1) by the first term, we obtain

where factors of *L*(*ω*)were retained only in the numerators, and the term involving *G _{ee}* was dropped, since it was previously found by numerical exploration to be of secondary importance for peak locations (Robinson et al., 2001b). Peaks occur where equation (5) is closest to being solved, frequencies depending on the relative strengths and signs of direct and indirect thalamocortical feedback. Generally,

*G*> 0,

_{ese}*G*< 0, and |

_{esre}*G*| > |

_{ese}*G*| reflecting the waking state with positive overall thalamocortical interactions and a relatively inactive thalamic reticular nucleus (Robinson et al., 2002, 2004; Van Albada et al., 2010). Minimums of the left-hand side of equation (5) then occur approximately where the complex argument of the

_{esre}*G*term is 2

_{ese}*πn*(

*n*= 1, 2, …). The strongest resonance or putative alpha rhythm corresponds to

*n*= 1, leading to the frequency estimate

Peaks in the low and high beta ranges, which we will term beta_{1} and beta_{2} peaks, correspond to *n* = 2 and *n* = 3, and are located around 2 and 3 times the alpha peak frequency. Due to the approximations made, equation (6) tends to underestimate peak frequencies; more precise estimates are made in Section 3.1.

In some cases, *G _{ese}* > 0,

*G*< 0, and |

_{esre}*G*| > |

_{esre}*G*|, so that thalamocortical resonances arise in an overall negative feedback loop. Peaks then occur where the argument of the

_{ese}*G*term is

_{esre}*π*+ 2

*πn*(

*n*= 0, 1, …). The first of these resonances is a putative theta rhythm with frequency (Robinson et al., 2002)

Note that this is close to half the alpha peak frequency in equation (6) if *t*_{0} + 2/α + 2/β + 2/γ ≫ 1/α + 1/β. In our sample, this was generally the case, the difference between the estimated theta frequency and half the estimated alpha frequency being of order 10%. Higher-order peaks are expected for *n* = 1, 2 with frequencies around 3 and 5 times *f _{θ}*, respectively.

Since no hard limit was imposed during fitting to force *G _{esre}* to be negative, there were also some cases with

*G*> 0,

_{ese}*G*> 0. For

_{esre}*G*>

_{ese}*G*, the frequency becomes

_{esre}in analogy with the previous derivations. Higher-order peaks are expected around integer multiples of this frequency.

We used the estimates equations (6–8) to label peaks in fitted model spectra as theta, alpha, etc., and to obtain more precise predictions of relationships between peak frequencies. For spectra with a theta peak, higher-order peaks are expected to lie between alpha and beta_{1}, and between beta_{1} and beta_{2}. Following the tradition of denoting EEG rhythms by Greek letters, we refer to these rhythms as iota and xi. The definitions of these bands are illustrated in Figure 7, where different sets of band limits were used depending on the location of the highest peak, as described in Section 2.4.1.

#### 2.2.3. Qualitative predictions of amplitude relationships

The thalamocortical model also predicts the amplitudes of the various peaks to covary. We here provide qualitative predictions of such relationships, while quantitative estimates are obtained from fitted model spectra in Section 3.1.

Since beta peaks arise as near-harmonics of alpha peaks in the model, the prediction of a positive association between their amplitudes is straightforward. Predicting the relationship between theta and alpha peaks is more complicated. Simultaneous theta and alpha peaks in empirical spectra could be due to activity in parallel thalamocortical pathways with different gains, or to temporal changes in gain in a single pathway. For instance, positive net feedback may exist in some regions, with negative feedback in others, especially in the drowsy state near the sleep-wake transition, thereby allowing theta and alpha peaks to coexist. Concurrent peaks in what are traditionally considered the theta and alpha ranges could also arise due to parallel thalamocortical loops with different delays, or due to spatial variations in loop delays (Robinson et al., 2001a, 2003b). The version of the model considered here does not account for concurrent theta and alpha peaks via these mechanisms, due to static gains and the lumping of possible parallel or spatially varying thalamocortical loops into a single loop. However, empirical theta peaks can be considered to be superposed on the model background and on troughs which the thalamocortical model also predicts in this range, as described in the next section.

Correlations between theta and alpha peak amplitudes are expected to have contributions from opposing mechanisms. In our model, positive and negative *G _{ese}* +

*G*generally lead to alpha and theta peaks, respectively, and their amplitudes tend to be large when |

_{esre}*G*+

_{ese}*G*| is large. The common dependence of

_{esre}*G*and

_{ese}*G*on the thalamocortical gain

_{esre}*G*will contribute positively to the correlation between empirical theta and alpha peak amplitudes. If concurrent peaks in what are traditionally labeled the theta and alpha ranges arise due to spatial variations in thalamocortical loop delays, a positive association between their amplitudes is also expected. Note that such peaks should actually be labeled by their generating mechanisms rather than by frequency ranges; however, this is difficult to do in practice based directly on empirical spectra.

_{es}In the following, we require that the frequency of theta peaks differ by more than 3 Hz from that of the alpha peak. Since it is possible in principle to have split alpha peaks with a larger frequency difference, alpha splitting may provide a small positive contribution to the relationship between empirical peaks in the nominal theta and alpha ranges determined in this paper.

On the other hand, substantial spatial or temporal variations in *G _{ese}* +

*G*are required to produce large alpha peaks at one time or location, and large theta peaks at another. Assuming that small variations are more likely, especially within the limited time window from which spectra are computed, this will provide a negative contribution to the association between alpha and theta peak amplitudes.

_{esre}### 2.3. Fitting of Gaussian Peaks

The fitting routine is illustrated in the flowchart in Figure 4. It differed in several respects from one previously used to identify alpha peaks in the subjects considered here plus 32 additional subjects (Chiang et al., 2008, 2011). First, the current fitting routine covers not just the alpha frequency range but the larger range 2–35 Hz. Another notable difference is that Chiang et al. (2008, 2011) considered spectra at multiple electrode sites to find clusters of alpha peaks with similar frequencies. Furthermore, in those papers, peaks were fitted with Gaussian functions of log power vs. *f*, whereas we use Gaussian functions of log power vs. log *f*. However, we compared our results with fits of log power vs. *f*, finding no strong differences. The methods also differ in the type of background used: the previous papers considered peaks superposed on a power-law background, while the current paper examines peaks and troughs that modulate a model-based background. This is done to assess spectral features due to thalamocortical interactions (cf., Section 2.2). The model-fitting routine has been validated and its properties analyzed in a number of publications (Rowe et al., 2004; Van Albada et al., 2007, 2010). We do not analyze troughs further in this paper, but including them in the fitting routine enables their future analysis and is relevant for the theta range, as further explained in the following. Finally, the papers cited used a single degree of spectral smoothing, whereas we compared moving averages with different ranges and selected the closest fit, which helped ensure that no peaks were missed, and yielded close agreement with visually identified peaks. The fitting method was developed independently of the results on relationships between peaks described here.

**Figure 4. Flowchart showing the steps involved in fitting peaks and troughs to empirical spectra**. The different degrees of smoothing refer to the step where extremums adjacent to the current peak or trough are found. A low degree of smoothing tends to yield narrow peaks/troughs, whereas a high degree of smoothing tends to yield broad ones. Different degrees of smoothing may be appropriate for different levels of noise. Such smoothing improved the agreement between fitted and visually identified peaks.

Before fitting, spectra were smoothed using a five-point moving average to reduce noise. Up to 12 peaks and troughs were then fitted to the difference of log spectra and log background in the range 2–35 Hz. This number of peaks was chosen since it was seen to adequately capture all visually identified peaks in the range considered. The theta range was fitted first, since background-subtracted empirical spectra suggested that overlapping peaks and troughs were present in this range, and therefore an adjusted method was used for theta. Peaks were first sought in that part of the range 2–9 Hz where the spectrum was below the background, corresponding naively to the theta range. Additional smoothing was then applied until at most a single peak was present in the range. If a peak was present and the distance between its adjacent minimums was ≥1 Hz (to avoid spurious sharp peaks), recursive fitting was performed of the overlapping peak and trough, as illustrated in Figure 5. This entailed the following steps: first, the peak was fitted on a possible constant baseline, and the fitted values were subtracted. Then the trough was fitted with zero baseline, the residual was calculated with only the trough subtracted, and the peak was again fitted. The latter sequence was repeated up to 10 times as long as this decreased the residual computed by subtracting both peak and trough. A further constraint was that the fitted trough was not more negative on average than the empirical one in the first and last quarters of its frequency range. This ensured that recursive fitting did not lead to very large peaks and troughs where these were not present in empirical spectra. Two examples of spectra with overlapping theta peak and trough are given in Figure 6.

**Figure 5. Illustration of recursive fitting of overlapping peak and trough**. **(A)** Shows a Gaussian peak at *f* = 5 and trough at *f* = 5.2 (dashed), and their sum in gray. A small peak remains visible over the range shown in **(B)**. Fitting this peak leaves the residual shown in gray in **(C)**. The dash-dotted line indicates the trough fitted to this residual. Subtracting this trough yields the residual shown in gray in **(D)**. After two more steps, shown in **(E)** and **(F)**, it is seen that the fitted peak and trough closely match the actual ones.

**Figure 6. Example spectra from two subjects with overlapping theta peak and trough**. The residuals and fitted peaks and troughs (in log power) are exponentiated for visual clarity on the logarithmic scale.

The remaining peaks and troughs were fitted in order of decreasing amplitude. Peaks were sought in those ranges where no peaks or troughs had yet been fitted, except for theta, where initially only the range of the peak (the closed range of frequencies between its adjacent minimums) was excluded to allow possible additional peaks in this range to be found. If a frequency range was not bounded by already-fitted peaks or troughs, it was bounded by the closest frequencies where the residual was of opposite sign to that of the extremum. To avoid fitting spurious narrow peaks, only those peaks or troughs were considered that extended over at least max(1, *f*/16) Hz, where *f* is the frequency of the extremum in Hz. Locations of extremums were identified using eight different degrees of smoothing (but the same degree for each peak/trough), and the fit with the lowest absolute residual was selected. The limit to the number of peaks and troughs fitted prevented under-smoothing, resulting in approximately uniformly distributed degrees of smoothing. Gaussian peak or trough values were subtracted over the range where they were ≥0.05 and ≤−0.05, respectively. This fitting algorithm yielded good agreement with visually identified peaks and troughs.

### 2.4. Peak Classification

Frequency bands for the analysis of peak parameters are defined in Section 2.4.1. Peak classification took into account putative split alpha and beta peaks, as explained in Section 2.4.2, but the detailed analysis of split peaks is left to future work considering multiple electrodes. Peak exclusion criteria are described in Section 2.4.3. These take into account the statistical nature of EEG spectra, eliminating peaks that may have occurred by chance.

#### 2.4.1. Band limits

We defined band limits based on the location of the largest peak in the range 2–13 Hz. To prevent influencing correlations by the choice of band limits, we assigned subjects to five groups with appropriate band definitions, and analyzed correlations for each group separately. Figure 7 gives example spectra with fits from each group and illustrates the corresponding bands. If the largest peak was in the range 2–5 Hz (Group 1, *N* = 62) it was considered to be a theta peak, and if it was in the range 5–13 Hz it was treated as an alpha peak. A further subdivision was made based on alpha peak frequency: 5–7 Hz (Group 2, *N* = 49), 7–9 Hz (Group 3, *N* = 461), 9–11 Hz (Group 4, *N* = 797), and 11–13 Hz (Group 5, N = 55). Symmetric bands were defined around these 2-Hz ranges, bandwidth increasing with alpha peak frequency (by 1 Hz for each consecutive group) to maximize coverage of the frequency space. Bands were then defined via the linear regression equations for peak frequencies derived from fitted model spectra (cf., Figure 8). For Group 1, iota and xi limits were calculated from theta limits, while for Groups 2–5, beta_{1} and beta_{2} limits were calculated based on the alpha band. As a result, the iota and xi bands were relatively wide for subjects whose main peak was a theta peak, while the beta_{1} and beta_{2} bands were relatively wide for subjects whose main peak was an alpha peak. This helped ensure that no peaks were missed in the relevant bands. The resulting bands are listed in Table 1. Correlations between peak parameters were determined using only the largest peak in each band.

**Figure 7. Example spectra of subjects from each group having a different set of band limits, as listed in Table 1**. Vertical lines indicate peak locations: magenta, primary peaks; blue, secondary alpha peak; green, secondary beta_{1} peak.

**Figure 8. Relationships between frequencies of peaks or spectral enhancements determined from fitted model spectra**. **(A)** Beta_{1} vs. alpha; **(B)** beta_{2} vs. alpha; **(C)** iota vs. theta; **(D)** xi vs. theta. Thick lines, fits of frequency ratios; thin lines, linear fits with intercept. Dashed lines represent 99% non-simultaneous confidence intervals for the linear trends, and the corresponding 99% confidence bounds for the slopes and intercepts are indicated. Note that beta_{2} frequencies can exceed 35 Hz (the maximum frequency of fitted Gaussian peaks), since model spectra were evaluated up to 50 Hz.

#### 2.4.2. Split peaks

Peaks in the range 5–13 Hz, differing from the primary alpha peak by no more than 3 Hz and less than a factor of two in height, were considered to be secondary alpha peaks. If there were several peaks fulfilling these criteria, the one with the smallest frequency difference with the primary peak was chosen. If such a peak was the highest in the theta or iota band, the next-highest peak in the relevant band was taken to be the primary theta or iota peak, if present.

Secondary beta_{1} peaks were considered to be those peaks lying within 6 Hz of the primary beta_{1} peak, at higher frequency than the highest-frequency alpha peak and not directly flanking it, and differing by less than a factor of two in height from the primary peak. As for alpha, if several such peaks were present, the one closest to the primary peak was selected. If the secondary beta_{1} peak fell outside the beta_{1} band, the next-highest peak in the relevant band was considered to be the primary peak for that band.

This classification of split peaks may be refined and further analyzed in future studies using data from multiple electrodes, as done for alpha peaks by Chiang et al. (2008, 2011).

#### 2.4.3. Rejection criteria

Peaks in the theta or iota bands were rejected if they immediately flanked alpha peaks (using the criterion that their frequency ranges had an overlap of at least two points, corresponding to a range of 0.25 Hz) and were more than four times smaller than the alpha peak, since such peaks usually appeared to result from non-Gaussianity of the alpha peak. The rejected peak was replaced by the next-highest peak in the same range, if present.

Another criterion for peak identification was a good signal-to-noise ratio. At each frequency, a nine-point root mean square deviation between log raw and log smoothed spectra was determined as an estimate of noise. The 10% of peaks with the lowest ratio of height to this RMS deviation at the nearest frequency point were rejected. It is of course possible that some spurious peaks were nevertheless fitted, but these will be randomly scattered and not influence the main trends.

In some cases model spectra did not closely fit empirical spectra. After classifying peaks into bands, we therefore determined the mean deviation between log empirical and log model spectra, and excluded peaks when the model fit was among the worst 15% for the given group and band. Visual inspection showed this to be a relatively conservative exclusion criterion, so that only peaks were considered where the model fit well. Mean deviations rather than mean absolute deviations were used because the reliability of the background depends mainly on whether the model fit is systematically above or below the empirical spectrum. This criterion was not applied to the theta band, since model fits did not yet adequately capture theta peaks. Instead, theta peaks were rejected if the fit in the alpha band was among the worst 15%. Note that empirical theta peaks could nevertheless be investigated, since the model background was fitted in this range, and Gaussian theta peaks and troughs were fitted on top of this background, as explained in Section 2.3.

The rejection criteria were chosen to obtain a maximal set of fitted peaks showing good correspondence with visually identified peaks. Thus, the criteria were independent of the results reported here. The qualitative results were robust to variations of rejection levels.

## 3. Results

Section 3.1 concerns relationships between peak frequencies and heights found from fitted model spectra. These may be regarded as theoretical predictions using physiologically realistic parameters, and as such are intermediate between theoretical predictions and empirical results. The limited number of model parameters prevents overfitting and ensures that relationships between model peaks do not simply reflect empirical ones. The results for Gaussian peaks fitted to empirical spectra are discussed in Section 3.2.

### 3.1. Peak Relationships Based on Fitted Model Spectra

The following two sections respectively describe the frequency and amplitude relationships of peaks in fitted model spectra.

#### 3.1.1. Frequency relationships

Figure 8 shows the dependences of beta_{1} and beta_{2} frequencies on alpha frequency, and of iota and xi on theta frequency, where peaks were labeled as described in Section 2.2. Spectra (including background) with *G _{ese}* +

*G*< 0 often showed a theta enhancement but no actual theta peak. Therefore, theta frequencies were determined from sign changes of the second derivative of the spectrum with respect to frequency. We excluded those cases from analysis where the spectrum was below the background at the theta frequency thus determined (9 out of 64). It is seen that theta peaks or shoulders in model spectra tend to occur much below half the normal alpha peak frequency. This may be an artifact due to the fact that the fitted version of the model has only a single set of gains and therefore does not account for concurrent theta and alpha peaks. Thus, cases with

_{esre}*G*+

_{ese}*G*< 0 have theta and iota peaks in fitted spectra, with a frequency ratio close to three. Since the fitting routine emphasizes the goodness of fit for the peak around 10 Hz, theta peaks are fitted around 3 Hz even when they empirically occur around 5 Hz. In contrast, the analysis in Section 2.2 showed that variations in thalamocortical gains will tend to cause theta and alpha peaks with a frequency ratio close to two. We will take this as our prediction of the relationship between theta and alpha peak frequencies.

_{esre}The mean ratios of beta_{1} and beta_{2} peak frequencies to the alpha peak frequency were 2.123 ± 0.008 and 3.32 ± 0.01, respectively, slightly above the ratios of 2 and 3 predicted based on the approximations in Section 2.2. Due to the non-zero intercepts of the linear trends, these ratios depend somewhat on the constitution of the sample. For instance, for the 50% of subjects with the lowest alpha peak frequencies, the mean ratios were closer to 2.2 and 3.4. The iota-to-theta frequency ratio was 3.4 ± 0.2, again somewhat above the approximate theoretical prediction of 3. Similarly, the xi-to-theta ratio was 6.3 ± 0.3, to be compared with the value of 5 based on the simplified equations in Section 2.2. We note that the reliability of the latter estimates is somewhat compromised by the possible fitting of iota peaks to actual alpha peaks, as mentioned above, which may have affected the parameter values and hence relative peak locations.

#### 3.1.2. Amplitude relationships

Relationships between alpha and beta peak heights are illustrated in Figure 9. The regressions were performed without an intercept term, since no beta peaks arise in the model without alpha peaks. The heights are more scattered than the frequencies, but clear positive trends remain. The slopes of the trend lines are slightly reduced by the few spectra with particularly strong alpha, even though the two spectra with the highest alpha peaks in fitted spectra were excluded, since fitted peaks did not accurately reflect empirical ones in these cases. For instance, excluding all cases with alpha height >3, the slopes become 0.278 and 0.079 for beta_{1} and beta_{2}, respectively.

**Figure 9. Relationships between peaks heights from fitted model spectra**. **(A)** Beta_{1} vs. alpha; **(B)** Beta_{2} vs. alpha. Dashed lines and the text indicate 99% confidence intervals for the fits.

No significant correlations were found between theta peak heights on the one hand and iota and xi peak heights on the other hand (*p* > 0.5). However, iota and xi peak heights did have a positive association (*r* = 0.63, *p* = 9.9*e* − 8), with the slope of the regression line for xi vs. iota height being 0.4 ± 0.1. A similar word of caution applies to these amplitude relationships as to the corresponding frequency relationships, since the model parameters for subjects with *G _{ese}* +

*G*< 0 may be affected by the simultaneous presence of theta and alpha peaks in empirical spectra.

_{esre}### 3.2. Empirical Peak Relationships

Here we respectively present the empirical findings on frequency and amplitude relationships between spectral peaks, and compare these with the model predictions.

#### 3.2.1. Frequency relationships

Figure 10 shows the average empirical spectrum and average fitted Gaussian peaks of log power vs. log frequency plotted against *f*/*f _{α}*, where

*f*is the individual alpha frequency; this permits frequency ratios to be explored. Averages consisted of mean spline-interpolated spectra across those subjects for which an alpha peak was fitted and not rejected based on the model fit.

_{α}**Figure 10. (A,B)** Mean empirical spectrum (black) and fitted peaks (gray) vs. *f*/*f _{α}*.

**(A)**All subjects for which alpha peak parameters were obtained.

**(B)**The 10% of subjects with the highest alpha peak amplitudes. Fitted peaks are scaled for clarity, with the same scale in

**(A,B)**.

**(C)**Ratios of all pairs of peak frequencies within spectra. Dotted lines are drawn at 1/3, 1/2, 2/3, 3/2, 2, and 3.

**(D)**Pairs of peak frequencies within spectra. The grayscales indicate the density of points. The empty diagonal band reflects the omission of the 1:1 points and the separation necessary for peak resolution.

It is clearly seen that beta_{1} peaks occurred on average close to twice the alpha peak frequency, while theta peaks occurred around half the alpha peak frequency. Third harmonics of alpha may have been too small or scattered to be visible in the overall average. Since we expect large alpha peaks to be concurrent with large beta peaks, we also plotted the average for the 10% of subjects with the largest alpha peaks separately. This average does seem to have a shoulder around three times the alpha peak frequency (Figure 10B). The average fits additionally show small peaks around 1.5 times the alpha peak frequency, which are however not clearly apparent in the mean empirical spectra. This effect might be explained by the presence of superposed positive and negative power modulations.

Figures 10C,D provide a visualization of these findings that avoids labeling peaks as “alpha” or otherwise, and does not depend on band limits. Figure 10C shows the frequencies of all peaks not rejected based on signal-to-noise ratio vs. their ratios to all other peaks in the same spectrum. Figures 10C,D confirm the association of peaks around 10 Hz with peaks at half, twice and three times that frequency. In particular, the horizontal stripes around (20, 1/2) and (30, 1/3) in Figure 10C clearly show the presence of second and third harmonics of alpha. The constant ratios indicate that these frequencies covary on an individual basis. The individual covariation of theta and alpha frequencies is somewhat less clear, but on average, theta peaks occurred close to half the alpha frequency. Pairs of peaks around 8 and 10 Hz are also seen, possibly representing split alpha. The finite width of the diagonal band arises because a certain minimum separation was necessary in order to resolve peaks; this does not imply a discontinuity in the frequencies of rhythms that can co-occur. The slopes of the frequency relationships are brought out in Figure 10D. The Hann-windowed spectra and the log-linear fits of the Welch-windowed spectra showed the same progression of peak frequencies as the log-log fits of the Welch-windowed spectra.

The relationships between peak frequencies are further illustrated in Figure 11, and Table 2 lists the corresponding correlation coefficients and fit parameters. The plots show only those subjects whose main peak was an alpha peak (Groups 2–5), in order to have only a single set of band limits for each range of alpha frequencies. Intercepts are included since these were found to be significant for both empirical and model peaks in many cases, and since fits without intercept would mainly reflect band limits.

**Figure 11. Relationships between peak frequencies**. Linear regressions were performed for each group separately to avoid spurious correlations induced by the adjustment of band limits to alpha peak frequencies. Correlation coefficients and parameters of the fits are indicated in Table 2. Dashed red lines indicate model predictions based on the approximations in Section 2.2 (*f _{θ}* =

*f*/2, ${f}_{{\beta}_{1}}=2{f}_{\alpha},{f}_{{\beta}_{2}}=3{f}_{\alpha}$); continuous red lines linear fits with intercepts based on fitted model spectra from Section 3.1. The dashed black line is a reminder that no peaks were fitted above 35 Hz. Significance levels: **0.01, ***0.001.

_{α}**Table 2**. **Correlation coefficients, the corresponding p-values, slopes, and intercepts for linear fits of theta, beta_{1}, and beta_{2} peak parameters vs. alpha peak parameters**.

Theta, beta_{1}, and beta_{2} frequencies of Group 4 (alpha frequencies in the range 9–11 Hz) had significant positive correlations (at the 0.05 level) with alpha peak frequencies. The theta trend line for this group is close to the theoretical prediction of *f _{θ}* = 0.5

*f*. Group 3 (alpha frequencies in the range 7–9 Hz) also showed a significant positive correlation between alpha and beta

_{α}_{1}frequencies. The same correlations were significant for the Hann-windowed spectra, apart from a positive theta-alpha correlation for Group 3 but not Group 4. Using log-linear instead of log-log fits of Welch-windowed spectra also yielded the same pattern of trends, except none of the theta-alpha correlations reached significance. However, all these correlations were positive.

The slopes of the beta_{1} trends for Groups 3 and 4 were 0.9 ± 0.5 and 1.2 ± 0.4, respectively. However, these slopes are affected by the rectangular sampling regions defined by the group-specific band limits, causing many points to lie to the top left and bottom right of a central region of higher density. The slope of this region is very close to 2, matching predictions based on the approximations in Section 2.2. The prediction based on peaks in fitted model spectra (cf., Section 3.1) yields beta_{1} frequencies slightly above the high-density region, and thus seems to be a somewhat poorer fit to the data.

For beta_{2} frequencies, we note that the model-based predictions may be better than they appear visually, since no empirical peaks were fitted above 35 Hz, producing a selection effect. Higher upper limits for the beta_{2} band might therefore have yielded additional points in the upper right-hand corner of the plot, giving a closer correspondence between empirical peak frequencies and model predictions. The Hann-windowed spectra and the log-linear fits of the Welch-windowed spectra showed the same significance levels as the log-log fits of the Welch-windowed spectra for both beta_{1}-alpha and beta_{2}-alpha frequency correlations.

As noted, the relationships between peaks in fitted model spectra are influenced by the empirical data themselves. The corresponding predictions may be considered as theoretical predictions with physiological parameter distributions, yet the findings should be interpreted with caution. To achieve a level of prediction intermediate between the parameter-independent ones from the approximations in Section 2.2, and the ones from fitted model spectra, we considered model spectra based on independent Gaussian distributions for the model parameters (e.g., gains and delays) with the empirical means and standard deviations, thus destroying any true correlations between the model parameters. This yielded an approximate mean alpha:beta_{1}:beta_{2} frequency ratio of 1:2.2:3.8, exceeding both the empirical ratios and the model predictions with correlated parameters. This implies that correlations between the parameters are important for the model to reproduce the empirical frequency relationships.

#### 3.2.2. Amplitude relationships

Figure 12 shows relationships between peak heights in the different bands, both differentiating between groups and for the sample as a whole. Note that for generality an intercept term was included in the regressions, in contrast to Figure 9. An *F*-test revealed that the intercept significantly improved each of the three whole-sample fits (*p* ≪ 0.001). However, for direct comparison with Figure 9, we also considered fits without intercept.

**Figure 12. Relationships between peak heights**. Linear least-squares fits are indicated by lines in the corresponding colors for each group, and in red for all subjects combined. No beta_{2} peaks were identified for Group 1. Correlation coefficients and parameters of the fits are listed in Table 2. Significance levels: n.s. non-significant, *0.05, ***0.001.

Alpha and theta peak heights of the combined groups lack a positive relationship. This matches the trend for Group 4 (alpha frequencies in the range 9–11 Hz), while Groups 3 and 5 have significant positive trends.

More convincing positive correlations are seen for beta_{1}, being significant for Groups 3–5 as well as for the sample as a whole. The overall slope is 0.11 ± 0.03. Discarding the intercept, the slope is 0.28 ± 0.01, consistent with the prediction of 0.266 ± 0.009 based on model fits.

The overall correlation between beta_{2} and alpha peak heights is 0.14 (*p* = 2.4*e* − 4). For beta_{2} peaks, the slopes are 0.04 ± 0.03 and 0.21 ± 0.01 with and without inclusion of the intercept, respectively, thus bracketing the predicted value of 0.076 ± 0.004. The beta_{2} trends are significantly positive for Groups 3 and 4, and similar in slope to each other and to the trend for Group 5.

The large variability of trends in theta peak height may be partly due to the requirement that theta peaks be higher than alpha peaks for Group 1 and vice versa for Groups 2–5. This constitutes a selection effect that may have increased the slopes of all trend lines, but that would have been strongest for Group 1, due to alpha peaks generally being higher than theta peaks. For Group 5 (alpha frequencies in the range 11–13 Hz), the positive trend may be partly explained by actual alpha and beta_{1} peaks being mislabeled respectively as theta and alpha peaks in a small proportion of cases. Thus, the definition of alpha as corresponding to the largest peak in the range 5–13 Hz may not be optimal, and it could for instance help to take subjects’ ages into account (Van Albada et al., 2010; Chiang et al., 2011). All in all, the relation between alpha and theta peak heights merits further investigation.

Peak height correlations for Hann-windowed spectra differed from those for Welch-windowed spectra for some groups and bands, but for the combined groups, the theta-alpha correlation was still insignificant, while beta_{1}-alpha and beta_{2}-alpha correlations were positive and highly significant. Moreover, for those cases where the significance levels differed greatly (theta height of Group 1 and beta_{2} height of Group 4), the linear trends were nevertheless quite similar. The same held for the log-linear fits of the Welch-windowed spectra.

We checked whether the positive overall associations between alpha and beta peak heights could be due to relationships between fit deviations in each band. The partial correlation between alpha and beta_{1} peak heights, corrected for deviations between log empirical and model spectra in both bands, is 0.33, close to the uncorrected correlation. However, the corrected correlation between alpha and beta_{2} peak heights is only 0.036. The positive correlation between fit deviations in these bands (*r* = 0.15) may itself be partly due to positively correlated peak heights, but this is impossible to verify without an objectively appropriate background subtraction.

Using independent Gaussian model parameter distributions with the empirical means and standard deviations, model spectra exhibited greater relative beta_{1} and beta_{2} amplitudes, the slopes of the fits without intercept being 0.35 for beta_{1} and 0.15 for beta_{2}. This provides a poorer match to the empirical results for beta_{1} but a better match for beta_{2}.

## 4. Discussion

Using a large sample (1424) of resting eyes-closed EEG spectra, we have shown clear interdependences between the frequencies and amplitudes of peaks in different bands in this condition, frequencies of many peaks following an approximately harmonic progression. These results strongly suggest that a common process contributes to the different rhythms.

Our main findings are: (i) a positive correlation between theta and alpha peak frequencies for subjects with alpha peak frequencies in the range 9–11 Hz, theta peaks occurring on average near half the alpha peak frequency for the sample as a whole; (ii) peaks in the low beta range tending to occur near twice the alpha peak frequency on an individual-subject basis, with positive correlations between frequencies of alpha and low beta peaks reaching significance for subjects with alpha peak frequencies in the range 7–11 Hz; (iii) peaks in the high beta range tending to occur near three times the alpha peak frequency on an individual-subject basis, with positive correlations between frequencies of alpha and high beta peaks reaching significance for subjects with alpha peak frequencies in the range 9–11 Hz; (iv) a lack of correlation between theta and alpha peak amplitudes for the sample as a whole; and (v) a positive, approximately linear, relationship between alpha and beta_{1} peak amplitudes for the sample as a whole. A positive correlation between alpha and beta_{2} peak amplitudes was also found, but further tests are needed to verify this result.

The harmonic progression of peak frequencies closely matches predictions based on an approximation of a linearized mean-field model of thalamocortical activity (Robinson et al., 2001b, 2005). It is not consistent with any of the following proposals: (i) a geometric progression with a peak spacing of Euler’s number (Penttonen and Buzsáki, 2003; Buzsáki and Draguhn, 2004) or the golden ratio (Roopun et al., 2008a,b; Pletzer et al., 2010); (ii) pacemakers that would not *a priori* be related in frequency or occurrence; (iii) Nunez’s theory of purely cortical eigenmodes, which predicts a non-harmonic sequence of peaks (Nunez, 1995). More generally, to our knowledge there is no model of purely cortical oscillations that predicts the observed peak relationships.

In view of the predicted relationships between peak frequencies, we adjusted band limits to the alpha peak frequency for peak classification. These limits could have been set separately for each subject, followed by a statistical analysis attempting to correct for this. However, to more strongly control for the effects of band limits, we defined only five sets of band limits and investigated trends in each group separately. This yielded group-specific frequency relationships that only reached significance for subjects with alpha peak frequencies in the range 7–11 Hz, probably related to the fact that these were the largest subject groups. It may be investigated in future studies whether larger sample sizes or different band definitions yield significance and similar slopes also in the other groups.

Relationships between peak amplitudes were in good agreement with predictions based on physiological considerations and model spectra. We identified possible contributions to both positive and negative associations between theta and alpha peak amplitudes, consistent with the overall lack of correlation found. Ratios between alpha, beta_{1}, and beta_{2} peak heights were close to those from fitted model spectra. Since the latter were partly influenced by the data themselves, we conclude that the empirical results match the model predictions at least semiquantitatively. Amplitude relationships, especially those between alpha and theta peaks, displayed variability between groups of subjects with different alpha peak frequencies. We discussed possible confounding effects that may account for this, suggesting a closer investigation of the theta band in particular.

According to the thalamocortical model, the observed peaks are largely explained by two scenarios: either the inhibitory thalamic reticular nucleus is weakly active, creating a positive thalamocortical feedback loop; or it is strongly active, creating a negative feedback loop. In the first case, the lowest-frequency resonance gives rise to an alpha peak corresponding to one pass through the loop, while in the latter case, it produces a theta peak corresponding to two passes. These basic rhythms are associated with near-harmonics around odd numbers of times the theta frequency, and integer numbers of times the alpha frequency. This agrees with the observed peaks around twice and three times the alpha peak frequency, and the hints of peaks around three times the theta peak frequency.

The covariation of peak frequencies suggests that band limits should be adjusted on an individual basis (at least for the resting-state condition considered here), as also proposed for instance by Doppelmayr et al. (1998) and Klimesch (1999), and consistent with age-associated changes in alpha peak frequency (Van Albada et al., 2010; Chiang et al., 2011) and various theoretical predictions (Nunez and Srinivasan, 1981; Nunez, 1995; Robinson et al., 2001b). It seems most expedient to base the limits on the alpha peak frequency, provided of course that an alpha peak is present. Consistent with the present study, Doppelmayr et al. (1998) argued for a positive association between bandwidth and alpha peak frequency. They measured task-related increases in theta and decreases in alpha peak power, and defined a transition frequency between ranges of increase and decrease. Individual alpha frequency was positively correlated not only with this transition frequency, but also with the difference between the alpha and transition frequencies. Thus, task-related activity confirms that higher alpha peak frequencies imply wider and higher-frequency EEG bands.

The present theoretical (cf., also Robinson et al., 2001b) and empirical results suggest that, for peak identification, band limits may be placed at approximately *n* + 1/4 and *n* + 3/4 times the individual alpha peak frequency (*n* = 0, 1, …), with theta and low beta peaks respectively being sought in the ranges 1/4–3/4 and 7/4–9/4 times the alpha peak frequency (see Figure 13). This follows especially from results where we avoided classifying peaks and examined pairs of peak frequencies within spectra. This clearly showed that peaks around 20 Hz often occurred together with peaks around half that frequency, and that peaks around 30 Hz often occurred together with peaks around one-third that frequency. These ratios were nearly constant despite individual variations in absolute frequencies. After peak classification, points also clustered around ${f}_{{\beta}_{1}}=2{f}_{\alpha}$ in those subjects with alpha peak frequencies in the range 7–11 Hz. The relationship between theta and alpha frequencies was slightly less clear, but theta peaks occurred on average very close to half the alpha peak frequency. Moreover, for the one group of subjects having a highly significant correlation between alpha and theta peak frequencies (those with alpha peak frequencies of 9–11 Hz), the trend line was close to *f _{θ}* =

*f*/2. Further research could ascertain whether the bands depicted in Figure 13 are appropriate for detecting task- or state-related changes.

_{α}The strictly individual adjustment of frequency bands is appropriate for within-subject comparisons where the alpha peak frequency is relatively stable. It may also be used for group comparisons when the distribution of band limits does not differ systematically between groups. However, depending on the questions asked, individual band adjustment may complicate analyses, since the band limits (and possibly associated filter characteristics, for instance when using wavelet analysis) affect spectral band power, peak characteristics, and the structure of the corresponding oscillations in the time domain. One option for dealing with this can be to define several subgroups, each with fixed frequency bands, as done in the present study. Alternatively, one could correct or account for differences in band definitions, for instance using analysis of variance where band limits constitute one of multiple factors.

Defining algorithms for peak fitting and classification naturally involves many choices that may influence the results. However, fitting Gaussian functions of frequency and of log frequency yielded qualitatively identical results. Moreover, we designed our methods to yield good agreement with visually identified peaks, and we consider it likely that any algorithm fulfilling this criterion will give similar results. This may be further investigated in future studies.

We allow for the possibility of contributions to the EEG which do not conform to the simple pattern of (sub)harmonics of alpha due to thalamocortical resonances. These may include cortically generated rhythms, rhythms originating in the hippocampus or amygdala, and intrathalamically generated rhythms such as sleep spindles (Robinson et al., 2002; Niedermeyer and Lopes da Silva, 2005). However, we argue that interpretations of EEG rhythms in terms of mechanisms, state dependence, and functional correlates should take into account their partially overlapping origins.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

We thank the individuals who gave their time to take part in the study. We are also grateful to Alan Chiang for elucidation of his fitting methods, to Eylem Kirlangic for several critical readings of the manuscript, and to Junji Ito for comments on the statistics. This work was supported by the Australian Research Council, the Westmead Millenium Institute, the National Health and Medical Research Council of Australia, and the Helmholtz Alliance on Systems Biology. We acknowledge the support of the Brain Resource International Database (under the auspices of Brain Resource, Ltd.; www.brainresource.com) for data acquisition and processing. Access to the database for scientific purposes is administered independently via the scientific network (BRAINnet; www.brainnet.net), which is coordinated by the non-profit BRAINnet Foundation.

## References

Barnett, T. P., Johnson, L. C., Naitoh, P., Hicks, N., and Nute, C. (1971). Bispectrum analysis of electroencephalogram signals during waking and sleeping. *Science* 172, 401–402.

Berger, H. (1933). Über das Elektrenkephalogramm des Menschen. *Arch. Psychiatr. Nervenkr.* 98, 231–254.

Braitenberg, V., and Schüz, A. (1998). *Cortex: Statistics and Geometry of Neuronal Connectivity*. Berlin: Springer.

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. *Science* 304, 1926–1929.

Carlqvist, H., Nikulin, V. V., Strömberg, J. O., and Brismar, T. (2005). Amplitude and phase relationship between alpha and beta oscillations in the human electroencephalogram. *Med. Biol. Eng. Comput.* 43, 599–607.

Chen, A. C. N., Feng, W., Zhao, H., Yin, Y., and Wang, P. (2008). EEG default mode network in the human brain: spectral regional field powers. *Neuroimage* 41, 561–574.

Chiang, A. K., Rennie, C. J., Robinson, P. A., Roberts, J. A., Rigozzi, M. K., Whitehouse, R. W., et al. (2008). Automated characterization of multiple alpha peaks in multi-site electroencephalograms. *J. Neurosci. Methods* 168, 396–411.

Chiang, A. K., Rennie, C. J., Robinson, P. A., van Albada, S. J., and Kerr, C. C. (2011). Age trends and sex differences of alpha rhythms including split alpha peaks. *Clin. Neurophysiol.* 122, 1505–1517.

David, O., and Friston, K. J. (2003). A neural mass model for MEG/EEG: coupling and neuronal dynamics. *Neuroimage* 20, 1743–1755.

Doppelmayr, M., Klimesch, W., Pachinger, T., and Ripper, B. (1998). Individual differences in brain dynamics: important implications for the calculation of event-related band power. *Biol. Cybern.* 79, 49–57.

Gordon, E., Cooper, N. J., Rennie, C. J., Hermens, D., and Williams, L. M. (2005). Integrative neuroscience: the role of a standardised database. *Clin. EEG Neurosci.* 36, 64–75.

Gratton, G., Coles, M. G. H., and Donchin, E. (1983). A new method for off-line removal of ocular artifact. *Electroencephalogr. Clin. Neurophysiol.* 55, 468–484.

Hill, S., and Tononi, G. (2005). Modeling sleep and wakefulness in the thalamocortical system. *J. Neurophysiol.* 93, 1671–1698.

Izhikevich, E. M., and Edelman, G. M. (2008). Large-scale model of mammalian thalamocortical systems. *Proc. Natl. Acad. Sci. U.S.A.* 105, 3593–3598.

Jirsa, V. K., Jantzen, K. J., Fuchs, A., and Kelso, J. A. S. (2002). Spatiotemporal forward solution of the EEG and MEG using network modeling. *IEEE Trans. Med. Imaging* 21, 493–504.

Klem, G., Lüders, H., Jasper, H., and Elger, C. (1999). “The ten-twenty electrode system of the international federation,” in *Recommendations for the Practice of Clinical Neurophysiology: Guidelines of the International Federatioin of Clinical Neurophysiology*, eds G. Deuschl and A. Eisen (Amsterdam: Elsevier), 3–7.

Klimesch, W. (1999). EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. *Brain Res. Rev.* 29, 169–195.

Liley, D. T. J., Alexander, D. M., Wright, J. J., and Aldous, M. D. (1999). Alpha rhythm emerges from large-scale networks of realistically coupled multicompartmental model cortical neurons. *Network* 10, 79–92.

Liley, D. T. J., Cadusch, P. J., and Dafilis, M. P. (2002). A spatially continuous mean field theory of electrocortical activity. *Network* 13, 67–113.

Lopes da Silva, F. H., van Lierop, T. H. M. T., Schrijer, C. F., and Storm van Leeuwen, W. (1973). Organization of thalamic and cortical alpha rhythms: spectra and coherences. *Electroencephalogr. Clin. Neurophysiol.* 35, 627–639.

Lopes da Silva, F. H., Vos, J. E., Mooibroek, J., and van Rotterdam, A. (1980). Relative contributions of intracortical and thalamo-cortical processes in the generation of alpha rhythms, revealed by partial coherence analysis. *Electroencephalogr. Clin. Neurophysiol.* 50, 449–456.

Lumer, E. D., Edelman, G. M., and Tononi, G. (1997). Neural dynamics in a model of the thalamocortical system I. Layers, loops and the emergence of fast synchronous rhythms. *Cereb. Cortex* 7, 207–227.

Marshall, P. J., Bar-Haim, Y., and Fox, N. A. (2002). Development of the EEG from 5 months to 4 years of age. *Clin. Neurophysiol.* 113, 1199–1208.

Niedermeyer, E. (1997). Alpha rhythms as physiological and abnormal phenomena. *Int. J. Psychophysiol.* 26, 31–49.

Niedermeyer, E., and Lopes da Silva, F. (eds). (2005). *Electroencephalography: Basic Principles, Clinical Applications, and Related Fields*. Philadelphia: Lippincott Williams & Wilkins.

Nunez, P. L. (1978). The relationship of head size to alpha frequency with implications to a brain wave model. *Electroencephalogr. Clin. Neurophysiol.* 44, 344–352.

Nunez, P. L., and Srinivasan, R. (1981). *Electric Fields of the Brain*. Oxford: Oxford University Press.

Palva, J. M., Palva, S., and Kaila, K. (2005). Phase synchrony among neuronal oscillations in the human cortex. *J. Neurosci.* 25, 3962–3972.

Penttonen, M., and Buzsáki, G. (2003). Natural logarithmic relationship between brain oscillators. *Thalamus Relat. Syst.* 2, 145–152.

Pletzer, B., Kerschbaum, H., and Klimesch, W. (2010). When frequencies never synchronize: the golden mean and the resting EEG. *Brain Res.* 1335, 91–102.

Posthuma, D., Neale, M. C., Boomsma, D. I., and de Geus, E. J. C. (2001). Are smarter brains running faster? Heritability of alpha peak frequency, IQ, and their interrelation. *Behav. Genet.* 31, 567–579.

Rennie, C. J., and Robinson, P. A. (2002). Unified neurophysical model of EEG spectra and evoked potentials. *Biol. Cybern.* 86, 457–471.

Roberts, J. A., and Robinson, P. A. (2008). Modeling absence seizure dynamics: implications for basic mechanisms and measurement of thalamocortical and corticothalamic latencies. *J. Theor. Biol.* 253, 189–201.

Robinson, P. A. (2003). Neurophysical theory of coherence and correlations of electroencephalographic and electrocorticographic signals. *J. Theor. Biol.* 222, 163–175.

Robinson, P. A., Loxley, P. N., O’Connor, S. C., and Rennie, C. J. (2001a). Modal analysis of corticothalamic dynamics, electroencephalographic spectra, and evoked potentials. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 63, 041909.

Robinson, P. A., Rennie, C. J., Wright, J. J., Bahramali, H., Gordon, E., and Rowe, D. L. (2001b). Prediction of electroencephalographic spectra from neurophysiology. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 63, 021903.

Robinson, P. A., and Rennie, C. J. (2010). *Brain Function Parameter Measurement System and Method*. U.S. patent 20100106043.

Robinson, P. A., Rennie, C. J., and Rowe, D. L. (2002). Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 65, 041924.

Robinson, P. A., Rennie, C. J., Rowe, D. L., and O’Connor, S. C. (2004). Estimation of multiscale neurophysiologic parameters by electroencephalographic means. *Hum. Brain Mapp.* 23, 53–72.

Robinson, P. A., Rennie, C. J., Rowe, D. L., O’Connor, S. C., and Gordon, E. (2005). Multiscale brain modelling. *Philos. Trans. R. Soc. Lond. B Biol. Sci.* 360, 1043–1050.

Robinson, P. A., Rennie, C. J., Rowe, D. L., O’Connor, S. C., Wright, J. J., Gordon, E., et al. (2003a). Neurophysical modeling of brain dynamics. *Neuropsychopharmacology* 28(Suppl. 1), 74–79.

Robinson, P. A., Whitehouse, R. J., and Rennie, C. J. (2003b). Nonuniform corticothalamic continuum model of electroencephalographic spectra with application to split-alpha peaks. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 68, 021922.

Robinson, P. A., Rennie, C. J., and Wright, J. J. (1997). Propagation and stability of waves of electrical activity in the cerebral cortex. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 56, 826–840.

Robinson, P. A., Rennie, C. J., Wright, J. J., and Bourke, P. D. (1998). Steady states and global dynamics of electrical activity in the cerebral cortex. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 58, 3557–3571.

Roopun, A. K., Kramer, M. A., Carracedo, L. M., Kaiser, M., Davies, C. H., Traub, R. D., et al. (2008a). Temporal interactions between cortical rhythms. *Front. Neurosci.* 2, 145–154. doi:10.3389/neuro.01.034.2008

Roopun, A. K., Kramer, M. A., Carracedo, L. M., Kaiser, M., Davies, C. H., Traub, R. D., et al. (2008b). Period concatenation underlies interactions between gamma and beta rhythms in neocortex. *Front. Cell Neurosci.* 2:1. doi:10.3389/neuro.03.001.2008

Rowe, D. L., Robinson, P. A., and Rennie, C. J. (2004). Estimation of neurophysiological parameters from the waking EEG using a biophysical model of brain dynamics. *J. Theor. Biol.* 231, 413–433.

Shackman, A. J., McMenamin, B. W., Greischar, L. L., and Davidson, R. J. (2010). Identifying robust and sensitive frequency bands for interrogating neural oscillations. *Neuroimage* 51, 1319–1333.

Snyder, A. C., and Foxe, J. J. (2010). Anticipatory attentional suppression of visual features indexed by oscillatory alpha-band power increases: a high-density electrical mapping study. *J. Neurosci.* 30, 4024–4032.

Steriade, M. (2000). Corticothalamic resonance, states of vigilance and mentation. *Neuroscience* 101, 243–276.

Steriade, M., Contreras, D., Amzica, F., and Timofeev, I. (1996). Synchronization of fast (30-40 Hz) spontaneous oscillations in intrathalamic and thalamocortical networks. *J. Neurosci.* 16, 2788–2808.

Steriade, M., McCormick, D. A., and Sejnowski, T. J. (1993). Thalamocortical oscillations in the sleeping and aroused brain. *Science* 262, 679–685.

Valdés-Hernández, P. A., Ojeda-González, A., Martínez-Montes, E., Lage-Castellanos, A., Virués-Alba, T., Valdés-Urrutia, L., et al. (2010). White matter architecture rather than cortical surface area correlates with the EEG alpha rhythm. *Neuroimage* 49, 2328–2339.

Van Albada, S. J., Gray, R. T., Drysdale, P. M., and Robinson, P. A. (2009). Mean-field modeling of the basal ganglia-thalamocortical system II: dynamics of Parkinsonian oscillations. *J. Theor. Biol.* 257, 664–688.

Van Albada, S. J., Kerr, C. C., Chiang, A. K. I., Rennie, C. J., and Robinson, P. A. (2010). Neurophysiological changes with age probed by inverse modeling of EEG spectra. *Clin. Neurophysiol.* 121, 21–38.

Van Albada, S. J., Rennie, C. J., and Robinson, P. A. (2007). Variability of model-free and model-based quantitative measures of EEG. *J. Integr. Neurosci.* 6, 279–307.

Van Albada, S. J., and Robinson, P. A. (2009). Mean-field modeling of the basal ganglia-thalamocortical system I: firing rates in healthy and parkinsonian states. *J. Theor. Biol.* 257, 642–663.

Van Rotterdam, A., and Lopes da Silva, F. H. (1982). A model of the spatial-temporal characteristics of the alpha rhythm. *Bull. Math. Biol.* 44, 283–305.

Keywords: EEG, spectra, peaks, mean-field model, neural field theory, harmonic frequencies

Citation: van Albada SJ and Robinson PA (2013) Relationships between electroencephalographic spectral peaks across frequency bands. *Front. Hum. Neurosci.* **7**:56. doi: 10.3389/fnhum.2013.00056

Received: 10 August 2012; Accepted: 11 February 2013;

Published online: 04 March 2013.

Edited by:

John J. Foxe, Albert Einstein College of Medicine, USACopyright: © 2013 van Albada and Robinson. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.

*Correspondence: S. J. van Albada, Institute of Neuroscience and Medicine (INM-6) and Institute for Advanced Simulation (IAS-6), Jülich Research Centre and JARA, Leo-Brandt-Straße, 52428 Jülich, Germany. e-mail: s.van.albada@fz-juelich.de