Cross-Frequency Power-Power Coupling Analysis: A Useful Cross-Frequency Measure to Classify ICA-Decomposed EEG

Magneto-/Electro-encephalography (M/EEG) commonly uses (fast) Fourier transformation to compute power spectral density (PSD). However, the resulting PSD plot lacks temporal information, making interpretation sometimes equivocal. For example, consider two different PSDs: a central parietal EEG PSD with twin peaks at 10 Hz and 20 Hz and a central parietal PSD with twin peaks at 10 Hz and 50 Hz. We can assume the first PSD shows a mu rhythm and the second harmonic; however, the latter PSD likely shows an alpha peak and an independent line noise. Without prior knowledge, however, the PSD alone cannot distinguish between the two cases. To address this limitation of PSD, we propose using cross-frequency power–power coupling (PPC) as a post-processing of independent component (IC) analysis (ICA) to distinguish brain components from muscle and environmental artifact sources. We conclude that post-ICA PPC analysis could serve as a new data-driven EEG classifier in M/EEG studies. For the reader’s convenience, we offer a brief literature overview on the disparate use of PPC. The proposed cross-frequency power–power coupling analysis toolbox (PowPowCAT) is a free, open-source toolbox, which works as an EEGLAB extension.


Background
Power spectral density (PSD) calculations using fast Fourier transform (FFT) are ubiquitous in magnetoencephalography (MEG) and electroencephalography (EEG) analyses. However, this approach's drawback is the lack of temporal resolution; hence, sliding window solutions such as short-term Fourier Transform and wavelet transform become necessary for temporal resolution.
Cross-frequency power-power coupling (PPC) addresses the temporal resolution dilemma by calculating spectral covariance. When the amplitudes of two oscillations in different frequencies covaries along with time, PPC occurs. The resulting spectral covariance is a square matrix plot (also called a comodulogram [1][2][3][4][5][6][7][8] or a comodugram [9][10][11][12][13][14][15]. The following examples illustrate two different types of cross-frequency coupling, the distinction of which is of common interest in EEG and MEG studies. Imagine two peaks occurring in a PSD plot at X Hz and at Y Hz. Power changes in the two oscillations causing the two peaks may or may not be correlated. Consider the first case where a PSD calculated from data recorded on or adjacent to the sensorimotor area shows a mu rhythm with the first peak within the alpha band and the second harmonic in the beta band (similar to the Greek alphabet

Goal of the Current Study
The study's goal is to demonstrate the usefulness of cross-frequency PPC measure within EEG data analysis, specifically as a post-independent component analysis (ICA) [26][27][28][29][30]. One of the post-ICA processing challenges is classifying the obtained independent components (ICs). Many applications aim for automated labeling of the ICs including MARA [31,32], ADJUST [33], FASTER [34], eyeCatch [35], IC MARC [36], SASICA [37], ICLabel [38], and others [39,40]. Multiple solutions exist because IC classification is a difficult task. The properties of ICA-decomposed EEG data may be represented by event-related potential (ERPs), event-related power-spectral perturbation, inter-trial phase coherence (ERSP and ITC [41]), scalp topography, and dipole locations. Adding different types of EEG property that carry unique information not represented by these cited measures improves the characterization of the obtained ICA-decomposed EEG data. Here, we demonstrate that brain, eye, and muscle classes of ICA-decomposed scalp-recorded signals have characteristic PPC readily visible in comodulograms. The visible features of PPC provide complimentary temporal information useful for classification and labeling of ICs. The proposed method is available as a plugin Power-Power Coupling Analysis Tool (PowPowCAT) (https://sccn.ucsd.edu/wiki/Power-Power_Coupling_Analysis_Tool). The plugin is for a free, open-source library EEGLAB [42].

Data Acquisition and Preprocessing
The EEG data used is from a previously published study on facial familiarity [43]. Twenty-four paid volunteers participated in the study (mean 21.1 ± 1.4 years old, range 20-26, 11 women). All participants gave written informed consent. All experimental procedures followed the Declaration of Helsinki. Visual stimuli with three levels of familiarity (self-relevant, familiar, and unfamiliar) and with two levels of categories (face and cup) were presented. The participant's task was to judge the level of object familiarity by pressing one of three buttons as quickly as possible. The stimulus was presented either unilaterally (to the left or right visual hemifield) or bilaterally (identical images presented to both visual hemifields). During the task, EEG signals were recorded using the Nihon Koden EEG 1100 (Tokyo, Japan) with 33 electrodes (AFp9, Fp1, Fp2, Fz, F3, F4, F7, F8, FC1, FC2, FC5, FC6, Cz, C3, C4, T7, T8, CP1, CP2, CP5, CP6, TP9, TP10, Pz, P3, P4, P7, P8, O1, O2, PO9, PO10, and Iz), positioned in correspondence with international 10-10 system, at the sampling frequency of 500 Hz, which was later downsampled to 200 Hz. Electrode impedance was kept typically below 5 kΩ. EEG electrodes were initially referenced to a nose tip, and later, a common average re-referencing was performed to gain better estimates of reference-independent potentials [44]. An FIR high-pass filter of 0.5 Hz @ −6dB (transient bandwidth 1 Hz, Hamming windowed) [45] was applied to the continuous EEG data. Finally, adaptive mixture ICA [46] was performed to obtain 32 ICs from each of the 24 datasets as it is most efficient among all decomposition methods in reducing mutual information in the component time courses [47]. For each scalp topography of ICs, which are the column of the mixing matrix (i.e., the inverse of the weight matrix), the estimated equivalent current dipole was fit using Dipfit2.2. Dipfit2.2 is a group of functions taken from Fieldtrip Toolbox [48]. It uses a four-shell spherical head model co-registered to each subject's electrode locations by warping the electrode locations to the template head from Montreal Neurological Institute (MNI).
We chose ICA as a preprocessing technique. ICA relies on these three assumptions: (1) Functional independence of cortical regions should be represented by temporal independence of the local field potentials; (2) volume conduction (i.e., signal paths from an EEG source to scalp sensors) and signal mixing process are linear and instantaneous; and (3) the functional mapping of the neocortex is stationary across the recording. If these assumptions hold, the ICA's spatial filter produces (1) effective sources of EEG as temporally maximally independent signals as well as (2) their linear and instantaneous mixing processes (i.e., volume conduction) that are (3) stationary. The beauty of ICA in application to scalp-recorded EEG is that ICA solves the spatial problem without requiring any spatial information.
Sensors 2020, 20, 7040 4 of 16 ICA does not require either electric forward model of the human head or knowledge of scalp channel locations. Onton and Makeig [30] established the physiological interpretation of ICA on EEG. The only confirmed case of violation of the nonstationary assumption is traveling wave, which is thought of as moving EEG sources in the brain. In [30], the authors simulated that a traveling wave spreading across 3 cm of distance with the speed of 2 m/s oscillating in 10 Hz produces 27 degrees of phase shift, which is only 7.5% of a single cycle of 10 Hz. They determined that this level of phase shift could be trivial and regarded as nearly synchronous for ICA. However, traveling waves with larger scales, such as the one reported by Muller and colleagues [49], would result in failure in applying ICA. Nonetheless, previous studies on EEG experiments during awake states reported the successful application of ICA [28][29][30][50][51][52].
For the group-level data integration, independent component (IC) clustering using the k-means algorithm, implemented in EEGLAB STUDY for grouping similar ICs across subjects, after PCA-dimension reduction was used. The equivalent current dipoles were normalized with moments and the weight factor was empirically determined as 5. PSD and scalp topography were also precomputed for each IC, and they were dimension-reduced to 10 for each using PCA. Next, k-means clustering was applied on the dimension-reduced data feature vectors to categorize 790 ICs into 31 clusters (the number of clusters was based on [53]). Each cluster was examined to label as brain, muscle, eye, and other noise by visual inspection on PSD, equivalent current dipole locations, and scalp topography. As a result, 5 of 31 clusters were identified as eye-movement IC clusters for their locations near the edge of ventral frontal brain regions and for their PSD being relatively high at low frequencies with no spectral peak; 7 of 31 clusters were identified as clusters associated with muscle artifacts for their locations outside of the brain or near lower head region and for their characteristic mean spectral plateau above 25 Hz. Thus, 16/31 clusters were identified as clusters associated with brain activity. They showed PSD feature of the 1/f characteristic typical in EEG signals sometimes with peaks in the alpha band (8)(9)(10)(11)(12)(13). Finally, 3/31 clusters were excluded from the analysis as outliers after manually cleaning the other clusters by gathering cluster-inconsistent ICs due to limitations of the algorithm in clustering accuracy. Based on this initial classification, we combined each class of clusters to form aggregated clusters of brain activities (consisting of 380 ICs: mean 15.83 ± 3.60 ICs from each subject), eye movements (consisting of 87 ICs: mean 3.63 ± 1.47 ICs from each subject), and muscle activities (consisting of 213 ICs: mean 8.88 ± 3.66 ICs from each subject), and outliers (consisting of 110 ICs: mean 4.58 ± 1.82 ICs from each subject). The cross-frequency PPC of each IC was computed and studied separately in each class.

Assessment of Cross-Frequency Power-Power Coupling
On continuous IC time-series data, Short-Term Fourier Transform (STFT) with a 50%-overlapping 1-s sliding window was applied to compute spectrograms using Matlab's spectrogram function to generate semi-logarithmically spaced 1-50 Hz divided into 100 frequency bins. Additionally, 20% of the individual PSD windows were excluded using distance from median values. This manner of cleaning the data preprocessing is justified in the Discussion section.
The correlation between the power spectrum S( f ) at a frequency f i and at another frequency f j can be calculated using the expression where S k ( f i ) is the PSD at frequency f i in time-window k, S( f i ) the averaged PSD at frequency f i overall sliding windows, σ i the standard deviation of the PSD at frequency f i , and k ranges over all sliding windows. To illustrate common patterns in comodulograms in each class, we generated a standardized matrix from correlation matrices of all ICs in the class. Given that a class contains m ICs, the standardized correlation at the frequency f i and at another frequency f j , namely z ij , is computed as where corr ij = m corr ij m is the averaged correlation across all m ICs at the frequency f i and at another frequency f j , and σ ij is the standard deviation of correlation. The coupling of any frequencies that is commonly found across ICs in the class is expected to be salient in the standardized matrix.

Visualizations by Power-Power Coupling Analysis Toolbox (PowPowCAT)
The matrix of correlation coefficients computed from a representative recording is shown in Figure 1. This figure is a screen capture from our in-house developed PowPowCAT (ver 1.00) graphical user interface (GUI). In the right large plot, the color scale shows that red indicates positive coupling (correlation), and the blue indicates negative coupling (anti-correlation). This panel also serves as an interactive interface to select two frequencies by moving the mouse cursor over the frequencies' another frequency , and is the standard deviation of correlation. The coupling of any frequencies that is commonly found across ICs in the class is expected to be salient in the standardized matrix.

Visualizations by Power-Power Coupling Analysis Toolbox (PowPowCAT)
The matrix of correlation coefficients computed from a representative recording is shown in Figure 1. This figure is a screen capture from our in-house developed PowPowCAT (ver 1.00) graphical user interface (GUI). In the right large plot, the color scale shows that red indicates positive coupling (correlation), and the blue indicates negative coupling (anti-correlation). This panel also serves as an interactive interface to select two frequencies by moving the mouse cursor over the frequencies' desired combination. The top left panel shows the mean PSD across all 1-s sliding windows, as well as its standard deviations. The middle left panel shows the vertical profile of the color map on the right large plot. The bottom left small panel shows the downsampled and normalized time-series of the PSD of the selected pair of the frequencies.
The bottom right small panel shows the scatter plot of the downsampled data points shown in the bottom left small panel.
To test the null hypothesis that the power spectra time series of two different frequencies, and , are not coupled in the data, non-parametric surrogate data method is supported as it makes minimum assumptions. It preserves the original data's statistical properties while generating time series that are randomized such that any possible nonlinear coupling is removed [54]. It has also been actively used to test coupling in nonlinear systems in many EEG studies [55,56]. Specifically, the surrogate data method randomizes time-window k differently for each frequency bin to build surrogate time-frequency time series and computes surrogate cross-frequency PPC. This process is repeated 5000 times (default) to produce distributions for the dataset in which the null hypothesis holds, i.e., the chance-level of the cross-frequency PPC. The original, non-permuted data are then compared to the surrogate distribution to obtain uncorrected p-values. The significance threshold can be selected to be either 0.05 or 0.01 after false discovery rate (FDR) correction across all ICs and all pairs of frequencies. Pixels with non-significant results are masked. To test the null hypothesis that the power spectra time series of two different frequencies, f i and f j , are not coupled in the data, non-parametric surrogate data method is supported as it makes minimum assumptions. It preserves the original data's statistical properties while generating time series that are randomized such that any possible nonlinear coupling is removed [54]. It has also been actively used to test coupling in nonlinear systems in many EEG studies [55,56]. Specifically, the surrogate data method randomizes time-window k differently for each frequency bin to build surrogate time-frequency time series and computes surrogate cross-frequency PPC. This process is repeated 5000 times (default) to produce distributions for the dataset in which the null hypothesis holds, i.e., the chance-level of the cross-frequency PPC. The original, non-permuted data are then compared to the surrogate distribution to obtain uncorrected p-values. The significance threshold can be selected to be either 0.05 or 0.01 after false discovery rate (FDR) correction across all ICs and all pairs of frequencies. Pixels with non-significant results are masked.

Clustering ICs Using Comodulogram as Clustering Criterion
To test whether comodulogram is useful for clustering the data, we re-clustered the ICs using their comodulograms as the clustering criterion (k-means algorithm without PCA-dimension reduction) to see if the clusters show physiologically meaningful source distributions. The three major classes of the scalp recorded signals, brain, eye, and muscle classes, have characteristic source distributions. Brain sources should be localized within the brain; eye sources within the peri-ocular regions; and muscle sources in the periphery of the anterior, temporal, and occipital regions.
After computing comodulograms of all 790 ICs, the non-redundant parts of the comodulograms were used as feature vectors. PCA dimension reduction was not applied for this analysis. In performing the same k-means algorithm, we set the number of clusters to be 8 to give it sufficient degrees of freedom to represent at least three classes of physiological signals (brain, eye, and muscle) plus noise. For each resulting cluster, the standardized correlation matrix was calculated from cluster members' ICs using Equation (2). In addition, the corresponding brain regions were illustrated by calculating the probabilistic density of the equivalent current dipoles using std_dipoleDensity plugin for EEGLAB that implemented 3-D Gaussian kernel transformation with full-width half maximum (FWHM) of 20 mm [47].
To test further if the comodulogram-based clustering produced reasonable classifications of data, we used ICLabel [38] to classify ICs into classes and then evaluated the proportion of ICs classes in each cluster resulted from comodulogram-based clustering. The ICLabel is an open-source, freely available plugin for EEGLAB that serves as a state-of-the-art IC classification tool. The classification model was trained on a large-scale ICs collection using a set of sophisticated features and labels given by human experts. It provides probabilistic class labels of Brain, Muscle, Eye, Heart, Line Noise, Channel Noise, and Others, and the highest class-label probability determines the class of that IC. We adopted this tool to generate the class labels for ICs. Thus, ICLabel served as the operational ground truth generator in the current study. We examined the percentage of the most prevailing class labels per cluster where higher probability indicates higher consistency with the operational ground truth. It should also be consistent with our visual evaluation on standardized comodulogram and probabilistic dipole density. The percentage of the most prevailing class label per cluster is called consistency level, indicating consistency between the results from ICLabel and our proposed comodulogram-based clustering. Moreover, among the total 790 ICs, we calculated the percentage of ICs classified by ICLabel into each class, regardless of the cluster, to deliver the distribution of classes. We found that the most prevalent class label is the brain class, and the percentage of this class was 42.5%. We adopted this number as the chance level of classification for comparison.

Standardized Comodulograms for Brain, Eye, and Muscle Classes
The grand comodulograms of standardized correlation for brain activity, eye movement, and muscle activity IC classes are plotted below (Figure 2) to confirm the general impressions. The strong couplings and the diagonal lines represent self-correlation and should be ignored. The visual examination of these comodulograms indicates that (1) anti-correlation (i.e., corr ij < 0) was rare.
(2) Differences across the classes are present: In the brain-activity IC class (left), the cross-frequency PPC was found between 6.5 and 30 Hz, and that in the off-diagonal area was weakly anti-correlated (see the blue color in the off-diagonal area); in the eye-movement IC class (Figure 2b), it showed characteristic block-diagonalization on the lower end of the frequencies in 1-9.4 Hz and the higher end of the frequencies in 23-50 Hz; lastly, in the muscle IC class (Figure 2c), characteristic parallel lines were observed between 30 and 50 Hz. In visual examination on PSD, such a high-frequency PPC is hard to find because they do not form a noticeable peak in the frequency domain. To sum, this highest-level result showed successful distinction at the level of visual inspection among the three classes of the EEG signals decomposed by ICA.   Figure 3 shows an example of a representative comodulogram from the brain-class IC (the equivalent current dipole coordinate was (−43, −4, 27) in Talairach coordinates, localized within left Brodmann Area 6), its normalized time series of the selected pair (by cursor position on the main comodulogram plot) of frequencies and , scatter plot of against , and the corresponding equivalent current dipole location. The IC in Figure 3 is the same as in Figure 1. The comodulogram showed alpha harmonics up to 4th order, though the corresponding PSD showed only 2nd harmonic. Such a cross-frequency coupled alpha oscillation was reported in the visual experiment known as resonance effects in the visual cortex [57]. However, as far as we know, this is the first report of 4th order harmonic of the human alpha-band oscillation. The results demonstrated that the proposed method allows us to reveal a cross-frequency coupling structure not available for standard PSD by FFT.  Figure 3 shows an example of a representative comodulogram from the brain-class IC (the equivalent current dipole coordinate was (−43, −4, 27) in Talairach coordinates, localized within left Brodmann Area 6), its normalized time series of the selected pair (by cursor position on the main comodulogram plot) of frequencies f i and f j , scatter plot of f i against f j , and the corresponding equivalent current dipole location. The IC in Figure 3 is the same as in Figure 1. The comodulogram showed alpha harmonics up to 4th order, though the corresponding PSD showed only 2nd harmonic. Such a cross-frequency coupled alpha oscillation was reported in the visual experiment known as resonance effects in the visual cortex [57]. However, as far as we know, this is the first report of 4th order harmonic of the human alpha-band oscillation. The results demonstrated that the proposed method allows us to reveal a cross-frequency coupling structure not available for standard PSD by FFT. equivalent current dipole location. The IC in Figure 3 is the same as in Figure 1. The comodulogram showed alpha harmonics up to 4th order, though the corresponding PSD showed only 2nd harmonic. Such a cross-frequency coupled alpha oscillation was reported in the visual experiment known as resonance effects in the visual cortex [57]. However, as far as we know, this is the first report of 4th order harmonic of the human alpha-band oscillation. The results demonstrated that the proposed method allows us to reveal a cross-frequency coupling structure not available for standard PSD by FFT.  Figure 4 shows an example of a representative comodulogram from the eye-class IC (the equivalent current dipole coordinate was (−1, 72, −5) in Talairach coordinates, localized outside the brain). The eye-class characteristic is the broad, within-low-frequency coupling (delta and theta bands, below 9.4 Hz on the plot) that composes a square on the bottom left of the diagonal line. This is in line with the report that blinks and eye movements generate large electrical activities-relative to brain-source EEG-within delta and theta frequency ranges but not in the alpha range and above [58]. However, this feature is often hard to distinguish from the EEG signal's 1/f distribution in PSD.   Figure 4 shows an example of a representative comodulogram from the eye-class IC (the equivalent current dipole coordinate was (−1, 72, −5) in Talairach coordinates, localized outside the brain). The eye-class characteristic is the broad, within-low-frequency coupling (delta and theta bands, below 9.4 Hz on the plot) that composes a square on the bottom left of the diagonal line. This is in line with the report that blinks and eye movements generate large electrical activities-relative to brain-source EEG-within delta and theta frequency ranges but not in the alpha range and above [58]. However, this feature is often hard to distinguish from the EEG signal's 1/f distribution in PSD.

Representative Comodulogram of Eye Class
It is hard to distinguish frontal brain component from eye component just by referring to PSD because of these overlapping features. Comparing the comodulograms of the brain class and the eye class in Figure 2 and between Figures 3 and 4, it is clear that PPC within the below-theta range is observed only in the eye class.  Figure 5 shows an example of a representative comodulogram from the muscle-class IC (the equivalent current dipole coordinate was (72, −10, −6) in Talairach coordinates, localized outside the brain). The comodulogram showed distinct parallel lines in the beta and gamma frequency ranges It is hard to distinguish frontal brain component from eye component just by referring to PSD because of these overlapping features. Comparing the comodulograms of the brain class and the eye class in Figure 2 and between Figures 3 and 4, it is clear that PPC within the below-theta range is observed only in the eye class. Figure 5 shows an example of a representative comodulogram from the muscle-class IC (the equivalent current dipole coordinate was (72, −10, −6) in Talairach coordinates, localized outside the brain). The comodulogram showed distinct parallel lines in the beta and gamma frequency ranges (20-40 Hz) as a second harmonic of 10-20 Hz, which is also visible in this class's standardized comodulogram in Figure 2. The known fact is that PSD of muscle activity lies below 40 Hz though it extends to 300 Hz [59]. However, it is the first time to reveal that muscle activity has characteristic cross-frequency PPC continuously from 10-20 Hz (coupled with 20-40Hz). Figure 5b shows the time-series of PSD changes at 16.0 Hz and 32.6 Hz, which confirmed a strong correlation (r = 0.504). To sum, we found a characteristic comodulogram pattern unique to the muscle class.   Figure 6 shows the standardized correlation matrices and the dipole density maps projected on anatomical images in axial view at different z coordinates in standard MNI space. Our visual inspection identified Cluster 2 as the brain class for the well-suppressed off-diagonal areas (near blue) and the frontocentral dipole density distribution. Cluster 5 was also identified as the brain class for the centroparietal dipole density distribution. However, Cluster 5's comodulogram showed an increase of correlations between 10 and 30 Hz, which differs from that of Cluster 2. We concluded that the comodulogram-based clustering identified (at least) two separate clear brain classes. Cluster 3 and Cluster 4 showed some brain class evidence due to their weak frontal dipole density distribution, but they also suggest dipole distribution within the frontal the right temporal areas. Their comodulograms showed that Cluster 3 is similar to Cluster 2, but Cluster 4's comodulogram is not. This qualitative evaluation concluded that Cluster 3 may consist of many brain sources, while Cluster 4 may consist of artifacts. Cluster 1 showed a block diagonal comodulogram centered on the low frequency (< 9.4 Hz). From what we learned above, we identified this cluster as the eye class. The dipole density distribution of the cluster supported this interpretation. Cluster 6, Cluster 7, and Cluster 8 showed block diagonal comodulograms centered on the higher frequency (above 13, 18,  Figure 6 shows the standardized correlation matrices and the dipole density maps projected on anatomical images in axial view at different z coordinates in standard MNI space. Our visual inspection identified Cluster 2 as the brain class for the well-suppressed off-diagonal areas (near blue) and the frontocentral dipole density distribution. Cluster 5 was also identified as the brain class for the centroparietal dipole density distribution. However, Cluster 5's comodulogram showed an increase of correlations between 10 and 30 Hz, which differs from that of Cluster 2. We concluded that the comodulogram-based clustering identified (at least) two separate clear brain classes. Cluster 3 and Cluster 4 showed some brain class evidence due to their weak frontal dipole density distribution, but they also suggest dipole distribution within the frontal the right temporal areas. Their comodulograms showed that Cluster 3 is similar to Cluster 2, but Cluster 4's comodulogram is not. This qualitative evaluation concluded that Cluster 3 may consist of many brain sources, while Cluster 4 may consist of artifacts. Cluster 1 showed a block diagonal comodulogram centered on the low frequency (<9.4 Hz). From what we learned above, we identified this cluster as the eye class. The dipole density distribution of the cluster supported this interpretation. Cluster 6, Cluster 7, and Cluster 8 showed block diagonal comodulograms centered on the higher frequency (above 13, 18, and 4 Hz, respectively). From our observations above, we identified these clusters as the muscle class. Dipole density distribution also showed peripheral distributions.

IC Clustering, According to Comodulogram
Further, there is a tendency that lower z coordinates are associated with temporo-occipital distribution, while higher z coordinates are related to frontal distributions. This makes sense in terms of the scalp electrode locations in which the edge channels tend to pick up muscle potentials. To sum, the comodulogram-based clustering produced reasonable classifications of the data supported by the estimated source distributions and their functions. The dipole density maps projected on anatomical images in axial view at different MNI z coordinates ranging from 0 mm (mm) to 80 mm; warmer color represents a higher density of equivalent dipoles. The same k-means clustering was used with the number of clusters specified to be 8. Cluster 2 and Cluster 5 and potentially Cluster 3 as well were identified as the brain class. Cluster 1 was identified as the eye class. Cluster 6, Cluster 7, and Cluster 8 were identified as the muscle class. Cluster 4 was not classified to any of these categories, thus left as uninterpretable noise.
In the test of suitability of using comodulogram-based clustering to classify ICs, we evaluated the consistency between our interpretation on each cluster as explained above and the majority class of ICs in that cluster as labeled by ICLabel. Table 1 shows the results. In Clusters 1, 5, and 6, more than 80% of ICs were classified into a single class, showing the high consistency level between comodulogram method and ICLabel. All ICs from these clusters are combined to represent 32.7% (258 ICs) of total ICs. The classes with maximal proportion are also consistent with our interpretation. Meanwhile, Clusters 2, 4, 7, and 8 show reasonable consistency levels between 0.5 and 0.8, showing moderate success. These clusters constitute 60.3% (476 ICs) of total ICs. The majority of ICs in Cluster 2 were classified into brain class. As we expected, clusters 7 and 8 had most ICs classified as muscle. However, we considered Cluster 4 as uninterpretable noise. This cluster had many eye ICs, demonstrating the limitation of comodulogram-based clustering. Nevertheless, this cluster with 30 ICs constituted to only 3.8% of the total ICs. The remaining ICs in the cluster were approximately The dipole density maps projected on anatomical images in axial view at different MNI z coordinates ranging from 0 mm (mm) to 80 mm; warmer color represents a higher density of equivalent dipoles. The same k-means clustering was used with the number of clusters specified to be 8. Cluster 2 and Cluster 5 and potentially Cluster 3 as well were identified as the brain class. Cluster 1 was identified as the eye class. Cluster 6, Cluster 7, and Cluster 8 were identified as the muscle class. Cluster 4 was not classified to any of these categories, thus left as uninterpretable noise.
In the test of suitability of using comodulogram-based clustering to classify ICs, we evaluated the consistency between our interpretation on each cluster as explained above and the majority class of ICs in that cluster as labeled by ICLabel. Table 1 shows the results. In Clusters 1, 5, and 6, more than 80% of ICs were classified into a single class, showing the high consistency level between comodulogram method and ICLabel. All ICs from these clusters are combined to represent 32.7% (258 ICs) of total ICs. The classes with maximal proportion are also consistent with our interpretation. Meanwhile, Clusters 2, 4, 7, and 8 show reasonable consistency levels between 0.5 and 0.8, showing moderate success. These clusters constitute 60.3% (476 ICs) of total ICs. The majority of ICs in Cluster 2 were classified into brain class. As we expected, clusters 7 and 8 had most ICs classified as muscle. However, we considered Cluster 4 as uninterpretable noise. This cluster had many eye ICs, demonstrating the limitation of comodulogram-based clustering. Nevertheless, this cluster with 30 ICs constituted to only 3.8% of the total ICs. The remaining ICs in the cluster were approximately evenly classified into brain, muscle, and eye class, making this cluster inconclusive. The consistent level was limited to 26.8%, but this cluster constitutes only 7.1% (56 ICs) of the total ICs. To sum up, even simple k-means can successfully categorize ICs with similar properties into reasonably separate subgroups. Around one-third of total ICs can be classified correctly with impressive consistency level (all > 82%) with ICLabel (an SCCN/EEGLAB plugin). The majority of remaining clusters are classified with a reasonable consistency level (all > 51.5%) than the chance level of 42.5%. We concluded that the results demonstrate the usefulness of comodulogram PPC for IC classification.

Discussion
In the present study, we demonstrated the usefulness of cross-frequency PPC analysis on ICA-decomposed EEG data. We showed that the brain, eye, and muscle classes have their features in comodulograms. The classification from comodulogram-based clustering was supported by physiological validity provided by their estimated source distributions and the reasonable consistency with the results from a standardized IC classification tool (ICLabel).
From the analysis of representative ICs, we found that the proposed method allows us to reveal a cross-frequency coupling structure unavailable for PSD by FFT. Besides, there exist many unique characteristics among brain, eye-movement, and muscle activity classes, which have been overlooked so far. PPC is a unique metric and is an added value for EEG IC classifications. We demonstrated that comodulogram-based clustering using a simple k-means algorithm yields satisfactory results. Using more sophisticated supervised/unsupervised learning with fine-tuned PPC calculation is an easy target toward better IC classification performance.
In this study, we evaluated our classification's consistency via three methods: first, by comparing our results with classification labels generated by ICLabel; second, by visual examination of PPC clustering results; and third, by examining the PPC clustering results using a standardized comodulogram and the probabilistic dipole density. Our comodulogram-based clustering was judged against the operational ground truths for the current purpose generated by ICLabel. While label probabilities calculated by ICLabel are typically very high (mostly > 80%), sometimes the highest label probability can be no more than 50%. When we examined individual ICs with low label probability, we often found that comodulogram can inform us more about the ICs with additional insight. Examples from representative ICs selected from four different subjects are illustrated in Figure 7. These plots show PSD, the class label suggested by ICLabel with label probability, and the corresponding comodulogram. In the case of successful prediction in Figure 7a, the Brain class shows convincing evidence of the first and second harmonics of alpha rhythm. In contrast, the Muscle class in Figure 7b has low label probability of 0.391, but the corresponding comodulogram shows clear parallel lines in the beta-gamma coupling, which is the common pattern we found for the case of muscle class. The ICLabel-generated label of "Brain" in Figure 7c seems challenged by high-frequency broadband coupling, which is a unique characteristic of muscle ICs. Finally, comodulogram in Figure 7d clearly shows the presence of the first alpha harmonic, despite a relatively low label probability of 0.519. These examples demonstrate that PPC comodulogram can deliver complementary insights about ICs. It might help improve the classification of ICs, especially when used together with traditional methods, thus opening new possibilities for better data-driven EEG data classification. A future systematic study may test this possibility. Note that any correlation measure is susceptible to the presence of extreme values; data with non-normal distributions require cautious use of correlational measures [60]. In the current application, muscle-class ICs showed higher powers and more sporadic activations over time compared with constant and stationary brain-class ICs. We indeed observed that the intermittent, high-power muscle activity calculated comodulogram was often unsuccessful. After investigation, we learned that extensive removal of outliers is essential for successful computation. Using robust statistics, we empirically determined trimming the top and bottom 10% from the data distribution for data cleaning resulted in satisfactory outcomes. However, we are dissatisfied this approach provides neither theoretical nor qualitative solutions to the problem. Future studies must address outliers in calculating comodulograms.
From our earlier review, we observed cross-frequency PPC had been used more in the animal electrophysiological studies using LFP. In contrast, we show that the comodulogram measure alone can classify effective EEG sources into representative physiological categories. This evidence strongly supports generalization of the PPC within EEG analyses across disparate fields. We await future studies that will clarify the generative mechanism of the cross-frequency PPC for the physiological classes (brain, eye, and muscle), sub-categories of the brain (i.e., perceptual, cognitive, behavioral functions), as well as anatomical regions. Such studies would make PPC electro-physiologically informed physiologically established measure.

Conclusions
We demonstrated that brain, eye, and muscle classes of ICA-decomposed scalp-recorded EEG signals have characteristic features in comodulograms of cross-frequency PPC analysis. The (c) data from Subject 3 led to the discordance between "Brain" label suggested by ICLabel and comodulogram pattern unique to the muscle class indicated by the current study; (d) the presence of the alpha harmonics in comodulogram from Subject 4 can boost the label probability of ICLabel-suggested class label as "Brain" despite the current low likelihood.
Note that any correlation measure is susceptible to the presence of extreme values; data with non-normal distributions require cautious use of correlational measures [60]. In the current application, muscle-class ICs showed higher powers and more sporadic activations over time compared with constant and stationary brain-class ICs. We indeed observed that the intermittent, high-power muscle activity calculated comodulogram was often unsuccessful. After investigation, we learned that extensive removal of outliers is essential for successful computation. Using robust statistics, we empirically determined trimming the top and bottom 10% from the data distribution for data cleaning resulted in satisfactory outcomes. However, we are dissatisfied this approach provides neither theoretical nor qualitative solutions to the problem. Future studies must address outliers in calculating comodulograms.
From our earlier review, we observed cross-frequency PPC had been used more in the animal electrophysiological studies using LFP. In contrast, we show that the comodulogram measure alone can classify effective EEG sources into representative physiological categories. This evidence strongly supports generalization of the PPC within EEG analyses across disparate fields. We await future studies that will clarify the generative mechanism of the cross-frequency PPC for the physiological classes (brain, eye, and muscle), sub-categories of the brain (i.e., perceptual, cognitive, behavioral functions), as well as anatomical regions. Such studies would make PPC electro-physiologically informed physiologically established measure.

Conclusions
We demonstrated that brain, eye, and muscle classes of ICA-decomposed scalp-recorded EEG signals have characteristic features in comodulograms of cross-frequency PPC analysis. The classification based on the comodulogram was supported by physiological validity provided by the estimated source distribution and the reasonable consistency with the results from a standardized IC classification tool. The cross-frequency PPC measure is a new, data-driven EEG-IC classifier capable of exploiting underrepresented univariate single-channel/component time-frequency structure information. The data processing methods herein and the interactive visualization searching environment are freely available as open-source Matlab library plugins in EEGLAB [42] to facilitate researchers' replication of our analysis.