Whole-brain analytic measures of network communication reveal increased structure-function correlation in right temporal lobe epilepsy

The in vivo structure-function relationship is key to understanding brain network reorganization due to pathologies. This relationship is likely to be particularly complex in brain network diseases such as temporal lobe epilepsy, in which disturbed large-scale systems are involved in both transient electrical events and long-lasting functional and structural impairments. Herein, we estimated this relationship by analyzing the correlation between structural connectivity and functional connectivity in terms of analytical network communication parameters. As such, we targeted the gradual topological structure-function reorganization caused by the pathology not only at the whole brain scale but also both in core and peripheral regions of the brain. We acquired diffusion (dMRI) and resting-state fMRI (rsfMRI) data in seven right-lateralized TLE (rTLE) patients and fourteen healthy controls and analyzed the structure-function relationship by using analytical network communication metrics derived from the structural connectome. In rTLE patients, we found a widespread hypercorrelated functional network. Network communication analysis revealed greater unspecific branching of the shortest path (search information) in the structural connectome and a higher global correlation between the structural and functional connectivity for the patient group. We also found evidence for a preserved structural rich-club in the patient group. In sum, global augmentation of structure-function correlation might be linked to a smaller functional repertoire in rTLE patients, while sparing the central core of the brain which may represent a pathway that facilitates the spread of seizures.


Introduction
Characterization of the structure-function relationship in brain networks is crucial to the understanding of both normal and diseased brain organization. This relationship can be defined by correlating functional and structural connectivity (Honey et al., 2009). The complete set of structural or functional links between brain areas has been referred to as the structural or functional connectome (Sporns et al., 2005) and can be analyzed by graph theoretical measures (Achard et al., 2006;Sporns, 2013). Graph theory holds promise as a method to reveal pathological reorganization of brain networks (Fornito et al., 2015;Guye et al., 2010). When the network is altered in pathologies, localized structural changes in the brain are likely to result in complex direct and indirect functional alterations at various different scales and levels of integration.
This propagation of local changes through the network renders the structure-function relationship highly non-trivial. Here we study the network changes in temporal lobe epilepsy (TLE), whilst taking into account the complexity of indirect structure-function correlations (Damoiseaux and Greicius, 2009) by means of a graph analytic network communication approach (Goñi et al., 2014). TLE is now regarded as a pathology involving areas not only participating in a localized epileptogenic network (producing transient electrical abnormalities) but also having long-lasting large scale network effects over the whole brain (Bartolomei et al., 2001;Laufs, 2012;Spencer, 2002). These effects have been revealed and quantified mostly by in vivo MRI studies, accessing both the structural and the functional organization of the whole brain (Bernhardt et al., 2013).
Structural alterations of a priori brain subnetworks have been assessed in epileptic patients using graph analysis, finding widespread subcortical alterations in epileptic subnetworks (Bonilha et al., 2012). In contrast, a more global approach comparing whole brain structural connectivity between patients with bilateral TLE and controls identified a network of structural connectivity alterations in the temporal lobe distributed to different networks for left and right TLE (Besson et al., 2014). Substrates of these structural alterations are still debated, though studies employing post mortem histology in TLE show the presence of white matter damage, neuronal loss and gliosis both within (Thom et al., 2000(Thom et al., , 2001 and beyond epileptogenic regions (Blanc et al., 2011), most likely caused by the propagation of recurrent epileptic seizures throughout the brain.
From a functional perspective, rsfMRI studies in epileptic patients reported various, widespread and complex patterns of increased and decreased functional connectivity, effecting key features of large scale networks such as the default mode network (Haneef et al., 2012;Liao et al., 2010;Pittau et al., 2012;Voets et al., 2012). Discrepancies in reported functional brain reorganization could be related to the broad range of methodological approaches used in these studies (e.g. whole brain network corrected analysis vs. a-priori ROI analysis), which make it difficult to reach a synthetic interpretation of these results (Centeno et al., 2014;Constable et al., 2013).
In epilepsy, modeling the brain as a graph can be of use in assessing the link between structural and functional connectivity (Bernhardt et al., 2013;Guye et al., 2008). Several recent studies have endeavored to combine structural and functional approaches in epilepsy (Chiang et al., 2015;Douw et al., 2015;Vaessen et al., 2014;Voets et al., 2012;Zhang et al., 2011). Function-structure relationship has been studied in the epileptic brain in generalized epilepsies (Zhang et al., 2011), frontal lobe epilepsy (Vaessen et al., 2014) and in TLE (Voets et al., 2012). While focusing only on limbic and pre-limbic structures Chiang et al. (2015) found decreased structural-functional relationship in TLE. Douw et al. (2015) observed higher anatomical between-module connectivity which correlated with default mode disintegration in TLE compared to a control group.
Structural parameters of brain networks such as Euclidian distance or the combination of more complex graph analytical metrics of network communication such as path length, search information and path transitivity have also provided good predictions of functional connectivity in the normal brain (Goñi et al., 2014). Considering a structural link in the context of its neighbors makes structure-function analysis more sensitive to accumulated small structural changes which trigger widespread functional connectivity changes across the network. As such, we hypothesize that these metrics could help to better characterize structural-functional reorganization in TLE, particularly in the context of epilepsy where dMRI tractography cannot grasp the full spectrum of structural change induced by the pathology. For example, gradual linear structural changes (secondary effects of both ictal and interictal epileptic phenomena) such as demyelination may not be visible to probabilistic tractography, but by their distributed effect on functional connectivity they may be discernable via changes in structurefunction relationship. The role of highly connected regions has also been stressed out recently. Indeed meta-analysis of graph theoretical studies has revealed that most of pathological lesions in the brain are linked to 'hubs' (Crossley et al., 2014). In the same line, recently, a "rich club" architecture has been proposed to describe the highly inter-connected core of hubs inside the brain (van den Heuvel and Sporns, 2011), which was found to be modified by pathology (Fornito et al., 2015;van den Heuvel et al., 2013). This structural core overlaps with propagation networks of epileptic activity in TLE, for example the precuneus (Arthuis et al., 2009) and as such might undergo structural and functional alterations in TLE patients, though this has yet to be demonstrated.
To our knowledge there is no study that directly evaluates the edgewise structure-function relationship at the whole brain scale, as well as within the rich-club and pathologically altered functional sub-networks in patients while taking into account network-based communication processes. In this context, we aimed to determine structure-function relationship within altered brain networks observed in rTLE. To do so we used a framework of several analytic approaches at different scales in a group of rTLE patients compared to a group of age and sex matched controls. First, structural and functional connectivity derived from dMRI/rs-fMRI data were used to examine whole brain reorganization in rTLE by using the Network Based Statistics (NBS) approach (Zalesky et al., 2010a), and by extracting the rich club structure of the brain. Second, to characterize whole brain functional-structural relationship in altered brains of rTLE patients, we sought to predict the functional connectivity from i) graph-analytical edgewise metrics derived from structural data (Goñi et al., 2014), and ii) Euclidian distance, and by comparing these results against current analytical and simulation models (Robinson et al., 2014). Third, we aimed to analyze structural-functional relationship within the altered functional subnetworks, richclub networks and peripheral networks outside the rich club.

Subjects
Seven patients diagnosed with drug-resistant epilepsy of the right temporal lobe (rTLE, 4 males, mean age 31.8, range 19-50, 6 right handed, 1 left handed, for heterogeneity we selected only right lateralized patients, for detailed information see Supplementary material 1) and fourteen healthy subjects with no history of neurological disease were recruited into the study. The participants signed an informed consent form according to the rules of the local ethics committee (Comité de Protection des Personnes (CPP) Marseille 2). One control was excluded due to excessive head motion resulting in diffusion artifacts (resulting in a final control group of 13, mean age 31.8, range 20-59, 7 males, 12 right handed, 1 ambidextrous).

MRI acquisition
The participants were scanned on a Siemens Magnetom Verio 3 T MR-Scanner (Siemens, Erlangen, Germany). 350 functional MRI images were acquired in a BOLD-sensitized EPI T2*-weighted sequence with a TR of 3.6 s (2.0 × 2.0 × 2.5 mm, TE = 27 ms, 50 slices, FA = 90°), resulting in a total fMRI time series of 20 min. During the resting state protocol the subjects were asked to keep their eyes closed and not to fall asleep. The dMRI-sequence was acquired with the following parameters: angular gradient set of 64 directions, TR of 10.7 s (2.0 × 2.0 × 2.0 mm, TE = 95 ms, 70 slices, b weighting of 1000 s/ mm 2 ). T1-weighted anatomical images were acquired with a MPRAGE-sequence (TR = 1900 ms, TE = 2.19 ms, 1.0 × 1.0 × 1.0 mm, 208 slices).

Diffusion MRI preprocessing
To correct for head motion effects, the gradient direction matrix was rotated using a customized in-house algorithm (Leemans and Jones, 2009;Raffelt et al., 2012). Next, to reduce spatial intensity inhomogeneities, bias correction was performed on the b0 image and subsequently applied to all other diffusion volumes (Sled et al., 1998). Lastly, a Higher Order Model Outlier Rejection model (Pannek et al., 2012) identified voxels with residual outliers in the diffusion-weighted signal.

Fiber tracking
The fiber orientation distribution (FOD) function was estimated within MRtrix software , by performing constrained spherical deconvolution (CSD, lmax = 8) of the diffusion signal (Tournier et al., 2008) within single-fiber (FA N 0.7) populations. The probabilistic streamline algorithm (iFOD2) (Tournier et al., 2010) was employed to reconstruct 2.5 million whole-brain fiber tracks of high-angular resolution. As streamlines were generated from random seeds throughout the brain volume, two tractography runs (i.e. 2 × 2.5 million) were performed to ensure comprehensive brain coverage. A set of plausible fiber trajectories are generated by iFOD2 by random sampling from the orientation likelihood inherent in each FOD along each candidate path. The default tracking parameters were employed for the acquisition parameters of the present study (step size = 1.25 mm, minimum length = 20 mm, max length = 200 mm, FOD termination threshold = 0.1, curvature constraint = 1 mm).

Structural connectome construction
Anatomical nodes of interest were generated by parcellating the standard AAL template (Tzourio-Mazoyer et al., 2002) into 512 cortical and sub-cortical regions (Supplementary material 2, see also (Perry et al., 2015)) of approximately uniform size (Zalesky et al., 2010b). In order to transform parcellations into subject-space the parcellated template was first co-registered to the MNI T1 2 mm brain template. The parcellated template (in MNI space) was subsequently transformed into subject-space by applying the transformation matrix generated from registering the MNI template to the subject's FA image. All co-registrations were affine linear registrations (to better conserve uniform regions size we did not use non-linear transformations) employed within the FSL software package (FSL 5.0 FLIRT, http://fsl.fmrib.ox.ac. uk/fsl/fslwiki/, (Smith et al., 2004)).
Within a weighted graph G w , a weighted connection w ij represents the number of streamlines from region i terminating within j. Because we want to use Euclidian distance as a covariable in our later analysis, here we do not divide the tract count by distance (Cheng et al., 2012). For the resulting number of streamlines of the two tractography runs we excluded every connection w ij with zero streamlines in at least one run or a difference of streamlines larger than three standard deviations (according to the variation across the whole individual connectome) between the two runs. This resulted in an individual structural connectome matrix with average sparsities of 0.46 (control group) and 0.44 (rTLE group, control-patient difference non-significant p N 0. 05, Wilcoxon rank-sum test). Our main interest is to keep a maximum of information for measurement of whole brain structural-functional relationship so further thresholding is only applied were absolutely needed (see below).

Functional connectome construction
We used wavelet analysis on the resulting region-averaged time series (Brainwaver toolbox, version 1.6, http://cran.r-project.org/web/ packages/brainwaver/index.html) keeping only the time series of wavelet coefficients of the second wavelet scale which represent the frequency band from 0.03 Hz to 0.07 Hz (Achard et al., 2006) for a TR of 3.6 s. Finally we created a functional connectivity matrix by calculating the Pearson-correlation between each region's wavelet coefficient time series.

NBS
To compare the functional connectivity matrices of the two groups, we used network based statistics (NBS), a statistical test to correct for multiple comparisons in a network which generates a p-value by comparing network sizes of permuted random group samplings (Zalesky et al., 2010a). Correlations and anticorrelations between regional BOLD signals can both be regarded as a functional connection (Achard et al., 2006). To measure the functional connectivity strength we thus took the absolute correlation value for statistical testing. According to literature the single-connection t-value threshold was arbitrarily set sufficiently high (t ≥ 5) to avoid problems raised by low t-thresholds (as discussed for t b 4 by Smith and Nichols (2009)). The corrected pvalue of NBS was calculated by finding the maximum network size of connectivity changes between 10,000 times randomized group labels (networks with p b 0.05 were considered to be significant, NBS implementation of the ConnectomeViewer toolbox (V2.0.0): http://cmtk. org/viewer/).

Graph analytical measures
To analyze the topological properties of structural connectivity we used the weighted path length (the inverse of streamline count); search information (a metric which penalizes the shortest path by the number of edges branching off in any of the nodes along the path); and path transitivity (a metric taking into account the weight of the shortest path and all paths having one supplementary step compared to the shortest path to reach the target region). Path length, search information and path transitivity metrics are illustrated in Fig. 1. Detailed formulae can be found in (Goñi et al., 2014).
In addition to graph analytical predictors of functional connectivity we defined the following metrics: Euclidian distance (distance in mm between regions centers), streamline count and an analytic prediction of functional connectivity (FCA) (Robinson et al., 2014). The FCA was defined by a linear model estimating the edgewise covariance Cov from the structural weight matrix W by with I being the identity matrix with the same size of the weight matrix W and c being the free coupling parameter. FCA of the covariance matrix is obtained by transforming the edgewise covariance estimates into correlation values: The coupling parameter c was defined in respect to the maximal correlation between FCA derived from the averaged structural connectivity and averaged functional connectivity of the control group (c = 0.32). This coupling parameter was also applied to calculate the FCA of the rTLE group.

Rich-club analysis
To compare whole brain differences the averages over the whole structural connectivity matrix of each metric were compared between the groups (Wilcoxon rank-sum test). The rich club is a high strength core of hubs which are more densely connected to each other than expected by chance. To extract the rich club of each group we used the group averaged structural connectome. To exclude effects from varying sparsity on rich club analyses (Ercsey-Ravasz et al., 2013) and to have a sparsity comparable to earlier rich club studies (van den Heuvel et al., 2013) we limited the connections of each subject to the 30% strongest connections in the structural connectome. The weighted rich club coefficient is defined as: where k refers to the nodal degree; W Nk the summed weights of all links greater than the degree threshold k; E Nk is the number of links in the node subset given by threshold k, and; w ranked the vector of the ranks of all weights in the observed network having a degree greater than k. Thus ϕ w (k) is the ratio between the summed weights and the sum of ranks of each link in the subset thresholded by k (see (Opsahl et al., 2008;van den Heuvel and Sporns, 2011) for details). A rich club coefficient for each averaged connectome was calculated and normalized according to the rich club coefficient distribution of 10,000 random networks with equivalent nodal degree distributions as implemented in the brain connectivity toolbox (https://sites.google. com/site/bctnet/, version: bct-cpp r396) (Rubinov and Sporns, 2010). We characterized the profiles of normalized rich club coefficients according to a sliding degree threshold (150 b k b 350) and extracted rich club networks at a minimum degree threshold of 250, 300 and 330 connections. The lower border was selected for comparability to previous descriptions of rich club regions (van den Heuvel and Sporns, 2011). The upper limit was defined by the maximal degree of the two connectomes.
To better understand the robustness of control and patient subnetworks at high thresholds we gradually increased the threshold k and compared the size of each subnetwork in patients and controls, from which it was possible to obtain a p value by permuting of group labels 10,000 times.

Role of lesions in structural and functional connectome
Previous work has shown that network organizations in TLE patients can differ depending on the presence of brain lesions for both structural and functional connectivity (Campos et al., 2015;Coan et al., 2014). To rule out the possibility that the network reorganization observed here are caused by the underlying lesion type we show that search information is augmented even in patients with normal MR (Supplementary material 3). Additionally we validated the observed functional NBS subnetwork by showing that almost all connections in the subnetwork in every subject represented higher correlations than the control mean (except five individual connections which were found to be smaller in different subjects), with no differences for lesional or non-lesional patients (Supplementary material 3).

Correlation of structural and functional networks
We defined structure-function relationship as the Pearson correlation between functional connectivity and the various metrics derived from the underlying structure. To analyze the structure-function relationship we correlated these edgewise metrics with edgewise functional connectivity. This network communication approach has been shown to describe the structure-function relationship through high correlations between structural and functional connectivity (Goñi et al., 2014). Correlation of the Euclidian distance enables the analysis of gradual change in functional connectivity as a function of topological distance representing structural change accumulating along the streamlines particularly in long-range highly myelinated structures not quantified by streamline count. Accounting for distance is also necessary to control for the strong relationship between tract length and tract count (Ercsey-Ravasz et al., 2013). Correlation with analytical functional connectivity (FCA) can be used to estimate functional connectivity between regions with no direct structural connections. This enables us to identify the main factors of structure-function relationship organization in rTLE. To better predict functional connectivity we averaged all edgewise metrics over each group and correlated them with the averaged functional data (Goñi et al., 2014;Honey et al., 2009). To test for differences in correlation values and the slopes of the linear regressions we permuted group labels and recalculated mean group connectomes 10,000 times. Additionally to these metrics we tested stepwise multi-parametric linear regression models based on combined streamline count, path length, path transitivity and search information with and without accounting for the Euclidian distance.
Finally, we systematically applied streamline count, search information and Euclidian distance predictors from our previous analysis to investigate structure-function relationship in four subnetworks: an NBSderived network of altered functional connectivity; and three structurally defined networks in the form of the rich club, feeder and peripheral/local networks. These structural networks were classified into three different classes according to its relationship to the rich club (at the selected degree threshold of k N 250, determined by best fit with previous rich club work (van den Heuvel and Sporns, 2011)). The first class consists of connections interconnecting rich club regions, the second class consists of feeder connections connecting rich club regions with non- Fig. 1. Illustration of graph metrics to characterize features of the shortest path from A to B: a) weighted path length (summed connection strengths); b) search information (weight path length by paths branching off in one of the shortest path's nodes) and c) path transitivity (path lengths weighted by additional triangle detours which can be used to arrive at B). rich club regions and the third class consists of peripheral local connections connecting non-rich club nodes with other non-rich club nodes.

Functional connectivity in rTLE
NBS analysis of the functional connectome revealed a widespread hypercorrelated network defined through absolute values of correlation coefficients which were significantly higher in the rTLE group than in controls (NBS parameters: t ≥ 5, p(corrected) ≤ 0.05, Fig. 2, Supplementary material 4). Main parts of the network consist of fronto-occipital and fronto-temporal connections with a dominance of connected nodes in the left hemisphere both connected inter and intra hemispheric.
Hubs of the altered NBS subnetwork (degree ≥ 5 considering only subnetwork connections) were bilateral middle frontal gyri, superior frontal gyri and bilateral supplementary motor area, left lingual gyrus, precuneus, superior medial gyrus, middle temporal gyrus and left calcarine sulcus. Connections between those hubs can be split into a short-ranged dense locally-connected frontal network and a widespread long range connected network with the strongest connection between bilateral frontal lobes to bilateral lingual gyri but also to the right superior temporal gyrus and left middle temporal gyrus. Additionally we observed a smaller altered network of long-range connections between bilateral hippocampi, the left precuneus, thalamus and left calcarine sulcus (Fig. 2).

Structural connectivity in rTLE
The mean values of the streamline counts, path length, FCA, path transitivity and mean structural nodal degree did not differ significantly between the two groups (Supplementary material 1). Only average search information differed significantly between the controls and rTLE patients (p = 0.007 Bonferroni corrected, Wilcoxon rank sum test).
At a degree threshold of 150-250 we observed mostly equivalent regions in patients and controls and a similar normalized rich club coefficients profile (see Fig. 3a, b. Fig. 3c shows the rich club subnetwork size as a function of the nodal degree threshold to define this subnetwork for both controls and rTLE patients. From a degree threshold of 300 connections rTLE patients showed a significantly larger rich club subnetwork than observed in controls. The rich club started to disintegrate at a degree threshold of 300 for the control group and at a threshold of 330 for rTLE patients. Fig. 3a shows rich club regions at different nodal degree thresholds of 250, 300 and 330.
At a degree threshold of 250, rich club regions for controls included bilateral precentral gyri, superior frontal gyri, superior orbital frontal gyri, insulas, superior parietal gyri, precuneus and bilateral putamen, right superior temporal pole and right support motor area, left inferior triangularis and left inferior orbital frontal gyrus. At the same degree threshold of 250 for rTLE patients, the rich club did not include the left superior frontal gyrus and the bilateral orbital frontal gyri, but rather the left hippocampus, left superior temporal pole, left rolandic operculum, left supplementary motor area and right calcarine were part of the rich club.
At a degree threshold of 300, the rich club of controls included bilateral precuneus, right thalamus and left inferior frontal triangularis. For rTLE patients, the rich club was composed by bilateral precuneus, putamen and bilateral thalami, left precentral gyrus, inferior frontal triangularis and insula and right superior parietal gyrus.
At a degree threshold of 330 no regions could be considered to belong to the rich club in controls while in rTLE patients, rich club properties were found for the right precuneus, putamen and thalamus. Thus, at Fig. 3. a) High degree rich club nodes and their intra-connections at a degree threshold of k = 250, 300 and 330 connections for controls and rTLE patients forming a rich club at different levels (for better visualization 512 regions are projected on the 90 regions of the AAL atlas); b) Rich club coefficient for averaged connectomes (connectome sparsity = 30%, rich club coefficient weighted by streamline counts); c) Subnetwork size (number of nodes) of the largest component as a function of the nodal degree threshold (connectome sparsity = 30%), grey background marks significant difference between the two groups (p b 0.05, two-sided t-test, 10000 permutations of group labels to build the averaged connectome).

Table 1
Pearson correlation between functional connectivity and structural derived metrics in the whole brain and in subnetworks (all predictors were significant with P b 10 −6 , group averaged connectivity matrices, NBS: altered NBS subnetwork of functional connectivity rTLE N controls, Rich Club, Feeder and Local: rich club, feeder and local connections for structural connectome sparsitiy 30%, and nodal degree N 250); FCA: analytic functional connectivity, ED: Euclidian distance, PL: Path length, PT: Path transitivity, SI: search information, PP path predictors (path length, path transitivity, search information).

Structure-function relationship in rTLE
Regression of using each structural connectivity metric as a single predictor (streamline count, FCA, path length, path transitivity, search information, Euclidian distance) showed significant correlations with functional connectivity (P b 10 −6 , Table 1). Euclidian distance had the highest correlation values both in control and in rTLE patients. Multiparametric stepwise regression including all graph analytical parameters (streamline count, path length, path transitivity, search information and Euclidian distance) showed that each of the parameters contributed significantly to the model (P b 10 −6 ).
Comparing the relative strength of the structure-function relationship between controls and rTLE patients using permutation testing showed significant differences (10,000 permutations, p b 0.05) in correlation for all metrics and significant differences for slope values of simple regression models for FCA, path transitivity, search information, Euclidian distance (Table 2, Fig. 4).
The altered functional NBS subnetwork comparing rTLE and controls showed higher structure-function correlation for all metrics (Table 1) compared to the whole brain analysis. The Euclidian distance was also the best predictor of structure-function relationship but contributions of search information and streamline count remained significant when controlling for Euclidian distance (Fig. 5). The rich club edges had showed an increased correlation for Euclidian distance and path length, while the other metrics decreased in comparison to the whole brain analysis. The streamline count was more highly correlated in the peripheral network may be because of the use of a fixed sparsity of 30% which excludes connections with streamline counts of zero, in contrast to the more inclusive whole brain analysis. Correlation of the other metrics was lower than in whole brain analysis.
Slopes of the linear fit differed between the two groups for the altered NBS subnetwork and the peripheral network. Neither the rich club nor feeder connections showed a different slope of structure-function correlation between controls and rTLE patients (Table 2). We did not estimate multiparametric models because both the altered NBS subnetwork and the rich club subnetwork lacked sufficient data to obtain significance of all regression coefficients in the model (P b 0.001).
Looking at structure-function correlation differences between rTLE patients and controls within the rich club, feeder and peripheral networks respectively the only significant difference was observed for peripheral connections (Fig. 5).

Discussion
In this study we systematically investigated the descriptive and explanatory power of graph theory and network communication measures to characterize the structure-function relationship in the normal and epileptic brain. By focusing on right lateralized TLE we emphasized homogeneity across patients to allow for a maximum of comparability of our findings. In rTLE patients, on a whole brain level we observed a general pattern of increased functional connectivity, increased global search information, a preserved rich club, as well as an increased structure-function correlation for all metrics compared to healthy controls. Notably the prediction of functional connectivity by structural parameters such as Euclidian distance and search information capture structure-function reorganization in rTLE. Compared to peripheral connections the structurally preserved rich club shows no differences in structure-function relationship.

Reorganization of the functional connectome in the rTLE group
NBS-analysis revealed a widespread hypercorrelated network for rTLE patients in comparison to controls mainly consisting of fronto-occipital and fronto-temporal connections. This altered NBS subnetwork overlapped with default mode network regions such as the precuneus and precentral gyrus, middle temporal gyrus and the temporal pole with more connections in the left hemisphere. This is in line with previous whole brain pairwise analysis of rTLE which also found connectivity increases in the left hemisphere (Su et al., 2015). Though concerned only with pairwise connectivity data the current work is in line with previous graph theoretical demonstrations of topological rewiring in temporal lobe epilepsy and widespread functional connectivity changes for ipsi-and contralateral temporal lobe and default mode network regions for rTLE patients (Liao et al., 2010;Voets et al., 2012). Our results differ regarding the previously found decreases of connectivity (Liao et al., 2010;Su et al., 2015). A possible explanation is that our results are not directly comparable to previous work because we took absolute connectivity values to treat account for positive and negative correlations (as these both represent an inter-regional coupling) in a homogeneous manner.
Previously we found increased functional connectivity outside the epileptogenic zone correlated to altered memory performance (Bettus et al., 2009) suggesting a global reorganization of functional connectivity in TLE. This increased synchronization of BOLD signals could be an expression of more homogenous brain functioning with a decrease in functional repertoire (Maccotta et al., 2013) caused by strong self-sustaining oscillatory states induced by an increased noise level (Deco et al., 2009) due to epileptic activity. In line with this, we previously found a homogenization of hubs by flattening the differences between hubs and non-hubs (Ridley et al., 2015). This might be another consequence of the here observed globally decreased functional repertoire. Of note, these increased correlations did not involve functional connectivity between regions commonly implicated in the individual epileptogenic network in rTLE patients (Bettus et al., 2009) but relates to edges terminating in core epileptic nodes such as the right hippocampus (See right temporal lobe, Fig. 2).

Alteration of the structural connectome in the rTLE group
Among the different parameters we used to characterize structural connectivity (streamline count, FCA, path length, search information, path transitivity), only the search information metric (averaged over the whole brain) was shown to be different between rTLE patients and controls. Previous whole brain structural studies found decreased path length in TLE patients (see (Haneef and Chiang, 2014) for review), by either exploring a range of arbitrary sparsity thresholds of the connectivity matrix (Bernhardt et al., 2011) or by restricting the sparsity artificially through building consensus networks of largely heterogeneous groups of temporal and frontal lobe epilepsy (Vaessen et al., 2012). Our approach was to retain only connections varying less than three standard deviations in the two tractography runs. This resulted in a sparsity of 0.45. In line with (Bernhardt et al., 2011) path length does not differ at this sparsity. Nevertheless our search information parameter was Table 2 P-values for permutation of group labels of the averaged connectome on the regression slopes of a linear fit (Controls vs. rTLE, n = 10.000, *p b 0.05); NBS: altered NBS subnetwork of functional connectivity rTLE N controls, Rich Club, Feeder and Local: rich club, feeder and local connections for structural connectome sparsitiy 30%, and nodal degree N 250), FCA: analytic functional connectivity, ED: Euclidian distance, PL: Path length, PT: Path transitivity, SI: search information, PP combined path predictors: path length, path transitivity, search information. shown to be sensitive to rTLE with no need to arbitrarily adjust the matrix sparsities.

Whole-brain NBS
Conversely NBS-corrected pairwise connections on a global level using streamline count were not different, concordant with findings at the whole brain scale in other studies in rTLE (Bonilha et al., 2012). This suggests search information as a plausible new parameter to measure global scale reorganization of the structural network in rTLE. These results are in line with recent findings stressing a global network involvement with loss in white matter integrity and rise in connections at the local scale (Bonilha et al., 2012) as well as at the whole brain scale (DeSalvo et al., 2014). Investigations of histology have provided evidence of abnormal myelinated cortical fiber tracts in the temporal lobe (Thom et al., 2000), altered white matter neuronal density (Thom et al., 2001) and widespread white matter alterations (Blanc et al., 2011) that might trigger the observed modified structural and functional connectivity in TLE. Structural brain lesions in many brain disorders such as TLE are more likely to occur in hub regions of the human brain (Crossley et al., 2014). By extracting the rich club coefficient to examine nodal degree changes in central interconnected hubs, we observed a rich club organization in controls and rTLE patients consistent with the high resolution rich club described in van den Heuvel and Sporns (2011). Comparison of rich club coefficients (van den Heuvel et al., 2013) between rTLE patients relative to controls did not show any differences.
Despite this strong overlap an increased nodal degree threshold reveals a differing disintegration of the rich club network showing different core components between rTLE patients and controls. When increasing the nodal degree thresholds, resilience of the rich club network was higher for rTLE patients relative to controls. Interestingly, among nodes of networks that are divergent between rTLE and controls at higher degree thresholds, the hubs constituting the rich club included bilateral thalami, which have been shown to be i) highly involved during seizures in TLE (Guye et al., 2006), ii) prone to atrophy (Bernhardt et al., 2012;Keller et al., 2014) and iii) prone to topological reorganization (Bonilha et al., 2012). In addition, thalamo-parietal connectivity has previously been found to play a central role in seizure propagation in TLE (Arthuis et al., 2009). Bonilha et al. (2012) previously reported increased nodal degree of insula and bilateral thalami. Structural changes in the putamen have been previously found in DTI studies (Keller et al., 2013) and furthermore the parietal lobe, thalamus, posterior cingulate gyrus/precuneus and insula were shown to be implicated during TLE seizures (Bartolomei et al., 2001;Gotman et al., 2005;Kahane and Bartolomei, 2010).
It is important to note that the observed increased streamline count does not automatically imply an increase in anatomical fibers when comparing patients to controls. The tractography algorithm employed here which uses a fixed number of streamlines cannot discern whether the more highly connected regions of the rich club in patients reflect genuinely augmented degree in those specific regions in comparison to controls, or are merely the regions with preserved degree despite a global fiber loss (Calamante et al., 2015). This point can be resolved by taking into account structure function alterations of these rich club regions which are not influenced by the limitation of constant streamlines (see Section 4.5). Given the role these regions are known to play in epilepsy, the notion of a preserved rich club despite widespread fiber loss could explain disease symptoms by providing a propagation channel of epileptic activity through this preserved network.

Structure-function relationship in the whole brain
In line with previous work in healthy controls, we observed significant correlation between functional connectivity and structural parameters such as Euclidian distance and streamline count (Honey et al., 2009). Furthermore, we demonstrated the high predictive value of the search information metric in relation to functional connectivity similar to previous reports in healthy brains (Goñi et al., 2014). The multivariate analysis showed the significant added value of Euclidian distance in the prediction of functional connectivity when added to the regression model of search information, streamline count, path length and path transitivity.
By relating structural to functional alterations, this graph analytical approach could help deciphering discrepancies in results between structural and functional low-level connectional disruptions in TLE (Bernhardt et al., 2013). In contrast to results obtained in generalized epilepsy showing a structure-function decorrelation (Zhang et al., 2011), we observed in this group of rTLE patients an increased structure-function correlation. This might be due to the different phenotypes of disease. Interestingly, idiopathic generalized epilepsy is less characterized by cognitive dysfunction and particularly social cognitive dysfunction than TLE Epilepsy (Realmuto et al., 2015), which suggests different underlying reorganizations of structure and function between the two phenotypes.
From a pathophysiological point of view, differences between rTLE patients and controls in the prediction strength of the Euclidian distance metric with regards to functional connectivity might be related to altered cumulative structural damage along the streamline (e.g. gliosis or demyelination (Thom et al., 2000)) that did not change the fiber orientation. This observed altered distance-function correlation might reflect changed time-lag structure in functional resting-state connectivity in rTLE (Mitra et al., 2015). In addition, while we observed increased search information in patients and its ability to better predict functional connectivity, the mean path transitivity does not differ between the two groups. This suggests a reorganization of networks in rTLE patients via additional, spurious connections acting as branches off the shortest path. Increased search information might be also linked to previously found decreased global fractional anisotropy (Bernhardt et al., 2013) which would cause the tractography algorithm to find more secondary streamlines branching off the main tract. Especially for directly connected nodes this can add useful information for prediction of functional connectivity in opposition to weighted path length which reflects only the inverse streamline count for directly connected nodes. Like the earlier discussed increased functional connectivity, augmented structure-function correlation suggests more homogeneous functional connectivity in rTLE approaching the underlying structural connectivity, which might represent a loss in the functional repertoire of the rTLE brain.

Structure-function relationship in the hypercorrelated functional network identified by NBS
Connections of the altered functional NBS subnetwork between controls and rTLE patients had higher structure-function correlations than observed for whole brain analysis (both controls and patients). This suggests that the functional changes of epilepsy predominantly manifest in brain subnetworks with pre-existing strong structure-function coupling in controls. This relationship was observed for both Euclidian distance and tractography-(derived) measures like streamline count and search information (Fig. 5). Notably the effect for streamline count and search information persists even when we control for Euclidian distance. Interestingly in contrast to whole brain analysis we found that after Fig. 5. Prediction of functional connectivity for controls and rTLE patients using a) Euclidian distance, b) search information and c) streamline count based on group averaged connectomes for connections in the altered NBS subnetwork of functional connectivity; Euclidian distance (middle row) and search information (lower row) vs. functional connectivity in three subnetworks at a degree threshold k N 250 splitting the brain in d) Rich club connections, e) Feeder connections and f) Local connections. For detailed correlation values and slope differences see Table 1 and 2. regressing out the distance, streamline count predicts functional connectivity better than search information (while both metrics remain significantly correlated to functional connectivity).

Structure function relationship in rich club and peripheral regions
Function-structure correlations observed in the rich club, feeder and peripheral connections (Fig. 5d-f) revealed correlation changes in peripheral connections only. This might be in relation to observations in macaque data attributing slower dynamics to the rich club compared to the decentralized periphery (Gollo et al., 2015), suggesting that variation due to dynamic pathological epileptic activity will also be more likely to be observed in the peripheral network given the relative scarcity of additional stabilizing connections as found in the rich club. Taking into account the rich club alterations (see Section 4.2) -which could be linked either to a preserved core in the presence of global fiber loss or a hyperconnected core in an unchanged connectome, this result is more in favor of the latter alternative. This is supported by previous studies which found localized but widespread structural degeneration in rTLE (Bernhardt et al., 2013;Besson et al., 2014). This is also in line with previous findings of rich club preservation in the elderly brain (Perry et al., 2015) and in Alzheimers disease (Daianu et al., 2015). Taking into account the general degenerative progression of epilepsy (Avanzini et al., 2013;Goddard et al., 1969) this preservation of the central core could possibly be a necessary mechanism to prevent a collapse of global brain functioning at the price of constant added noise by the epileptogenic zone. Future work should validate if this preserved structural functional network indeed eases spreading of seizures, for example through hub nodes such as thalamus, insula or precuneus. This might provide new insight on the characterization of epilepsy as a functional disease of uncontrolled information spreading in a closed network with secondary degenerative effect on the structural connectome. In the light of these results it would be especially interesting to better understand the possibility of brain adaptation to decrease potential seizure propagation through this central core limited by preserving a general functioning of the brain.

Structure-function relationship and non-stationary brain dynamics
It should be noted that in the current paper we have considered structure function correlation as averaged over a relatively large acquisition period (20 min). Previous literature (Honey et al., 2007) would lead us to expect that under this paradigm, variation through time as the brain oscillates through various states will be averaged out and thus the functional connectome approaches the underlying structural connectome. The importance of the current work is that structure-function correlation is augmented in epileptic patients. However, recent studies pointed out that the time averaged correlation of the fMRI is an oversimplification of resting state fMRI data (Allen et al., 2014;Hansen et al., 2015;Tagliazucchi and Laufs, 2015;Zalesky et al., 2014). In particular, time averaged data may suppress non-stationarities in the dynamic behavior of networks. In this context, the approach taken by Deco et al. (2013) and incorporated into the FCA methodology used here is insufficient to capture the vagaries of structure-function correlation as the brain transitions between states over time. A range of possible brain organization states has been identified with the brain's dynamic repertoire (Deco et al., 2011). The tighter coupling observed here in rTLE patients between analytical structural models and functional connectivity may reflect a reduced dynamic repertoire of possible brain states which could be observed via Functional Connectivity Dynamics (FCD) (Hansen et al., 2015). Future time-resolved dynamic analyses are needed to examine the extent of structure-function augmentation in brain states resolved over smaller analysis windows.

Limitations
We acknowledge the limits of the sample size, chosen in order to have a relatively homogeneous group in terms of epileptogenic network location, including only right lateralized epileptic patients. Nevertheless we have shown that our data is in line with both controls (Goñi et al., 2014) and rTLE findings (Bernhardt et al., 2013). Furthermore, the use of connectomes averaged over the whole groupas per the approach of Goñi et al. (2014) is not optimal especially in terms of the potential clinical applicability of the results. Future work should focus on deriving reliable brain network descriptors for single subject analysis of structure-function correlations. Lateralization of TLE is an important factor in both structural (Besson et al., 2014) and functional connectivity (Ridley et al., 2015) so our findings might not be generalizable to TLE.
This study cannot answer the question whether search information is a good predictor to classify rTLE structure-function relationship in general or if it is related to hub alterations generally observed in pathologic brain states (Crossley et al., 2014) and it might be interesting to apply this approach to different pathologies.

Conclusion
In this study we have demonstrated that an edge-based graph analytical approach can elucidate the impact of pathological alterations of structural connectivity on function by linking widespread hypercorrelation of networks in rTLE to structural changes. We maximized the structure-function correlation by using search information, a metric which represents the random reorganization of structural connections, revealing a significantly higher correlation in rTLE patients. Additionally we highlighted the importance of considering Euclidian distancewhich is both highly correlated with structural and functional connectivity data -as a major confounding variable of functional-structural correlations in rTLE, possibly explaining structural alterations in TLE which cannot be captured by dMRI based tractography. Alterations seem to be more likely to appear in the local periphery of the brain while core connections show preserved structural and functional stability. Increased functional connectivity and increased structure-function correlation suggest a fallback to structure and loss of functional variability in rTLE patients.