Mild traumatic brain injury (mTBI) is a significant cause of brain insult (Len and Neary, 2011; Huang et al., 2014) representing close to 90% of all brain injuries (Len and Neary, 2011). Approximately 5–20% of the irremediable patients (Bharath et al., 2015) still suffer from post-concussion symptoms several months after the initial injury (Huang et al., 2014). These symptoms are often characterized by physical, emotional, cognitive, and sleep disturbances and may take many months to return to the baseline (Huang et al., 2014). Several neuropsychological studies have reported reduced cognitive efficiency, especially in tests measuring executive function, processing speed, attention, connectivity, and memory in mTBI patients with persistent symptoms (Huang et al., 2014; Pang et al., 2016). Management of mTBI is crucial due to its deleterious effects on certain brain functions, including attention (De Monte et al., 2006), working memory (Vanderploeg et al., 2005), and verbal learning (De Monte et al., 2006).

The human brain can be viewed as a neurophysiological network of brain areas that are structurally and functionally interconnected. These distinct networks are temporally and spatially distributed and exist in a range of spatiotemporal scales. Spatially, they extend from microscopic networks of single neurons and local synaptic interactions to large-scale brain networks interconnected via long white-matter tracts (Eierud et al., 2014). The time domain scales vary from milliseconds to seconds (Dimitriadis et al., 2013a, 2015a; Hansen et al., 2015; Betzel et al., 2016). Considering that interactions among large-scale networks are significant for high-level cognitive functions, most recent functional connectivity (FC) studies of mTBI focus on large-scale intrinsic connectivity networks (ICNs) aiming at identifying the changes they undergo as a result of injury. FC is reflected in the neurophysiological activity of neural populations that mediate cortical communication and information integration (Wang, 2010). Clinically, FC has been shown useful in the study of several neurological and neuropsychiatric disorders and their symptoms (Tewarie et al., 2013, for a review see Eierud et al., 2014).

During the past several years, numerous studies, including ours, have attempted to develop reliable biomarkers of mTBI based on resting state MEG. While other imaging modalities measure brain activation indirectly, such as functional Magnetic Resonance Imaging (fMRI) which relies on hemodynamic events, MEG measures neuronal activity directly. MEG possesses most of the highly desired neuroimaging properties, including high sensitivity and efficient handling of environmental noise, and combines very good spatial details with excellent temporal resolution (Dimitriadis et al., 2013b; Antonakakis et al., 2015, 2016). Numerous analysis techniques applied to EEG and MEG recordings have clearly demonstrated altered functional connectivity in TBI that is closely correlated with disease severity (Castellanos et al., 2010). Zouridakis et al. (2012) using Granger causality showed that mTBI patients exhibited a sparsely distributed network of long-range connections compared to controls during the first few weeks after mTBI. Tarapore et al. (2013) found that resting state MEG could detect abnormal connectivity in TBI, while Da Costa et al. (2015) suggested that MEG could detect subtle neural changes associated with cognitive dysfunction in mTBI. Huang et al. (2014) uncovered anatomical and functional correlations between abnormal slow waves and mild axonal injury. Dunkley et al. (2015) found increased connectivity in mTBI that was limited to slow wave amplitude coupling. Dimitriadis et al. (2015b) used phase-locking value estimates to quantify intra-frequency couplings and showed that this pattern was mostly seen in the delta band, whereas Antonakakis et al. (2015, 2016), in a follow-up analysis based on inter-frequency couplings, showed that controls formed a dense network of stronger local and global connections in agreement with other studies (Rapp et al., 2015).

Among the various procedures to estimate FC networks, also known as FC graphs (FCG), intra-frequency measures, such as mutual information (MI) (Bullmore and Bassett, 2011; Tsiaras et al., 2011), and inter-frequency measures, such as cross-frequency coupling (CFC) (Antonakakis et al., 2015, 2016; Dimitriadis et al., 2015c; Florin and Baillet, 2015), are of special interest in the topological analysis of FCGs.

Several brain connectivity studies using fMRI and electro- (EEG) and magneto- (MEG) encephalography have suggested that FCGs exhibit properties of a “small world” (SW) network organization (or simple organization) (Micheloyannis et al., 2006; Palva and Palva, 2011; Dimitriadis et al., 2015a; Vértes and Bullmore, 2015). SW is a special type of mathematical graph in which the majority of the nodes are not direct neighbors; yet, most nodes can be reached from any other node within a small number of steps. Thus, SW networks are simultaneously highly clustered and highly efficient. In particular, when considered as nodes that are connected to each other, SW networks are likely to have many first degree neighbors in common, and the average path length between a pair of nodes is short (Vértes and Bullmore, 2015).

An additional property often seen in SW organization is the formation of certain nodes (termed “hubs”) that are more densely connected to each other than the rest of the nodes (Van den Heuvel and Sporns, 2011). Hubs with the highest interconnectivity values follow what is known as “rich club” (RC) organization and appear to be the most relevant nodes in a network in terms of global information processing (Van den Heuvel and Sporns, 2011; Van den Heuvel et al., 2013; Schroeter et al., 2015; Vértes and Bullmore, 2015). Thus, the RC organization can be seen as a variation of the SW network organization, with different topological properties (Mišić et al., 2014) featuring disproportionately dense interconnections (Van den Heuvel and Sporns, 2011; Bullmore and Sporns, 2012; Mišić et al., 2014). These hubs support more traffic than the ordinary SW nodes (Mišić et al., 2014), and effectively define the top-level structure of a network, its hierarchical ordering, and node specialization (Van den Heuvel and Sporns, 2011).

Typically, the largest amount of information flow between pairs of nodes in human brain networks passes through RC and SW hubs (Palva and Palva, 2011; Van den Heuvel and Sporns, 2011; Vértes and Bullmore, 2015). As a consequence, the SW organization has been studied in several brain disorders, including schizophrenia (Micheloyannis et al., 2006), Alzheimer's disease (Stam et al., 2007), autism (Rubinov and Sporns, 2010; Tsiaras et al., 2011), and epilepsy (Bharath et al., 2016). The RC networks have been explored in brain ischemia (Van den Heuvel et al., 2013; Watanabe and Rees, 2015; Crossley et al., 2016), in healthy subjects using diffusion tensor imaging MRI (Van den Heuvel and Sporns, 2011) and recently in studies of brain activation (Van den Heuvel and Sporns, 2011; Bullmore and Sporns, 2012; Mišić et al., 2014; Antonakakis et al., 2015) whereby RC hubs were shown to play a pivotal role in information integration.

Given our previous use of connectivity analysis to study mTBI using Granger causality (Zouridakis et al., 2012), phase synchronization (Dimitriadis et al., 2015b), and cross-frequency coupling (Antonakakis et al., 2015, 2016) of spontaneous MEG, as well as brain activation patterns of both EEG and MEG at the sensor (Li et al., 2015) and source (Zouridakis et al., 2016; Li et al., 2017) levels, an obvious question is whether the possible presence of an RC organization could provide some complementary features to the SW organization that is typically seen in mTBI FCGs. Thus, in the present study, we hypothesize that exploring brain connectivity network models derived from spontaneous MEG activity using estimators for both intra and cross-frequency couplings (Buzsáki and Watson, 2012) would help identify meaningful network topological features in compromised mTBI brain networks that could be used as guideline biomarkers for validating the recovery from mTBI (Bharath et al., 2015). For a better understanding of network topologies linked to mTBI, we followed an attack strategy to reveal the key network model, either small-world or rich-club, that best describes mTBI functional brain networks. Our network analysis was based on intra- inter-frequency, and integrated functional brain networks to cover the different aspects of the multiplexity of brain rhythms via brain connectivity.

The SW and RC organizations, however, are not mutually exclusive, considering that an RC network may also present SW characteristics in a subnetwork (Bullmore and Sporns, 2012). In particular, subareas of an RC network can simultaneously be part of a SW network featuring SW behavior. Thus, to identify the network organization that explains best the network topology of intra-frequency and inter-frequency FCGs, for both normal controls and mTBI patients, we implemented an “attack strategy” (Van den Heuvel and Sporns, 2011; Antonakakis et al., 2015) on SW and RC nodes to reveal their relative importance in information transfer and neural communication within the entire brain network. The present study is an extension of our recent short report (Antonakakis et al., 2015) which compared the RC and SW organizations only in the δ-β cross-frequency pair and found that resting state MEG FCGs followed a rich-club organization.

In the present study, we compared the RC and SW organizations under three different FCG representations, namely IFCG, CFCG, and ICFCG. Our aim has been to develop reliable biomarkers for the accurate detection of network abnormalities caused by mTBI by analyzing functional brain connectivity profiles that combine the RC organization with prominent MEG intrinsic coupling modes. We assumed that the highest reduction in %GE would reveal the organization that best described the topology of the “true” underlying brain networks. Our results showed that separation of the mTBI patients from the normal control group was feasible with ~91% accuracy based on the ICFCG representation (Table 1).

The importance of the RC organization was first investigated in neuroimaging brain studies (Van den Heuvel and Sporns, 2011) using diffusion tensor imaging (DTI) that emphasized the role of RC organization in information integration and in conferring robustness to its structural core. Since then, many human connectomic studies have followed (Van den Heuvel et al., 2013; Senden et al., 2014; Vértes and Bullmore, 2015; Crossley et al., 2016) mostly using DTI or fMRI on patients and normal control subjects. Recently, we investigated a prominent RC organization obtained from MEG recordings of spontaneous brain activity of mTBI patients for a specific cross-frequency pair (Antonakakis et al., 2015). In another recent study, we presented a promising imaging method for detecting network abnormalities caused by mTBI during task-free spontaneous MEG activity based on intra- (Dimitriadis et al., 2015b) and cross-frequency coupling (Antonakakis et al., 2015, 2016). In this study, we performed a thorough investigation on these issues based on different interaction metrics, network formation approaches, and evaluation scenarios, over all frequency bands and band pairs. We observed a significant direct dominance (Figure 2) of low frequency bands and band pairs regarding the distribution of RC nodes (Figure 2) for both groups. More specifically, temporal areas played a crucial role in the definition of the CFCG-RC topology in both groups. The RC nodes for the mTBI group were mainly seen in right temporal areas compared to controls (Figure 2B).

In a previous analysis of the same dataset, we classified correctly the two groups with more than 90% accuracy (Antonakakis et al., 2015, 2016) using only the CFCG frequency pairs. A recent study using resting state MEG (Dunkley et al., 2015) showed slow abnormal brain activity in mTBI patients. The trend was an increase of low frequency amplitude in patients with mTBI, and this spectral alteration appeared most prominently in temporal regions. In numerous studies, these findings have been linked to brain injury (Lewine et al., 2007; Huang et al., 2014). We found reliable evidence (Figure 2) that the highest SW and RC distributions were mainly located over temporal areas. The absence of significant differences in the distribution of the CFCG-SW (Figure 2B, upper part) was also reported by a recent study (Athanasiou et al., in press) using the complementary modalities of EEG/MEG and multivariate functional connectivity. Over-representation in right frontal areas has also been found by previous studies, suggesting that the effects of mTBI are more prevalent in frontal regions which are more vulnerable to brain injury (Eierud et al., 2014; Antonakakis et al., 2016). A recent MEG study with a mixed group of mild, moderate, and severe TBI patients showed reduced FC in frontal areas bilaterally and increased FC in left temporo-parieto-occipital regions and in the right thalamus (Tarapore et al., 2013).

Revealing a consistent over-representation of hubs in right frontal areas in mTBI subjects across all frequency bands provides further support that structural alterations cause these frequency changes over the whole frequency spectrum (Dunkley et al., 2015). The significant over-representation of hubs in left occipital brain areas (Figure 2B, lower part) of the mTBI group may be indicative that, unconsciously, subjects were experiencing mental images caused by the trauma as spontaneous brain activity. It would be interesting to explore the activity of mTBI patients in the fusiform gyrus and compare it to controls using neuromagnetic source reconstruction (Zhan et al., 2016). Further analysis on the source level could reveal if there is increased connectivity between visual system and default mode network.

Furthermore, we found that, compared to SW, the RC organization underwent significantly higher GE changes in both groups, with a much larger reduction in the mTBI group. These findings indicate that the RC organization can encode important topological features of the underlying brain networks (Figure 6). These findings in particular extend our previous results from comparing the SW and RC models only on CFCGs and in the δ-β frequency pair (Antonakakis et al., 2015). Through the adopted attack strategy, the RC organization demonstrated higher levels of %GE reduction and confirmed the damage inflicted on a node by attack, which was simulated as a decrease in the weights of the node's connections (Van den Heuvel and Sporns, 2011). Mišić et al. (2014) also demonstrated that a large part of information flows in the RC organization, suggesting that it is better modeled by an RC organization than an SW. The reduction of GE in mTBI could be linked to the significant reduction of functional connectivity within the default mode network revealed by a recent mTBI study working with MEG resting-state on the source level (Alhourani et al., 2016). In a next study working at the source level, we will try to reproduce our results.

Additionally, to uncover the different role and definition of hubs in a brain network, we explored the mean degree and mean strength of RC and SW hubs in two conditions, namely within their subnetworks and with the rest of the network. Our analysis revealed significantly higher mean values in degree and strength in the control group compared to the mTBI in RC topology for both conditions in the δ and θ bands (Figure 4A). Based on these findings, we conclude that the RC definition can uncover the subset of nodes from a network that plays a pivotal role in global information integration. A collapse of such highly interconnected hub regions can cause communication aberrations between different parts of the brain (Van den Heuvel and Sporns, 2011). Thus, these brain hubs should be further studied both structurally and functionally in various brain diseases and disorders.

Through the examination of the strength ratio of RC nodes (Figure 7), we found hyper synchronization in mTBI patients compared to normal controls within RC subnetworks based on intra- and inter-frequency intrinsic couplings in the δ, δ-γ₁, θ-γ₁, and β-γ₂ frequencies and frequency pairs. These results agree with the basic findings of Hillary et al. (2014) who showed that brain connectivity increased after TBI— they examined the mean degree of 52 nodes and showed significantly greater connectivity in TBI patients compared to normal controls. Instead, the highest degree nodes (i.e., RC nodes in our case) were selectively observed in several core subnetworks (either TRC or HRRC). Our analysis further explains hyper-connectivity in mTBI patients under the RC model and the use of both intra- and inter-frequency coupling.

In conclusion, this study demonstrated that the RC organization of graphs reflecting the local distribution of activity from resting-state brain networks can encode characteristic aspects of two types of FCG for both mTBI patients and control subjects. RC hubs in mTBI subjects were overrepresented in right temporal areas from θ to γ₁ frequencies and underrepresented in left occipital areas in δ-β, δ-γ₁, θ-β, and β-γ₂ frequency pairs. Therefore, our analysis does not only support the use of resting state MEG for the extraction of meaningful features that describe abnormal brain connectivity after mTBI (Buzsáki and Watson, 2012; Zouridakis et al., 2012; Huang et al., 2014; Dimitriadis et al., 2015b; Dunkley et al., 2015; Antonakakis et al., 2016) it also indicates the need to explore the RC organization under different types of interactions (intra- and cross-frequency) for the development of complementary connectomic biomarkers of recovery from mTBI that can be useful in clinical research (Buzsáki and Watson, 2012). The potential functional implications of RC organization of MEG intrinsic coupling modes, considering its role in network integration and its vulnerability in various disorders like mTBI, seem to deserve further investigation for diagnostic and clinical purposes. Furthermore, our approach is suitable for accessing the recovery process following mTBI using resting state MEG (Zouridakis et al., 2016; Li et al., 2017) and focusing not only on the strength of the couplings but also on the dominant type of interactions (Bharath et al., 2015; Losoi et al., 2015).

In order to provide a robust spatial mapping for brain functions in both control subjects and mTBI patients, it is necessary to adopt a dynamic functional connectivity approach (Dimitriadis et al., 2012, 2013a, 2014, 2015a,c, 2016) through the definition of functional connectivity microstates (Dimitriadis et al., 2013a) and/or network microstates (Dimitriadis et al., 2015c). Our current studies focus on exploring a dynamic combination of IFCGs and CFCGs into ICFCGs (Dimitriadis et al., 2016), their related microstates, and their symbolic dynamical signature in a possible combination with compression methods (Luo et al., 2013).

MA is currently a Ph.D. student in the Institute of Biomagnetism and Biosignal analysis. MA has performed the whole data analysis and written the first draft of manuscript. SD significantly contributed his ideas for the structure of the pipeline and co-supervised together with MZ the analysis of the data. AP offered the data-set and his old group was responsible for the collection of it. GZ developed with his experience the main body of the manuscript and made the final check of the whole manuscript.