Mutual Information Better Quantifies Brain Network Architecture in Children with Epilepsy

Purpose Metrics of the brain network architecture derived from resting-state fMRI have been shown to provide physiologically meaningful markers of IQ in children with epilepsy. However, traditional measures of functional connectivity (FC), specifically the Pearson correlation, assume a dominant linear relationship between BOLD time courses; this assumption may not be valid. Mutual information is an alternative measure of FC which has shown promise in the study of complex networks due to its ability to flexibly capture association of diverse forms. We aimed to compare network metrics derived from mutual information-defined FC to those derived from traditional correlation in terms of their capacity to predict patient-level IQ. Materials and Methods Patients were retrospectively identified with the following: (1) focal epilepsy; (2) resting-state fMRI; and (3) full-scale IQ by a neuropsychologist. Brain network nodes were defined by anatomic parcellation. Parcellation was performed at the size threshold of 350 mm2, resulting in networks containing 780 nodes. Whole-brain, weighted graphs were then constructed according to the pairwise connectivity between nodes. In the traditional condition, edges (connections) between each pair of nodes were defined as the absolute value of the Pearson correlation coefficient between their BOLD time courses. In the mutual information condition, edges were defined as the mutual information between time courses. The following metrics were then calculated for each weighted graph: clustering coefficient, modularity, characteristic path length, and global efficiency. A machine learning algorithm was used to predict the IQ of each individual based on their network metrics. Prediction accuracy was assessed as the fractional variation explained for each condition. Results Twenty-four patients met the inclusion criteria (age: 8–18 years). All brain networks demonstrated expected small-world properties. Network metrics derived from mutual information-defined FC significantly outperformed the use of the Pearson correlation. Specifically, fractional variation explained was 49% (95% CI: 46%, 51%) for the mutual information method; the Pearson correlation demonstrated a variation of 17% (95% CI: 13%, 19%). Conclusion Mutual information-defined functional connectivity captures physiologically relevant features of the brain network better than correlation. Clinical Relevance Optimizing the capacity to predict cognitive phenotypes at the patient level is a necessary step toward the clinical utility of network-based biomarkers.


Introduction
Computational methods now have the capacity to model the cerebral network at the whole-brain scale [1]. In this context, the brain is represented as a collection of anatomical elements, or nodes; connections between pairs of nodes, referred to as edges, are then measured noninvasively. Once network can be quantified. Elements of the brain that interact to support a given function continue to exhibit similar BOLD fluctuations at rest [4]. Hence, the strength of a connection between each pair of nodes can be inferred from the similarity of their BOLD signal time courses. As this sequence is task free, it offers the potential to measure the functional status of children who are too young or too impaired to cooperate with traditional functional imaging.
ese attributes point to the potential for resting-state approaches to deliver new clinical tools, especially in disorders of the brain that emerge from reorganization of the cerebral network such as epilepsy [3]. Recent work has demonstrated the potential of network metrics derived from rs-fMRI to provide clinically meaningful markers of cognitive function in adults [5,6], in healthy children [7], and in children with focal epilepsy [5][6][7][8]. Despite this promise, exactly how neuronal interaction across the cerebrum is reflected by these spontaneous fluctuations in the BOLD signal-and therefore how to best measure similarity in BOLD time courses-is yet to be determined. e most commonly used measure of functional connectivity in resting-state studies is the Pearson correlation coefficient, defined as the linear covariance of two variables divided by the product of their standard deviations. e Pearson correlation coefficient is simple to calculate and facilitates communication among researchers of diverse disciplines. However, a critical assumption inherent to the use of correlation in the resting state-that the physiologically relevant information about interactions between two discrete brain regions is reflected by a linear relationship between the values of their respective BOLD signals at the same time-may not be valid. In particular, recent studies have shown that nonlinearities inherent to resting-state acquisitions, predominantly hemodynamic in origin, affect both the timing and the amplitude of the measured BOLD signal [9]. As a result, relationships between time series are influenced by the profile of temporal interactions rather than by zero-lag interactions alone [10,11]. Furthermore, recent work has suggested that nonlinear relationships may play an even more prominent role in the connectivity of pathologic tissues [12]. Beyond issues of linearity, there is a great deal of uncertainty in terms of how the true neuronal interactions we hope to measure are represented by fluctuations in blood flow (BOLD); this challenge highlights the importance of generality [13]. Mutual information is an alternative measure of similarity that quantifies in a very general way how much one random variable tells us about another. It is a dimensionless quantity and can be thought of as the reduction in uncertainty about one variable given knowledge of another. Mutual information has been shown to outperform other methods for characterizing association between time series in simulated networks, in part for its generality and equitability [10,13]. It has also been shown to provide a repeatable estimate of network connectivity in the brains of normal subjects [10]. Little data exist, however, regarding the importance of nonlinear association within resting-state networks with regard to the emergence of cognitive dysfunction. We therefore sought to compare brain networks constructed from mutual information to those based on correlation in terms of their capacity to support patient-level inferences on the relationship between brain network architecture and brain function.

Patients.
e HIPAA-compliant study was approved by a local institutional review board. Informed consent was waived. Patient medical records were retrospectively reviewed to identify patients with the following inclusion criteria: (1) pediatric age group (21 years of age or younger); (2) a clinical diagnosis of focal epilepsy; (3) an available 3 tesla MR imaging of the brain, including an rs-fMRI sequence; and (4) full-scale intelligence quotient (IQ) using an age-appropriate Wechsler Intelligence Scale measured by a pediatric neuropsychologist within 3 months of the MR imaging. e above-defined cohort was refined by applying the following exclusion criteria: (1) any brain operations performed prior to the MR imaging or (2) having poor image quality due to either motion or other artifacts.
Imaging was performed from January 2013 to June 2015. irty-four patients met the inclusion criteria. Ten were excluded on the basis of prior brain surgery. Twentyfour patients (age range: 8-18 years; median: 13.4; 12 (46%) females) made up the final cohort. Of this cohort, 5 patients had structurally normal brains and 19 patients had demonstrable structural abnormalities at MRI, including focal cortical dysplasia (n � 8), mesial temporal sclerosis (n � 5), low-grade tumor (n � 4), and a single epileptogenic tuber in the setting of tuberous sclerosis (n � 2). An ageappropriate version of the Wechsler intelligence test was successfully administered in all patients; full-scale intelligence quotient in the cohort ranged from 52 to 129 (median: 91).

Image Processing and Analysis.
e processing pipeline was implemented using MATLAB scripts (version 7.13; MathWorks, Inc.) in which adapter functions were embedded to execute FreeSurfer reconstruction (version 5.3.0; http://surfer.nmr.mgh.harvard.edu) and several FMRIB software library (FSL) suite tools [14]. Details regarding this pathway have been previously described [3,8]. A brief summary is provided here.

Network Node Definition.
e reference space was created from images of one patient in our database, who had no visible abnormality and with optimal registration to the MNI space [15]. Structural imaging data for each patient were aligned to a standard reference template (MNI152) using the FSL's nonlinear registration algorithm [14,16]. Nodes in the network were defined on the template according to parcellation of whole-brain gray matter. First, FreeSurfer reconstruction of cerebral cortical surfaces was performed on the T1 structural image.
is processing stream includes motion correction, skull stripping, intensity normalization, segmentation of white matter and gray matter structures, parcellation of the gray matter and white matter boundary, and surface deformation following intensity gradients which optimally place the gray matter/white matter and gray matter/cerebrospinal fluid borders [17,18]. e pial and gray white surfaces were visually inspected using the Freeview software for accurate placement.
Next, a self-developed MATLAB program was applied to the FreeSurfer output to further subdivide the 75 standard gray matter parcels according to their surface area. During this process, each parcel was iteratively divided into two new parcels of equal size until the surface area of each parcel (as defined on the FreeSurfer gray-white surface mesh) was less than a predetermined threshold value. Networks were constructed with a size threshold of 350 mm 2 . e final parcellation contained 780 nodes ( Figure 1). Each surface parcel was then converted into a volume mask of gray matter at that region to form a node on the network. All nodes defined in the reference space were transformed into each individual patient's space by applying the nonlinear transformation matrix (12 degrees of freedom) obtained during registration.

FMRI Data Preprocessing.
e first 5 volumes in each resting-state functional datum were removed to allow magnetization to reach equilibrium. Standard preprocessing and independent component analysis (ICA) of the functional datasets were performed using FSL MELODIC [14], consisting of motion correction, interleaved slice timing correction, brain extraction, spatial smoothing with a Gaussian kernel full width at half maximum of 5 mm, and high-pass temporal filtering equivalent to 100 seconds (0.01 Hz). Noise related to motion and other physiologic nuisance was addressed according to an independent component analysis technique [19]. Nonsignal components were removed manually by an expert operator with 6 years of experience using independent component analysis in this patient population. Although the optimal strategies for noise removal in fMRI are debatable [20,21], an independent component analysis was selected because it has been shown to minimize the impact of motion on network metrics while, at the same time, decreasing the loss of temporal degree of freedom and preserving the signal of interest across a variety of resting-state datasets [21]. Affine boundary-based registration as implemented in FSL FLIRT was then used to align the preprocessed functional image volumes for each patient to that individual's structural T1 dataset using linear registration. e inverse transformation matrix was calculated in this step and subsequently used to transform all masks from structural to functional space. Mean BOLD signal time series were then computed for each node.

Network Edge Definition.
e strength of an edge (connection) between 2 nodes was defined in two ways: (1) the absolute value of the Pearson correlation coefficient between their BOLD time series and (2) the mutual information calculated based on the following method.
For two discrete random variables X and Y, their mutual information takes the following form: where S x and S y are possible values of X and Y, p(x, y) is the probability that the pair (X, Y) takes values x in S x and y in S y , and p(x) and p(y) are two marginal probabilities of X and Y. For a pair of time series taking small number of values, the probability functions p(x), p(y), and p(x, y) are estimated by frequency counts of the values x and y appeared in the time series. Applying this formula to continuous time series seen in most studies requires a grid to discretize the continuous space into small boxes. e probability functions p(x), p(y), and p(x, y) then are the frequency counts of values within the boundaries of a box centered at x and y. e boundaries and the resolution of the grid affect the value of mutual information. To avoid the ambiguity of choosing a grid, we take the largest mutual information of all possible grids of a predetermined resolution.
is maximization applied in the mutual information calculation shares the same principal as the maximal information coefficient (MIC), where the maximization is taken over all grids up to a maximal resolution [13]. It is computationally impossible to search all grids for over 300 thousand pair time series of a single patient in our study however. erefore, the resolution of the grids had to be predetermined. By testing data of several randomly selected patients, the 3-by-3 grids were chosen in the study as they provided similar mutual information compared to finer grids but required much shorter computation time. e boundaries of 3 bins on x-or y-axis of 3-by-3 grids were determined by 4 values. e two ends are min and max of a time series. e middle two values were determined by mean ± a multiple of standard deviation of the time series. We chose 5 values for the multiple, which yielded 25 possible choices for a pair of the middle two values, 25 choices for 3 bins on x-or y-axis, and 125 choices for 3-by-3 grids for a pair of time series.

Graph Construction and Network Metric Calculation.
Two weighted, undirected connection matrices of each patient were constructed, named as "Pearson and mutual Computational and Mathematical Methods in Medicine 3 information graphs," consisting, respectively, of the pairwise Pearson correlations and the mutual information between BOLD time series over all network nodes. e following topologic properties were calculated by using MATLAB scripts provided in the Brain Connectivity Toolbox (https://sites.google.com/site/bctnet/): clustering coefficient, modularity, characteristic path length, and global efficiency. A short description of each metric is provided in Table 1.
Clustering coefficient and modularity are metrics that measure the brain's tendency to segregate into relatively independent, local neighborhoods. In other words, these measures reflect the ability of the brain to process specialized functions within highly interconnected functional subnetworks. Characteristic path length and global efficiency measure the global integration of the brain. A short characteristic path length or a high global efficiency indicates that information can be integrated easily across the brain.

Statistical Analyses.
Statistical testing was performed using SAS version 9.3 and R language version 3.4.0 (R Foundation for Statistical Computing, Vienna, Austria). e primary endpoint was the predictive value of the 4 output metrics of the brain network architecture (derived from either Pearson or mutual information graphs) with respect to individual intelligence. is multivariate analysis was accomplished using a random forest approach, which has been previously described in detail in [22]. In short, this ensemble learning method operates by constructing a multitude of decision trees during training and outputting the mean of predictions from individual trees. It is based on bootstrap aggregating, or bagging, in which numerous models are fitted during individual bootstrap sampling and then combined by averaging. During training, approximately one-third of the cohort is omitted at random from the training set-this omitted portion of the dataset is considered "out of bag." e IQ of each individual held out of bag is then predicted based on the "learned" model. Prediction accuracy for the out-of-bag cohort was quantified in two ways: (1) mean absolute error and (2) fractional variation explained [23]. To be specific, the random forest algorithm was given access to only the four network metrics and no other patient information during this analysis. e absolute errors of predictions from Pearson graphs were compared to those from mutual information graphs using the Wilcoxon signed rank test. All random forest models were run 500 times to obtain the 95% confidence interval (CI) for fractional variation explained. e random forest algorithm was also used to measure the independent contribution of individual network metrics to the prediction of IQ. In other words, it measures the association of each variable after accounting for all other variables. is contribution is estimated for each variable  by measuring the error for IQ prediction in the out-of-bag cohort compared with the error that results when that particular variable is negated during bagging. Connections from a mutual information graph were compared to those from the corresponding Pearson graph through a scatter plot of each patient. Differences in network metrics computed on Pearson versus mutual information graphs were assessed using the Wilcoxon signed rank test. Relationships between Pearson-derived network metrics and mutual information-derived network metrics were also quantified using the Pearson correlation coefficient. e Pearson graph was chosen to measure this association since monotonic relationships were expected. Finally, the univariate association of each metric with IQ was measured in a univariate analysis by the Spearman correlation coefficient.

Association between Pearson and Mutual Information Graphs.
A representative example of a scatter plot of Pearson versus mutual information connections of one patient is provided in Figure 2. e reference line on the graph is −1/2 log (1−r 2 ), which is the relationship between mutual information and Pearson correlation if the joint distribution is Gaussian [24]. Deviation of our data from the reference line, reflecting nonlinear relationships between resting-state time series, was observed in all patients (data not shown).
Network metrics derived from Pearson graphs versus those from the mutual information graphs for each patient are presented in Figure 3. Although association between the Pearson and mutual information graph metrics was generally high-correlation coefficients ranging from 0.84 to 0.88, differences were apparent (Table 2). On average, clustering coefficient and global efficiency became smaller when computed on the mutual information graph. By contrast, path length tended toward higher values under the mutual information graph. Notably, modularity was not statistically different between the graphs (Table 2). Table 3. For most metrics, the association with patient IQ was greater when computed on the mutual information graph. Using a multivariate approach, mutual information graph metrics made the dominant contribution to subject IQ prediction by the random forest model (Figure 4).

Network Architecture and Intelligence. Univariate correlation of the mutual information and Pearson-derived graph metrics with subject IQ is presented in
Accuracy of the machine learning algorithm's prediction of IQ based on network metrics is presented in Table 4. Metrics derived from mutual information graphs demonstrated a significantly higher predictive value compared to that of the Pearson graph. e relationship between the magnitudes of prediction error for the two methods is demonstrated graphically in Figure 5.

Discussion
We evaluated two measures of association-the Pearson correlation and mutual information-that are commonly used to infer connectedness in brain networks constructed from resting-state functional MRI. We specifically assessed the impact of these measures on output metrics of the global brain architecture in terms of their capacity to support the prediction of global intelligence in children with focal epilepsy. We report that measuring brain network edges using mutual information significantly outperformed the use of the Pearson correlation in this setting.
Higher-order functions of the human brain are not accomplished by individual functional centers compartmentalized to a particular region of cortex. Rather, they emerge from parallel processing within subspecialized, but distributed, functional systems. e ability to decode these neuronal interactions, particularly as they relate to the emergence of brain function, has become a major focus in current neuroimaging research. Resting-state functional MRI is one modality that has been used extensively as a surrogate for connectedness in the human brain. A significant body of work now exists in support of its capacity to probe physiologically meaningful features of the human brain in a diversity of settings [4]. For example, studies have demonstrated an abnormal network architecture in a variety of disease states, including those with prominent cognitive dysfunction [25][26][27][28][29][30]; network reorganization has also been observed in adults [31][32][33] and children [29,34] with focal epilepsies. Beyond group-level comparisons, a relationship of brain network features quantified by graph theory with intelligence has been demonstrated in many populations, including healthy adults [5,6,35], healthy children [7], normal aging [36], Alzheimer's disease [37,38], autism [39], and epilepsy [8]. Given this capacity of resting-state networks to capture interindividual phenotypic variance in brain function, there is great interest in the development of subject-level markers that could be used to guide patient care [3,8]. Despite this promise, exactly how neuronal interaction across the cerebrum is reflected by spontaneous fluctuations in the BOLD signal-and therefore how to best measure similarity in BOLD time courses-is yet to be determined. We observed that metrics of the network architecture computed on mutual information graphs outperformed network metrics based on the Pearson correlation in terms of the ability to predict subject-level intelligence in a cohort of children with epilepsy.
We observed that graph-based metrics from Pearson and mutual information graphs were relatively similar and demonstrated high linear correlation. Although this finding is consistent with prior work demonstrating only a small contribution of nonlinear associations to the rs-fMRI time   series [24], we found that this "small amount" made a significant difference in terms of subject-level prediction in this population. Interestingly, we also observed larger non-Gaussian dependencies among the time series than what has been reported in a healthy adult population [24]. is idea is consistent with the work by Rummel et al. who used a uniform surrogate-based approach to study interrelations that significantly exceed linear correlation in EEG data of epilepsy patients [12]. ey observed that nonlinearity occurred predominantly for epileptogenic tissue as well as during epileptic seizures [12]. Our results align with these studies, suggesting that the dynamics of the abnormal brain may be more complex than those of normal brains and that nonlinear associations may be more prevalent. erefore, a general measure of brain interactions may be more important when analyzing a disease population. Our results are in line with the previous work that has used mutual information to quantify network edges. Reshef et al. calculated the maximal information coefficient (MIC), based on the same concept as our calculation of mutual information; MIC was shown to be superior to linear association measures in terms of discovering important relationships [13]. It allows one to capture a wide range of interesting associations, not limited to specific function types, or even to all functional relationships. is generality is very crucial as many important relationships are not well modeled by a function. It was also shown that MIC was equitable in the sense of being able to retain the discovery of various types of associations even with increased noise in the simulation data [13]. ese attributes may explain the superior prediction of subject intelligence using network metrics observed in our study. Along similar lines, a study in    a cohort of patients with schizophrenia [28] demonstrated that nonlinear functional connectivity provided useful discriminative power toward making the diagnosis in each patient.
is study has several limitations. First, our cohort was a selected population of pediatric patients with focal epilepsy. e results may not be generalizable to other patient populations, or to normal subjects. Second, our sample size was small, which did not allow a study stratified by disease severity or a study on characteristics of patients who benefit more from using mutual information. Nevertheless, our goal was to show the general advantages of a nonlinear method used to quantify functional connections. Finally, the extensive computation time required to generate whole-brain networks under a range of nonlinear methods precludes comparison of an exhaustive list of available methods. Nevertheless, mutual information has been proven to be an effective measure in various disciplines for its generality and equitability.

Conclusion
Brain networks constructed using edges defined by mutual information significantly outperformed the use of the Pearson correlation for predicting global intelligence in a pediatric cohort with focal epilepsy. Network methodologies specifically optimized to make predictions about individuals will be critical to the development and implementation of clinical tools based on resting-state constructs.
Data Availability e patient characteristics and image data used to support the findings of this study are restricted by the HIPAA in order to protect patient privacy. Data are available from Wei Zhang (wxzhang1@texaschildrens.org) for researchers who meet the criteria for access to confidential data.