traumatic brain injury is a major cause of sustained morbidity and disability both in the military and civilian populations ( tarapore et al . , 2013 ) .
in the case of mild traumatic brain injury ( mtbi ) , an initial brief change in mental state or consciousness is typically followed by postconcussion symptoms ( cassidy et al . , 2004 ) , such as headaches , fatigue , and dizziness , which usually emerge on the day of injury and persist for at least the first few days thereafter ( boccaletti et al . , 2006 ) .
in most patients cognition recovers and symptoms resolve within 3  months ; however , up to 25% of patients ( sigurdardottir et al . , 2009 ) suffer residual symptoms , long - term impairment , and sometimes disability ( levin , 2009 ) . in a large percentage of patients
mtbi is difficult to diagnose because of the absence of apparent external injuries and obvious focal brain lesions in conventional computed tomography ( ct ) or magnetic resonance imaging ( mri ) scans ( tarapore et al . , 2013 ) , although mtbi - related morphological changes at a microscopic level have been recently reported using mri ( pasternak et al . , 2014 ;
a common practice is to classify tbi as mild , moderate , or severe based on the level of consciousness assessed using the well - known glasgow coma scale ( gcs , teasdale and jennet , 1974 ) .
some studies have attempted to classify tbi and predict long - term outcome based on ct findings ( zhu et al . , 2009 ) or the amount of exosomes released in the peripheral circulation due to injury ( taylor and taylor , 2014 ) .
one mechanism by which tbi is thought to affect cognition and behavior is through changes in functional connectivity ( castellanos et al .
, 2010 ; martino et al . , 2011 ) . functional brain connectivity refers to temporally correlated activity in spatially distinct brain regions and is believed to mediate neural communication among these regions ( aertsen et al . , 1989 ; salvador et al . , 2005
several imaging studies have demonstrated the existence of a network of functionally connected brain regions that support a default mode of brain function ( guskiewicz et al . , 2005 ) .
changes in brain activity caused by focal lesions result in abnormal interactions among local and remote brain regions that are functionally connected ( quigley et al . , 2001 ;
van cappellen van walsum et al . , 2003 ) , and these changes are reflected on the configuration of the default brain network .
thus , functional connectivity analysis offers a unique opportunity to quantify changes in neurophysiological processes that correlate well with clinical symptoms ( huang et al .
2009 , huang et al . 2012 ; lewine et al . , 2007 ) .
the literature shows several examples of pathological alterations , including increase ( van dellen et al . , 2012 ) and decrease ( sharp et al
2010 , rubinov , and bullmore 2013 ) , which provide evidence that a balance in the synchronization level in healthy controls is required for optimal brain function .
several brain connectivity studies in clinical populations focus on the description of topological differences with control group employing various synchronization measures , network metrics and statistics ( sporns , tononi , and ktter 2005 , guggisberg et al .
stam et al . , 2007 ; bassett et al . , 2009 ; bullmore and sporns , 2009 ; castellano et al . , 2010 ; dimitriadis et al . , 2013b ; bigler , 2013 ) . in the area of brain injury in particular
, during the last two decades there has been a growing effort to characterize the structural and functional effects of mtbi using the full range of neuroimaging methods ( see the excellent review by eierud et al . ,
2014 ) , and various biomarkers have been proposed based on functional mri ( fmri ) and diffusion tensor imaging ( dti ) ( kou et al . , 2010 ; hunter et al . , 2012
, meg plays a special role because it measures neuromagnetic fields resulting from primary neuronal currents that are not subject to the distortion induced by the variable conductance of brain tissues ( tarapore et al . , 2013 ) .
2011 , qian et al . 2013 ) because they require no training or experience with cognitive tasks , and they impose minimal demands on a patient , which is especially important after brain injury .
thus , resting state meg is particularly suitable for detecting changes in functional connectivity networks in patients with mtbi compared to matched controls .
one of the very first studies using resting - state meg as a possible neurophysiological biomarker for mtbi provided clear evidence that analysis of functional connectivity patterns could be a valuable tool for early detection of mtbi ( zouridakis et al . , 2012 ) .
other studies showed abnormal slowing in brain areas affected by tbi ( huang et al . , 2014 ) and reduced overall functional connectivity in tbi patients compared to controls ( tarapore et al . , 2013 ) . in particular ,
using resting - state meg and low - frequency source imaging , huang et al .
( 2012 ) showed 87% and 100% accuracy in detecting abnormalities in mild and moderate tbi , respectively .
furthermore , lewine et al . ( 1999 ) used mri and meg to examine postconcussive symptomatology and were able to discriminate between healthy adults and individuals with resolved mtbi . another meg study attempting to characterize system complexity and potential neural network damage found reduced complexity in multiple brain areas in mtbi subjects relative to healthy controls ( luo et al . , 2013 ) .
other studies further demonstrated that decreased connectivity in resting - state meg may persist for years after mtbi , even in mild tbi cases ( castellanos et al . , 2011 ) , but the abnormally reduced connectivity might improve over time across serial meg scans , suggesting neuroplasticity during recovery from tbi ( tarapore et al . , 2013 ) .
continuing our earlier attempt ( zouridakis et al . , 2012 ) to understand how mtbi affects communication in the human brain network , the present study extends our previous work in several ways : first , we employ a much larger sample of patients and controls , 31 mtbi patients and 50 normal controls , compared to only 10 and 10 , respectively ; second , we use the phase locking value estimator and graph theory to compute functional connectivity graphs ( fcgs ) , whereas in the earlier study we used granger causality as the estimator ; and third , we follow a tensor - based instead of vector - based approach for subject classification . here
we consider fcgs as second - order tensors , which , after dimensionality reduction , are used for subject classification , whereas the earlier study considered the granger connectivity matrices as one - dimensional vectors and used them for subject classification ; and fourth , we performed the analysis in several frequency bands , instead of only one , based on recent studies that found differences in fcgs across frequency bands ( illebrand et al . , 2012 ; engel et al . , 2013 ) .
the remainder of this paper is structured as follows : section 2 describes the data and the detailed analysis procedure , section 3 presents the classification results , section 4 is devoted to discussing our findings , and section 5 summarizes some concluding remarks .
the overall methodology employed in this study is outlined in fig . 1 and includes the following steps : data collection , preprocessing , wavelet decomposition , computation of functional connectivity graphs ( fcgs ) , selection of significant links , feature extraction , optimization of the number of features , designing a classifier , evaluating its performance , and assessing the topological properties of the fcgs .
the present project is part of a larger study of mtbi , supported by the department of defense ( dod ) .
the definition of mtbi used in this study followed the guidelines of dod ( assistant secretary , 2007 ) and the american congress of rehabilitation medicine ( kay t. , 1993 ) .
this work was approved by the institutional review boards ( irbs ) at the participating institutions and the human research protection official 's ( hrpo ) review of research protocols for dod .
all procedures were compliant with the health insurance portability and accountability act ( hipaa ) .
subjects included in this analysis included a group of mtbi subjects from the dod project and a group of age - matched normal controls drawn from a database that was being assembled as a normative data repository at uthsc - houston .
the mtbi group consisted of 31 right - handed mtbi participants ( 29.33    9.2  years of age ) .
subjects were recruited from the emergency departments ( eds ) of two level 1 trauma centers and one level iii community hospital in a large ethnically diverse southwestern metropolitan area .
subjects were recruited by healthcare professionals ( rn , md , emt - p ) who had clinical experience with brain injury patients , knowledge of research , and excellent interpersonal and problem - solving skills .
screening occurred through review of data in the eds ' electronic healthcare system ( ehs ) , consultation with ed staff , and subject interviews .
special permission was obtained from the institutional irbs to administer the galveston orientation and amnesia test ( goat ) ( levin et al . , 2008 ) prior to obtaining informed consent to identify cognitive impairment that would preclude provision of informed consent .
all subjects showed goat scores of 75 or greater and so have provided informed consent .
the inclusion criteria for all subjects in the dod project included age 1850  years , injury occurring within the preceding 24  h , and no requirement for hospitalization for the injury for which the participant was enrolled . for mtbi subjects , the inclusion criteria also required the presence of a head injury ( documented in medical records and/or verified by witnesses ) , glasgow coma scale ( gcs ) ( teasdale and jennett , 1974 ) score 1315 , loss of consciousness < 30  min including 0  min , post - traumatic amnesia < 24  h including 0  min , and a negative head computed tomography ( ct ) scan .
the exclusion criteria included a score on the abbreviated injury scale ( ais ) > 3 for any body part , history of significant pre - existing disease ( e.g. , psychotic disorder , bipolar disorder , post - traumatic stress disorder ( ptsd ) diagnosed by a psychiatrist or psychologist , past treatment for alcohol dependence or substance abuse ) , blood alcohol level > 80  mg / dl at the time of consent , documentation of intoxication , left - handedness , and contraindications for mri ( including claustrophobia and pregnancy ) .
the demographics of mtbi subjects and the location of injury are given in table 1 .
the control group ( group n ) consisted of 55 right - handed age - matched normal subjects ( 29.25    9.1  years of age ) .
previous hospitalization for mtbi , history of neurologic disorder , schizophrenia , bipolar disorder , substance abuse , and extensive dental work and implants incompatible with meg were also the exclusion criteria for the control subjects .
the research protocol for the control group also received irb approval , and all participants provided written informed consent .
subjects were asked to lie down and remain as still as possible during the recording procedure with eyes closed .
data collection in this large multicenter study involved two different meg scanners , a 248-channel whole - head magnes wh3600 system ( 4d neuroimaging inc . , san diego , ca ) featuring axial gradiometer sensors , and a 306-channel elekta neuromag system ( elekta ab , stockholm , sweden ) featuring 204 planar gradiometers and 102 magnetometers .
to make the recordings from various projects of the larger study comparable , all data collected with axial gradiometers were transformed to planar gradiometer field approximations using the sincos method implemented in fieldtrip ( oostenveld et al . , 2011 ) .
planar gradiometers have maximum sensitivity to superficial cortical sources directly under them ( hmlinen , 1995 ) , which makes them less sensitive to artifacts and distant disturbances .
natick , ma , usa ) and fieldtrip ( oostenveld et al . , 2011 ) .
axial gradiometer data were originally collected at a sampling rate of 1017.25  hz and bandpassed between 0.1200  hz using hardware filters , whereas the planar gradiometer data were collected at 1000  hz and bandpass filtered via software between 0.1200  hz .
artifact reduction was performed based on independent component analysis ( ica ) as implemented in the eeglab toolbox ( delorme and makeig , 2004 ) . the extended informax algorithm ( delorme et al . , 2007a ; onton et al .
, 2006 ; romero et al . , 2008 ) was used to remove components associated with eye , muscle , and cardiac artifactual activity ( dimitriadis et al .
in addition to visually inspecting the time course and topographical layout of all ics , we estimated their kurtosis and skewness values . components reflecting cardiac activity were recognized from their typical rhythmic pattern in the time domain and their widespread topography .
ics related to muscle activity were identified based on statistical measures ( kurtosis higher than a predefined threshold , kurtthr  =  12 , determined empirically ) , spectral characteristics ( increased energy in the frequency range between 2060  hz ) , and topographies encompassing temporal brain areas ( delorme et al . , 2007b ; dimitriadis et al .
artifactual ics were zeroed , whereas the remaining ics were used to generate artifact - free signals that were back - projected to the original meg sensor space via the inverse ica transformation .
occasionally , activity from a bad sensor was replaced with the interpolated activity of its immediate neighbors after applying ica to avoid introducing dependencies in the recordings .
this interpolation procedure was necessary for five control and two mtbi subjects for at most four meg sensors .
meg signal oscillations in specific frequency bands were studied using wavelet transform analysis ( mallat , 2008 ) .
after preprocessing , all data were further lowpass filtered and decimated by a factor of six , giving an effective sampling rate of 166  hz , which constrained the frequency bands of the wavelet transform to roughly conform to the classical frequency bands ( assett et al . , 2009 ) .
decimation was used to reduce the processing time , while the new sampling rate was higher than twice the minimum rate required by the nyquist theorem to avoid distortion of the meg signals in the effective range of frequencies of interest ( 160  hz ) .
each frequency band was isolated by applying the maximum overlap discrete wavelet transform to each time series ( daubechies , 1992 ) using the daubechies 5 wavelet in line with previous work ( assett et al .
thus , wavelet scale 1 ( 3060  hz ) roughly corresponded to the classical low- band , scale 2 ( 1530  hz ) to  , scale 3 ( 815  hz ) to  , scale 4 ( 48  hz ) to  , and scale 5 ( 14  hz ) to . we constructed fcg by quantifying the coupling between pairs of sensors using the phase locking value ( plv ) as an estimator ( lachaux et al . , 1999 ) , separately for each frequency band .
plv is defined as(1)plv=1nt=1nej(t)where (t ) is the difference in phase between the two time series and n is the number of samples .
if the phase difference varies by a little , plv is close to one ; otherwise it is close to zero .
plv can assess phase relationships between two signals but it can not detect any causal relationships ( lags or delays ) between the signals . to construct intra - frequency functional connectivity matrices
, we calculated the plv between the time series in each frequency band for all possible pairs of sensors to create the weights matrix fcg . through this process
, we created 5 such connectivity matrices (  ,  ,  ,  , and  ) per subject . for each subject , the aforementioned procedures yielded a fully connected , weighted , symmetric fcg representing the mutual influences among all cortical regions .
the maximum number of possible undirected connections in a network with k nodes is nmax  =  k(k    1 )  /  2 . when k  =  248 , the corresponding nmax  =  30,628 , and the fcg is extremely dense .
therefore , it must be filtered so that the pattern of the most significant connections can emerge .
we performed two kinds of filtering based on statistical analysis and graph theory principles . for statistical filtering a surrogate data method ( lachaux et al . , 1999 ) with multiple realizations was used to select the most significant fcg values at the 99% confidence level . generating 200 surrogates from the original time series was enough to reveal significant interactions ( lachaux et al . , 1999 ) .
surrogates were obtained by randomizing the phase of the original signals while keeping their spectral profile intact to remove all temporal correlations ( theiler , 1992 ) .
after computing the surrogate plvs , we assigned a p - value to each connection in the fcg matrix by estimating the proportion of surrogate plv values that were higher than the observed values ( theiler et al . , 1992 ) .
to correct the effects of multiple comparisons , p - values were adjusted using the false discovery rate ( fdr ) method ( benjamini and hochberg , 1995 ) . a threshold of significance q
was set such that the expected fraction of false positives was restricted to q    0.01 ( dimitriadis et al . , 2012 ; ioannides et al . , 2012 ) .
graph theory - based analysis ( bullmore and sporns , 2009 ; bassett and bullmore , 2009 ; he and evans , 2010 ; stam , 2010 ) was used to capture the structure of the neural system under investigation and the relationship between separation and integration of neural populations .
small - word structures are characterized by a dense network of local connections and a limited number of long - range connections that provide efficient communication between distant nodes .
efficiency in information transmission between nodes is measured as the inverse of the shortest distance between the nodes , while the average of all pair - wise efficiencies represents the global efficiency of the graph .
the function cost relates to the energy expenditure needed for a network to maintain its efficiency , and it is given by the ratio of existing connections divided by the total number of possible pairwise connections in a network .
global cost efficiency is defined as the global efficiency e at a given cost c minus the cost ( e    c ) , which typically has a positive maximum value at some cost cmax , for an economical small - world network .
importantly , this metric of network topology is independent of arbitrary , investigator - specified thresholds . instead
, the cost efficiency curve is estimated over a wide range of thresholds , and the behavior of the curve is summarized by its maximum value , which occurs at a data - driven connection density or cost c ( bassett et al . , 2009)1 the physical distance among the meg sensors was used to separate the local from the long - range connections .
for each sensor , we computed its neighbors of order one , two , three , and four , i.e. , neighbors in all directions that were one , two , three , and four steps away , based on the topographical layout of the meg sensors ( zouridakis et al . , 2012 ) . to optimize the performance and generalizability of our methodology ( hastie et al . , 2003 )
, we implemented a true external cross - validation procedure by splitting the data into two subsets : a training set that included 80% of the subjects ( 44
n and 25 mtbi ) and was used for designing the classifier and a test set that included the remaining 20% of the subjects ( 11 n and 6 mtbi ) and was used for assessing the classifier performance . as a metric
, we used the correct classification rate , which was defined as the proportion of subjects in the test set for which the correct label was predicted .
the entire cross - validation procedure was repeated 100 times to avoid any adaptation effects , and the average value of the correct classification rate was used to compute the overall accuracy of the procedure .
the optimum number of feature to represent the data was determined using only the training dataset in a five - fold cross - validation procedure as described in fig .
test subjects were removed from the analysis at the very beginning of the procedure , prior to feature selection and five - fold cross validation .
most of the previous studies in this area represent fcgs as vectors in a high - dimensional space ( shen et al . , 2010 ; pollolini et al . , 2010 ; richiardi et al . ,
the main limitation of this approach ( and the related feature extraction algorithms ) is that it overlooks the inherent form of fcgs .
each fcg has a natural tabular representation and hence it can be thought of as a second - order tensor .
the relationship between the row and column vectors of the associated matrix might be important for deriving a suitable low - dimensional representation ( projection ) , especially when the number of available connectivity patterns is small .
to overcome the limitations of previous approaches , in this study we treated fcgs as tensors and resorted to tensor subspace analysis ( tsa ) ( he and cai , 2005 ) as the most appropriate feature extraction algorithm . in our formulation
, the tensor form was given by ( subjects    sensors    sensors ) ( dimitriadis et al .
fcgs from individual subjects were represented based on a recent methodology that blends ideas from multi - linear algebra and manifold data learning ( he and cai , 2005 ) .
briefly , given a set of fcgs sampled from the space of functional connectivity patterns , we approximated the underlying manifold by first building an adjacency graph that captured the proximity relationships among the connectivity patterns , and then derived a tensor subspace that faithfully represented these relationships .
mathematically , an fcg xn1n2 of size n1  n2 , where n1 and n2 denote the dimensions of a second order tensor ( in our case , n1 and n2 are equal to the number of meg sensors ) , can be thought of as a second order tensor in the tensor space     . then , the generic problem of linear dimensionality reduction in the second order space is expressed as follows : given a set of m tensors x1, ... ,xmn1n2 , find two transformation matrices u and v of size ( n1l1 ) and ( n2l2 ) , respectively , that map these m tensors to a new set of tensors y1, ...
,yml1    l2 , where yi = uxiv , l1 <  n1 , and l2 <  n2 , such that each new tensor yi represents uniquely a tensor xi in the original set
. the method is of particular interest in the special case when x1 , x2 ,  , xm    m , where m is a nonlinear sub - manifold embedded in n1n2 .
traditional linear dimensionality reduction algorithms , such as principal component analysis and linear discriminant analysis , find first - order mappings from  to   with l  <  n , but tsa finds a second - order mapping from n1    n2 to l1    l2 with l1 <  n1 and
domain of fcgs most probably forms a nonlinear sub - manifold embedded in the ambient space of second order tensors ( lu et al . ,
, we attempted to find a linear subspace approximation to the sub - manifold in the sense of local isometry using a technique that is the tensorial counterpart of the locality preserving projection ( he et al . , 2005 ) .
given a set of m tensors { xi}i  =  1:m , where each tensor is the tabular version of a single fcg having the group membership as class label ( n or mtbi ) , tsa starts by building an ( m    m ) weight matrix s that represents the nearest neighbor graph g among the tensors .
the value of m is given by the total number of subjects in the training dataset .
  xj||2/r )  condition  1  0  otherwise}where the exponential is known as the heat kernel and is employed with the frobenius norm ( golub and van loan , 1996 ) , r is a control parameter usually referred to as radius of influence , and condition 1 states that tensors xi and xj share the same class label and each of them is among the k - nearest neighbors of the other , while the indices i and j vary across the m tensors ( m equals the number of subjects in the training dataset ) .
then tsa seeks two transformation matrices u and v which , when applied to each tensor , result in a mapping that preserves the neighborhood relationship encoded in graph g. mathematically , this is formulated in the following optimization problem,(3)minu , vijutxivutxjv2sij . denoting with d a diagonal matrix with elements dii=jsij
, the above problem is reformulated as two coupled problems of eigenvector analysis ( he et al .
, 2005):(4)(du    su)v=  duvdu=idiixituutxi ,  su=ijsijxituutxj  ( 5)(dv    sv)u=  dvudv=idiixivvtxit  sv=ijsijxivvtxjt .
the optimal matrices u and v can be obtained by iteratively computing the generalized eigenvectors of eq .
( 4 ) and eq . ( 5 ) after initializing u with the identity matrix .
tsa detects the intrinsic local geometric structure of the tensor space by learning the topology of a new tensor subspace of lower dimension . the objective function in eq .
( 3 ) incurs a heavy penalty if neighboring points xi and xj are mapped far apart .
( 3 ) ensures that when the original tensors xi and xj are close , the new tensors yi  =  u xiv and yj  =  u xjv will be close as well .
the dimensionality of the reduced tensors yi , i.e. , the number of eigenvectors used for mapping yi  =  u xi v , was optimized using a five - fold cross - validation procedure that achieved maximum separation between the n and mtbi groups . specifically , we randomly partitioned the training dataset into five equal sized subsets and used four subsets for building a k - nearest neighbor ( k - nn ) classifier ( altman , 1992 ) and one subset for testing its performance .
the number of neighbors k and the heat parameter were set in a similar way .
our methodology employed the so - called extreme learning machine ( elm ) classifier ( huang et al . , 2006 ) , an emerging learning technique that provides efficient unified solutions to generalized feed - forward networks , including single- and multi - layer neural networks .
we chose this relatively new classification procedure because of its computational elegance and fast - learning capabilities , which lead to improved performance compared to other learning algorithms , like back propagation neural networks , radial basis function networks , and support vector machines ( kim et al .
, 2009 ) . the proposed scheme combined tensor subspace analysis with extreme learning machine classification and it was thus denoted as tsa+elm2 . for comparison purposes , we used two additional popular classifiers as a baseline : the k - nn algorithm with the frobenius norm ( horn and mathias , 1990 ) as a measure of similarity , and a support vector machine ( svm ) with a linear kernel ( cortes and vapnik , 1995 ) . furthermore , to demonstrate the superiority of the proposed tensorial approach
the fcg tensors were first vectorized , i.e. , converted to high dimensional vectors by traversing the corresponding matrices in a systematic way , and then their dimensionality was reduced using linear discriminant analysis ( lda ; fisher , 1936 ) for feature extraction .
the dimensions kept by lda were selected to optimize classification performance of the two groups .
therefore , our analysis employed various combinations of approach ( tsa , lda , and vectorial ) and classifier ( elm , k - nn , and svms ) . after selecting a classifier / approach combination , and optimizing it using the training data , the classifier was validated using the test data as follows : each fcg in the test set was compared against fcgs of known group membership ( normal or mtbi ) developed from the training dataset .
the fcg was finally assigned to the group with which it had the highest similarity . following the significant link selection procedure described previously ,
the resulting fcgs were described using the well - known topological metrics of global and local efficiency defined for weighted graphs .
the global efficiency ge of a network is given by(6)ge=1nin1n1jn , ji(dij)1where n represents the total number of nodes in the network and dij is the shortest path length between nodes i and j. ge is the inverse of the harmonic mean of paths between nodes and reflects the overall efficiency of parallel information transfer in the network ( achard and bullmore , 2007 ; latora and marchiori , 2001 ) . on the other hand ,
local efficiency le is a measure of fault tolerance in the network and is defined as(7)le=1n  innodal_lei  = 1ninj , hgi , j , hi(djh)1ki(ki1)where ki corresponds to the functional neighbors of the ith node , n represents the total number of nodes in the network , and dij is the shortest path length between nodes i and j. le indicates how well subgraphs exchange information when a particular node is eliminated ( achard and bullmore , 2007 ) . to identify cost - efficient functional networks ,
the same link selection procedure was applied to the average fcgs , independently for each frequency band and group .
statistical analyses were performed using the wilcoxon rank sum test ( p  <  0.001 ) .
table 2 summarizes the classification performance obtained for various combinations of approach ( vectorial , lda , and tsa ) and classifier ( k - nn , elm , and svms ) . for tsa , the parameters used were as follows , weight mode : eat kernel , neighbor mode : 10 , supervised learning , and number of dimensions : 6 ( dimitriadis et al .
2013a , dimitriadis et al . 2015 ) . for the k - nn algorithm we used majority voting . for each of the three classifiers ( k - nn , elm , and linear svm ) , the corresponding dimensions kept by lda in each frequency band were as follows ,  : ( 5 , 6 , 5 ) ,  : ( 6 , 6 , 5 ) ,  : ( 6 - 6 - 5 ) ,  : ( 7 - 6 - 7 ) , and  : ( 6 - 6 - 5 ) .
as it can be seen in table 2 , the proposed classification scheme of tsa+elm provided clear separation of the two groups , in all frequency bands , with high sensitivity and specificity .
in particular , in the  band normal controls and mtbi patients could be separated with 100% classification accuracy .
the combination of lda  +  k - nn showed reduced classification accuracy ( approximately 80% across all frequency bands ) compared to tsa  +  knn ( table 2a ) . on the other hand ,
the elm classifier showed improved performance compared to k - nn ( table 2b ) in all cases . using the svm classifier with tsa resulted in decreased performance compared to elm classifier ( table 2c ) .
however , the use of the svm classifier improved performance for both the lda and vectorial approaches compared to the elm classifier ( table 2c ) .
to demonstrate visually the ability of the proposed tsa+elm scheme to absolutely separate mtbi patients from normal controls in the  frequency band , we computed a three - dimensional representation for the fcgs , each of which was originally described by a vector of 36 dimensions ( in tsa each dimension of the second order tensor was 6 , thus the total number of dimensions needed for the vector representation was 6    6  =  36 ) .
we estimated the pair - wise similarity between all possible pairs of vectors a and b across all subjects using cosine similarity ( cos ) , which was given by(8)cosab
= abab=i=1naibii=1nai2i=1nbi2 cos(a , b ) measures the cosine of the angle between vectors a and b , and thus it is a judgment on their orientation , not on their magnitude .
the final cos estimates were tabulated in a matrix with dimensions 86    86 ( since 86 is the total number of subjects in the two groups ) . by employing multidimensional scaling ( borg and groenen , 2003 ) , a well - known dimensionality reduction technique
, we were able to project the original multidimensional matrix onto a set of three - dimensional points .
finally , we estimated the convex hull ( borst et al . , 1999 ) of the 3d points corresponding to each group .
2 , which shows complete separation of the two groups and also illustrates their variance .
3 illustrates the average functional connectivity patterns for each frequency band and group after applying the cost  efficiency thresholding scheme . as it can be seen in fig .  3 , fig .  4 , the functional connectivity profile of the control group exhibited a consistent network with stronger local connections ( level 1 and 2 neighbors ) whereas the mtbi group showed a variable network with stronger distant connections ( level 3 and 4 neighbors ) . across most frequencies , with the exception of the  band , in the n group more than 80% of the total number of connections was distributed locally ( level 1 and 2 neighbors ) , whereas in the mtbi group , more than 60% of the connections were shared between distant sensors ( level 3 and 4 neighbors ) .
additionally , from  to  , the n group exhibited significantly stronger connections locally , whereas the mean strength of distant connections was significantly higher in the mtbi group ( see fig .
5 shows , when long - range connections were compared across subjects , the n group exhibited mostly central - to - central or central - to - peripheral connections , especially in the  to  frequency bands , whereas the mtbi group showed connections involving primarily peripheral locations , practically bypassing fronto - central and centro - parietal brain regions . in all frequency bands , global efficiency ( ge ) values were in general higher in normal controls compared to mtbi patients , but without reaching significance ( p  =  0.61 , p  =  0.57 , p  =  0.61 , p  =  0.83 , and p  =  0.76 ) .
in contrast , local efficiency ( le ) values were significantly higher ( p  <  0.001 ) in the n group compared to the mtbi in all frequency bands , except for  ( fig . 6 ) .
in this study we used plv estimation and graph theory to investigate topological differences in the functional connectivity networks of mtbi patients and normal controls computed from resting - state meg .
configuration , i.e. , a dense network of local connections and a limited number of long - range connections . level 1 and 2 neighbors accounted for more than 80% of the total number of connections and , on average , they were significantly stronger compared to the ones seen in mtbi subjects .
in contrast , mtbi subjects showed sparse arrangements of weak local connections and a large number of stronger distant connections .
the average strength was significantly higher compared to the control group ( fig .  4 , fig .  5 )
however , these long - range connections in the mtbi group involved primarily peripheral locations , whereas in the controls , long - range links connected frontal - to - central or central - to - peripheral locations .
these findings echo the results of our earlier study ( zouridakis et al . , 2012 ) .
marked differences between the two groups were seen in all frequency bands , except for  , where the connection strength and sensor location were not significantly different between the two groups , at any neighborhood level . as an extension of our earlier work , in the present study we proposed a novel approach that features two improvements : estimation of functional connectivity graphs based on the plv estimator and a multivariate tensorial treatment of these graphs .
as table 2 shows , the elm classifier provided 100% classification accuracy in the  frequency band and more than 96% accuracy in the rest of the frequency bands studied .
biomarker to characterize mtbi , instead of the typical univariate estimates , like relative power .
this can be clearly seen in table s1 in the supplementary material , where we provide a detailed analysis based on relative power .
our findings of frequency - specific differences between the two groups match current research on spontaneous activity . a number of studies used meg source localization ( de pasquale et al .
, 2010 ; brookes et al . , 2012 ; hipp , engel , and siegel 2011 , hipp et al .
2012 ; hillebrand et al . , 2012 ; bardouille and boe , 2012 ; marzetti et al . , 2013 ) to explore two different modes of intrinsic coupling , namely phase coupling between band - limited oscillations and coupling between aperiodic fluctuations of signal envelopes ( hipp et al . , 2012 ) .
there is now evidence that the former is more sensitive to functional changes due to injury even in the absence of structural damage , as in the case of mtbi , while the latter can quantify functional abnormalities when structural damage predominates ( ipp et al . , 2012 ) .
these studies used various connectivity estimators , including power correlations ( de pasquale et al .
, 2010 ; brookes , woolrich , and barnes 2012 , brookes , woolrich , and barnes 2012 ; hipp et al . , 2012 ) , phase coherence ( hipp et al . , 2011 ; bardouille and boe , 2012 ) , phase lag index ( hillebrand et al . , 2012 ) , and imaginary coherence ( marzetti et al . , 2013 ) .
thus , our findings of  band differences between the two groups using the plv estimator are in line with the current literature . in terms of topology , even though our analysis was performed at the sensor level , our findings follow closely the results of studies that explored frequency - dependent networks at the source level ( hillebrand et al . , 2012 ) . in general , the brain areas that are part of the default - mode network showed a dense network of connections extended over the entire brain for low frequencies and a more anatomically constrained network for higher frequencies ( buzsaki , g. and draguhn , 2004 ; fox and raichle , 2007 ; florin and baillet , 2015 ) .
we found that , in the range of high frequencies , both the number ( 80% of total connections ) and strength of local connections were higher in the control group , whereas in the  band there were fewer ( only 20% ) and weaker connections .
additionally , we showed local topological differences with higher local efficiency for controls compared to mtbi subjects in the  band , which resulted in 96.48% classification accuracy and partially matched previous findings of abnormal brain activity in mtbi patients in the  band ( huang et al .
2012 ) . nevertheless , despite the very interesting results obtained , the present study has some limitations .
for instance , we estimated fcgs at the sensor level , although some studies have applied source reconstruction methods to identify activation patterns at the source level ( huang et al .
2009 , huang et al . 2012 ; hillebrand et al . , 2012 ) . even though we can not interpret the activation patterns in terms of intracranial cortical locations , we are still able to describe statistically significant differences in the connectivity patterns of mtbi and control subjects using surface brain regions .
another potential criticism of the present study could be attributed to volume conduction effects , i.e. , the transmission of electromagnetic fields from a primary intracranial source through biological tissue towards the measurement sensors resulting in distorted recordings . although these effects can be reduced using connectivity estimators that are robust to volume conduction , such as the phase lag index ( hillebrand et al .
, 2012 ) and imaginary coherence ( marzetti et al . , 2013 ) , or through orthogonalization preprocessing ( hipp et al . , 2012 ) , we decided not to employ any correction , since the absence of correlation between power and connectivity profiles in all frequency bands ( see the supplementary material ) was a strong argument that a common source did not alter the functional connectivity patterns ; otherwise , a common source would have changed also the signal power .
furthermore , any possible volume conduction effects would be common to both groups , but we were looking only for differences between the two groups . finally , some details on subject demographics are incomplete , as table 1 shows , making it difficult to establish correlations between lesion locations and connectivity patterns for every participant in the study . the present study has focused on the development of noninvasive biomarkers that can be used to uniquely characterize mtbi patients and normal controls .
it follows a general trend in current research in the area of neuroinformatics , where the big challenge for the next decade is the development and deployment of a large database with pre - computed brain connectivity profiles and derived network metrics ( e.g. , a  connectome database  ) in both healthy and mtbi subjects ( sporns , 2012 ) .
this would allow recording resting - state meg in the laboratory , estimating the related connectivity maps in both sensor and source space , querying the connectome database , and visualizing the results in a laboratory pc ( akil et al .
understanding brain function and visualizing the human connectome are slowly becoming a reality ( human brain project , europe and usa)3 . in the near future ,
computational approaches applied to connectomics will undergo rapid expansion , especially with shared neuroimaging resources , and complex network analysis and modeling ( jirsa et al . , 2010 ; sporns , 2012 ) , and may be helpful in designing optimal recovery strategies following traumatic brain injury ( irimia et al .
the proposed methodology can estimate patterns of brain activation in the form of connectivity networks that are highly consistent and can uniquely characterize subject groups using resting state meg , a noninvasive modality with superior temporal resolution over other neuroimaging techniques .
furthermore , resting state recordings are easy to obtain , as they do not require engagement from the subject , which makes such tests suitable for all sorts of patient types and age groups , including children .
these encouraging findings support the hypothesis that functional connectivity patterns may be used as general biomarkers for predicting a subject 's brain state ( dimitriadis et al .
2013a , dimitriadis et al . 2015 ) , for detecting various brain diseases ( dimitriadis , 2015 ) , and for developing image - based tools that can provide more accurate diagnosis , help guide treatment , and monitor the effectiveness of intervention in mtbi patients .
 mild traumatic brain injury ( mtbi ) may affect normal cognition and behavior by disrupting the functional connectivity networks that mediate efficient communication among brain regions . in this study , we analyzed brain connectivity profiles from resting state magnetoencephalographic ( meg ) recordings obtained from 31 mtbi patients and 55 normal controls . 
 we used phase - locking value estimates to compute functional connectivity graphs to quantify frequency - specific couplings between sensors at various frequency bands . 
 overall , normal controls showed a dense network of strong local connections and a limited number of long - range connections that accounted for approximately 20% of all connections , whereas mtbi patients showed networks characterized by weak local connections and strong long - range connections that accounted for more than 60% of all connections . 
 comparison of the two distinct general patterns at different frequencies using a tensor representation for the connectivity graphs and tensor subspace analysis for optimal feature extraction showed that mtbi patients could be separated from normal controls with 100% classification accuracy in the alpha band . 
 these encouraging findings support the hypothesis that meg - based functional connectivity patterns may be used as biomarkers that can provide more accurate diagnoses , help guide treatment , and monitor effectiveness of intervention in mtbi . 
Introduction
Methods
Results
Discussion
Conclusions