Components of cross-frequency modulation in health and disease

The cognitive deficits associated with schizophrenia are commonly believed to arise from the abnormal temporal integration of information, however a quantitative approach to assess network coordination is lacking. Here, we propose to use cross-frequency modulation (cfM), the dependence of local high-frequency activity on the phase of widespread low-frequency oscillations, as an indicator of network coordination and functional integration. In an exploratory analysis based on pre-existing data, we measured cfM from multi-channel EEG recordings acquired while schizophrenia patients (n = 47) and healthy controls (n = 130) performed an auditory oddball task. Novel application of independent component analysis (ICA) to modulation data delineated components with specific spatial and spectral profiles, the weights of which showed covariation with diagnosis. Global cfM was significantly greater in healthy controls (F_1,175 = 9.25, P < 0.005), while modulation at fronto-temporal electrodes was greater in patients (F_1,175 = 17.5, P < 0.0001). We further found that the weights of schizophrenia-relevant components were associated with genetic polymorphisms at previously identified risk loci. Global cfM decreased with copies of 957C allele in the gene for the dopamine D2 receptor (r = −0.20, P < 0.01) across all subjects. Additionally, greater “aberrant” fronto-temporal modulation in schizophrenia patients was correlated with several polymorphisms in the gene for the α2-subunit of the GABA_A receptor (GABRA2) as well as the total number of risk alleles in GABRA2 (r = 0.45, P < 0.01). Overall, our results indicate great promise for this approach in establishing patterns of cfM in health and disease and elucidating the roles of oscillatory interactions in functional connectivity.

Introduction

The clinical presentation of schizophrenia includes diverse symptoms involving perceptual, cognitive and emotional impairments. These abnormalities implicate extended brain networks, suggesting that the pathophysiology of schizophrenia involves impaired coordination between regions, rather than localized deficits (Andreasen et al., 1998; Friston, 1998). However, the neural mechanisms underlying normal network integration, and their possible deficits in schizophrenia, remain uncertain. Past studies have focused on the role of neuronal oscillations, particularly in the gamma band (∼30–100 Hz), that reflect the paced synchronous activity of local ensembles and are linked with many cognitive and sensory processes (Lee et al., 2003; Herrmann and Demiralp, 2005; Fries et al., 2007; Gonzalez-Burgos and Lewis, 2008; Uhlhaas et al., 2008). Reports of reduced gamma oscillations in schizophrenia have lead to speculations of altered cortical circuitry less able to synchronize at higher rhythms (Lee et al., 2003; Symond et al., 2005; Cho et al., 2006), however some studies report no group differences (Spencer et al., 2008) or even patient increases (Herrmann and Demiralp, 2005; Tekell et al., 2005; Flynn et al., 2008), suggesting that gamma power alterations do not reflect an inability to engage in higher frequency oscillations, but rather an inability to appropriately modulate this activation. Furthermore, the local nature of gamma oscillations is discrepant with the scale of network integration. Gamma synchronization is typically observed in intra-area assemblies where cells are separated by few synaptic connections, making high-frequency oscillations unlikely to mediate integration between distal structures (von Stein and Sarnthein, 2000). In contrast, low-frequency oscillations (∼1–20 Hz), are well-suited to coordinate long-distance communication as they are robust to conduction delays in polysynaptic pathways and are observed in more global, inter-area interactions, e.g., top–down processing (von Stein and Sarnthein, 2000) and attentional modulation (Lakatos et al., 2008).

It is hypothesized that lower-frequency oscillations coordinate information between areas by modulating excitability of local ensembles (von Stein and Sarnthein, 2000; Fries, 2005; Lisman and Buzsáki, 2008). Modulation of high-frequency amplitude by low-frequency phase, known as cross-frequency modulation (cfM), permits greater gamma synchronization during specific phases of the slower rhythm, theoretically allowing effective interactions between neurons with similar phase preferences (von Stein and Sarnthein, 2000; Fries, 2005; Lisman and Buzsáki, 2008). The phase of low-frequency oscillations has been noted to affect not only rhythmic gamma activity, but also the spiking rates of individual neurons (Canolty et al., 2010) and non-rhythmic broadband power (He et al., 2010; Miller et al., 2010), which has been associated with multi-unit activity (Ray and Maunsell, 2011). CfM has been observed in numerous brain regions, is associated with sensory perception and task performance, and is emerging as a robust organizational principal of neural signals (Lakatos et al., 2005, 2008; Mormann et al., 2005; Canolty et al., 2006; Demiralp et al., 2007; Cohen et al., 2009; Maris et al., 2011) and complex systems in general (He et al., 2010).

Given the symptomology of schizophrenia and the proposed role of cfM in inter-areal communication, we hypothesized that patterns of cfM would be significantly altered in schizophrenia patients. Departing from previous studies where investigation scope is limited to particular regions or phase–amplitude relationships due to the combinatorial nature of oscillatory interactions and complexity of data analysis, we propose a new approach for the thorough exploration of cfM that is also amenable to comparative analysis. We utilize independent component analysis (ICA) to discover fundamental patterns of cfM with no specific assumptions regarding their spatial or spectral composition. We apply this data-driven method to EEG recordings acquired from a large, heterogeneous group of participants and explore (1) cfM components identified with ICA and their variation over task conditions, (2) differences in cfM strength between healthy controls and schizophrenia patients, and (3) possible genetic contributions to cfM. Our analyses identify plausible components of cfM, some of which are different between subject groups, and implicate genes involved in dopaminergic and GABAergic transmission as contributors to cfM. Overall, these findings suggest great promise for this approach in advancing the investigation of oscillatory interactions in health and disease.

Materials and Methods

Participants

The current work combines existing data from several related studies performed at the Olin Neuropsychiatry Research Center. We consider data from 189 participants, including 47 schizophrenia patients (SZ), 130 healthy controls (HC) and 12 subjects who were unaffected first-degree relatives of individuals with schizophrenia (REL). Subjects were recruited by newspaper, word of mouth, and outpatient units at the Institute of Living in Hartford, CT, USA. First-degree relatives were specifically recruited as part of a large-scale study looking at probands and their unaffected relatives. Healthy controls were free from any history of DSM IV Axis-I diagnosis, as assessed with the structured clinical interview for DSM IV axis-I disorders (SCID) and were also interviewed to determine that there was no history of psychosis in any first-degree relatives. Patients met DSM IV TR diagnostic criteria for schizophrenia on the basis of the SCID. SCID diagnoses were made first by the clinician performing the interview, which is videotaped, and second by an independent clinician who rated the videotape. If diagnoses disagreed, a consensus meeting was called to review all available information including the patient’s case file. All participants were substance-abuse free as assessed by urine toxicology and provided written, IRB-approved consent by Hartford Hospital.

To include as many datasets as possible in our study, we decided to use all participants instead of the relatively balanced subset, while assessing the effect of age, gender, and other factors on our measurements. Subject demographics are displayed in Table 1. Control and patient groups showed no statistical differences on race, sex or handedness, but did show a significant age difference (HC: 26.0 ± 11.7, SZ: 37.6 ± 11.6; REL: 36.5 ± 13.0, F_{2, 186} = 18.89, P < 10⁻⁷). Therefore, age was regressed from all performance measures (i.e., IQ scores, reaction times, etc.) before comparing groups. IQ scores, estimated from the national adult reading test (NART) were available for 32 HC and 34 SZ. Positive and negative syndrome scale (PANSS) scores were available for 35 patients, however we restricted our analysis to symptom scores obtained within 2 weeks of EEG data acquisition (25 patients).

TABLE 1

Table 1. Demographic, clinical, and behavioral data for subject groups.

Of 47 schizophrenia patients, 35 were taking antipsychotic medication at the time of assessment; 6 took one or more typical antipsychotics, 24 took one or more atypical antipsychotics, and 5 patients took a combination. In addition, 14 patients took one or more anti-depressant medications (including both SSRI and non-SSRIs), 4 patients were prescribed anti-anxiety medication, and 4 patients were taking both anti-depressant and anti-anxiety medications. Based on the available PANSS scores, we found no difference in symptom severity for patients taking (19/25) and not taking (6/25) antipsychotics (positive scores, SZ with medication: 16.5 ± 4.5; without medication: 15.3 ± 9.1, t₂₃ = 0.41, P = 0.68, P_permute = 0.69; negative scores, SZ with medication: 19.2 ± 7.0; without medication: 15.2 ± 2.7, t₂₃ = 1.35, P = 0.19, P_permute = 0.19; general scores, SZ with medication: 33.9 ± 11.2; without medication: 31.6 ± 14.0, t₂₃ = 0.41, P = 0.68, P_permute = 0.68; two-tailed t-test and permutation test with n = 10,000 permutations).

Auditory Oddball Discrimination Task

Subjects were presented with a series of standard tones (500 Hz), target tones (1000 Hz) and novel stimuli (e.g., tone sweeps, whistles). Stimuli were 200 ms in duration with an inter-stimulus interval varying between 1000 and 2000 ms. Standard tones occurred with a probability of 0.80 and target and novel stimuli each occurred with a probability of 0.10. Subjects were asked to respond as quickly and accurately as possible with their right index finger to target stimuli and to ignore standard and novel stimuli. The task was performed in two 8-min runs, with approximately 250 stimuli per run. Age-corrected performance measures are provided in Table 1.

EEG Acquisition

Multi-channel EEG data were acquired using a Bioelectric Amplifier System (SA Instrumentation Co., San Diego, CA, USA). Scalp potentials were recorded at 62 electrode sites in accordance with the International 10–20 System. Electrooculogram data were recorded from electrodes located on the lateral and supra-orbital ridges of the right eye. All electrodes were referenced to the nose. Electrical impedances were maintained below 10 kΩ. EEG signals were amplified (20,000 gain) and digitized at 500 Hz.

Genotyping

Genotype data at 26 polymorphic loci putatively involved with schizophrenia were available for 170 of 189 participants. Polymorphisms, detailed in Table 2, were defined previously as part of several ongoing studies investigating the genetic basis of schizophrenia. Genotyping was performed at Yale using a fluorogenic 5′ nuclease assay or, for VNTR polymorphisms, by PCR amplification followed by agarose gel size fractionation. The SLC6A4 triallelic system was genotyped as reported previously (Stein et al., 2006). Samples analyzed with the taqman method were genotyped in duplicate for quality control; for gel-based genotyping, 8% of genotypes were repeated.

TABLE 2

Table 2. Genotyped polymorphic loci.

Data Analysis

All analyses were conducted using MATLAB (MathWorks, Inc, Natick, MA, USA). Statistical comparisons were thresholded at P < 0.05 following Bonferroni correction for the number of comparisons within each analysis. Continuous EEG recordings were band-pass filtered from 1 to 200 Hz, and 60 Hz line noise was removed by subtracting the best-fit sine wave from the raw data in overlapping windows (http://chronux.org). Eye blink artifacts were removed using the ICA algorithm within in the EEGLAB toolbox as described previously (Delorme and Makeig, 2004; Liu et al., 2009a). EEG signals were segmented around the arrival of the tones, from 500-ms pre-stimulus to 1200-ms post-stimulus, and data epochs were baseline corrected. Only trials with correct responses (targets detected within 1000 ms; no false alarms) were included in the analysis. Data were excluded on the basis of noisy or dead channels (<1% occurrence), trial variance across time > 3 mean ± SD, or voltage > 100 μV amplitude safeguard. On average, 20% of correct trials were rejected from each subject with no significant difference in the number of trials retained between groups (Table 1).

Cross-frequency modulation was quantified using the modulation index, as described previously (Canolty et al., 2006; Cohen et al., 2009; Figure 1A). Briefly, the preprocessed voltage trace from a single trial, x(t), was filtered into low-frequency (f_P) and high-frequency (f_A) bands, generating x_fP(t) and x_fA(t), respectively. Using the Hilbert transform, the analytic phase, φ_fP(t), and amplitude envelope, A_fA(t), were extracted from x_fP(t) and x_fA(t), respectively, to form a composite signal z(t) = A(t)e^iφ(t), which projects high-frequency power onto low-frequency phase. The raw modulation index was then computed as the mean vector length of z(t), . Intuitively, coupling between low-frequency phase and high-frequency amplitude is present if the probability distribution of z(t) is circularly non-uniform (see Figure 1A), equivalent to larger values of m_raw. Because the vector length is dependent on signal amplitude, raw modulation indices were z-scored using the mean, μ, and standard deviation, σ, of 50 surrogate values: m = (m_raw − μ)/σ, where surrogate values were generated under the null hypothesis of no modulation by randomly shuffling the amplitude time series to disrupt the phase–amplitude relationship (Cohen et al., 2009). Previous studies have used more surrogate values to generate the null distribution (e.g., 200 in Cohen et al., 2009), however this was not possible in the current work given the computational time required to generate distributions for our large dataset. In preliminary testing, we confirmed that 50 surrogates yielded m values similar to and unbiased from those obtained with 200. We further note that 50 surrogates is twice the number suggested by Tukey’s rule of thumb for estimating moments of a distribution, which recommends that calculation of the kth moment is based on at least 5^k observations (Sachs, 1992). Because cfM metrics can be affected by sharp edges in time series data (Kramer et al., 2008), we examined high-frequency peak-locked averages of data as well as bicoherence between bands to detect edges effects; we found little evidence for spurious coupling. We also examined spectrograms of trough-locked data segments and found clear modulation of high-frequency power by low-frequency phase (see Figure 3).

FIGURE 1

Figure 1. Schematic of cross-frequency modulation analysis. (A) Steps to compute the cfM index (m). AOD trials are pre-processed and filtered into low-frequency bands [e.g., f_P = 12–16 Hz, forming x_fP(t)] and high-frequency bands [e.g., f_A = 110–120 Hz, forming x_fA(t)]. Analytic phase, φ(t), and amplitude envelope, A(t), are extracted from x_fP(t) and x_fA(t), respectively, to form a composite signal: z(t) = A(t)e^iφ(t). Coupling between low-frequency phase and high-frequency amplitude is present if the probability distribution of z(t) is circularly non-uniform, equivalently, if the length of z(t), , is different from zero. Raw modulation indices are transformed into z-scores based on a null distribution from surrogate datasets. (B) Steps in (A) are repeated for all f_P and f_A combinations to produce the comodulogram. This is repeated over trials, conditions (target, novel or standard stimuli), EEG channels, and subjects to generate the full dataset. (C) For each subject, comodulograms are averaged over trials and data from all channels are concatenated to form a single row for each condition. Vertical concatenation of these rows forms the data matrix (X_M). ICA estimates a mixing matrix (A_M) and set of independent components (S_M) from X_M. The columns of A_M indicate the loading parameters or weights of a particular component for each subject and condition. The rows of S_M correspond to the component topography and spectral composition.

For each trial of EEG data, cfM was calculated between six low-frequency phase bands (f_P = 1–24 Hz, width = 4 Hz) and 15 high-frequency amplitude bands (f_A = 30–200 Hz, width = 11.3 Hz) to generate a comodulogram (e.g., Figure 1B). This procedure was performed at 62 scalp channels and comodulograms were averaged over trials within each condition (target, novel or standard stimuli), creating a data set of 90 modulation indices × 62 locations × 3 conditions for each subject.

We used ICA to decompose patterns of cfM from the large set of comodulograms (Hyvärinen et al., 2001). To apply ICA to the modulation data, comodulograms were flattened into row vectors and the means were subtracted. Because modulation indices have been z-scored to form appropriate metrics of coupling (see above), no further amplitude standardization between frequencies, channels, or subjects is necessary. Rows were then vertically concatenated over conditions and subjects to form the data matrix (X_M). The number of components in X_M was estimated following the method of Liu et al. (Liu et al., 2009a): we used the Akaike information criterion (AIC) as an upper bound of dimensionality since it tends to overestimate component number (De Ridder et al., 2005) and reduced the number of components until they were stable between 500 bootstrap subject resamplings. AIC estimated the number of components at 9, which we reduced to 6 based using a stopping criterion of cluster quality index of greater than 0.8 for all components (Himberg et al., 2004). ICA was performed using the Infomax algorithm (Bell and Sejnowski, 1995), which linearly decomposed X_M into six sources (rows of S_M) after PCA compression, and their associated loading parameters (columns of A_M). From the full decomposition, we observed that modulation loading parameters showed association with subject age (IC_M 1: r = −0.19, P < 0.01; IC_M 4: r = 0.27, P < 0.001; IC_M 5: r = −0.13, P = 0.07) or sex (IC_M 1: t₁₈₇ = 2.04, P = 0.04; IC_M 2: t₁₈₇ = 2.99, P = 0.05). As a conservative correction, we linearly regressed these factors from the columns of A_M before proceeding with statistical comparisons between groups or across genotypes.

Along with the bootstrap subject resamplings, we ensured that the sources were representative of features present in all subjects and trials by visually comparing (1) separate decompositions of HC and SZ groups, (2) separate decompositions using only target, novel or standard conditions, and (3) decompositions using randomly selected trials to match the number of trials across all subjects. In all cases, the identified sources were similar to those found from the full data set.

Trough-locked spectrograms of gamma power (see Figure 3) were computed as outlined in (Canolty et al., 2006). From data at electrodes FZ and OZ, we identified theta (4–8 Hz) or beta (12–16 Hz) troughs as local minima in the filtered phase time series, then extracted the raw data centered on each trough in 600 ms (beta) or 1200 ms (theta) windows. Normalized gamma power time series of these data segments were then calculated by (1) band-pass filtering the extracted signal into 30 linearly spaced bands between 30 and 200 Hz (narrower bands than used to generate the comodulograms), (2) z-scoring each band-passed signal, (3) applying the Hilbert transform to obtain the amplitude time series of each band, and (4) point-wise squaring the amplitude time series to obtain the power time series. Spectrograms of the normalized gamma power show the average over data segments from all trials and conditions. The number of segments included for single subjects were 2,442 (beta-locked, Figure 3A, left) and 1,274 (theta-locked, Figure 3B, left); the number of segments for the group average were 449,631 (Figure 3A, right) and 221,697 (Figure 3B, right).

For analyses of genetic data, genotypes at biallelic loci were coded as −1, 0, or 1, corresponding to minor homozygous, heterozygous, and major homozygous. At loci with more than two alleles (SLC6A3 VNTR, DRD4 exon 3 VNTR, and SLC6A4 promoter region L_A/L_G/S) we initially sorted genotypes in order of possible functional effects (Asghari et al., 1995; Fuke et al., 2001; Hu et al., 2006). For comparison, we also re-classified genotypes to conform to a biallelic ([A|B]) system (SLC6A3: [≥10|≤9 repeats]; DRD4: [≥7|≤6 repeats]; SLC6A4: [L_A|L_G or S]) and found no difference in the results. Thus we retained biallelic classification at all loci for simplicity.

Before performing ICA on the genetic data, we imputed missing values (<5% of all genotypes). For loci in linkage disequilibrium (LD; r² > 0.8), we used linked polymorphisms as a reference. For remaining loci, missing values were imputed using the genotype from the individual with the most similar set of genetic data. We compared this method with filling missing values with the most common, least common, or a randomly assigned genotype. None of these methods produced significantly different results in the sources identified by ICA, signifying little impact of imputation. For direct associations of genotypes with IC_M loading parameters (Figures 6C–F) we used raw (un-imputed) data, thus the number genotyped individuals varies slightly between polymorphisms.

Following imputation, subject rows of genetic data were concatenated and the means were removed to form X_G. AIC estimated 25 genetic sources, which we reduced to 16 based on stability in n = 1000 bootstrap resamplings (Liu et al., 2009a,b). Relative contributions of polymorphisms to components were determined from the variability between bootstrap resamplings. Homologous components between runs were identified using the minimum distance (d) between components, where d = 1−|r_ij|, and r_ij is the correlation coefficient between component S_M(i) and bootstrapped component (Himberg et al., 2004).

Results

We examined patterns of phase-to-amplitude cfM in EEG data recorded from 47 schizophrenia patients (SZ) and 130 healthy controls (HC), and 12 participants who were unaffected first-degree relatives of individuals with schizophrenia (REL) during performance of an auditory oddball paradigm. This paradigm asks subjects to detect and respond to rare target tones presented in a series of frequent standard tones and infrequent novel sounds, and has been repeatedly demonstrated to evoke differential scalp potentials in HC and SZ (Braff et al., 2007; Allen et al., 2009). Table 1 provides demographic and task performance data.

Cross-frequency modulation was quantified using a modulation index based on the distribution of high-frequency amplitude as a function of low-frequency phase (Canolty et al., 2006; Figure 1A). Modulation indices were calculated between pairs of low-frequency phase bands (f_P = 1–24 Hz, 6 bands) and high-frequency amplitude bands (f_A = 30–200 Hz, 15 bands). The pattern of cfM is displayed in a comodulogram, with phase frequency (f_P) along the abscissa and amplitude frequency (f_A) on the ordinate (Figure 1B). Comodulograms over 62 channels and 3 conditions (target, novel and standard stimuli) from the 189 subjects were decomposed into six independent sources using ICA (Hyvärinen et al., 2001; Figure 1C). Applied here, ICA decomposes the comodulograms into linear factors based on the covariation of cfM patterns across subjects and the maximal independence of factors over spatial topography (EEG channels) and frequency composition (profile over f_P and f_A). This approach allows for different patterns of cfM to have overlapping topographies and makes no specific assumptions regarding the frequency composition of each factor.

Cross-Frequency Modulation Components

Six cfM components (IC_M) accounting for 96.1% of the data variance are shown in Figure 2. Each IC_M has a scalp topography (Figure 2A), identified from the frequency bin with greatest modulation amplitude, and a comodulogram (Figure 2B) which is displayed for a representative channel. Across channels, comodulograms are quite consistent, with differences largely captured by changes in scale. To aid interpretation, the IC_M are scaled to the original data units (z-scores) using multiple regression. As indicated by the color bars adjacent to each panel in Figure 2, the magnitudes of cfM vary greatly between components, with IC_M 1 capturing the large majority of cfM variance. Because IC_M 1 exhibits such widespread spatial and spectral activation, other components can be considered as representing more localized addition or subtraction of modulation in particular frequency bands.

FIGURE 2

Figure 2. Cross-frequency modulation components. Component scalp topographies (A), spectral comodulograms (B), and loading parameters (C). The scalp distribution of each component is identified using the amplitude of the modulation index from the frequency bin with the greatest average amplitude over channels. Comodulograms in (B) are displayed for a single representative channel, indicated by the white dots in (A). As indicated by the color scales adjacent to each panel, the magnitudes of cfM vary greatly between components; components are listed in decreasing order of their variance. (C) Component loading parameters for each condition stratified by SZ (red) and HC (blue) groups (REL subjects not shown). Error bars in this and all subsequent plots denote ± 1 SEM. T* and G* indicate significant differences over conditions or groups, respectively, as determined with a repeated measures one-way ANOVA (P < 0.05, Bonferroni corrected for 6 tests).

Several features are apparent in the spatial and spectral patterns of cfM components. First, some components exhibit great specificity. For example, IC_M 2 shows significant modulation almost exclusively between theta phase (f_P = 4–8 Hz) and gamma amplitude (f_A = 30–200 Hz), spatially localized to posterior/occipital electrodes. IC_M 5 shows a similar pattern, with delta phase (f_P = 1–4 Hz) modulating the amplitudes of high-frequency oscillations (f_A = 30–120 Hz) over central/posterior electrodes. A second notable feature is a general antagonism/opposition between cfM at lower and higher phase frequencies. For example, IC_M 3, which is localized to fronto-temporal electrodes, shows mild positive modulation between phases at f_P ≤ 8 Hz and gamma amplitude, with simultaneous negative modulation between phases at f_P ≥ 12 Hz and high-frequency amplitude. IC_M 4, which also peaks at fronto-temporal electrodes but extends through central and more posterior electrodes, shows a similar pattern, though with opposite polarity. Third, there is a diagonal trend across the comodulograms, most apparent in IC_M 1. Although IC_M 1 exhibits nearly global activation, modulation is notably absent or slightly negative below a diagonal representing roughly f_A/f_P ≤ 4 (i.e., right, lower corner). The diagonal trend, which is also seen in IC_M 3, IC_M 4, and the comodulogram grand average (Figure 1B), is consistent with a theorized optimal relationship between the phase and amplitude bands for cfM to occur (Roopun et al., 2008).

To ensure that our ICA approach accurately captured cfM patterns, we examined fluctuations in gamma power at individual channels, shown in Figure 3 (see Materials and Methods). We verified robust oscillations in FZ gamma power locked to beta phase (f_P = 12–16 Hz, Figure 3A), as well as OZ gamma power modulated by theta phase (f_P = 4–8 Hz, Figure 3B), validating patterns of cfM reflected in components IC_M 1 and IC_M 2, respectively. We also ensured that the cfM components represented features found in all subjects by performing separate ICA decompositions on the HC and SZ groups separately. All six components were present in the separate decompositions (not shown), though were noticeably noisier for the SZ group which is expected given the smaller sample size.

FIGURE 3

Figure 3. Examples of phase-to-amplitude modulation. (A,B) Spectrograms of mean normalized gamma power, time-locked to the beta [f_P = 12–16 Hz, (A)] or theta [f_P = 4–8 Hz, (B)] trough for a single subject (top left) and all subjects (top right). Spectrograms include data segments from all trials and conditions. Gamma power is normalized within narrow bands (∼6 Hz) to permit comparison across frequencies. (bottom) Plot of the averaged beta (A) or theta (B) trough-locked signal (black) and normalized gamma power, averaged over frequency (f_A = 30–200 Hz, gray). Examples of beta modulation (A) are recorded from frontal electrode FZ (see Figure 2A), while examples of theta modulation (B) are from posterior electrode OZ.

We next examined condition-specific loading parameters, which indicate the contribution of each component to the modulation data for an individual subject. Component loading parameters for HC and SZ, corrected for age and sex via regression, are displayed in Figure 2C. Weights between conditions were highly correlated within subjects (Figures 4A,B), thus we evaluated differences between groups and conditions with a repeated measures ANOVA. After Bonferroni correcting for tests over 6 components, the ANOVAs indicated significant group differences for the loading parameters of IC_M 1 and IC_M 4 (F_1,175 = 9.25, P < 0.005 and F_1,175 = 17.5, P < 0.0001, respectively). These differences were in opposite directions: IC_M 1 weights were larger in HC, while IC_M 4 weights were larger in SZ. Though not included in the ANOVA, we observed that the loading parameters of the REL subjects were more similar to those of SZ for IC_M 1 (mean ± SD; HC: 0.48 ± 0.03; SZ: 0.45 ± 0.05; REL: 0.45 ± 0.05), and were intermediate for IC_M 4 (REL: 0.08 ± 0.05; SZ: 0.12 ± 0.08; REL: 0.10 ± 0.09). Despite the small number of REL subjects, the average loading parameters of IC_M 1 were significantly different between HC and REL (t₁₄₀ = 2.27, P < 0.05) and showed a trend for IC_M 4 (t₁₄₀ = 1.68, P = 0.10). Loading parameters of components 1 and 4 were also significantly negatively correlated (r = −0.56, P < 10⁻¹⁵, Figure 4C), and this relationship was particularly strong for schizophrenia patients (HC: r = −0.37, P < 0.0001; SZ: r = −0.64, P < 0.00001; REL: r = −0.68, P < 0.05). These findings suggest that subjects with less global cfM (IC_M 1), have altered modulation in fronto-temporal regions (IC_M 4), and that this differential pattern of modulation is most pronounced in SZ patients.

FIGURE 4

Figure 4. Relationships between cross-frequency modulation loading parameters. (A) Full correlation matrix between loading parameters (A_M) for each cfM component and each condition (target, novel and standard stimuli). In general, component weights are well correlated between conditions (values along the diagonal). (B) Example of correlation between loading parameters of target tones and standard tones for IC_M 1. Solid gray line indicates unity. (C) Negative correlation between the loading parameters of IC_M 1 and IC_M 4. Here, the loading parameters have been averaged over condition. Dashed black line shows the least-squares linear fit to the data.

Given the group difference, we then examined whether cfM patterns covaried with symptoms in SZ. We found no relationships between PANSS scores and loading parameters of components 1 or 4, but found a significant correlation between negative symptom scores and the weights of IC_M 6 after correcting for 18 tests (r = −0.69, P < 0.005, Figures 5A,B). This indicates that more severe negative symptoms are associated with less cfM between low-frequency phase and high-gamma power.

FIGURE 5

Figure 5. Relationship between cross-frequency modulation and patient symptoms. (A) Plot of correlation coefficients between IC_M loading parameters (averaged over condition) and positive (red), negative (blue), and general (black) symptom scores. Only patients with PANSS scores collected within 2 weeks EEG acquisition were included in this analysis (n = 25). Error bars ( ± 1 SEM) estimated with 1000 bootstrap resamplings. A single significant correlation (indicated by the asterisk, after Bonferroni correction for 18 tests) was found between negative symptoms and IC_M 6 (r = −0.69, P < 0.005). (B) Scatter plot of negative symptom scores and the loading parameters of IC_M 6. Dashed black line represents the least-squares linear fit.

The ANOVAs also identified significant differences over the target, novel and standard conditions for IC_M 1, IC_M 2, and IC_M 5 (F_2,175 = 10.3, P < 0.00005; F_2,175 = 5.23, P < 0.01; and F_2,175 = 6.73, P < 0.005, respectively; Bonferroni corrected for 6 tests; see Figure 2C). For IC_M 5, subject weights increased with the salience or significance of the tone (target > novel > standard). This trend was reversed for IC_M 1 and IC_M 2, which showed the least activation during target trials. Note that we found no relationship between IC_M weights and various measures of task performance (mean reaction time, percent correct responses and percent false alarms; P > 0.05 for all tests), thus it is unlikely that any of the components represent neuronal activity directly underlying target detection. Rather, trends likely reflect general changes in the patterns of activity that accompany different conditions. An increase in delta phase modulation (IC_M 5) may reflect sensory selection (Lakatos et al., 2008; Monto et al., 2008; Händel and Haarmeier, 2009), or enhanced coordination between auditory and motor systems in preparation for a response. The cfM in theta and other bands (IC_M 2 and IC_M 1) may indicate ongoing internal processes that are disrupted or dampened by the arrival of attention-demanding stimuli.

Genetic Associations

Having identified schizophrenia-relevant components of cfM, we explored whether these were associated with genetic factors. Genotypes from 26 polymorphic loci in putative SZ risk genes were available for 170/189 subjects (36/47 SZ and 134/142 HC, including 9/12 REL). Genetic data included single nucleotide polymorphisms (SNPs) and variable number tandem repeats (VNTRs) from 16 different genes (see Table 2). Several of the polymorphisms were in close physical proximity and showed some degree of LD (Figure A1 in Appendix). To avoid redundant tests for association with cfM components, we again employed ICA to identify independent genetic factors (Liu et al., 2009a).

Independent component analysis decomposed the genetic data into 16 components, capturing 94.9% of the original variance. We evaluated associations between cfM and genetic components (IC_G) by performing linear regression between the loading parameters of IC_M 1 and IC_M 4 and the loading parameters of the IC_Gs, with subject diagnosis and genetic × diagnosis interaction terms included in the regression models (Begleiter and Porjesz, 2006). Associations between (1) IC_M 1 and IC_G 6 (Figure 6A, top; F_3,166 = 6.25, P < 0.0005), (2) IC_M 4 and IC_G 6 (Figure 6A, bottom; F_3,166 = 8.50, P < 0.00005), and (3) IC_M 4 and IC_G 4 (Figure 6A, bottom; F_3,166 = 7.65, P < 0.0001) were significant at P < 0.05 after correction for 32 tests. Examination of the individual beta weights revealed a main effect of genetic component 6 on IC_M 1 (i.e., independent of diagnosis; Figure 6A, top; t₁₆₆ = 2.92, P < 0.005). Associations with IC_M 4 were both due to interactions between the genetic components and the diagnosis (Figure 6A, bottom; t₁₆₆ = 3.59, P < 0.0005 and t₁₆₆ = 3.35, P < 0.005 for IC_G 6 and IC_G 4, respectively). Including subject race in the regression model yielded nearly identical results (Figure A2A in Appendix), suggesting that correlations were not due to population stratification. We also verified that the loading parameters IC_G 4 and IC_G 6 showed no significant association with age or sex (P ≥ 0.10 for all comparisons).

FIGURE 6

Figure 6. Genetic contributions to cross-frequency modulation. (A) Plot of the negative natural logarithm of the P-values for the regression models between each of the 16 genetic components and cfM components 1 (top) and 4 (bottom). Pink shaded region indicates the P-values for the regression model; black and white regions show the P-values for the genetic and genetic × diagnosis terms, respectively. Asterisks denote significant models at P < 0.05 after Bonferroni correction for 32 tests. (B) Genetic component 6 (top) and component 4 (bottom). Contributions of individual polymorphisms are z-scored to facilitate interpretation. Error bars were determined with 1000 bootstrap resamplings. Asterisks denote polymorphisms contributing weights significantly different from zero (P < 0.05, corrected for 26 tests). (C–F) Plots of modulation component loading parameters as a function of genotype. In (D) and (E), data is stratified by SZ (red) and HC (blue) to highlight group × diagnosis interactions. In (F), the number of risk alleles is determined from the genotypes of rs279869, rs279858, rs279837, and rs567926; for clarity, only the SZ group data is shown. Dashed line shows the least-squares linear fit to the data. The number of individuals with each genotype is indicated adjacent to the data marker.

Examining the genetic components, IC_G 6 showed a large contribution from rs6277 (see Figure 6B, top), a polymorphism within the gene for the dopamine D2 receptor (DRD2), and a smaller contribution from rs1799732, located in the promoter region of DRD2. IC_G 4 (Figure 6B, bottom) showed contributions from rs279869, rs279858, rs279837, and rs567926, all found within or flanking the gene for the α2-subunit of the GABA_A receptor (GABRA2). Using bootstrap resamplings to assess the component stability (n = 1000) we found that for IC_G 6, only SNP rs6277 contributed a weight significantly different from zero (P < 0.05, corrected for 26 polymorphisms). For IC_G 4, SNPs rs279869, rs279858, and rs567926, all in relatively high LD with one another (r² > 0.5, see Figure A1), contributed significant weights.

Given that the genetic components were largely composed of single (or a few linked) polymorphisms, we performed post hoc association tests between cfM components and individual SNPs. Shown in Figure 6C, the DRD2 rs6277 genotype was significantly correlated with the loading parameters of IC_M 1 (r = −0.20, P < 0.01). As expected from the regression model, this trend was present within both HC and SZ groups (r = −0.20, P < 0.05 and r = −0.30, P = 0.08, respectively), with copies of the C allele associated with weaker cfM. For IC_M 4, the regression model indicated an interaction between genetic components and diagnosis, thus we examined the effect of DRD2 rs6277 genotype while stratifying by group (Figure 6D). We found a strong association for the SZ subjects, (r = 0.64, P < 0.00005) and no trend for the HC group (r = 0.02, P = 0.82). We confirmed that these relationships were not due to heterogeneous population structure, as correlations persisted in the subset of subjects identifying as Caucasian (IC_M 1, all subjects, n = 115: r = −0.24, P < 0.05; IC_M 4, SZ, n = 26: r = 0.61, P < 0.001). We additionally verified that DRD2 rs6277 genotypes were not significantly different across subject age (r = 0.01, P = 0.87), sex or (χ²_2,166 = 1.79, P = 0.40).

Relationships between the loading parameters of IC_M 4 and GABRA2 genotypes also varied as a function of diagnosis. SZ exhibited a significant correlation with GABRA2 rs279869 genotype (Figure 6E; r = 0.42, P < 0.05), which we verified within the Caucasian subset (r = 0.48, P < 0.05). In contrast, HC showed no association between rs279869 genotype and IC_M 4 loading parameters (r = −0.07, P = 0.46). Similar trends were observed for the additional three GABRA2-associated polymorphisms (Figure A2B in Appendix). As with DRD2 rs6277, GABRA2 genotypes showed no statistical distinction across age (e.g., GABRA2 rs279869: r = −0.06, P = 0.45), or sex (e.g., GABRA2 rs279869: χ²_2,162 = 3.50, P = 0.17; with similar findings for rs279858, rs279837 and rs567926). To putatively combine information over the different loci, we computed the number of risk alleles from the four GABRA2 SNPs. The risk allele was determined as the allele with the quantitative trait furthest from the HC group (e.g., allele A in Figure 6E). We found a significant positive correlation between the number of risk alleles in SZ and the magnitude of modulation component 4 (Figure 6F, r = 0.45, P < 0.01), indicating a roughly additive effect of the individual GABRA2 polymorphisms. No such relationship was found for the HC group (r = −0.01, P = 0.91), suggesting that these GABRA2 genotypes interact with other factors (genetic or environmental) specific to the schizophrenia background.

Discussion

In this study, we have used data-driven methods to examine patterns of cfM in a large (n = 189), heterogeneous group of subjects. This approach has yielded three main contributions: (1) the identification of spatially distinct and concurrent interactions between low-frequency phase and high-frequency power, (2) the demonstration of robust differences in cfM patterns between controls and schizophrenia patients, and (3) the demonstration of associations between cfM and genetic polymorphisms. We discuss the significance of each of these findings and their limitations in turn.

Phase-to-amplitude coupling between low- and high-frequency signals has emerged as a fundamental and ubiquitous feature of neural oscillations. First described between the phase of the robust hippocampal theta rhythm and firing of individual place cells (O’Keefe and Recce, 1993), cfM has since been identified within numerous cortical and sub-cortical areas (Lakatos et al., 2005; Mormann et al., 2005; Canolty et al., 2006; Demiralp et al., 2007; Osipova et al., 2008; Cohen et al., 2009), as well as between distinct regions (Bruns and Eckhorn, 2004; de Lange et al., 2008; Sirota et al., 2008; Maris et al., 2011). The spectral characteristics of cross-frequency coupling have also been expanded: high-frequency amplitudes are modulated not only by the theta rhythm (Canolty et al., 2006) but also by oscillations in the delta (Lakatos et al., 2005; Monto et al., 2008; Händel and Haarmeier, 2009), alpha (Osipova et al., 2008; Cohen et al., 2009; Voytek et al., 2010), and beta bands (de Lange et al., 2008). Furthermore, there is evidence for a “nested hierarchy” of interacting oscillations (Lakatos et al., 2005, 2008; Monto et al., 2008).

Our examination of cfM used ICA to identify modulation features at the group level. A principle advantage of this approach is the reduction of a large dataset that is both spatially and spectrally redundant into a manageable number of components well-suited for visual and statistical comparisons. We note that the combinatorial nature of this dataset (i.e., 90 modulation indices × 62 locations × 3 conditions) would have prohibited effective comparisons of cfM between groups or task conditions without the benefit of dimension reduction/feature identification. Such approaches will become increasingly important as studies investigate “inter-areal” coupling between oscillations at different spatial locations (Bruns and Eckhorn, 2004; de Lange et al., 2008), as demonstrated recently by Maris et al. (2011). In an innovative and systematic investigation, Maris et al. (2011) computed phase-to-amplitude coupling between all pairs of channels and used tensor decomposition, a generalization of singular value decomposition for arrays with more than two dimensions, to identify sources of cfM and characterize the spatial and spectral properties of low-frequency coupling oscillations and the corresponding high-frequency phase-locked activity. While tensor decomposition makes assumptions regarding the spatio-spectral separability of cfM components (i.e., that coupling can be fully described as the product of spatial maps and frequency profiles), such assumptions are quite plausible for phase-amplitude coupling, and we expect that this and other decomposition methods will be used increasingly to describe and understand oscillatory interactions.

A second benefit of ICA is the discovery of multiple independent features of modulation, without any constraints regarding their spatial or spectral composition. In agreement with previous studies, we found evidence for the presence of specific theta-to-gamma (IC_M 2; Canolty et al., 2006; Demiralp et al., 2007) and delta-to-gamma coupling (IC_M 5; Lakatos et al., 2005; Monto et al., 2008; Händel and Haarmeier, 2009). However, ICA additionally revealed the occurrence of numerous, simultaneous, cross-frequency interactions with specific spatial profiles (see Figure 2), suggesting a role for cfM in the concurrent coordination of multiple neural ensembles (Canolty et al., 2010; Maris et al., 2011). The differential modulation of cross-frequency components by condition supports this view: component loading parameters were significantly increased (IC_M 5) or decreased (IC_M 1, IC_M 2) as the tone saliency and behavioral relevance increased. Further support comes from recent work by Voytek and colleagues, who demonstrate localized differences in the relative strength of theta-to-gamma coupling and alpha-to-gamma coupling as a function of task demands (Voytek et al., 2010), putatively reflecting dynamic alterations in functional connectivity and regional engagement.

While ICA affords advantages in terms of analytic feasibility and scope, it is important to mention the interpretive challenges of such an approach. Because components have no inherent meaning ascribed to them beyond statistical independence, researchers must make their own interpretations of topographic and frequency patterns. Importantly, researches must view the components as a whole, keeping in mind that they simply represent a linear decomposition of the data. As an example, consider the pattern of cfM represented by IC_M 4. One notices subtle “gaps” and “bumps” in the modulation indices around f_A = 140 Hz and f_A = 170 Hz. These “gaps” are puzzling: it is commonly believed that high-gamma power (e.g., f_A > 80 Hz) represents a single localized broadband signal that should be uniformly modulated by low-frequency phase (He et al., 2010; Miller et al., 2010), thus one expects a somewhat continuous pattern of cfM extending over the high-gamma range (e.g., see Figure 3). However, we note that IC_M 1 shows the inverse pattern of “gaps,” i.e., greater modulation in frequency bands that are reduced for IC_M 4. Because these components share some similar spatial and spectral features (e.g., greater weights over fronto/temporal electrodes) as well as covariance in loadings (see Figure 4A), the separation of these sources with ICA may be imperfect and result in a general “give and take” in modulation values between the components. The discontinuous pattern of cfM in IC_M 4, and well as the cfM specific to very high frequencies (f_A > 150 Hz) in IC_M 6, may truly represent modulation within specific high-gamma sub-bands, however they may also represent an incomplete “un-mixing” of sources particular to this dataset and decomposition. Only with replication in independent datasets can we assess the robustness and functional relevance (if any) of these components. We note that ICA has proven to be of great value for delineating meaningful EEG timeseries and imaging spatial maps (for review see Calhoun et al., 2009), however, as with any decomposition technique, there is no guarantee that identified components will have physiological relevance or that they will represent specific neural processes.

Consistent with our predictions, patterns of cfM were significantly altered in patients, who showed decreased global cfM (IC_M 1) and increased coupling at fronto-temporal electrodes (IC_M 4). Intriguingly, the small number of unaffected SZ relatives (n = 12) exhibited intermediate magnitudes of IC_M 1 and IC_M 4. We speculate that specific cfM patterns could constitute an endophenotype or disease liability marker (Braff et al., 2007; Allen et al., 2009), as supported by our genetic association results (see below). While it is tempting to relate the altered cfM found in schizophrenia patients to the etiology of the disease, we feel that such inference is premature. Cross-frequency interactions were assessed while subjects engaged in a single behavioral paradigm known to evoke different event-related responses between subject groups (Braff et al., 2007; Allen et al., 2009) and it is quite possible that the observed differences are related to the evoked potentials or task-specific behaviors. We noted significant differences in task performance between the subject groups (see Table 1), and while the strength of cfM did not directly correlate with any recorded measures of task performance (see Results), we cannot rule out the possibility that group effects are related to differences in attentional engagement or behavioral responses to the AOD task. Future studies should to investigate patterns of cfM while subjects are resting or engaged in other behaviors. It will also be important to determine whether distinctions in cfM are specific to schizophrenia patients, or whether they represent a more general marker of altered neural coordination that may common to other neuropsychiatric disorders (Herrmann and Demiralp, 2005).

Genetic association analysis identified several SNPs co-varying with schizophrenia-relevant cfM components. It is interesting to note that none showed a direct association with diagnosis in the same set of subjects (see Table 2), supporting the biomarker approach toward elucidating the genetics of complex disorders (Braff et al., 2007). While this study is the first to identify genetic contributions to cross-frequency coupling, it complements a large body of work demonstrating genetic influences on cortical oscillations (Begleiter and Porjesz, 2006). Investigations of event-related potentials and resting-state dynamics implicate many genes within dopaminergic, adrenergic, cholinergic, GABAergic and glutamatergic pathways (Begleiter and Porjesz, 2006; Liu et al., 2009a).

In the current study, we observed an association between the magnitude of global cfM (IC_M 1) and the rs6277 genotype. The rs6277 polymorphism, also known as C957T, is a synonymous mutation in exon 7 of the DRD2 gene. In both patients and controls, the strength of cross-frequency coupling decreased with copies of the 957C allele. Consistent with the suggested role of cfM in neural coordination, presence of the 957C allele has been associated with poorer performance on executive function (Rodriguez-Jimenez et al., 2006), working memory (Jacobsen et al., 2006), and avoidance learning (Frank and Hutchison, 2009) in healthy controls. Furthermore, several case-control association studies have supported an over-representation of the 957C allele in schizophrenia (Lawford et al., 2005; Hänninen et al., 2006; Betcheva et al., 2009).

We also found an association between the weights of cfM component 4 and several SNPs within/near the GABRA2 gene; these were intronic (rs279869, rs279837), synonymous (rs279858), or in flanking regions of DNA (rs567926). Though their functional effects are unknown, these SNPs (or others in high LD) may affect the expression of the GABA_A receptor (Ittiwut et al., 2007). Because correlations between GABRA2 genotypes and IC_M 4 were observed only in SZ (see Figure 6E, Figure A2B), the GABRA2 genotype may only affect cross-frequency coupling only when neuronal circuitry is altered. Evidence of diminished GABAergic transmission in schizophrenia patients (Lewis et al., 2005; Bullock et al., 2008) supports this conjecture. Differences include decreased expression of enzymes regulating GABA synthesis and re-uptake, and compensatory increased expression of post-synaptic GABA_A, α2 subunit receptors (Lewis et al., 2005). Recently, the use of α2 selective GABA_A agonists to normalize GABAergic tone has been proposed, with preliminary evidence of improved cognitive function during treatment (Lewis et al., 2008). Our findings make important predictions regarding the success of these clinical trials, as they imply that the therapeutic effects of GABA_A agonists will be stratified over GABRA2 genotypes. Theoretically, patients with the most “risk” GABRA2 alleles (see Figure 6F) would receive the greatest benefit from this type of pharmacological intervention.

Limitations and Future Work

The results presented here must be interpreted in the context of several methodological limitations. First, because our analysis combined data from a number of related studies and constitutes a “convenience sample,” demographics were not well matched between HC and SZ groups, particularly with regard to age. We attempted to correct for this difference by linearly regressing age from dependent variables prior to group comparisons, however distinctions in modulation patterns between HC and SZ may still be impacted by age. To see if this was the case, we compared the original IC_M loading parameters in a subset of age-matched SZ (n = 42) and HC (n = 42; HC age: 36.1 ± 11.4; SZ age: 36.2 ± 11.6, t₈₂ = 0.03, P = 0.97), and found that the group differences persisted (IC_M 1: F = 4.16, P < 0.05; IC_M 4: F = 8.46, P < 0.005). While this suggests that age is not the sole factor contributing to the group difference, future investigations should use better-matched groups to address this issue and other possible confounds.

A second limitation is the possible effect of medication. As with any study of a clinical population, our findings may be influenced by the medication status of patients at the time of evaluation as well as life-long medication history. Here, there is a particular concern that prescribed antipsychotic medications often target the dopamine receptors identified in our genetic association analysis. While our patient sample size is too small (and available medication information too sparse) to fully assess these possible confounds, we have attempted to address concerns with what information we do have. Regarding the main effect of cfM group differences, the loading parameters of schizophrenia-relevant cfM components IC_M 1 and IC_M 4 are not significantly different between SZ patients taking (35/47) and those not taking antipsychotics (12/47) at the time of assessment (two-tailed t-test, IC_M 1: t₄₅ = 0.95, P = 0.34; IC_M 4: t₄₅ = −0.31, P = 0.76; permutation test, IC_M 1: P = 0.65; IC_M 4: P = 0.24; 10,000 permutations). Also, as described above, subjects with first-degree SZ relatives (REL) exhibited patterns of cfM that were intermediate between HC and SZ. Regarding the genetic association findings, both HC and SZ show the same trend between DRD2 genotype and the loading parameters of IC_M 1 (Figure 6C and related text), implying that medication does not affect this relationship. Additionally, for SZ patients with genotype data (36/47), we found no difference in cfM-associated genetic components between patients taking (28/36) and not taking antipsychotics (8/36). This was true whether we examined the loading parameters of genetic components (two-tailed t-test, IC_G 6: t₃₄ = 1.01, P = 0.32; IC_G 4: t₃₄ = −0.54, P = 0.52; permutation test, IC_G 6: P = 0.69, IC_G 4: P = 0.48) or the genotype distributions directly (rs6277: χ² = 1.46, P = 0.48; rs279869: χ² = 0.56, P = 0.76). These analyses support an interpretation that findings are not completely due to medication status, though we again note that with such small samples we lack sufficient statistical power to perform a conclusive investigation and the influence of medication must be addressed in future work.

A third limitation is the ability to detect and quantify cfM. Most previous studies have used intra-cranial recordings (though see e.g., Demiralp et al., 2007; Osipova et al., 2008) to assess interactions in oscillatory patters because these methods are much less susceptible to artifacts and are not limited by the severe attenuation of high-frequency power at the scalp. Though EEG and MEG recordings have the obvious benefit of being non-invasive, true patterns of cfM may be confounded by sources related to eye-movement (e.g., Yuval-Greenberg et al., 2008), scalp musculature (e.g., Whitham et al., 2008), or signal artifacts (e.g., Kramer et al., 2008), or simply may be impossible to detect at the scalp. Additional preprocessing with ICA can be used to mitigate the influence of these effects; however they remain a concern, especially in high-frequency bands where contributions are very large relative to neural sources. Additionally, we must consider the exact approach used to calculate cfM. As reviewed by Canolty and Knight (2010), there are numerous proposed methods to assess relationships between phase and amplitude, and as yet there is no consensus on a “gold standard” that performs best in all scenarios. Because the approach used here (modulation index based on the mean vector length, see Figure 1) is sensitive to the strength of modulation rather than just the presence or absence, we believe it to be appropriate for our purpose of comparing between groups and conditions, though other methods may have different sensitivities based on their assumptions of phase consistency (Penny et al., 2008; Tort et al., 2010).

Conclusion

Exploratory in nature, this study presents a data-driven approach to examine and compare cross-frequency interactions. We propose the application of ICA to discover independent features of cfM with no specific assumptions regarding their spatial or spectral composition. With this approach, we have successfully identified modulation features that are (1) similar to patterns observed in previous studies, (2) present with different strengths in healthy controls and schizophrenia patients, and (3) associated with genetic polymorphisms. Though there is much more work to be done, e.g., evaluating patterns of modulation during different mental states and validating group differences, we consider this a first and exciting step in the global examination of oscillatory interactions and network coordination in the human brain.

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

This work was supported by the National Institutes of Health under grants R01 EB005846, R01 EB006841, and P20 RR021938 to Vince D. Calhoun, grants R01 MH43775, R01 MH52886 and R01 MH074797 to Godfrey D. Pearlson, and grant R01 MH072681 to Kent A. Kiehl. We thank Mark Kramer of Boston University for providing the MATLAB code to calculate bicoherence and the reviewers for providing helpful comments that improved the manuscript.

References

Allen, A., Griss, M., Folley, B., Hawkins, K., and Pearlson, G. (2009). Endophenotypes in schizophrenia: a selective review. Schizophr. Res. 109, 24–37.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Andreasen, N., Paradiso, S., and O’Leary, D. (1998). “Cognitive dysmetria” as an integrative theory of schizophrenia: a dysfunction in cortical-subcortical-cerebellar circuitry? Schizophr. Bull. 24, 203.

Pubmed Abstract | Pubmed Full Text

Asghari, V., Sanyal, S., Buchwaldt, S., Paterson, A., Jovanovic, V., and Van Tol, H. (1995). Modulation of intracellular cyclic AMP levels by different human dopamine D4 receptor variants. J. Neurochem. 65, 1157–1165.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Begleiter, H., and Porjesz, B. (2006). Genetics of human brain oscillations. Int. J. Psychophysiol. 60, 162–171.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bell, A. J., and Sejnowski, T. J. (1995). An information-maximization approach to blind separation and blind deconvolution. Neural. Comput. 7, 1129–1159.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Betcheva, E. T., Mushiroda, T., Takahashi, A., Kubo, M., Karachanak, S. K., Zaharieva, I. T., Vazharova, R. V., Dimova, I. I., Milanova, V. K., and Tolev, T. (2009). Case–control association study of 59 candidate genes reveals the DRD2 SNP rs6277 (C957T) as the only susceptibility factor for schizophrenia in the Bulgarian population. J. Hum. Genet. 54, 98–107.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Braff, D., Freedman, R., Schork, N., and Gottesman, I. (2007). Deconstructing schizophrenia: an overview of the use of endophenotypes in order to understand a complex disorder. Schizophr. Bull. 33, 21.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bruns, A., and Eckhorn, R. (2004). Task-related coupling from high-to low-frequency signals among visual cortical areas in human subdural recordings. Int. J. Psychophysiol. 51, 97–116.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bullock, W., Cardon, K., Bustillo, J., Roberts, R., and Perrone-Bizzozero, N. (2008). Altered expression of genes involved in GABAergic transmission and neuromodulation of granule cell activity in the cerebellum of schizophrenia patients. Am. J. Psychiatry 165, 1594.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Calhoun, V. D., Liu, J., and AdalI, T. (2009). A review of group ICA for fMRI data and ICA for joint inference of imaging, genetic, and ERP data. Neuroimage 45, S163–S172.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Canolty, R., Edwards, E., Dalal, S., Soltani, M., Nagarajan, S., Kirsch, H., Berger, M., Barbaro, N., and Knight, R. (2006). High gamma power is phase-locked to theta oscillations in human neocortex. Science 313, 1626–1628.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Canolty, R. T., Ganguly, K., Kennerley, S. W., Cadieu, C. F., Koepsell, K., Wallis, J. D., and Carmena, J. M. (2010). Oscillatory phase coupling coordinates anatomically dispersed functional cell assemblies. Proc. Natl. Acad. Sci. 107, 17356–17361.

CrossRef Full Text

Canolty, R. T., and Knight, R. T. (2010). The functional role of cross-frequency coupling. Trends Cogn. Sci. 14, 506–515.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cho, R., Konecky, R., and Carter, C. (2006). Impairments in frontal cortical synchrony and cognitive control in schizophrenia. Proc. Natl. Acad. Sci. 103, 19878.

CrossRef Full Text

Cohen, M., Axmacher, N., Lenartz, D., Elger, C., Sturm, V., and Schlaepfer, T. (2009). Good vibrations: cross-frequency coupling in the human nucleus accumbens during reward processing. J. Cogn. Neurosci. 21, 875–889.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

de Lange, F., Jensen, O., Bauer, M., and Toni, I. (2008). Interactions between posterior gamma and frontal alpha/beta oscillations during imagined actions. Front. Hum. Neurosci. 2:7. doi: 10.3389/neuro.09.007.2008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

De Ridder, F., Pintelon, R., Schoukens, J., and Gillikin, D. P. (2005). Modified AIC and MDL model selection criteria for short data records. IEEE Trans. Instrum. Meas. 54, 144–150.

CrossRef Full Text

Delorme, A., and Makeig, S. (2004). EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Demiralp, T., Bayraktaroglu, Z., Lenz, D., Junge, S., Busch, N., Maess, B., Ergen, M., and Herrmann, C. (2007). Gamma amplitudes are coupled to theta phase in human EEG during visual perception. Int. J. Psychophysiol. 64, 24–30.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Flynn, G., Alexander, D., Harris, A., Whitford, T., Wong, W., Galletly, C., Silverstein, S., Gordon, E., and Williams, L. (2008). Increased absolute magnitude of gamma synchrony in first-episode psychosis. Schizophr. Res. 105, 262–271.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Frank, M., and Hutchison, K. (2009). Genetic contributions to avoidance-based decisions: striatal D2 receptor polymorphisms. Neuroscience 164, 131–140.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fries, P. (2005). A mechanism for cognitive dynamics: neuronal communication through neuronal coherence. Trends Cogn. Sci. 9, 474–480.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fries, P., Nikoli, D., and Singer, W. (2007). The gamma cycle. Trends Neurosci. 30, 309–316.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Friston, K. (1998). The disconnection hypothesis. Schizophr. Res. 30, 115–125.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fuke, S., Suo, S., Takahashi, N., Koike, H., Sasagawa, N., and Ishiura, S. (2001). The VNTR polymorphism of the human dopamine transporter (DAT1) gene affects gene expression. Pharmacogenomics J. 1, 152.

Pubmed Abstract | Pubmed Full Text

Gonzalez-Burgos, G., and Lewis, D. A. (2008). GABA neurons and the mechanisms of network oscillations: implications for understanding cortical dysfunction in schizophrenia. Schizophr. Bull. 34, 944.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Händel, B., and Haarmeier, T. (2009). Cross-frequency coupling of brain oscillations indicates the success in visual motion discrimination. Neuroimage 45, 1040–1046.

CrossRef Full Text

Hänninen, K., Katila, H., Kampman, O., Anttila, S., Illi, A., Rontu, R., Mattila, K., Hietala, J., Hurme, M., and Leinonen, E. (2006). Association between the C957T polymorphism of the dopamine D2 receptor gene and schizophrenia. Neurosci. Lett. 407, 195–198.

CrossRef Full Text

He, B. J., Zempel, J. M., Snyder, A. Z., and Raichle, M. E. (2010). The temporal structures and functional significance of scale-free brain activity. Neuron 66, 353–369.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Herrmann, C., and Demiralp, T. (2005). Human EEG gamma oscillations in neuropsychiatric disorders. Clin. Neurophysiol. 116, 2719–2733.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Himberg, J., Hyvärinen, A., and Esposito, F. (2004). Validating the independent components of neuroimaging time series via clustering and visualization. Neuroimage 22, 1214–1222.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hu, X., Lipsky, R., Zhu, G., Akhtar, L., Taubman, J., Greenberg, B., Xu, K., Arnold, P., Richter, M., and Kennedy, J. (2006). Serotonin transporter promoter gain-of-function genotypes are linked to obsessive-compulsive disorder. Am. J. Hum. Genet. 78, 815–826.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hyvärinen, A., Karhunen, J., and Oja, E. (2001). Independent Component Analysis. New York: John Wiley & Sons.

Ittiwut, C., Listman, J., Mutirangura, A., Malison, R., Covault, J., Kranzler, H., Sughondhabirom, A., Thavichachart, N., and Gelernter, J. (2007). Interpopulation linkage disequilibrium patterns of GABRA2 and GABRG1 genes at the GABA cluster locus on human chromosome 4. Genomics 91, 61–69.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jacobsen, L., Pugh, K., Mencl, W., and Gelernter, J. (2006). C957T polymorphism of the dopamine D2 receptor gene modulates the effect of nicotine on working memory performance and cortical processing efficiency. Psychopharmacology 188, 530–540.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kramer, M., Tort, A., and Kopell, N. (2008). Sharp edge artifacts and spurious coupling in EEG frequency comodulation measures. J. Neurosci. Methods 170, 352–357.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lakatos, P., Karmos, G., Mehta, A., Ulbert, I., and Schroeder, C. (2008). Entrainment of neuronal oscillations as a mechanism of attentional selection. Science 320, 110.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lakatos, P., Shah, A., Knuth, K., Ulbert, I., Karmos, G., and Schroeder, C. (2005). An oscillatory hierarchy controlling neuronal excitability and stimulus processing in the auditory cortex. J. Neurophysiol. 94, 1904–1911.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lawford, B., Young, R., Swagell, C., Barnes, M., Burton, S., Ward, W., Heslop, K., Shadforth, S., van Daal, A., and Morris, C. (2005). The C/C genotype of the C957T polymorphism of the dopamine D2 receptor is associated with schizophrenia. Schizophr. Res. 73, 31–37.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lee, K., Williams, L., Breakspear, M., and Gordon, E. (2003). Synchronous gamma activity: a review and contribution to an integrative neuroscience model of schizophrenia. Brain Res. Rev. 41, 57–78.

CrossRef Full Text

Lewis, D., Cho, R., Carter, C., Eklund, K., Forster, S., Kelly, M., and Montrose, D. (2008). Subunit-selective modulation of GABA type A receptor neurotransmission and cognition in schizophrenia. Am. J. Psychiatry 165, 1585.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lewis, D., Hashimoto, T., and Volk, D. (2005). Cortical inhibitory neurons and schizophrenia. Nat. Rev. Neurosci. 6, 312–324.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lisman, J., and Buzsáki, G. (2008). A Neural Coding Scheme Formed by the Combined Function of Gamma and Theta Oscillations. Schizophr. Bull. 34, 974.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, J., Kiehl, K., Pearlson, G., Perrone-Bizzozero, N., Eichele, T., and Calhoun, V. (2009a). Genetic determinants of target and novelty-related event-related potentials in the auditory oddball response. Neuroimage 46, 809–816.

CrossRef Full Text

Liu, J., Pearlson, G., Windemuth, A., Ruano, G., Perrone Bizzozero, N. I., and Calhoun, V. (2009b). Combining fMRI and SNP data to investigate connections between brain function and genetics using parallel ICA. Hum. Brain Mapp. 30, 241–255.

CrossRef Full Text

Maris, E., van Vugt, M., and Kahana, M. (2011). Spatially distributed patterns of oscillatory coupling between high-frequency amplitudes and low-frequency phases in human iEEG. Neuroimage 54, 836–850.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Miller, K. J., Hermes, D., Honey, C. J., Sharma, M., Rao, R. P. N., Den Nijs, M., Fetz, E. E., Sejnowski, T. J., Hebb, A. O., and Ojemann, J. G. (2010). Dynamic modulation of local population activity by rhythm phase in human occipital cortex during a visual search task. Front. Hum. Neurosci. 4:197. doi: 10.3389/fnhum.2010.00197

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Monto, S., Palva, S., Voipio, J., and Palva, J. M. (2008). Very slow EEG fluctuations predict the dynamics of stimulus detection and oscillation amplitudes in humans. J. Neurosci. 28, 8268–8272.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mormann, F., Fell, J., Axmacher, N., Weber, B., Lehnertz, K., Elger, C., and Fernández, G. (2005). Phase/amplitude reset and theta-gamma interaction in the human medial temporal lobe during a continuous word recognition memory task. Hippocampus 15, 890–900.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

O’Keefe, J., and Recce, M. (1993). Phase relationship between hippocampal place units and the EEG theta rhythm. Hippocampus 3, 317–330.

CrossRef Full Text

Osipova, D., Hermes, D., and Jensen, O. (2008). Gamma power is phase-locked to posterior alpha activity. PLoS ONE 3, 3990. doi: 10.1371/journal.pone.0003990

CrossRef Full Text

Penny, W., Duzel, E., Miller, K., and Ojemann, J. (2008). Testing for nested oscillation. J. Neurosci. Methods 174, 50–61.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ray, S., and Maunsell, J. H. R. (2011). Different origins of gamma rhythm and high-gamma activity in macaque visual cortex. PLoS Biol. 9, e1000610. doi: 10.1371/journal.pbio.1000610

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rodriguez-Jimenez, R., Hoenicka, J., Jimenez-Arriero, M., Ponce, G., Bagney, A., Aragues, M., and Palomo, T. (2006). Performance in the Wisconsin Card Sorting Test and the C957T polymorphism of the DRD2 gene in healthy volunteers. Neuropsychobiology 54, 166–170.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Roopun, A., Kramer, M., Carracedo, L., Kaiser, M., Davies, C., Traub, R., Kopell, N., and Whittington, M. (2008). Temporal interactions between cortical rhythms. Front. Neurosci. 2:2. doi: 10.3389/neuro.01.034.2008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sachs, L. (1992). Angewandte Statistik: Anwendung statistischer Methoden, 7th Edn. Berlin: Springer.

Sirota, A., Montgomery, S., Fujisawa, S., Isomura, Y., Zugaro, M., and Buzsáki, G. (2008). Entrainment of neocortical neurons and gamma oscillations by the hippocampal theta rhythm. Neuron 60, 683–697.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Spencer, K. M., Niznikiewicz, M. A., Shenton, M. E., and McCarley, R. W. (2008). Sensory-evoked gamma oscillations in chronic schizophrenia. Biol. Psychiatry 63, 744–747.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stein, M., Seedat, S., and Gelernter, J. (2006). Serotonin transporter gene promoter polymorphism predicts SSRI response in generalized social anxiety disorder. Psychopharmacology 187, 68–72.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Symond, M., Harris, A., Gordon, E., and Williams, L. (2005). Gamma synchrony in first-episode schizophrenia: a disorder of temporal connectivity? Am. J. Psychiatry 162, 459–465.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Tekell, J., Hoffmann, R., Hendrickse, W., Greene, R., Rush, A., and Armitage, R. (2005). High frequency EEG activity during sleep: characteristics in schizophrenia and depression. Clin. EEG Neurosci. 36, 25–35.

Pubmed Abstract | Pubmed Full Text

Tort, A., Komorowski, R., Eichenbaum, H., and Kopell, N. (2010). Measuring phase-amplitude coupling between neuronal oscillations of different frequencies. J. Neurophysiol. 104, 1195–1210.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Uhlhaas, P. J., Haenschel, C., Nikolic, D., and Singer, W. (2008). The role of oscillations and synchrony in cortical networks and their putative relevance for the pathophysiology of schizophrenia. Schizophr. Bull. 34, 927–943.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

von Stein, A., and Sarnthein, J. (2000). Different frequencies for different scales of cortical integration: from local gamma to long range alpha/theta synchronization. Int. J. Psychophysiol. 38, 301–313.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Voytek, B., Canolty, R. T., Shestyuk, A., Crone, N. E., Parvizi, J., and Knight, R. T. (2010). Shifts in gamma phase–amplitude coupling frequency from theta to alpha over posterior cortex during visual tasks. Front. Hum. Neurosci. 4:191. doi: 10.3389/fnhum.2010.00191

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Whitham, E. M., Lewis, T., Pope, K. J., Fitzgibbon, S. P., Clark, C. R., Loveless, S., DeLosAngeles, D., Wallace, A. K., Broberg, M., and Willoughby, J. O. (2008). Thinking activates EMG in scalp electrical recordings. Clin. Neurophysiol. 119, 1166–1175.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yuval-Greenberg, S., Tomer, O., Keren, A. S., Nelken, I., and Deouell, L. Y. (2008). Transient induced gamma-band response in EEG as a manifestation of miniature saccades. Neuron 58, 429–441.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text