Building an EEG-fMRI Multi-Modal Brain Graph: A Concurrent EEG-fMRI Study

The topological architecture of brain connectivity has been well-characterized by graph theory based analysis. However, previous studies have primarily built brain graphs based on a single modality of brain imaging data. Here we develop a framework to construct multi-modal brain graphs using concurrent EEG-fMRI data which are simultaneously collected during eyes open (EO) and eyes closed (EC) resting states. FMRI data are decomposed into independent components with associated time courses by group independent component analysis (ICA). EEG time series are segmented, and then spectral power time courses are computed and averaged within 5 frequency bands (delta; theta; alpha; beta; low gamma). EEG-fMRI brain graphs, with EEG electrodes and fMRI brain components serving as nodes, are built by computing correlations within and between fMRI ICA time courses and EEG spectral power time courses. Dynamic EEG-fMRI graphs are built using a sliding window method, versus static ones treating the entire time course as stationary. In global level, static graph measures and properties of dynamic graph measures are different across frequency bands and are mainly showing higher values in eyes closed than eyes open. Nodal level graph measures of a few brain components are also showing higher values during eyes closed in specific frequency bands. Overall, these findings incorporate fMRI spatial localization and EEG frequency information which could not be obtained by examining only one modality. This work provides a new approach to examine EEG-fMRI associations within a graph theoretic framework with potential application to many topics.


INTRODUCTION
Graph theory-based analysis is a powerful technique to characterize the architecture of human brain networks (Avena-Koenigsberger et al., 2014;Pessoa, 2014). Graph metrics can quantitatively describe the topological properties of brain connectivity (Klimm et al., 2014). Previous studies that applied network science and graph theory based analysis to brain imaging data have reported "economical" small-world organization of connectivity, which reflects an economic balance between network cost and network efficiency (Bullmore and Sporns, 2012;Avena-Koenigsberger et al., 2014). The brain connectome shows modular and rich club organization with sets of hub regions that are crucial for efficient neuronal signaling and communication (van den Heuvel and Sporns, 2013b;Senden et al., 2014). Further studies show that the graph metrics and network structures of brain connectivity are altered in brain disorders (Collin et al., 2014;Crossley et al., 2014;Deco and Kringelbach, 2014;Hong et al., 2014;Korgaonkar et al., 2014;van den Heuvel and Fornito, 2014;Fornito and Bullmore, 2015a,b;Fornito et al., 2015;Gong and He, 2015;Wheeler et al., 2015). However, most of these studies have built the brain graphs with a single modality of brain imaging data. Additional insights of brain connectivity may thus be obtained by combining information from multiple modalities (Sui et al., 2012;Reid et al., 2016).
Different imaging techniques are sensitive to different aspects of brain dynamics. For example, functional magnetic resonance imaging (fMRI) measures the highly localized hemodynamic response throughout the brain, with a good spatial resolution (about 2-3 mm) but relatively poor temporal resolution. Electroencephalography (EEG) measures cortical electrical activity with a much higher temporal resolution, but its poor spatial resolution precludes precise anatomical identification of underlying neural sources. FMRI and EEG therefore represent complementary imaging signals, and combining concurrently collected data is a particularly useful way to examine brain dynamics over a broad range of spatial and temporal scales Herrmann and Debener, 2008;Eichele et al., 2009;Rosenkranz and Lemieux, 2010;Wu et al., 2010;Lei et al., 2011;Laufs, 2012;Bridwell et al., 2013;Mulert, 2013).
For coupling concurrent EEG-fMRI data, a popular approach is to analyze correlations between fMRI voxel time-series and EEG spectral power fluctuations (Valdes-Sosa et al., 2009;Rosa et al., 2010;Jorge et al., 2014). For example, it has been observed that low frequency EEG connectivity appears to best resemble fMRI connectivity (Deligianni et al., 2014). Brain connectivity between different regions detected by fMRI is associated with activity within different frequency bands of the EEG signal (Tagliazucchi et al., 2012;Chang et al., 2013;Liu et al., 2014). These findings provide electrophysiological signatures of functional brain connectivity identified in fMRI data (Mantini et al., 2007;Meir-Hasson et al., 2014). In this study, we make a further step to investigate topological organization of multimodal EEG-fMRI brain connectivity. Dynamic connectivity over time is an important feature in functional brain networks (Hutchison et al., 2013a,b;Li et al., 2013Li et al., , 2014Calhoun et al., 2014;Stephen et al., 2014;Yang et al., 2014), and dynamic connectivity patterns have been widely studied with fMRI Rashid et al., 2014;Yu et al., 2015). Graph theory based analysis has been successfully implemented to assess dynamics during cognitive tasks or during rest (Doron et al., 2012;Bassett et al., 2013;Cocchi et al., 2013;Betzel et al., 2014;Cole et al., 2014;Dwyer et al., 2014;Hermundstad et al., 2014;Zalesky et al., 2014;Davison et al., 2015;Liang et al., 2016). A few recent studies even examined the relationships between dynamic fMRI connectivity and EEG signals (Tagliazucchi et al., 2012;Chang et al., 2013). However, the graph properties of dynamic brain networks with multi-modal nodes are largely unknown.
The aim of this study is to explore graph properties of concurrent EEG-fMRI multi-modal brain connectivity. Static and dynamic EEG-fMRI brain graphs are built using concurrently collected data from 25 healthy subjects during eyes open and eyes closed. Graph nodes are represented by EEG electrodes and fMRI components identified using group independent component analysis (ICA; Calhoun et al., 2008). Graph edges are represented by correlations between fMRI time courses and/or the EEG spectral power time courses of five frequency bands (delta, theta, alpha, beta, low gamma; for details, see the Methods Section below). The findings characterize changes within graphical properties of connectivity across different states (eyes open vs. eyes closed), while incorporating the high spatial resolution of fMRI (by estimating nodal graph measures of specific brain regions) and the high temporal resolution of EEG (i.e., we integrate fluctuations of fMRI brain regions with EEG frequencies within a graphtheoretic framework). In addition to the different spatial and temporal resolutions, different aspects of neural activity are also integrated between the two modalities, since EEG is sensitive to synchronous cortical synaptic potentials, and BOLD fMRI is sensitive to BOLD oxygenation changes that follow increased post-synaptic metabolism (see Bayram et al., 2011;Bridwell and Calhoun, 2014). This work develops a novel framework to build multi-modal brain graphs, demonstrating associations between EEG and fMRI within a graphical theoretical framework.

Participants
Twenty-five healthy subjects (age: 29 ± 8; 8 females) were recruited via advertisements at the University of New Mexico and by word-of-mouth. Each individual had normal or corrected to normal vision and hearing. Prior to inclusion in the study, participants were screened to ensure they were free from DSM-IV Axis I or Axis II psychopathology [assessed using the SCID (First et al., 1995)] and to ensure that there was no history of neurological disease. All participants provided informed written consent at the Mind Research Network, and were compensated for their participation. The experiment design and simultaneous acquisition details were described in our previous study (Wu et al., 2010).

Experimental Design
Simultaneous EEG-fMRI data were recorded while individuals rested first with their eyes closed (8.5 min), and then with their eyes open (8.5 min). Individuals were instructed to relax, lie still, and remain awake for the duration of each recording.

EEG Acquisition
EEG was recorded with a 32-channel BrainAmp MR-compatible system (Brainproducts, Munich, Germany) and a BrainCap electrode cap (Falk Minow Services, Herrsching-Breitbrunn, Germany). The Ag/AgCI electrodes were placed according to the international 10-20 system. Electrocardiogram (ECG) and eye movement (EOG) signals were recorded in separate channels, reducing the number of scalp electrodes to 30. The reference channel was placed at FCz. The impedance of each electrode was kept lower than 5 K using conductive and abrasive electrode paste. The EEG signals were sampled at 5 KHz. To avoid temporal jitter, the EEG amplifier and fMRI were synchronized using an in-house device.

fMRI Acquisition
Functional MRI brain images were acquired with a Siemens Sonata (Siemens, Malvern PA) scanner at 1.5 T by means of a T2 * -weighted echo planar imaging sequence with the following parameters: repeat time (TR) = 2 s, echo time (TE) = 39 ms, field of view = 224 mm, acquisition matrix = 64 × 64, flip angle = 80 • , voxel size = 3.5 × 3.5 × 3 mm, gap = 1 mm, 27 slices, ascending acquisition. FMRI scans consisted of 256 volumes for each condition (eyes open and eyes close).

EEG Processing
EEG data were preprocessed in Matlab (http://www.mathworks.com) using custom and built-in functions and the EEGLAB toolbox (http://sccn.ucsd.edu/ eeglab). The EPI gradient artifact was attenuated by calculating the average artifact template (across 2 s epochs) separately for each channel, and subtracting the template from the individual epochs within that channel. The EEG data were down-sampled to 1 kHz, band-pass filtered (0.01 to 50 Hz), and average referenced. Additional artifacts (e.g., BCG, eye movement, and residual EPI gradients) were attenuated by conducting a temporal ICA decomposition on the individual recordings (Srivastava et al., 2005). Thirty components were estimated using the extended Infomax algorithm implemented in EEGLAB (Bell and Sejnowski, 1995;Lee et al., 1999). Artifactual components were identified by visual inspection of the component time-course, topographic distribution, and frequency spectrum and removed from the back reconstructed time-course. Seventeen components were eliminated on average (min 11, max 23). In this work, ballistocardiac artifacts were corrected only by ICA. Previous studies showed that using ICA to attenuate the BCG artifacts in EEG data collected in a low magnetic field of 1.5 T is acceptable rather than in higher field scanners like 3 T or 7 T . We visually inspected the EEG waves and found the BCG artifacts were indeed largely removed. The experiment design, simultaneous acquisition details and data preprocessing were described in our previous studies (Wu et al., 2010;Bridwell et al., 2013). The data used in this work is the same as in these two studies.
The preprocessed EEG data were variance normalized, segmented into 2 s epochs (resulting in 256 epochs, i.e., an epoch that corresponds to each concurrently recorded fMRI volume) and converted to the frequency domain by the fast Fourier transform (FFT). The spectral power was averaged within 5 frequency bands (delta:  . The gamma band was restricted to lower frequencies (30)(31)(32)(33)(34)(35) in order to avoid the pump and ventilation artifacts, which dominate above 40 Hz (Bridwell et al., 2013;Nierhaus et al., 2013). We noted that some other EEG literature define gamma as >40 Hz and call the band 30-35 Hz as high beta. However, in this study we name the band 30-35 Hz as low gamma.

Building EEG-fMRI Brain Graphs
A correlation matrix R was constructed with elements (r ij ) representing Pearson correlation coefficients computed using the 30 EEG electrodes' spectral time-courses and the 54 fMRI ICs' time-courses. This process was repeated for the five EEG frequency bands and the two conditions (EO and EC). When computing the correlation between EEG and fMRI signals, following previous studies (Goldman et al., 2002;Laufs et al., 2003a;Moosmann et al., 2003), EEG power time courses were convolved with a canonical hemodynamic response function (HRF) to account for the delayed hemodynamic response.
Consequently, undirected static connectivity EEG-fMRI graphs were built from each of the N × N (N = 84 in this study, including 30 EEG electrodes and 54 fMRI brain components) correlation matrices R. In order to preserve the information for both positive and negative correlations, weighted positive (W + ) and negative (W − ) connection graphs were built based on R.
In positive connection graphs, negative correlations in R were replaced by 0 and positive correlation values were maintained. In negative connection graphs, positive correlations were replaced by 0 and absolute values of negative correlations in R were maintained.
Dynamic EEG-fMRI graph analysis was performed by calculating correlation matrices along successive sliding windows of the matrix EF (256 × 84; width, L = 20 TRs, in steps of 1 TR; Allen et al., 2014;Yu et al., 2015). As with the static analysis, the first 30 columns correspond to EEG electrodes and the following 54 columns correspond to fMRI ICs. Two hundred thirty-seven EEG-fMRI correlation matrices were computed for 237 (237 = 256 − 20 + 1) windows. Positive and negative connection graphs were analyzed separately, as in the static analysis. See Figure 1 for the framework of building static and dynamic concurrent EEG-fMRI multi-modal brain graphs.
Here we use a window width of 20 TRs (40 s) based on a previous study indicating that cognitive states may be correctly identified with as little as 30-60 s of data (Shirer et al., 2012). A recent study has shown that non-stationary fluctuations in functional connectivity can in theory be detected with window length of 40 s . Also, others have demonstrated that changes of brain connectivity are not particularly sensitive to the specific time-window length (in the range of 10-20 2 s TRs, 20-40 s; Li et al., 2014). Previous works of our group and others consistently show that shorter time windows result in a lower number of statistically significant correlations in brain connectivity and greater variability of correlation values (Chang and Glover, 2010;Hutchison et al., 2013b;Allen et al., 2014;Yu et al., 2015), whereas a sliding window size of about 22 TRs (44 s) provides a good tradeoff between the ability to resolve dynamics and the quality of connectivity estimation Yu et al., 2015).
Connectivity strength (CS), clustering coefficient (CC), and global efficiency (GE) are three basic and important graph metrics which measure the functional segregation and integration of brain networks (Rubinov and Sporns, 2010). To quantitatively assess the topological properties of the brain connectivity, both global level and nodal level of these three graph metrics were derived for static and dynamic graphs in this study using the brain connectivity toolbox (http://www.brain-connectivity-toolbox.net/). For mathematical definitions and equations for computing the graph measures see the below Section Equations for Computing Graph Measures. More details are available in (Rubinov and Sporns, 2010)]. The variances of 237 dynamic graph metrics and the amplitude of low frequency (0-0.025 Hz) fluctuations of the time series of dynamic graph measures were computed in each subject. For statistical analysis, 5 (frequency band) × 2 (eyes condition) compound symmetry repeated measure analysis of variance (ANOVA) and paired t-tests were performed on static and dynamic measures.

Equations for Computing Graph Measures
We denote G as the set of all nodes in the weighted graph W, and N (N = 84 in this study) is the number of nodes. Connectivity strength of node i is defined as below: The connectivity strength of whole graph (global level) is the average of the connectivity strength of all the nodes in the graph: Nodal level Clustering coefficient is computed using the below equation: The clustering coefficient of whole graph is the average of the clustering coefficient of all the nodes in the graph: Global efficiency of node i is defined as: In which Where f is a map (here is an inverse) from weight to length and g i w ↔j is the shortest weighted path between i and j. Global efficiency of the whole graph is the average of the global efficiency of all the nodes in the graph.

Detecting Connectivity States
Recent fMRI studies showed that fluctuations of time-varying functional brain connectivity gives rise to discrete highlyorganized patterns that may emerge or dissolve over time, which are called connectivity states (Cribben et al., 2012; FIGURE 1 | Pipeline for building concurrent EEG-fMRI multi-modal brain graphs. x Segment EEG signal into 2 s time windows, and compute the average spectral power within a selected frequency window. y Perform group ICA on fMRI data. z Compute the correlation coefficient within and across the EEG spectral power's and fMRI ICA's full time courses, generating one EEG-fMRI static connectivity matrix for each frequency band. { Segment EEG spectral power and fMRI ICA time courses into time windows, then compute the correlation between each pair of time windowed time courses to get dynamic EEG-fMRI brain graphs. Both positive and negative correlations in the correlation matrix (R) are shown in this figure. These steps are repeated for each of the 5 frequency bands and during eyes open and eyes closed conditions. Allen et al., 2014;Yang et al., 2014;Yu et al., 2015). Here we performed the method developed in one of our previous studies (Yu et al., 2015) to detect connectivity states of the dynamic EEG-fMRI graphs in each individual. Firstly, nodal level connectivity strength of each time-varying EEG-fMRI graph was computed. To estimate how the EEG-fMRI network patterns of different time-windows were associated to each other, a new correlation matrix CCS (237 × 237; 237 is the number of time windows) was then computed based on correlations of the nodal connectivity strength between each pair of time windows across 84 nodes. Modular community structure is one of the most ubiquitous properties of complex networks (Newman, 2006;Bullmore and Sporns, 2009). Modularity is a function that measures the quality of a division of nodes into groups or communities, and modules of the matrix CCS may correspond to sets of time windows with similar brain connectivity patterns. Thus, the modular organization of CCS was analyzed with the modularity algorithm of Newman (2006) implemented in the brain connectivity toolbox. The number of modules of CCS is the number of connectivity states for the dynamic EEG-fMRI graph. Finally, the EEG-fMRI brain graphs from different time windows that categorized to the same module were averaged to get the graph of that connectivity state. More details of this approach were introduced in (Yu et al., 2015).  Allen et al., 2014;Yu et al., 2015), and a subset have been associated with cognitive functions in meta-analytic studies (Rottschy et al., 2012;Balsters et al., 2014;Kohn et al., 2014;Amft et al., 2015). Figure 2B displays the structure of stationary connectivity between graph nodes (ICNs and EEG channels), computed over the entire fMRI time courses and EEG spectra power time courses for the five frequency bands (delta, theta, alpha, beta, low gamma) and averaged over 25 subjects in each condition (eyes open and eyes closed). Patterns of connectivity within fMRI ICNs are consistent with prior literature, showing modular organization within sensory systems and default mode regions, as well as anticorrelation between these regions (Fox et al., 2005;Shirer et al., 2012;Allen et al., 2014;Yu et al., 2015).

Static EEG-fMRI Graph
For the global level graph metrics of positive connection networks, a five (frequency band: delta, theta, alpha, beta, low gamma) × two (eyes condition: open, closed) compound symmetry repeated measures ANOVA shows that the main effect of frequency band is significant (P < 0.001) for all three metrics, and the main effect of eyes condition is significant (P < 0.01) on clustering coefficient only. Post-hoc paired t-tests reveal significantly (P < 0.01) higher clustering coefficient during eyes closed than eyes open in beta band.
For the global level graph metrics of negative connection networks, a five (frequency band: delta, theta, alpha, beta, low gamma) × two (eyes condition: open, closed) compound symmetry repeated measures ANOVA shows that the main effects of frequency band and eyes condition are significant (P < 0.05) for all three metrics (connectivity strength, clustering coefficient, global efficiency). Post-hoc paired t-tests reveal significantly (P < 0.01) higher connectivity strength during eyes closed than eyes open in alpha and beta bands, and significantly higher global efficiency during eyes closed in the beta band. Figures S2, S3 show the group means of the graph metrics in the eyes closed and eyes open conditions computed within the delta, theta, alpha, beta and low gamma frequency bands of positive and negative connection graphs, respectively.
For nodal level graph metrics of positive connection graphs, the main effect of eyes condition is significant (FDR correction, q < 0.001) on all three graph metrics for 3 brain components which belong to somatomotor, visual, and auditory components, respectively. Graph measures are higher during eyes closed than during eyes open. For the spatial maps of the 3 ICNs see Figure S4.
For nodal level graph metrics of negative connection graphs, the main effect of eyes condition is significant (FDR correction, q < 0.001) on all three graph metrics for only one visual brain components. Graph measures are higher during eyes closed than during eyes open. For the spatial maps of that ICN see Figure S5.
To visually display the difference between eyes conditions of the brain network, as an example, we show the values of the nodal graph measures and the pattern of the connections from a visual component node to all of the other graph nodes in the five frequency bands during eyes open and eyes closed conditions in positive and negative connection graphs in Figures 3, 4, respectively.  Figures S17, S18, B1, B2) shows that low-frequency CS oscillations peak between 0.001 and 0.02 Hz.

Dynamic EEG-fMRI Graph
Variance (VAR) and amplitude of low frequency (LFA) [0-0.025 Hz] oscillations of the time varying global level graph metrics are computed. For time-varying positive connection graphs, five (frequency band: delta, theta, alpha, beta, low gamma) × two (eyes condition: open, close) compound symmetry repeated measure ANOVA shows that the main effect of eyes condition is not significant on VAR and LFA of all dynamic graph measures. The main effect of frequency band is significant (P < 0.01) on VAR and LFA of CS and GE (see Figure 5), indicating that the graph metrics computed using delta and theta EEG frequencies demonstrate low frequency patterns of time varying connectivity. For time-varying negative connection graphs, the main effect of eyes condition is significant (P < 0.01) on VAR of CS and GE, and on LFA of CS (P < 0.05). The main effect of frequency band is significant (P < 0.001) on VAR and LFA of CS and GE (see Figure 6). In general, we found greater variance (VAR) and LFA in CS and GE in eyes closed alpha compared to eyes open alpha.
For nodal level dynamic graph metrics of positive connection networks, the main effect of eyes condition on the VAR and LFA of all three dynamic measures is significant (FDR correction, q < 0.001) at two cognitive control brain components. For the spatial map of these two components see Figure S19. One visual component shows significant (FDR q < 0.001) main effect of eyes condition on VAR and LFA of two metrics (CS and GE). See Figure S20 for its spatial map. The same component is shown in (Figures 3, 4, 9, 10). To demonstrate the dynamic properties of nodal level graph metrics, we show VAR and LFA of the visual component as an example in Figures 7, 8.

Connectivity States
Consistent with previous dynamic fMRI connectivity studies Damaraju et al., 2014;Rashid et al., 2014;Yang et al., 2014), some connectivity patterns of the dynamic EEG-fMRI graphs reoccur over time. Connectivity states are detected using the method developed in our previous study (Yu et al., 2015). Specifically, 2-6 states are detected during each eyes condition in each subject for both positive and negative. See Tables S1, S2 for the details about how many connectivity states are detected in each subject. For a visual view of the structure of different connectivity states in an example subject, see Figures  S21-S24, and Figures 9, 10.

DISCUSSION
In the present study, concurrent EEG-fMRI resting state data collected during eyes open and eyes closed conditions are used to build multi-modal brain graphs. FMRI data are decomposed with group ICA into ICNs and corresponding time courses. EEG signals are segmented into 2 s-epochs and the spectral power is computed and averaged within five frequency bands (delta, theta, alpha, beta, and low gamma) for each segment. EEG-fMRI brain graphs are built by computing the correlations between and among fMRI ICA time courses and EEG spectral power time courses. Connectivity strength, local efficiency, and global efficiency are calculated for both static graphs, which are estimated using the full length of time courses, and dynamic graphs, which are estimated using a sliding window method. Five (frequency band: delta, theta, alpha, beta, low gamma) × two (eyes condition: open, close) compound symmetry repeated measure ANOVA and paired t-tests are performed to identify significant differences across frequencies and eyes conditions. For static graphs, in positive connection networks, graph metrics are higher during the eyes closed condition compared to eyes open mainly in delta and theta bands (Figure 3). In negative connection networks, graph metrics are higher during eyes closed compared to eyes open primarily in alpha and beta bands (Figure 4, Figure S3). For time varying graphs, in positive connection networks, the LFA and VAR of dynamic graph measures (in nodal level of specific nodes) are higher in eyes closed than eyes open mainly in delta and theta bands (Figure 7). In both positive and negative connection networks, generally, we found greater variance (VAR) and LFA in CS and GE in eyes closed alpha compared to eyes open alpha (Figures 5, 6) which is in line with previous studies (Wu et al., 2010;Bridwell et al., 2013). Consistent with previous single modality fMRI studies, dynamic EEG-fMRI connectivity shows some connectivity states which re-occur over time. This work provides an important first step in fusing EEG and fMRI using a graph theoretical framework.
In early studies which combine EEG and fMRI, EEG signals are traditionally separated into five frequency bands: delta, theta, alpha, beta, and gamma (Laufs et al., 2003b;Mantini et al., 2007;Keilholz, 2014). Different frequencies have been linked to different functional properties (Buzsaki, 2006). For example, alpha power increases at rest with eyes closed compared to eyes open (Pfurtscheller et al., 1996;de Munck et al., 2007), and increases when a greater number of items are held in working memory (Klimesch et al., 1997(Klimesch et al., , 1999Jensen et al., 2002). Experimental results have related power increases and synchronization in the gamma frequency band to the performance of perceptual and cognitive operations, including attention (Womelsdorf and Fries, 2007), conscious perception (Melloni et al., 2007), and decision making (Donner et al., 2009). Slower frequencies, such as delta, arise during sleep, and are hypothesized to reflect diminished temporal complexity underlying loss of conscious awareness (Tononi et al., 1994;Tononi and Edelman, 1998;Tononi, 2004Tononi, , 2008. Low frequency oscillations also appear to correspond to the cognitive events which primarily contribute to evoked potential's (Demiralp et al., 1999).
The majority of previous studies which combine EEG and fMRI data use correlation or general linear modeling (GLM) to link fluctuations between multiple EEG frequency bands and fMRI voxels (Bridwell and Calhoun, 2014). Correlations between EEG power variations of delta, theta, alpha, beta, gamma rhythms and BOLD activity of specific brain regions (ICNs; Laufs et al., 2003b), or BOLD connectivity between brain regions (ICNs) have been estimated (Tagliazucchi et al., 2012). These studies suggest that each functional brain ICN is characterized by a specific electrophysiological signature, and that BOLD fMRI fluctuations have a neurophysiological origin (Mantini et al., 2007).
In this work, we separate the EEG data into five frequency bands as in previous studies. However, in addition to computing the correlations between EEG signals and fMRI BOLD signals, we compute EEG-fMRI multi-modal brain graphs in which EEG nodes provide high temporal resolution information and fMRI nodes provide high spatial resolution. The finding that graph metrics show differences across frequency bands (the main effect of frequency band is significant) is consistent with the hypothesis that different EEG frequencies are associated with different BOLD activities. Within this study, we characterize different graph properties between eyes open and eyes closed conditions. In static positive connection EEG-fMRI graphs, nodal level graph metrics are higher during the eyes closed condition in three brain components which belong to somatomotor, visual, and auditory areas. In negative connection graphs, a visual component shows different nodal level graph measures (for all three metrics) between eyes conditions. The LFA and VAR of dynamic nodal graph measures of two cognitive control components are higher during eyes closed than eyes open for positive connection networks. In general, these findings are consistent with and add to previous studies demonstrating differences in BOLD amplitudes and functional connectivity across the two conditions (McAvoy et al., 2008;Zou et al., 2009Zou et al., , 2015. Importantly, the present findings provide new insights by incorporating fMRI spatial locations and EEG frequency bands within graph theoretic measures of brain connectivity between eyes closed and eyes open conditions. For example, the differences of graph metrics between eyes conditions are mainly in delta and theta bands for positive connection networks and mainly in alpha and beta bands in negative connection networks. Multiple recent brain imaging studies suggest that the functional brain connectivity is not stationary but changes over minute-to-minute intervals (Hutchison et al., 2013a;  et al., 2014). Here we assess dynamic properties (such as LFA and VAR) of the time-varying EEG-fMRI brain graphs and their associated connectivity states. Our results characterize dynamic measures of multi-modal functional brain organization by combining concurrent EEG-fMRI signals, and support the hypothesis that variability of brain connectivity emerges from structured connectivity patterns that emerge and dissolve over time .
Notably, the findings that alteration of graph measures of specific fMRI nodes across eyes conditions occurs in particular EEG frequency bands provide new electrophysiological signatures of functional brain connectivity examined in fMRI data, and imply that the graph-theory based analysis is powerful to assess the associations between EEG and fMRI. However, a few potential methodological limitations need to be discussed. Graph metrics may depend in part on the methods used to identify nodes. Thus, it is worth considering the difference between ICA-based and anatomically-based approaches. Within fMRI, brain graph nodes are often formed using predefined anatomical templates such as automated anatomical labeling (AAL; Tzourio-Mazoyer et al., 2002;Liu et al., 2008;Lynall et al., 2010), randomly generated templates (Hagmann et al., 2008; Fornito et al., 2010), and voxel-based divisions (Eguíluz et al., 2005;Buckner et al., 2009;Yu et al., 2011cYu et al., , 2013. Different approaches may significantly modulate the quantitative measures of graph metrics of brain connectivity de Reus and van den Heuvel, 2013). Also, prior work has shown a detriment to network estimation when using atlas-based regions of interest (ROIs) as graph nodes (Smith et al., 2011;Craddock et al., 2012;Shirer et al., 2012). Moreover, the ROIs provide an imperfect segregation of the functional boundaries of the human brain. However, ICA, which is adopted in this study, provides a data-driven approach to identify spatial brain components as functionally homogeneous nodes (Yu et al., 2011a,b;. We choose a relatively high model order (i.e., 100 ICs) ICA, because previous studies have demonstrated that such models yield refined components which correspond to known anatomical and functional segmentations (Kiviniemi et al., 2009;Abou-Elseoud et al., 2010). Importantly, previous work has shown that fMRI graph measures are relatively insensitive to high model orders (Yu et al., 2011b(Yu et al., , 2015. However, a limitation of this study is related to the differences in node number and edge weight from different modalities. Positive correlations within EEG signals are much higher than within fMRI signals FIGURE 9 | Patterns of connections from a visual component node to all other graph nodes for the connectivity states identified in an example subject in the alpha band for time-varying positive connection graphs. Color dots inside the brain map indicate the approximate centroid location of fMRI brain components. Color dots outside the brain map indicate EEG electrodes. (FE, frontal electrodes; CE, central electrodes; PE, parietal electrodes; OE, occipital electrodes; TE, temporal electrodes).
FIGURE 10 | Patterns of connections from a visual component node to all other graph nodes for the connectivity states identified in an example subject in the alpha band for time-varying negative connection graphs. Color dots inside the brain map indicate the approximate centroid location of fMRI brain components. Color dots outside the brain map indicate EEG electrodes. (FE, frontal electrodes; CE, central electrodes; PE, parietal electrodes; OE, occipital electrodes; TE, temporal electrodes). and between EEG-fMRI signals (see Figure 2B). Also, the Pearson correlation may contain some redundant information, though the ICA performed in data processing, which fits all the interacting network information in a single model with multiple components, may control some of them.
It's important to note that graph measures were computed based on the formula defined for a single-modal (classical) graph. The graph metrics may be interpreted in the same way as traditional graphs. But it is unclear how global level graph measures within a multi-modal graph would be affected by the distribution of edges and nodes from different modalities. This limitation is shared by each of the two conditions examined within the present study (eyes open and eyes closed). Thus, the observation of graph metric differences here motivates further studies which extend the single-modal (classical) graph formula for EEG-fMRI (Zhang et al., 2008). In addition, future work is needed to develop criteria for determining the number of nodes and edges in the context of a multi-modal brain graph. Also, future work may build EEG and fMRI graphs separately and evaluate the correlation between graph metrics between the two, or develop new methods for defining brain regions (graph nodes) with both EEG and fMRI information available as previous studies which estimated brain graph using multimodality imaging data Hermundstad et al., 2013;Liang et al., 2013;van den Heuvel and Sporns, 2013a;Tewarie et al., 2014a,b).

CONCLUSIONS
We believe that this work provides an important beginning step in characterizing EEG-fMRI associations within a graph theoretical framework. Both static and dynamic EEG-fMRI graphs are built in five EEG frequency bands on concurrently collected EEG-fMRI data while individuals rested with eyes open and eyes closed. Differences in global and nodal level static graph metrics including connectivity strength, local efficiency, and global efficiency, are revealed among frequency bands and between eyes conditions. Dynamic properties of the graph metrics also show differences between eyes conditions. These findings incorporate spatial location (provided by fMRI) information and frequency (delta, theta, alpha, beta, and gamma bands provided by EEG) information in identifying graph properties that differ between brain states (i.e., eyes open vs. eyes closed) by linking electro-hemodynamic responses. This paper proposes a novel approach for assessing associations among concurrent EEG and fMRI measures which couples electoral and hemodynamic BOLD signals in the brain at a network level.

AUTHOR CONTRIBUTIONS
QY designed the study; analyzed and interpreted the data; drafted and revised the manuscript; gave final approval. LW designed the study; collected, analyzed, and interpreted the data; revised the manuscript and gave final approval. DB analyzed and interpreted the data; revised the manuscript and gave final approval. EE analyzed and interpreted the data; revised the manuscript and gave final approval. YD analyzed and interpreted the data; revised the manuscript and gave final approval. HH analyzed and interpreted the data; revised the manuscript and gave final approval. JC collected, analyzed, and interpreted the data; revised the manuscript. PL analyzed and interpreted the data; revised the manuscript and gave final approval. JS collected, analyzed, and interpreted the data; revised the manuscript. GP designed the study; interpreted the data; revised the manuscript and gave final approval. VC designed the study; interpreted the data; revised the manuscript and gave final approval.

ACKNOWLEDGMENTS
This work is supported by the National Institutes of Health (NIH) grants including a COBRE grant (P20GM103472), R01 grants (R01EB005846, 1R01EB006841, 1R01DA040487, REB020407, EB000840), and other grants (5P20RR021938, R37 MH43775 PI: Pearlson). This work is partly supported by the "100 Talents Plan" of Chinese Academy of Sciences, the state high-tech development plan of China (863)

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnhum. 2016.00476