Multiplex Networks to Characterize Seizure Development in Traumatic Brain Injury Patients

Traumatic brain injury (TBI) may cause secondary debilitating problems, such as post-traumatic epilepsy (PTE), which occurs with unprovoked recurrent seizures, months or even years after TBI. Currently, the Epilepsy Bioinformatics Study for Antiepileptogenic Therapy (EpiBioS4Rx) has been enrolling moderate-severe TBI patients with the goal to identify biomarkers of epileptogenesis that may help to prevent seizure occurrence and better understand the mechanism underlying PTE. In this work, we used a novel complex network approach based on segmenting T1-weighted Magnetic Resonance Imaging (MRI) scans in patches of the same dimension (network nodes) and measured pairwise patch similarities using Pearson's correlation (network connections). This network model allowed us to obtain a series of single and multiplex network metrics to comprehensively analyze the different interactions between brain components and capture structural MRI alterations related to seizure development. We used these complex network features to train a Random Forest (RF) classifier and predict, with an accuracy of 70 and a 95% confidence interval of [67, 73%], which subjects from EpiBioS4Rx have had at least one seizure after a TBI. This complex network approach also allowed the identification of the most informative scales and brain areas for the discrimination between the two clinical groups: seizure-free and seizure-affected subjects, demonstrating to be a promising pilot study which, in the future, may serve to identify and validate biomarkers of PTE.


INTRODUCTION
Traumatic brain injury (TBI) is the third most common cause of death and debilitating secondary problems in adults and children worldwide. One common consequence of TBI that causes significant disability amongst patient populations is post-traumatic epilepsy (PTE) (Humphreys et al., 2013). This condition develops in up to 50% of patients with TBI. Post-traumatic epilepsy (PTE) is diagnosed if two or more unprovoked seizures occur at least 1 week after a TBI (Diaz-Arrastia et al., 2009). Recent investigations suggest that injury severity and especially epileptic activity are high risk factors of PTE, although the mechanisms by which trauma to the brain tissue leads to recurrent seizures is not known. Therefore, studying if specific structural Magnetic Resonance Imaging (sMRI) changes can be related to seizures after a TBI is of fundamental importance to carry out the first steps toward the discovery of early biomarkers of PTE (Kim et al., 2018). PTE is not a homogeneous condition and can appear weeks or several years after a TBI. As a consequence, the precise percentage of TBI patients who develop PTE is not known (Verellen and Cavazos, 2010). Currently, growing attention has been devoted to investigate PTE. In this regard, the Epilepsy Bioinformatics Study for Antiepileptogenic Therapy (EpiBioS4Rx) is an international, multi-center project conceived to identify biomarkers of epileptogenesis after a TBI in order to evaluate treatments that could prevent the development of PTE and design clinical trials of antiepileptogenic therapies on an extensive patient population. With this project, the scientific community can be granted access to a large amount of high quality, multi-modal data, including imaging, electrophysiology, and clinical data from both humans and animals.
Changes in gray matter and white matter related to epilepsy have been widely observed by using structural MRI (Immonen et al., 2018;Shah et al., 2019;Lutkenhoff et al., 2020). Many recent studies have shown that machine learning techniques and multiplex networks applied to completely non-invasive neuroimaging techniques, such as structural MRI, can be useful and efficient to detect pathological alterations in several neurological diseases, such as Alzheimer's disease, Parkinson's disease, and epilepsy (Amoroso et al., 2018c;La Rocca et al., 2018;Bharath et al., 2019). Multiplex networks overcome the limit of the existing complex network standard approaches not to be able to collectively study what happens to the same nodes as their interactions change. In our previous work , we used different machine learning strategies to identify alterations in functional brain connectivity that are related to seizure outcome following TBI. However, the present study is the first which uses the combination of multiplex networks of structural MRIs and machine learning techniques to distinguish patients who have developed at least one seizure after a TBI from those who have not experienced any seizures. This study is of paramount importance, because it offers an opportunity to observe alterations in TBI brain networks that may reflect structural MRI changes related to seizure development.
This paper provides three main results: (i) the implementation of a pipeline which combines complex network and machine learning models for the identification of TBI patients who have developed epilepsy; (ii) the investigation of the most appropriate scale or patch size to study seizure development in TBI patients; (iii) the implementation, on a TBI cohort, of a promising complex network model based on segmenting the brain in patches to obtain comprehensive clinical information on the whole brain. In the future, this pilot study may help clinicians localize the epileptogenic focus more precisely, relate brain lesions to seizure occurrence and understand the relationship between neuronal activity abnormalities and structural damage.  Age and GCS were provided in terms of mean and standard deviation. No statistically significant differences between the two classes were found with respect to age, GCS score, and gender. Statistical evaluations were performed with a Kruskal-Wallis statistic test except for the gender, for which a Chi-square test was used.

Dataset
In this work, we used 53 structural MRI scans of TBI subjects recruited in EpiBioS4Rx according to specific inclusion and exclusion available online 1 . Sixteen of these subjects have experienced at least one seizure within 6 months of a TBI and 37 have not experienced any seizures. As part of their clinical care, 14 subjects required a craniectomy (3 seizure and 11 non-seizure, p > 0.05). Additional clinical and demographic information are reported in Tables 1, 2. 3D T1-weighted volumes were acquired within 32 days (median 8 and interquartile range of [2, 15]) after the TBI using 3T Siemens, Philips, and GE scanners according to a magnetization-prepared rapid acquisition gradient echo (MPRAGE) sequence with the following parameters: 256 mm field of view (FOV); 1 mm slice thickness; 1,500-2,500 ms repetition time (TR); minimum echo time (TE); 1,100-1,500 ms inversion time (TI); 8-15 degree flip angle; 256 phase-encoding steps, number of excitations (NEX) >1 and 256 Hz frequency.

MRI Processing
Protocol compliance and quality control (QC) were undertaken using the Laboratory Of Neuro Imaging (LONI) QC System 2 . Structural MRI scans were processed with the Oxford FMRIB Software Library (FSL) (Jenkinson et al., 2012). Firstly, image skull-stripping was obtained with the optimized brain extraction script for patient brain (optiBET) (Lutkenhoff et al., 2014). For those few cases where the automatic brain extraction was particularly challenging due to the significant brain deformation caused by the trauma, we performed the skull stripping by manually adjusting the brain extraction threshold with FSL Brain Extraction Tool (BET). Then MRI scan intensity differences, yielded by bias field, were normalized. After intensity normalization and brain extraction, a spatial normalization was performed to co-register the different images into a common coordinate space by using an affine transformation. The MNI152 was adopted as the reference template, and registration was performed with the FSL Linear Registration Tool (FLIRT) with a standard parameter configuration. Even though a deformable registration would have given a better overlap among TBI subjects, we purposely used an affine registration for three main reasons: (i) proving the robustness of the method also in case of roughly overlap between the anatomical regions of different TBI subjects; (ii) avoiding misregistration issues due to the particularly challenging process to apply a non linear transformation to a cohort with huge brain deformations; (iii) registering all the subjects to a common reference space keeping as much as possible the individual differences of the subjects and the relative lesions. Besides these initial steps, the analysis pipeline includes two principal sub-pipelines: a complex network pipeline and a machine learning pipeline that are schematized in Figure 1 and are described in detail in the next two sections.

Multiplex Network Pipeline
After image processing, each scan was parceled in homologous non-overlapping parallelepipeds or patches of V voxels (where 1 voxel is 1 cubic millimeter) in order to obtain a 3D grid. These patches represent nodes of a brain network, and the absolute values of Pearson's correlation between patch pairs were considered the links between the nodes. In other words, each network link is obtained by computing the correlation voxelby-voxel between the T1 intensities of two patches. Therefore, for each scan, we obtained a weighted brain network using a patch-based segmentation. To remove links due to the noise, we neglected all the correlations lower than 0.3. This threshold has not been chosen arbitrarily but has been demonstrated in our previous works (Amoroso et al., 2018b,c) as a threshold that maximizes classification performance and is the best trade-off between minimizing noise and maintaining effective network information in this multiplex network methodology. This threshold choice is also confirmed by other works in literature, for example, Mukaka (2012) suggests that correlations lower than 0.3 are negligible in his guide about the appropriate use of correlation coefficients in medical research. To further avoid false positive links in the networks, we also excluded 2 https://qc.loni.usc.edu the patches with a non-brain number of voxels exceeding 10% of their volume. The idea behind this study is that seizure development in TBI patients may be related to injury severity which, as many study demonstrate, results in diffuse cerebral edemas, hemorrhages, contusions, and distortions of brain tissue localized in multiple brain regions both close to and distant from the lesion area. Patch-based approach is aimed to detect this alterations in terms of correlation variation between regions with and without tissue damage over the TBI cohort (Kurland et al., 2012). A patch-based approach is a beneficial trade-off between a voxel-based approach and an ROI-based approach and has already been found to be beneficial in the field of other neurodegenerative diseases (Suk et al., 2014). It has three main advantages: (i) it overcomes the problem of the "curse of dimensionality, " (ii) it does not depend on segmentation accuracy, and (iii) it is robust to misregistration errors (Amoroso et al., 2018a). Therefore, the patch representation can be very useful for TBI patients for whom ROI segmentation and spatial registration are particularly challenging tasks due to the large and irregular brain deformations caused by TBI lesions. To give a sense of the injury severity and the related processing challenges that were faced for this cohort, in Figure 2, for some of the TBI subjects with the most severe imaging findings, axial and coronal planes of brain scans after the processing (brain extraction and registration) are represented alongside the template that the subjects' scans are registered to. The size V of each patch was varied from 1,000 to 8,000 voxels in steps of 1,000 to investigate the most appropriate scales to study seizure development in TBI patients. At different scales, the patch number is not constant but is determined by the patch size. The grid's origin is fixed for all the scales because we start to segment each image from the medial sagittal plane which separates the two brain hemispheres in order to uniformly cover each hemisphere with an equal number of rectangular boxes. For each scale, we built a multiplex network G = {G 1 , G 2 , ..., G i , ..., G M } that is a collection of single subjects' weighted networks G i = (N, E, W) sharing the same nodes N, while the set of links E and weights W change depending on the subject's brain networks connections or layer connections. In other words, in each multiplex network, the same number of nodes (patches that each scan is segmented into) can be connected in different ways depending on the specific correlation coefficient values that characterize the network connections of a certain layer. Then, given N, the number of network nodes, we obtained 8N features for each subject, 4N features of single layer and 4N features of multiplex networks. The single layer features used in this work are strength and inverse participation ratio, given by the Equations (1) and (2), and their conditional means over the nodes with the same degree k, thus having the same connection number. Conditional means, given by the Equations (4) and (5), can be useful to examine whether, on average, the weights of central nodes and less connected nodes are identically distributed.
Frontiers in Neuroscience | www.frontiersin.org FIGURE 1 | Description of the complex network pipeline and machine learning pipeline. First each subject scan is preprocessed, segmented into patches, then for each subject a weighted undirected network was built and some complex network features were computed. Finally, the feature representation (subject × network features), obtained after the removal of null mean and variance features and highly correlated features, was used as input for the classification pipeline. The machine learning pipeline includes 1,000 rounds of cross-validation (CV). In each round the following steps are performed: (i) dataset was stratified; (ii) 80% of the stratified dataset was used as training set and 20% as validation set; (iii) training set was used, through a first nested Random Forest (RF) classifier, to select the most important features for the discrimination of seizure-free and seizure affected subjects: (iv) these selected features were used in turn to train a second RF classifier; (v) the important features and the classification models obtained on the training test were used to classify the subjects of the validation set; (vi) Averaging the classification performances over the 1,000 CV rounds, we obtained the final accuracy sensitivity, specificity, area under the receiver operating characteristics curve (AUC) and confidence interval on the validation test.
FIGURE 2 | Examples of coronal plane (Top) and axial plane (Bottom) of TBI subjects with the most severe imaging findings (Left) and the MNI152 reference template (Right) which the subjects' images are aligned to.
α = 1, .., M indicates the network layer, w i j is the correlation between gray level intensities of the nodes i, j = 1, .., N with M subject number, N is the node number, N k is the number of nodes having degree k, and δ is the Kronecker Delta function. Strength and inverse participation ratio indicate, respectively, the importance of a node and how evenly distributed the connections between nodes are. Specifically, (y α i ) −1 ∈ (1, k α i ) has value k α i if the weights of the links of node i are distributed uniformly and, it has value 1 if the weight of one link is much larger than the other weights (Bianconi, 2018). The multiplex network features were obtained by weighing the previous quantities on the multiplex network degree k multi , given by (5), indicating the number of connections of a node in the multiplex network. ( 5) where a multi ij is 1 if there is a layer with at least one link between nodes i and j and zero otherwise as described in Amoroso et al. (2018b). From now on, we will refer to these multiplex quantities as multi-strength, multi-inverse participation ratio, conditional strength, and conditional multi-inverse participation ratio. Overall, we obtained for each scale V, a M (subject network number) X 8N feature representation to analyze with the machine learning pipeline described in the next section.

Machine Learning Pipeline
For each scale V, multiplex network features were used to train a Random Forest (RF) classifier and obtain reliable classification models to identify, on the validation set, which TBI patients have developed seizures and which have not. This classification process, preceded by the removal of null mean and variance features, and highly correlated features (> 0.95), was carried out within a machine learning pipeline that includes 1,000 rounds of stratified cross-validation (Vabalas et al., 2019). For each round, we randomly picked the same percentage of seizurefree subjects and seizure-affected subjects in order to examine balanced datasets. After the stratification, a training set (80% of the stratified set) and validation set (20% of the stratified set) were defined. Subsequently, we first used a nested RF classifier, on the training set, to select and record features exceeding the third quartile of the importance distribution computed in terms of mean accuracy decrease. Then, we used those features to train a second RF classifier and obtain the classification models. Feature selection and training phases were nested within each crossvalidation round and were blind to the validation set to avoid the "double dipping" problem (Kriegeskorte et al., 2009). Finally, we used the classification models and the important features retrieved during the training phase to classify the two clinical classes. Classification performances for each scale were evaluated in terms of accuracy, specificity, sensitivity, and Area Under the receiver-operating-characteristic Curve (AUC) averaged over all the cross-validation rounds. For the average accuracy, we also reported the 95% confidence interval computed according to the Wilson score interval (Wilson, 1927). We chose to use RF model because it is a robust and easy-to-tune model, it does not overfit thanks to internal bagging and it is particularly appropriate for analyses with high-dimensional feature spaces and small sample sizes (even < 100) (Biau and Scornet, 2016;Floares et al., 2017). Each forest was grown with 500 trees, a number large enough for the out-of-bag error to reach the typical training plateau (Breiman, 1996). Therefore, in the internal bagging, given the training set, 500 bootstraps are formed obtaining 500 new subjects sets used to grow 500 trees. Each tree is grown by randomly choosing a subset of features equal to the square root of the feature number. The learning model built in this way can FIGURE 3 | The bar plot shows the mean and standard deviations of accuracy (gold) and AUC (cyan) over 1,000 rounds of cross-validation as patch volume changes (from 1,000 to 8,000 voxels). The best classification performances were obtained at 1,000, 3,000, and 5,000 voxels.
be then used to compute the out-of-bag error and the accuracy on the data left out of the training set.

Important Feature Assessment
After having found for each round the important features, we evaluated, for the best scales, which feature occurrences had not happened by chance over the 1,000 rounds by using the statistical test of equal or given proportions (Newcombe, 1998). Therefore, the most important features over all the cross-validation rounds were found by considering, after the Bonferroni correction for multiple comparison, a p − value < α * N −1 with α = 0.05. From the most important nodal features, it was possible to find the most important network nodes or patches, and thus the most important anatomical regions. We considered an anatomical region significantly related to the seizure development only if it occupied an important patch with a volume greater than 10% of the patch voxels. To identify the most important anatomical regions, we used Talairach labels projected in MNI 152 space (Lancaster et al., 2000). Other details on the reliability of the feature selection methods used in this work are discussed in the Supplementary Material. Figure 3 shows the mean and standard deviations of accuracy and AUC over all the cross-validation rounds as a function of the patch size. The best classification performances were found at three patch volumes: 1, 000, 3, 000, and 5, 000 voxels.

Classification Performance and Feature Evaluation
Accuracy, specificity, sensitivity, and AUC with the corresponding standard deviations obtained for these three optimal scales are reported in Table 3.
Even though we found the best accuracy for a patch size of 1,000 and 5,000 voxels (with a 95% confidence interval of [67, 73%]), and the best AUC for 3,000 voxels, the classification performances at these three scales are statistically comparable. Another important aspect is to examine which network properties are more important to discriminate the two clinical groups. In this regard, we evaluated the mean percentage of features associated with a certain network metric that are selected as important in a cross validation-round. In Figure 4, the mean percentage of features selected over the cross-validation rounds and relative to each of the eight network metrics, is reported for the three most informative scales. This experiment was TABLE 3 | Accuracy, specificity, sensitivity, and AUC with the relative standard deviations obtained at the scales of 1,000, 3,000, and 5,000 voxels to which the best classification performances were reached.

Patch volume
Accuracy Specificity Sensitivity AUC performed without excluding from the classification the highly correlated features. Even though, for each scale, all metrics extracted contribute to the discrimination of the two clinical groups, we can notice that the nodal metrics have a greater relevance compared with the conditional quantities.

ROI vs. Patch-Based Network Approach
We also compared the patch-based network approach with a standard ROI-based approach to evaluate the efficacy of the proposed complex network methodology to predict seizure development in TBI patients. We used the publicly available brain segmentation package, FreeSurfer (FS) v.6.0 (Fischl, 2012), which automatically performs: brain extraction, intensity normalization, spatial registration, volume labeling, segmentation, and all steps necessary to compute morphological features from each image. This tool allowed us to obtain 182 features, for each MRI scan, including subcortical and cortical gray matter parcellations, white matter parcellations, total gray and white matter volumes, and intracranial volume. These FS features were then used to distinguish TBI subjects who have developed epilepsy from those who have not by adopting the same machine learning pipeline used for the complex FIGURE 4 | Mean percentage of features, selected over the cross-validation rounds, relative to strength (S), inverse participation ratio (Y), multi-strength (multiS), inverse-participation ratio (multiY), and their conditional means(Sc, Yc, multiSc, multiYc) are reported for the scales of 1,000 voxels (yellow barplot), 3,000 voxels (blue barplot), and 5,000 voxels (dark green barplot). network features. In Figure 5, receiver operating characteristics (ROC) curve and the related area under the curve (AUC) are reported for the three best scales of the complex networks and for FS.

Anatomical Regions Related to Seizure Development
For the three scales proved to be more appropriate to identify brain network alterations related to seizure development, we reported, in Figure 6, the brain areas (highlighted in green) corresponding to the most significant complex network features for the classification of seizure-free subjects and subjects with one seizure.
At the scale of 1, 000 voxels, the significant patches (p < 1.563 * 10 −5 after the Bonferroni correction) identify anatomical regions mostly located in the right and left superior temporal gyrus lobe, but there are significant patches also in the left middle temporal gyrus, left inferior frontal and precentral gyrus, and in the right cerebellum within posterior lobe. At the scale of 3, 000, the most important brain area (p < 3.962 * 10 −5 after the Bonferroni correction) for the two group discrimination corresponds to the cingulate gyrus in the left parietal lobe and in the right and left limbic lobe, sub-gyral in left and right frontal and parietal lobe, right and left precuneus, right postcentral gyrus, left inferior parietal lobule, angular gyrus, medial frontal gyrus, and superior occipital gyrus. Finally, the important areas (p < 6.361 * 10 −5 after the Bonferroni correction) at the scale of 5, 000 voxels were the left and right cerebellum in the posterior lobe, right parahippocampal gyrus, right subcallosal gyrus, left inferior middle, and superior frontal gyrus, sub-gyral in the left frontal lobe, cingulate gyrus in the left limbic lobe, right and left extra-nuclear white matter, and left insula.
To make more understandable the relationship between patches and complex network features, in Figure 7, the distribution of the reciprocal of the inverse participation ratio of a patch, located in the left frontal lobe, is reported for both clinical groups. In the same figure, the representation of such a patch in a seizure-affected patient who has significant abnormalities in that area and in a seizure-free patient who does not have visible anomalies in that area is shown. The inverse participation ratio relative to the patch represented in Figure 7 is an example of a network feature which is important for the discrimination of the two clinical groups. Indeed, from the box plot, we can notice that the median of the distribution for the seizure-affected subjects is significantly greater than the median of the distribution for the seizure-free subjects.

DISCUSSION
In this work, an innovative multiplex network approach was used to find informative complex network features that can be used from machine learning systems for the identification of patients who have developed a seizure after a TBI. To the best FIGURE 5 | Classification performances in terms of area under the receiver operating characteristics curve (AUC). Performaces obtained with the complex network features at the the best three scales (red, orange, and green curves) are significantly greater than those obtained with FreeSurfer (FS) features (blue curves). of our knowledge, this is the first study to distinguish seizurefree subjects and seizure-affected subjects with an accuracy of 70% and an AUC of 76% by using T1-weighted MRI data. In Messori et al. (2005), PTE prediction using human MRI is based only on statistical evaluations. Classification performances obtained in this work are comparable and even higher than those found in La  and Garner et al. (2019) that examine functional and structural alterations related to seizure onset. All the network properties are proved to be useful to the classification, however, features relative to conditional metrics were selected less frequently in each cross-validation round. Classification performances were computed as a function of the patch volume, which was varied in an intermediate range (from 1,000 to 8,000 voxels) in order to avoid that the analysis was affected by a low sensitivity to subtle pathological changes in case of too large patches or by the 'curse of dimensionality' and misregistration in case of too small patches. The best classification performances were obtained at three different scales: 1,000, 3,000, and 5,000 voxels, proving that the study of seizure development in TBI patients requires multivariate analyses. Although some of the important anatomical regions, such as cingulate gyrus, sub-gyral, inferior central gyrus, and cerebellar tonsil in the right cerebellum are accordant for the three scales, we can notice that some morphological changes between the two clinical groups can be detected only at specific scales. This suggests that seizure development in TBI patients cannot be studied considering a unique scale, which is reasonable given the heterogeneity of the epileptogenesis process after a FIGURE 6 | Patches corresponding to the most important nodal complex network features for the discrimination of the two clinical groups (seizure-free and seizure-affected) are underlined in green along the axial planes of the MNI 152 template. (A-C) Display the most significant patches relative to the scales of 5,000, 3,000, and 1,000 voxels, respectively.
FIGURE 7 | As an example, (A,B) show the details of a patch in two different TBI subjects. As shown in the green box plots on the right, this patch has an inverse participation ratio that is significantly different (p < 2.2 * 10 −16 ) in the two clinical groups (seizure-free and seizure-affected patients). (B) Shows the patch pinpointing an area where voxels intensity is altered by a lesion and the surrounding edema and (A) shows the patch covering a brain area that does not have evident alterations.
TBI and the fact that TBIs affect the brain in different areas and at different scales. As a consequence, analyzing multiple scales can give a more exhaustive detection of the MRI changes related to seizure development after a TBI that a unique scale is not able to provide. Therefore, our methodology can be very useful to perform a multivariate analysis that take into account multi-scale features. Conventional volumetric analyses are based on manual segmentation of the Region-Of-Interest (ROI) that is time-consuming and affected by personal bias. The modern automated algorithms that allow the determination of volumes, thickness, and shape of anatomical structures often fail because of large lesion size and extensive tissue damage in TBI patients.
In this regard, we demonstrated that our approach (AUC of 76%) is more effectiveness than a ROI-based approach like FS (AUC of 67%). EpiBioS4Rx is an ongoing study that will enroll 300 patients, therefore in upcoming years, with a larger and more representative training sample, we will be able to fully exploit machine learning potentialities and obtain conclusive results about the generalization power of the model in predicting seizure development in TBI patients (Figueroa et al., 2012). Once more subjects will be enrolled and longitudinally examined, it will be also interesting to see if the proposed methodology is able to distinguish among immediate, early, and late seizures in order to take into account also the temporal aspect of the epileptogenic process. Besides, the completely automated complex network approach used in this work can be really beneficial because it allows an unsupervised identification of the important brain areas without being affected by ROI segmentation mistakes and time-consuming procedures. This methodology offers also other two main advantages: (i) it allows the identification of MRI changes which differentiate seizure-free and seizure-affected patients and which cannot be underlined using only CT MRI findings (see Table 1); (ii) it allows the identification of the brain scales at which the pathological changes related to seizure development occur. Indeed, as might be expected, given TBI variability, epileptogenesis mechanism will depend on alterations that happen at different scales (Cloots et al., 2013). It is interesting to notice that at the scale of 1,000 voxels, most of the patches are located at the periphery of the brain. This may be due to brain surface deformations or to subdural and epidural hematomas that are reported among the risk factors to develop epilepsy and are present in many subjects of this cohort as reported in Table 1 (Agrawal et al., 2006). An example of subdural hemotoma is shown in Figure 7B in the left posterior part of the brain. Most of the clinical results are in line with recent studies about seizure development. Norden and Blumenfeld (2002) states the increased likelihood of cerebellar alterations in patients with epilepsy. In Shultz et al. (2013), MRI alterations were found in hippocampus subfields of rodents with epilepsy after a lateral fluid percussion injury. Tubi et al. (2019) showed that subjects with lesions in the temporal lobe are at high risk to develop epilepsy, suggesting that morphological alterations in the temporal lobe may play a strategic role in seizure occurrence. Hippocampus, cingulate gyrus, precentral gyrus, postcentral gyrus, and middle and inferior frontal gyrus were proved to be regions related to the epileptogenesis process also in studies that apply machine learning techniques to fMRI and sMRI to characterize patients with epilepsy (Zhang et al., 2012;Garner et al., 2019;La Rocca et al., 2019). It is worthwhile to notice that the highly correlated (> 0.95) features that were excluded in the classification correspond to complex network metrics related to the same patch and thus, to the same brain area. This ensures that we did not exclude any important region in the clinical validation and suggests that for some patches more complex network metrics are accordant with each other.

CONCLUSION
We have demonstrated that the combined use of complex networks and machine learning techniques can be useful to study seizure development in TBI patients. Multiplex networks were able to provide network features that allow us to distinguish TBI patients who have developed epilepsy from those who have not with an accuracy of 70%. In addition, a patch-based approach used to build the multiplex networks made it possible to identify, in an unsupervised way, the brain areas important for the discrimination of the two clinical groups, even though a perfect solution of optimum features is a challenging and still open matter. EpiBioS4Rx is an ongoing study that will enroll 300 patients, thus in the near future, a larger dataset will be available, and we will be able to obtain more conclusive results.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/restrictions: access to data must be requested and approved by the EpiBioS4Rx steering committee.
Requests to access these datasets should be directed to epibiossteeringcommittee@loni.usc.edu.

ETHICS STATEMENT
This work was approved by the UCLA Institutional Review Board (IRB# 16-001 576) and the local review boards at each EpiBioS4Rx Study Group institution. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin. EpiBioS4Rx is the result of efforts of many investigators from a broad range of academic institutions and private corporations, and subjects have been recruited from over 30 sites across the world. EpiBioS4Rx data are disseminated by the Laboratory of Neuro Imaging at the University of Southern California. A complete listing of EpiBioS4Rx investigators can be found at: https://epibios.loni.usc.edu.