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 (F1,175 = 9.25, P < 0.005), while modulation at fronto-temporal electrodes was greater in patients (F1,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 GABAA 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 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  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(Lakatos et al., , 2008Mormann 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 datadriven 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.

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 largescale 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. fication 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.

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; see 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), m zt raw = | ( )| . 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 , we examined high-frequency peak-locked averages of 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 −7 ). 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).

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 ampli-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 lowfrequency phase bands (f P = 1-24 Hz, width = 4 Hz) and 15 highfrequency 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 187 = 2.04, P = 0.04; IC M 2: t 187 = 2.99,

GenBank accession number. c Tri-allelic locus with long variants (L A and L G ) and short variant (S).
Allen et al.
Cross-frequency modulation in health and disease
Before performing ICA on the genetic data, we imputed missing values (<5% of all genotypes). For loci in linkage disequilibrium (LD; r 2 > 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 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).

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: (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), m zt raw = | ( ) |, 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. 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 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 Ŝ M (j) (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 highfrequency 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 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).

Allen et al.
Cross-frequency modulation in health and disease Frontiers in Systems Neuroscience www.frontiersin.org small number of REL subjects, the average loading parameters of IC M 1 were significantly different between HC and REL (t 140 = 2.27, P < 0.05) and showed a trend for IC M 4 (t 140 = 1.68, P = 0.10). Loading parameters of components 1 and 4 were also significantly negatively correlated (r = −0.56, P < 10 −15 , 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.

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 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 frontotemporal 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 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 IC M 1 and IC M 4 and the loading parameters of the IC G s, 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 for 26 polymorphisms). For IC G 4, SNPs rs279869, rs279858, and rs567926, all in relatively high LD with one another (r 2 > 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 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 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.

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

FiGuRE 4 | Relationships between cross-frequency modulation loading parameters. (A) Full correlation matrix between loading parameters (A
Allen et al. Cross-frequency modulation in health and disease Frontiers in Systems Neuroscience www.frontiersin.org ing 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 subbands, 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 eventrelated 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 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(Lakatos et al., , 2008Monto 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 phaseamplitude 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 Maris et al., 2011). The differential modulation of cross-frequency components by condition supports this view: component load-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 82 = 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 45 = 0.95, P = 0.34; IC M 4: t 45 = −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 34 = 1.01, P = 0.32; IC G 4: t 34 = −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: χ 2 = 1.46, P = 0.48; rs279869: χ 2 = 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 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 . 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 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.

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. 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 noninvasive, 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). aPPendIx FiGuRE A1 | Correlation structure for the genotyped polymorphisms. (A,B) r 2 values computed between all pairs of loci. High correlations indicate data redundancy, which is reduced through ICA. Polymorphic loci in close physical proximity are bounded with black boxes. Correlation structure is displayed for all subjects in (A), and in (B) is computed separately for the 3 largest subsets of subjects identifying as Caucasian (left, n = 118), African American (middle, n = 21), and Hispanic (right, n = 16). Correlation structure is similar between groups, though differences in magnitude are evident.
Allen et al.
Cross-frequency modulation in health and disease Frontiers in Systems Neuroscience www.frontiersin.org

FiGuRE A2 | Genetic contributions to cross-frequency modulation. (A)
Genetic component association results, including race in the regression model. 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 (Left) and 4 (Right). The pink shaded region indicates the uncorrected P-values for the full model, while the black, white, and blue regions show the P-values for the genetic, genetic × diagnosis, and race terms, respectively. Asterisks denote significant models between (1) IC M 1 and IC G 6 (Left, pink; F 4,165 = 4.88, P < 0.001), (2) IC M 4 and IC G 6 (Right, pink; F 4,165 = 7.18, P < 0.00005), and (3) IC M 4 and IC G 4 (Right, pink; F 4,165 = 6.23, P < 0.0005). Associations were all significant at P < 0.05 after Bonferroni correction for 32 tests. Examination of the individual beta weights revealed a main effect of genetic component 6 on IC M 1 (Left, black; t 165 = 3.05 , P < 0.005). Associations with IC M 4 were both due to interactions between the genetic components and the diagnosis (Right, white; t 165 = 3.30 , P < 0.005 and t 165 = 3.18 , P < 0.005 for IC G 6 and IC G 4, respectively). No significant associations with race were found. (B) Plots of modulation component loading parameters as a function of GABRA2 genotype. Data is stratified by SZ (red) and HC (blue) to highlight group × diagnosis interactions. Error bars denote ± 1 SEM. 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.