EEG functional connectivity is partially predicted by underlying white matter connectivity

Over the past decade, networks have become a leading model to illustrate both the anatomical relationships (structural networks) and the coupling of dynamic physiology (functional networks) linking separate brain regions. The relationship between these two levels of description remains incompletely understood and an area of intense research interest. In particular, it is unclear how cortical currents relate to underlying brain structural architecture. In addition, although theory suggests that brain communication is highly frequency dependent, how structural connections influence overlying functional connectivity in different frequency bands has not been previously explored. Here we relate functional networks inferred from statistical associations between source imaging of EEG activity and underlying cortico-cortical structural brain connectivity determined by probabilistic white matter tractography. We evaluate spontaneous fluctuating cortical brain activity over a long time scale (minutes) and relate inferred functional networks to underlying structural connectivity for broadband signals, as well as in seven distinct frequency bands. We find that cortical networks derived from source EEG estimates partially reflect both direct and indirect underlying white matter connectivity in all frequency bands evaluated. In addition, we find that when structural support is absent, functional connectivity is significantly reduced for high frequency bands compared to low frequency bands. The association between cortical currents and underlying white matter connectivity highlights the obligatory interdependence of functional and structural networks in the human brain. The increased dependence on structural support for the coupling of higher frequency brain rhythms provides new evidence for how underlying anatomy directly shapes emergent brain dynamics at fast time scales.


Introduction
Network science provides an intuitive framework to study brain organization. Over the past decade, networks have become a leading model to illustrate both the anatomical relationships (structural networks) and the coupling of dynamic physiology (functional networks) linking separate brain regions. Alterations in functional and structural brain networks have been reported in normal cognitive processes (Uhlhaas et al., 2009;Hipp et al., 2011), across development (Hermoye et al., 2006;Micheloyannis et al., 2009;Power et al., 2010;Boersma et al., 2011;Smit et al., 2012;Chu et al., 2014), and in a wide range of neurological diseases (de Haan et al., 2009;Uhlhaas and Singer, 2010;Kramer et al., 2010). Although brain functional networks are embedded in anatomical space, a complete understanding of the associations between functional and structural networks remains an active research area (Rubinov et al., 2009;Honey et al., 2009;Honey et al., 2010). A better understanding of how anatomical scaffolds support complex, temporally organized brain activity is necessary to study normal cognitive processes, as well as to improve identification and prediction of alterations in brain function in disease states.
Brain structural networks are believed to conserve energy by minimizing longer axonal projections while simultaneously maximizing integration between local brain regions (Bullmore and Sporns, 2009;Chklovskii et al., 2002;Klyachko et al., 2003). Tissue based dissection and tracing methods in non-human primates have found that functionally related cortical areas are connected by shorter anatomical path length (Fellemen and van Essen, 1991;Hilgetag et al., 2000, Averbeck andSeo, 2008). New analysis techniques applied to modern neuroimaging data have enabled in vivo evaluation of white matter tracts to map the anatomical wiring between brain regions as well as the correlations between hemodynamic blood-oxygen level dependent (BOLD) signal fluctuations in different brain regions. Work in humans has consistently reported a strong association between structural and functional brain networks obtained using co-registered white matter tractography with low frequency BOLD oscillations Damoiseaux and Greicius 2009;Rubinov et al., 2009;Baria et al., 2011;Goni et al., 2014). Computational models simulating dynamic activity suggest that structural networks constrain functional networks across many times scales (Ponten et al., 2010;Honey et al., 2007;Hlinka and Coombes, 2012). However, in human studies and simulations, brain anatomy does not fully predict the spontaneous functional associations observed.
Prior work relating functional and structural brain networks has focused on functional networks derived from hemodynamic oscillations in the brain, a surrogate marker for neuronal activity. How and whether the associations between direct measures of cortical currents relate to underlying brain structural architecture remains unclear. Recent work has shown that fMRI BOLD signal fluctuations may co-localize with slow fluctuations in EEG gamma power (He and Raichle, 2009;Ko et al., 2011), suggesting that functional networks inferred from cortical currents should correspond to underlying structural networks. But, the relationship between the hemodynamic signals measured by fMRI and brain electrophysiology remains controversial (Bartels et al., 2008;Britz et al., 2010;Hlinka et al., 2010;Christen et al., 2014). In addition, although theory and observation suggest that higher frequencies tend to integrate more focal regions and lower frequencies broader cortical regions (Singer, 1999;Kopell, et al., 2000;Miller et al., 2007;He and Raichle, 2009;Tallon-Baudry, 2009; Baria et al., 2011), whether structural connections influence overlying functional connectivity in a similar manner across different frequency bands remains unknown.
To directly address these questions, we evaluated functional networks inferred from statistical associations between localized cortical currents using electrical source imaging (ESI) techniques, and corresponding cortico-cortical structural brain connectivity determined by probabilistic white matter tractography. Sophisticated techniques in ESI provide accurate localization of cortical source activity from high-density scalp EEG signals with fine temporal and spatial resolution (Hämäläinen and Sarvas, 1989;Michel et al., 2014). We evaluated the spontaneous fluctuating cortical brain activity over a long time scale (minutes) and related the inferred functional networks to underlying structural connectivity for broadband signals, as well as in seven distinct frequency bands. We found that cortical networks derived from ESI partially reflect both direct and indirect underlying white matter connectivity in all frequency bands evaluated. In addition, we found that although long-range functional connectivity is reduced in high frequency bands compared to low frequency bands, this reduction is significantly less when structural support is present. These results provide direct evidence for a link between structural anatomy and cortical functional connectivity in the human brain.

Material and Methods
Patients with high density EEG (70 electrodes), digitized electrode coordinates, and high resolution diffusion tensor imaging (DTI; 60 diffusion-encoding directions, 1.85 mm isotropic voxels) were retrospectively identified from clinical evaluations performed at the Massachusetts General Hospital Athinoula A. Martinos Center for Biomedical Imaging between 1/2009 and 12/2012. Patients identified for inclusion were 10-19 years of age, all female, and undergoing evaluation due to epilepsy from a variety of etiologies. Clinical information for these subjects is listed in Table 1. All EEGs were recorded in the interictal state. Analysis of the data from these patients was performed under protocols approved and monitored by the local Institutional Review Board according to National Institutes of Health guidelines.
EEGs were recorded with a 70-channel electrode cap, based on the 10-10 electrodeplacement system (Easycap, Vectorview, Elekta-Neuromag, Helsinki, Finland) in the quiet resting or sleep state. The positions of the EEG sensors were determined prior to data acquisition with a 3D digitizer (Fastrak, Polhemus Inc., Colchester, VA). The sampling rate was 600 Hz and the data were filtered with high-and low-pass filters (third-order Butterworth, zero-phase shift digital filtering) from 1-50 Hz for analysis using the MATLAB Signal Processing Toolbox and custom software. Data were visually inspected and large movement, muscle artifacts and electrodes with poor recording quality removed. A minimum of two minutes (range 122-199 seconds) of artifact-free recording from each patient was used for analysis. This epoch size has been previously demonstrated to be sufficient to produce high similarity between inferred functional networks within patients across different states of consciousness (Chu et al., 2012).

Electrical Source Imaging
Source analysis of EEG data was performed using the MNE software package (Hämäläinen and Sarvas, 1989, Sharon et al., 2007, Gramfort et al., 2013 with anatomical surfaces reconstructed using Freesurfer (Fischl, 2012). MNE provides a distributed source estimate of cortical currents incorporating constraints from the patients' MRI, transforming the data to brain space without requiring heuristic choices or strong assumptions about the sources.
Digitized electrode coordinates were aligned using the nasion and auricular points as fiducial markers. Each patient's cortical surface was reconstructed from T1-weighted magnetization-prepared rapid acquisition gradient-echo (MPRAGE) data (Fischl, 2012). Head modeling utilized a three-layer boundary element method (BEM) model that was generated using the reconstructed cortical surface and fast low-angle shot (FLASH) MRI data, composed of the scalp, skull and brain with electrical conductivities of 0.33 S/m, 0.0042 S/m and 0.33 S/m, respectively (Hämäläinen and Sarvas, 1989). A three dimensional grid with 5 mm spacing was used to form the solution space. The forward solution was calculated by using the BEM. The inverse operator was computed from the forward solution with a loose orientation constraint of 0.6 to eliminate implausible sources and 2 μV as the estimate of EEG noise. The closest gray/white junction point corresponding proximally to the digitized location of each scalp electrode was found for each patient (Fischl, 2012) and the source activities were extracted from the midpoints of these 70 regions of interest (ROIs).
MPRAGE data were co-registered to the DTI data using an affine transformation for each subject. The electrode positions were then transformed into DTI space based on the coregistration matrix, and the alignment was verified visually. ROIs with a sphere radius of 16 mm were defined as the nearest main vertex in the white matter surface to each digitized electrode coordinate. This ROI size was determined based on the expected surface resolution obtained using high density EEG recordings (Nunez and Srinivasan, 2005). In this way, ROIs for structural network analysis were chosen to overlap with ROIs used for constructing the functional networks. Overlapping volumes between neighboring ROI pairs were manually removed prior to tractography analysis.
For white matter connectivity analysis, probabilistic tractography (Probtrackx2 through FSL 5.0.4/FDT -FMRIB's Diffusion Toolbox 3.0; FMRIB's Software Library) was used to process the DTI data, as described in detail in (Behrens et al., 2003;Behrens et al., 2007). Briefly, bedpostx was applied to create an estimation of diffusion parameters (e.g., distribution of principal diffusion directions) on a voxel-wise basis. Probtrackx was then used to repeatedly sample from the distributions, each time computing a streamline (inferred fiber tract) through these local samples, thus generating a connectivity distribution on the location of the true streamline (Behrens et al., 2007). By accounting for the inherent uncertainty in the distribution of diffusion vectors within each voxel, the probabilistic tractography method provides greater reliability of quantitative connectivity measurements than deterministic tractography methods.
User-specified ROI-ROI pair masks were used to calculate the connectivity distribution between a seed ROI and a target ROI bidirectionally for each ROI pair based on the EEG ESI location ROIs. A termination mask was applied to the target ROI because we were interested in quantifying the number of streamlines that reached the target ROI and not the number of voxels traversed by streamlines within the target ROI. In this way, individual streamlines were counted only once when they reached the edge of the target ROI. To normalize for target ROI volume, the mean number of streamlines per voxel were then computed. 1000 streamlines were sampled with a 0.2 curvature threshold with distance correction on to produce connectivity distributions representing the mean length of the streamlines multiplied by the number of streamlines that originated from the seed ROI and reach the target ROI. We note that all of our reported results were qualitatively similar when using distance correction off, but the number of long distance structural connections identified was negligible without distance correction. Because proximal nodes are known to be highly connected in both structural and functional networks due to true anatomical and physiological connectivity as well as spatial bias in measurement techniques (Chu et al., 2012;Li et al., 2012), we report our results with distance correction on in order to highlight correlations between long distance structural and functional connections. Using the number of streamlines launched from the seed ROI and the number of voxels in the target ROI allows for calculation of a Connectivity Index (adapted from McNab et al., 2013): The CI values for each ROI pair were placed in a connectivity matrix to generate the structural network for each patient. The same procedure was performed switching the seed and target ROIs and a symmetric connectivity matrix generated for each patient. The structural network processing stream is outlined in Figure 1 A-D.

Construction of functional networks
To assess the associations between activities recorded at two ROIs, we use two measures of linear coupling: cross correlation and coherence. For the correlation networks, the cortical current activity obtained using MNE were segmented into contiguous 1s intervals. We chose this window size to balance signal stationarity and accurate assessment of the coupling measure. Within each window, the data were first normalized from each ROI to have zero mean and unit variance before coupling analysis. The maximal cross correlation between all ROI pairs was then calculated, allowing a lag of +/− 100 ms. The choice of lag time was selected to encompass the duration of known neurophysiological processes and crosscortical conduction times (Varela et al., 2001;Premoli et al., 2014). This resulted in a single value (the maximum value of the cross correlation computer over lags of +/− 100 ms) assigned to each ROI pair. For the coherence networks, the same segments of source activity were analyzed. The coherence was computed using a Hann taper and treating each segment of source activity as a separate trial. In this way, the coherence measures phase consistency across time. The coherence was computed in the frequency bands: delta (1-4 Hz), theta (4-8 Hz), alpha (8-12 Hz), sigma (12-15 Hz), beta 1 (15-20 Hz), beta 2 (20-30 Hz), and gamma (30-50 Hz).
For inference of significant coupling, for each patient we performed a bootstrap procedure to create 10,000 surrogate coupling values for each ROI pair by shuffling the time intervals and ROI labels. More specifically, to generate a realization of surrogate data, 1s of data at each ROI was chosen (randomly, with replacement) from all 1s intervals and ROIs for that patient. This procedure preserved properties of the signals (e.g., the power spectrum) while disrupting the temporal and spatial ordering. For each interval of surrogate data created in this way, we computed a functional network using the same procedures as applied to the observed data (described in detail below). The resulting distribution of 10,000 coupling values was used to compute a p-value for the observed coupling of each ROI pair. Those ROI pairs with p-values below a threshold were assigned an edge to create a functional network (Kramer et al 2009). The threshold was determined by correcting for multiple comparisons using a linear step-up false detection rate controlling procedure with q = 0.05. For this choice of q, 5% of the edges in a functional network are expected to be falsely declared (Benjamini and Hochberg, 1995).
Functional networks were then computed for each patient using this statistical threshold to identify significant coupling and remove spurious coupling due to the inherent properties of the dynamic data. For the correlation measure, a binary network was then inferred for each 1s segment. In some analyses and for visual inspection, these binary networks were averaged across segments to create a representative weighted network reflecting the average properties of the data over time, where the edge weight or strength reveals the consistency of an edge appearance across time. We have previously demonstrated that this technique identifies core functional networks that remain consistent across days (Kramer et al., 2011, Chu et al., 2012. The calculation of the coherence differs in that this measure requires a notion of "trials". Here, we chose to let each 1s interval represent a "trial", and therefore the coherence assessed the phase relationship between two signals across the entire ensemble of 1s segments for a patient. The result is a single coherence network that represents the entire duration of data for a patient. A high coherence value in the resulting network represents a consistent phase relationship between two signals over time. Edges declared significant were assigned a weight or strength equal to the coherence value. We note that both measures (the averaged correlation network and the coherence) produce a single weighted network that identifies the most consistent functional relationships across segments for a patient. The functional network processing stream is outlined in Figure 1 A, E-H.

Network Measures
To evaluate direct and indirect anatomical connectivity, we selected a fundamental measure of graph structure -the path length -defined as the smallest number of edges traversed in traveling between two nodes, and here computed using the function reachdist.m from the Brain Connectivity Toolbox (Rubinov and Sporns, 2010). We note that all nodes are reachable in all networks considered here. The average path length is defined as the average of the path length over all node pairs. To compute the path length, we first converted the weighted structural and functional networks for each patient to a binary network by including an edge between nodes whose weight exceeded zero.
To analyze the similarity between two networks, we computed the normalized twodimensional (2D) cross correlation with zero shift between the two networks, each represented by its adjacency matrix. The two-dimensional cross correlation is a template matching procedure used commonly to compare similarities between two images. When applied to two images of equal size with zero lag (as used here), it provides a measure of the point-by-point similarity between the two images. This method is widely used to compare pixel time courses in fMRI imaging analysis (see for example, Hyde and Jesmanowicz, 2012) and thus was selected for this analysis. For each adjacency matrix, the scale, s, was calculated as the sum of its elements squared. The 2D cross correlation is then normalized by the square root of the product of s for each matrix. We note that, in the normalization, the diagonal elements of each adjacency matrix and the edges identified at zero lag are fixed to zero.
In order to interpret the similarity between the observed functional and structural networks, we compared them to random networks. Because network structure can be dramatically influenced by a network's degree distribution (Faust, 2007, Newman et al., 2001, random networks (configuration models) were generated by shuffling each patient's functional or structural network while preserving the degree distribution (Newman, 2010).

Statistical tests
In order to evaluate the relationship between functional connectivity, and path length and frequency band, main effects between were evaluated using a one-way ANOVA test. For all pair-wise comparisons, t-tests were performed.

Generalized linear models
To clarify the relationship between functional connectivity strength, structural connectivity strength, a network measure (anatomical path length), and inter-node linear distance, we implemented a generalized linear model (GLM). Generalized linear models (McCullagh and Nelder, 1989) have been utilized in many neuroscience contexts, especially in the analysis and characterization of spike train data (Brown et al., 2004;Truccolo et al., 2005;Czanner et al., 2008;MacDonald et al., 2011; and provide a means to evaluate the independent contribution of different predictor variables to the value of a response variable. Here, we related the response variable -the functional connectivity -to three predictor variables and their combinations: the structural connectivity strength (scaled to maximum of 1), the internode linear distance (scaled to a maximum of 1), and the structural path length between nodes (measured in integer units). For the correlation networks, we counted the number of times an edge occurred between two nodes across the networks inferred from each 1s segment. We chose a binomial distribution for the response variables (the integer number of edges detected across segments), and the logistic link, where μ is the expected value of the functional connectivity between two nodes, the design matrix X is a function of the predictors (the structural connectivities, the distances, and the path lengths), and β are the unknown coefficients to determine.
We constructed three GLMs to fit the functional connectivity as a function of the predictors. In the first, we assumed that the number of edges observed in the correlation network depends only on the distance between the sources; we label this the distance (D) model, which represents the null hypothesis that the functional connectivity depends only on distance. This is the simplest model we consider, and is consistent with the notion that functional coupling depends strongly on distance. To this baseline model we add additional features and examine whether these additions improve the model accuracy. In the second model, which we label the distance+structure (D+S) model, we assumed that the functional connectivity depends on both distance and the structural connectivity between nodes. Again, this model is consistent with the notion that structural connectivity will impact functional connectivity. Finally, in the third model, we incorporated a measure of the structural network -the path length (see Network Measures) -as a third predictor; we label this the distance+structure+path length (D+S+PL) model. We note that many other models -with additional predictors -may be considered. However, we focus here on these three models, with the goal of characterizing the relationship between the physical distances between sources and the structural connectivity on the inferred functional networks.
To compare the three models, we compute the Akaike information criterion (AIC), defined as where Δ is the deviance of the model, and n is the number of coefficients in the model. We note that the last term ("2n") causes the I to increase as the model complexity (i.e., number of predictors) increases. Therefore, a model with many predictors may have a higher AIC (and therefore worse performance) than a model with fewer predictors. For each model, we computed the coefficients β and deviance using the combined results for all patients, and chose the model with the smallest AIC. To compare AIC, we report the change in AIC, which is the AIC value for a given model minus the AIC of the best performing model. We note that in itself, the value of the AIC for a given data set has no meaning. But, when compared among a series of models, the model with the lowest AIC explains the data set best. For this best performing model, we confirm the model fit by checking that the Pearson residuals are approximately normally distributed (Kass et al., 2014).
For the correlation networks, the model performance was consistently poor when all observed functional connections were used; the Pearson residuals deviated strongly from the normal distribution when the weakest functional network edges were included (Supplementary Figure 1, Left). To focus the model and improve performance, we refined the fitting procedure by selecting only the most common functional connections for each patient, here defined as edges that appear in more than 10 percent of the 1s networks for a patient. Modeling this more restricted set of strong functional connections, the Pearson residuals better approximated the normal distribution. We note that the functional connections are not completely described (as measured by the proportion of deviance explained, see Results) with the three predictor variables evaluated here. This is consistent with existing experimental and modeling work which has shown that functional connectivity cannot be completely explained by underlying structural connectivity Damoiseaux and Greicius 2009;Rubinov et al., 2009;Baria et al., 2011;Goni et al., 2014;Ponten et al., 2010;Honey et al., 2007;Hlinka and Coombes, 2012). Higher resolution structural connectivity data and introduction of additional covariates (e.g., intracortical and subcortical connectivity patterns, etc.) may further improve model performance.
For the coherence networks, we followed a similar approach to implement a generalized linear model. For these networks, the edge values are continuous numbers from 0 to 1. To model these edge values, we chose in the GLM an inverse Gaussian distribution for the response variables, and the identity link. The inverse Gaussian distribution is consistent with the expectation of many small coherence values, and few large values. As for the correlation networks, the predictors are the distance between sources, the structural connectivity, and the path length, and we estimate the unknown coefficients β for the different models: distance (D), distance+structure (D+S), and distance+structure+path length (D+S+PL).

Functional network connectivity partially reflects the underlying structural network
On direct visual inspection, each patient's weighted functional (correlation and coherence) and structural adjacency matrices were non-random and had overlapping features ( Figure  2A). To evaluate the similarity of the network connections, we computed the twodimensional cross-correlation between the weighted functional correlation and coherence networks and the underlying structural network for each patient. Functional and structural networks were significantly more similar within patients than when compared to configuration models, which randomize the edge weights but preserve the degree distribution of the functional and structural networks (correlation and coherence networks for each frequency band evaluated p<0.00001; Figure 2B). For coherence networks, we found that the gamma band networks were significantly more similar to the structural networks than the delta (p=0.026) and theta band networks (p=0.047, Figure 2B). There was no detectable main effect of age on functional and structural network similarity (ANOVA, p=0.99) and no difference detected between data obtained in the wake versus sleep state (p=0.63). To further evaluate the relationship between each patient's functional and structural networks, we computed the average edge weight in the functional connectivity correlation networks for structurally connected and for structurally unconnected nodes. We found that within each patient, the functional connectivity between structurally connected nodes was significantly stronger than between structurally unconnected nodes (p<0.00001, Figure 2C). This finding held for the coherence networks in each frequency band evaluated (p<0.00001). In addition, we evaluated the similarity of functional and structural networks obtained across patients. We found some consistency in both functional and structural networks beyond random; however, structural networks were significantly more similar between patients than functional networks ( Figure 2B, p<0.00001), consistent with the notion that all patients share some common features of anatomical connectivity.

Connections are more common between nodes that are spatially close
For each patient, we found that both structural and functional edges were more common between physically neighboring nodes ( Figure 3). The data possess a broad distribution of distances between nodes, ranging from 0.6 cm to 11.9 cm, with a mean inter-node distance of 5.9 cm ( Figure 3A). Both the functional connectivity correlation network strength and the structural connectivity strength tended to decrease as the inter-node distance increased, with the decrease in the structural connectivity strength occurring more rapidly ( Figure 3B,C). Qualitatively similar results were found for coherence networks in all frequency bands evaluated. These results are consistent with the notion that brain connectivity tends to decrease with distance. However, we note that weak functional and structural connectivity appear across all inter-node distances; therefore, although spatial proximity strongly impacts network connectivity, distance alone is not sufficient to deduce the functional or structural association between two nodes.

White matter connectivity predicts functional connectivity beyond inter-node distance alone
Given the strong relationship between inter-node distance and both functional and structural connectivity (Figure 3), we implemented a generalized linear model to evaluate the relationship between structural and functional connectivity while accounting for distance (see Methods -Generalized Linear Models).
For the correlation networks, we evaluated 3 models that included the predictor variables: distance, structural connectivity, and path length. We found that the model that contains all three predictors (i.e., the inter-node distance, structural connectivity strength, and inter-node structural path length, the D+S+PL model) outperformed a model consisting of distance alone (the D model) and the distance and structural connectivity (the D+S model, change in AIC > 1100 in both cases). We note that an additional model that includes three predictorsthe distance, structural connectivity, and an interaction term between distance and structural connectivity -also did not improve performance. Compared to a model consisting of only a constant predictor, the proportion of deviance (or variability) explained by the addition of three variables (D+S+PL) is 0.21, and is highly significant (p < 0.00001, by a likelihood ratio test using the chi-squared distribution with three degrees of freedom). The D+S+PL model predicts that, as the (scaled) distance between two nodes increases by 0.1, the odds of a functional connection decreases by a factor of 0.79 (95% CI [0.77, 0.81]). However, as the scaled structural connectivity between two nodes increases by 0.1, the odds of a functional connection increases by a factor of 1.14 (95% CI [1.11,1.18]). Finally, an increase of path length by one edge decreases the odds of a functional connection by a factor of 0.81 (95% CI ([0.75, 0.87]). We note that distance and structural connectivity values are scaled from 0 to 1 and path length is measured in integer units ranging from 1-4. In summary, these GLM results show that the odds of a functional connection between two nodes decreases with distance, as expected. While accounting for this distance dependence, we also find that the presence of a structural connection increases the odds of a functional connection. Moreover, the longer the structural path length between two nodes, the lower the odds of a functional connection. Comparing the subsets of patients in the wake and sleep states, we find that the proportion of deviance explained by the addition of three variables (D+S+PL) to a constant predictor model is highly significant (p <0.00001) in both cases, and there was no significant difference in the deviance explained between states (p=0.11).
Finally, to evaluate the impact of indirect structural connectivity on functional connectivity strength, we evaluated the correlation network functional connectivity of structurally connected and unconnected node-pairs grouped by anatomical path length. We found a significant main effect between functional connectivity strength and path length (ANOVA, p<0.00001). Node-pairs supported by direct white matter connections (i.e., separated by a path length of one) had the highest functional connectivity values, and a step-wise decrease in connectivity was seen for each increase in anatomical path length (path length 1 versus 2, p<0.00001, path length 2 versus 3, p<0.005; path length 3 versus 4, p=0.18; Figure 4B). Thus, cortical functional connectivity was found to reflect both direct and indirect white matter connectivity. Qualitatively similar results were found for coherence networks at each frequency band (delta: p<0. 00001, p=0.032, p=0.22; theta: p<0.00001, p=0.039, p=0.36; alpha: p<0.0001, p=0.0057, p=0.42; beta 1: p<0.00001, p=0.026, p=0.34; beta 2: p<0.00001; p=0.0023, p=0.083, sigma: p<0.00001, p=0.071, p=0.31; gamma: p<0.00001, p=0.021, p=0.28).

Structural networks influence coupling of high frequency oscillations more than low frequency oscillations
To evaluate the impact of structural connectivity and inter-node distance on functional connectivity in differing frequency bands, we implemented a generalized linear model for the coherence networks. We evaluated the predictor variables: distance, structural network strength, path length, and combinations of these variables. The response variables were the functional network edge weights in the coherence networks observed for each frequency band: delta (1-4 Hz), theta (4-8 Hz), alpha (8-12 Hz), sigma (12-15 Hz), beta 1 (15-20 Hz), beta 2 (20-30 Hz), and gamma (30-50 Hz). We found that the distance+structure +path length (D+S+PL) model performed the best in all frequency bands (change in AIC > 130; range 130-1620). We show in Table 2 the results for the coefficient estimates for the D +S+PL model in each frequency band. Compared to a model consisting of only a constant predictor, the proportion of deviance explained by the addition of three variables (D+S+PL) ranges from 0.13 to 0.38 (Table 2), and is highly significant for all frequency bands (p <0.00001, by a likelihood ratio test using the chi-squared distribution with three degrees of freedom). The deviance explained by the (D+S+PL) model is significantly greater for the gamma networks than the correlation networks or the lower frequency coherence networks (correlation: p=0.0002; delta: p<0.000001, theta: p=0.0004; alpha: p=0.0007; sigma: p=0.0045; beta 1: p=0.0137; beta 2: p=0.17).
The GLM results also show that, as the frequency increases, the functional connectivity strength decreases. We note that the constant term in the GLM represents the expected value of the functional connectivity (here, the coherence), excluding the impact of the other covariates (distance, structural connectivity, and path length). We find that the values of the constant term in the GLM model (i.e., the term with a constant predictor of 1) are smaller in the higher frequency bands (e.g., gamma and beta 2) than in the lower frequency bands (e.g., delta, theta, and alpha). This result is consistent with notion that higher frequency bands exhibit decreased long-distance coupling (Singer, 1999;Miller et al., 2007;Tallon-Baudry, 2009;Chu et al., 2012). These results also show that, for all frequency bands, as the structural connectivity between two nodes increases, the functional connectivity increases. Notice that for the structural connectivity predictor (second column of Table 2), the 95% confidence intervals exceed 0 in all frequency bands. Finally, for all frequency bands, as the distance between two nodes increases, or the path length increases, the functional connectivity decreases. We note that, as the frequency increases, the impact of the structural connectivity becomes stronger and the impact of inter-node distance and path length become weaker (i.e., closer to zero). Comparing the subsets of patients in the wake and sleep states, we find that the proportion of deviance explained by the addition of three variables (D+S +PL) to a constant predictor model is highly significant (p <0.00001) for both subsets of patients in all frequency bands. The relationship between the predictors (distance, structural connectivity, path length) and the functional coherence is the same in the wake and sleep states, and there was no significant difference in the deviance explained by the model between wake and sleep states in any frequency band evaluated (p>0.15, range 0.15-0.91).
Consistent with the GLM results, we found that higher frequency networks had lower functional connectivity. A significant main effect was found between mean edge weight of the functional networks and frequency, with lower coherence seen in higher frequency bands (ANOVA, p<0.00001, Figure 5A).
In order to further evaluate whether the decrease in functional connectivity strength with higher frequencies was related to underlying structural connectivity, we computed the average functional connectivity strength for structurally connected and unconnected nodepairs, normalized by the average functional network strength present in the delta band. We found no difference in the average functional connectivity strength of structurally connected versus structurally unconnected node-pairs within each lower frequency band: theta, alpha, and sigma, compared to delta. However, in higher frequency bands (beta 1, beta 2, and gamma), the average functional connectivity strength decreased significantly compared to the slower frequencies when structural connectivity was absent (p<0.00001 for each comparison, Figure 5B). Consistent with our finding that gamma band networks have highest similarity to underlying structural networks ( Figure 2B), these findings suggest that the decrease in functional connectivity strength evident in higher frequency bands is primarily due to the loss of coherence between structurally unconnected nodes.

Discussion
Here we evaluated the relationship between functional and structural human brain networks using principled measures to infer networks from source estimates of scalp EEG and probabilistic tractography, enabling the comparison of cortical activity with high temporal resolution and across multiple frequency bands to underlying white matter connectivity. We have found that the coupling of brain activity in each frequency band is shaped by the observed underlying structural connectivity, including both direct and indirect paths. These associations highlight the obligatory interdependence of structural and functional networks in brain function. In addition, we found that high frequency oscillations exhibit sparse functional connectivity that is highly dependent on the existence of underlying structural connections.
Multiple studies have demonstrated a robust relationship between structural brain networks and spontaneous coupling dynamics in brain BOLD fluctuations (Baria et al., 2011;Damoiseaux and Greicius 2009;Rubinov et al., 2009;Goni et al., 2014). We have extended this work to show that cortical brain currents are likewise yoked to underlying anatomy. We have previously shown that long-term human EEG reveals persistent functional networks across days and states of consciousness (Chu et al., 2012). Here we show that persistent structures in functional brain networks reflect not only the impact of spatial proximity, but also the direct structural connectivity between brain regions, as well as indirect structural topology as measured through path length.
Our finding of increased functional and structural connectivity between neighboring brain regions is consistent with prior work in animals and humans. Increased connectivity between spatially proximal and functionally related units has been observed in structural networks across spatial scales ranging from microns to centimeters (Bullmore and Sporns, 2009;Braitenberg and Schuz, 1998;Hellwig, 2000;Hagmann et al., 2008). Similarly, evaluation of human brain networks with fMRI has demonstrated increased functional connectivity between anatomically proximal regions (Bullmore and Sporns, 2009;Salvador et al., 2005;Achard et al., 2006). This modular network architecture has been proposed to both increase computational efficiency and decrease energy costs by organizing highly clustered modules that can be linked by few long distance connections (Cherniak, 1994;Sporns, 2011).
How the brain's relatively static structural scaffolding molds dynamic information processing and sophisticated cognitive processes remains an active question in neuroscience. Many studies have posited that coupling of complex signals between and across different frequency bands support and encode the integration of simultaneous brain processes (Buzsáki, 2004;Buzsáki and Draguhn, 2004;Fries, 2005;Singer, 1999;Fingelkurts and Fingelkurts, 2004;Laufs, 2008;Bullmore and Sporns, 2009;Deco et al., 2011). Here we show that long-range integration of cortical activity in each frequency band is biased toward anatomically linked regions. These observations suggest that cortico-cortical connections provide a reliable physical substrate for the long-range transmission of composite signals with high dimensional frequency content (Lakatos, et al., 2005;Klimesch, 1996;Friston, 1997;von Stein, et. al., 2000;Varela et al., 2001;Laufs, 2008;He and Raichle., 2009).
Here, we applied two measures of functional connectivity: the correlation (a broadband measure) and coherence (a narrowband measure). Although both measures produced qualitatively similar results, consistent with our finding that gamma coherence networks are most similar to underlying structural networks, the best performing model (D+S+PL) accounted for the most deviance in the gamma coherence networks. Increasing evidence suggests a spatial gradient in which slower oscillations couple more distributed brain regions, while faster oscillations are more focally distributed (Singer, 1999;Miller et al., 2007;He and Raichle, 2009;Tallon-Baudry, 2009;Baria et al., 2011). Our observations that faster oscillations are dependent on underlying structural connectivity for long-range integration complement this growing body of literature. In contrast, long-range cortical integration between lower frequency oscillations could be preferentially facilitated by alternate processes, such as high amplitude traveling waves (Ermentrout and Kleinfield., 2001;Sato et al., 2012), or shared subcortical pathways that synchronize neocortical delta activity (Steriade et al., 1993).
The observed relationship between cortical currents and underlying structural connectivity reported here likely underrepresents the true relationship between brain structural and functional connectivity, due to the limitations of current imaging and measurement techniques. For example, unmyelinated axons are not routinely reconstructed using current tractography techniques, including intra-cortical axonal and dendritic processes. Furthermore, inter-regional coupling mediated via shared subcortical sources were not evaluated. We focus here on observed long-distance connectivity patterns between superficial cortical regions in order to emphasize brain regions that can be anatomically approximated with surface EEG. In addition, current tractography techniques are known to underestimate inter-hemispheric connectivity, which is substantial between homotopic brain regions; this limitation was mitigated by using probabilistic tractography with multiple fiber orientations (Behrens et al., 2007). Furthermore, the b-value utilized for DTI in this study was relatively low, and diffusion data with higher angular resolution would improve accurate tract identification, especially in regions of crossing fibers (Setsompop et al., 2013). Finally, we evaluated young patients across a wide range of ages with epilepsy due to varying etiologies, who may have age or disease dependent alterations in local structural (Liu et al., 2014;Ji et al., 2014;Barkovich et al., 1988;Ashtari et al., 2007) and functional (O'Muircheartaigh et al., 2012;Bartolomei et. al., 2013;Chu et al., 2014) networks at baseline. We note that in spite of any individual variations expected in this population, our finding that large-scale brain structural networks support overlying brain functional connectivity in a frequency-dependent manner was evident in each patient, suggesting a robust finding. However, future work should evaluate whether these findings persist in a healthy population, across the lifespan, and in other disease states. In addition, detailed evaluation of these relationships at the site of known structural or functional lesions could identify subtle variations that were not investigated in this global analysis.

Conclusion
This work shows that coupling between human cortical brain dynamics partially reflects observed white matter connectivity across multiple brain rhythms. The increased dependence on structural support for the coherence of faster brain rhythms provides new evidence for how underlying anatomy directly shapes emergent brain dynamics at fast time scales. Although brain structure has not been shown to fully predict overlying cortical physiology, capturing the influence of brain structure on spontaneous brain activity provides new opportunities to manipulate these relationships to alter brain function and treat neurological disease.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.

Highlights
• Both structural and functional edges are more common between physically neighboring nodes.
• Direct and indirect white matter connectivity predicts functional connectivity beyond inter-node distance alone.
• Functional connectivity strength is higher in structurally connected node pairs in all frequency bands evaluated.
• Lack of structural connectivity disproportionately reduces functional connectivity in high frequency bands.  Examples of structural and functional adjacency matrices from one patient. Similarity between structural network architecture and cross-correlation and coherence functional networks is visually evident. B. The network similarity as measured by the two-dimensional cross correlation (+/− standard error) is plotted for configuration models with degree distribution preserved (black), each patient's structural and functional correlation network (light gray), all patient's structural and functional coherence networks (colored bars), and between all patients functional networks (medium gray), and all patients structural networks (dark gray). Patients' structural and functional (correlation and coherence) networks are significantly more correlated than random (p<0.00001). Coherence networks in the gamma band were more similar to structural networks than coherence networks in the delta (p=0.026) and theta (p=0.047) bands. Between patients, both structural and functional correlation networks were more similar than random networks (p<0.00001). Structural networks were significantly more similar across patients than functional networks (p<0.00001). C. Average correlation network functional connectivity edge strength (+/− standard error) is plotted for each patient for structurally connected (blue) and structurally unconnected (red) node pairs. Mean functional connectivity between structurally connected nodes is significantly higher in each patient (p<0.00001). Similar findings occur for coherence networks in each frequency band evaluated (p<0.00001). Inter-node distance ranged from 0.2 cm to 9.4 cm with a grossly normal distribution. B, C. Functional correlation network (B) and structural network (C) connectivity strength decreases with internode distance, with strong connections in both structural and functional networks tending to occur between nodes that are physically close (< 3 cm). In these figures, each dot indicates a single edge (or equivalently, single node-pair) for each patient (in color, see legend). Average correlation network functional connectivity strength (+/− standard error) is plotted for structurally connected (blue) and structurally unconnected (red) node pairs for three inter-node distances: <3 cm, 3-6 cm, >6 cm. Mean functional connectivity strength between structurally connected nodes is significantly higher at short distances (p<0.045) and medium distances (p=0.022), and tends to be higher at long distances (p=0.057). The relative paucity of data points available in longer distance bins may contribute to the lack of significant effect. Similar results were found for coherence networks (not shown). B. Average correlation network functional connectivity strength (+/− standard error) is plotted for node pairs by anatomical path length. Node pairs with direct white matter connectivity (path length 1) had significantly higher functional connectivity than node pairs connected by longer path lengths (p<0.00001 for all comparisons). This relationship tended to persist for each incremental increase in path length (path length 1 versus 2, p<0.00001, path length 2 versus 3, p<0.005; path length 3 versus 4, p=0.18). Qualitatively similar results were found for coherence networks in all frequency bands (not shown). Average coherence network functional connectivity edge strength (+/− standard error) decreases at higher frequencies (ANOVA, p<0.00001). B. Average functional connectivity strength (+/− standard error) normalized to the delta band is plotted for structurally connected (blue) and structurally unconnected (red) node-pairs for coherence networks in 7 frequency bands. Functional connectivity decreases more rapidly between structurally unconnected node pairs in the higher frequency bands. Structurally unconnected nodes (red) have significantly lower functional connectivity compared to structurally connected nodes (blue) in beta and gamma frequencies (B1, B2, and G, p<0.00001) but not lower frequency bands. In both subfigures: delta (D), theta (T), alpha (A), sigma (S), low beta (B1), high beta (B2), and gamma (G).