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).

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).

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.

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.

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.

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), mraw = |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).

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 S^M(j) (Himberg et al., 2004).

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.

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.

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.

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.

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.

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).

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.