Resting-state magnetoencephalography source magnitude imaging with deep-learning neural network for classification of symptomatic combat-related mild traumatic brain injury.

Combat-related mild traumatic brain injury (cmTBI) is a leading cause of sustained physical, cognitive, emotional, and behavioral disabilities in Veterans and active-duty military personnel. Accurate diagnosis of cmTBI is challenging since the symptom spectrum is broad and conventional neuroimaging techniques are insensitive to the underlying neuropathology. The present study developed a novel deep-learning neural network method, 3D-MEGNET, and applied it to resting-state magnetoencephalography (rs-MEG) source-magnitude imaging data from 59 symptomatic cmTBI individuals and 42 combat-deployed healthy controls (HCs). Analytic models of individual frequency bands and all bands together were tested. The All-frequency model, which combined delta-theta (1-7 Hz), alpha (8-12 Hz), beta (15-30 Hz), and gamma (30-80 Hz) frequency bands, outperformed models based on individual bands. The optimized 3D-MEGNET method distinguished cmTBI individuals from HCs with excellent sensitivity (99.9 ± 0.38%) and specificity (98.9 ± 1.54%). Receiver-operator-characteristic curve analysis showed that diagnostic accuracy was 0.99. The gamma and delta-theta band models outperformed alpha and beta band models. Among cmTBI individuals, but not controls, hyper delta-theta and gamma-band activity correlated with lower performance on neuropsychological tests, whereas hypo alpha and beta-band activity also correlated with lower neuropsychological test performance. This study provides an integrated framework for condensing large source-imaging variable sets into optimal combinations of regions and frequencies with high diagnostic accuracy and cognitive relevance in cmTBI. The all-frequency model offered more discriminative power than each frequency-band model alone. This approach offers an effective path for optimal characterization of behaviorally relevant neuroimaging features in neurological and psychiatric disorders.

INTRODUCTION

Combat‐related mild traumatic brain injury (cmTBI) mainly due to blast exposure is a leading cause of sustained physical, cognitive, emotional, and behavioral deficits in Veterans and active‐duty military personnel. Eighty‐nine percent of blast‐related head injuries experienced in Iraq and Afghanistan conflicts were mild TBIs (mTBI; MacGregor, Dougherty, & Galarneau, 2011). Although mild head injuries frequently resolve within days post‐injury (Bigler, 2008), post‐concussive symptoms (PCS) persist for 3 months or longer in 7.5–40% of Veterans (Cooper et al., 2015; Morissette et al., 2011; Schneiderman, Braver, & Kang, 2008; Terrio et al., 2009). Persistent PCS in cmTBI typically are related to problems in attention, working memory, and executive functioning (McInnes, Friesen, MacKenzie, Westwood, & Boe, 2017), yet these injuries frequently go undiagnosed and untreated despite their adverse effect on quality of life. Conventional magnetic resonance imaging (MRI) and computed tomography (CT) are typically insensitive to physiological alterations caused by mild and some moderate TBIs (Bigler & Orrison, 2004; Johnston, Ptito, Chankowsky, & Chen, 2001; Kirkwood, Yeates, & Wilson, 2006), even in individuals with persistent PCS and cognitive deficits. This underscores the need for neuroimaging measures that are sensitive to mild brain injuries and can serve as biomarkers for evaluating interventions aimed at improving functional capacity.

Although the pathophysiology of blast mTBI is not completely understood, diffuse axonal injury (DAI) is known to play a major role in brain dysfunction (Garman et al., 2011), producing an imbalance in excitatory/inhibitory neural activity after mTBI. It is widely assumed that white matter tracts are primarily vulnerable to DAI, which causes cortical network disconnection (see reviews [Asken, DeKosky, Clugston, Jaffee, & Bauer, 2018; Hannawi & Stevens, 2016]). As such, diffusion tensor imaging (DTI) techniques are commonly used to evaluate axonal injury in TBI. However, DTI has limited sensitivity in distinguishing individual patients with mTBI from healthy controls (HCs) (Davenport, Lim, Armstrong, & Sponheim, 2012; Mac Donald et al., 2011; Shenton et al., 2012), rendering clinical applications challenging at this time for measurements solely dependent on detecting white‐matter pathology (Asken et al., 2018; Douglas et al., 2015).

The possibility that other pathologies underlying mTBI are at play was recently suggested by an animal study showing that DAI substantially alters the integrity of gray matter (GM) (Vascak, Jin, Jacobs, & Povlishock, 2018). In this regard, magnetoencephalography (MEG) source imaging can detect subtle pathology that often goes undetected in individuals with mTBI when using structural neuroimaging techniques (Huang et al., 2012; Huang et al., 2014; Robb Swan et al., 2015). MEG directly measures the magnetic signal due to neuronal activation in GM with high spatial localization accuracy (2–3 mm in cortex; Leahy, Mosher, Spencer, Huang, & Lewine, 1998) and high temporal resolution (< 1 ms), which translates into excellent frequency specificity at different frequency bands (Hamalainen, Hari, Ilmoniemi, Knuutila, & Lounasmaa, 1993). Regional resting state (rs‐MEG) slow‐wave (delta band 1–4 Hz, extending to theta band 5–7 Hz) markers are highly sensitive (~85% sensitivity) in distinguishing chronic and sub‐acute mTBI patients with persistent PCS on a single‐subject basis from neurologically intact individuals (Huang et al., 2012; Huang, Nichols, et al., 2014; Lewine et al., 2007; Lewine, Davis, Sloan, Kodituwakku, & Orrison Jr., 1999; Robb Swan et al., 2015). Abnormal slow‐waves have also been reported in a cohort consisting of acute, sub‐acute, and chronic mTBI subjects (Kaltiainen, Helle, Liljeström, Renvall, & Forss, 2018). We recently discovered that rs‐MEG gamma‐band (30–80 Hz) markers showed striking hyperactivity in cmTBI, possibly due to injury of GABA‐ergic parvalbumin‐positive (PV+) interneurons (Huang et al., 2020). In addition, task evoked (te‐MEG) recordings during working memory detected abnormal alpha, beta, and gamma signals throughout the brain in individuals with cmTBI, which correlated with poorer cognitive functioning (Huang et al., 2019). However, the vast number of aberrant markers of neuronal dysfunction throughout the brain poses a significant challenge for clinical applications. In this regard, analytic approaches are needed to discover optimal combinations of aberrant features across different frequency bands that best distinguish cmTBI and also predict neurobehavioral deficits.

Machine learning is a data‐driven approach that optimally integrates high‐dimensional features in large datasets. It has been used to classify mTBI and HC subjects based on rs‐MEG phase‐synchronization analyses of functional connectivity. For example, Vakorin and colleagues applied a support vector machine (SVM) classifier to rs‐MEG source‐based measures of simple functional connectivity, reporting that reduced delta and gamma connectivity together with increased alpha‐band connectivity distinguished civilian mTBI patients from HCs with 88% accuracy (Vakorin et al., 2016). Dimitriadis and colleagues applied an extreme learning machine classifier to rs‐MEG sensor‐based functional connectivity measures derived from graph theory to measure the overall efficiency of information transfer across the brain and at local levels. They reported that alpha‐band local efficiency distinguished civilian mTBI patients from HCs with 100% accuracy (Dimitriadis, Zouridakis, Rezaie, Babajani‐Feremi, & Papanicolaou, 2015). These studies illustrate the great potential of machine learning approaches for uncovering optimal combinations of discriminating features, which can streamline the interpretation of analyses from high‐dimensional data. Despite the very good to excellent classification accuracies reported in these studies, however, functional connectivity metrics do not always easily pinpoint the injured brain regions due to the involvement of multiple regions in functional connectivity analysis and the many interconnecting pathways that connect these regions.

To our knowledge, machine learning has not been used to investigate rs‐MEG optimal features of regional activity across different frequency bands that distinguish cmTBI from HC in participants with similar combat experiences. The main goal of the present study was to develop a deep‐learning (DL) neural‐network algorithm, called 3D‐MEGNET, that could classify cmTBI and HC participants with at least 95% accuracy (sensitivity and specificity). Machine learning is a branch of artificial intelligence that includes methods, or algorithms, for automatically creating models from data. Unlike a system that performs a task by following explicit rules, a machine learning system learns from experience. DL is a subset of machine learning algorithms that use artificial neural network (ANN). DL performs the learning by adjusting the weights of the artificial neurons that are distributed within a large number of “hidden” layers to identify features from the data. Hidden layers come between the input and output layers. Each layer is made up of artificial neurons, often with some activation functions. DL approach has been the most promising for high‐dimensional structural medical imaging processing and classification (Eo et al., 2018; Hammernik et al., 2018; Hyun, Kim, Lee, Lee, & Seo, 2018; Jun et al., 2019; Kwon, Kim, & Park, 2017; Yang et al., 2018).

We investigated the performances of 3D‐MEGNET when modeling individual frequency bands and all bands combined (i.e., all‐Band model) to determine if regional changes in one or more frequency bands optimally distinguished cmTBI and HCs. The classification accuracies of DL based 3D‐MEGNET were also compared with those of SVM, a more a traditional machine learning approach. We also evaluated the behavioral significance of regional rs‐MEG features or classifiers by correlating them with measures of cognitive function.

Our focus was the performance of 3D‐MEGNET with Fast‐VESTAL as the MEG source magnitude imaging method (Huang et al., 2014; Huang et al., 2020) used in the input layer. However, for comparison, performance of 3D‐MEGNET with Beamformer as the source analysis approach for the input layer was also evaluated. Beamformer is a popular spatial filter approach that assumes source time‐courses are uncorrelated (Barnes & Hillebrand, 2003; Gross et al., 2001; Gross & Ioannides, 1999; Hillebrand & Barnes, 2003; Robinson & Vrba, 1999; Sekihara, Nagarajan, Poeppel, Marantz, & Miyashita, 2001; Van Veen, van Drongelen, Yuchtman, & Suzuki, 1997).

MATERIAL AND METHODS

Research participants

All study participants were 101 male U.S. active‐duty military service members or Operation Enduring Freedom / Operation Iraqi Freedom Veterans. Table 1 details demographic characteristics of the participants. The cmTBI group included 59 participants who met the inclusion criteria for a history of combat‐related mTBI with chronic sequelae and persistent PCS, and the age‐matched HC group included 42 individuals with combat experience but without self‐reported combat‐related head injury. Table S.1 lists additional characteristics of research participants that include the length of time since injury and injury type information for the participants with cmTBI. The cmTBI diagnoses were based on Veterans Affairs and Department of Defense diagnostic criteria (The Management of Concussion/mTBI Working Group, 2009): (a) loss of consciousness <30 min or transient confusion, disorientation, or impaired consciousness immediately after the combat‐related trauma; (b) post‐traumatic amnesia <24 hr; and (c) an initial Glasgow Coma Scale (Teasdale & Jennett, 1974) between 13 and 15 if available. Since the Glasgow Coma assessment was not accessible for most individuals who received their injury in theater, volunteers missing the Glasgow Coma Scale assessment, but who met other diagnostic criteria, were also enrolled. The research study clinical interview also assessed 21 enduring PCS (Table 1), modified slightly from the Head Injury Symptom Checklist (McLean, Dikmen, Temkin, Wyler, & Gale, 1984); requirement for inclusion in cmTBI was three or more PCS (McLean et al., 1984). The diagnosis of cmTBI or HC was corroborated using medical records, and research clinical interviews that assessed the nature of patients' injuries. The Supporting Information details exclusion criteria and procedures to prevent medication‐related confounds to rs‐MEG exam.

TABLE 1

Demographic characteristics, neuropsychological scores, and symptoms in the healthy control and cmTBI groups

	HC (n = 42)		cmTBI (n = 59)		Significance
	Mean	SD	Mean	SD
Age	32.00	8.53	29.86	6.31	n.s.
a Years of education	13.71	2.00	12.98	1.41	b n.s.
D‐KEFS
Number‐letter sequencing	11.05	2.00	10.17	2.31	p = .044
Letter fluency	11.14	3.66	10.09	2.82	n.s.
Category fluency	11.90	3.45	11.39	2.86	n.s.
Category switching	11.55	2.43	11.01	2.93	n.s.
WAIS
Digit symbol coding	10.50	2.70	9.77	2.94	n.s.
Processing speed index	105.43	16.29	102.95	14.53	n.s.

Note: D‐KEFS refers to the Delis‐Kaplan Executive Function System. WAIS refers to the Wechsler Adult Intelligence Scale‐Third Edition. Neuropsychological measures are scaled scores (mean = 10, SD = 3), except for The WAIS Processing Speed Index, which is a standard score.

Years of education: High school = 12; AA = 14; bachelor's = 16; master's = 18; JD = 19; MD, DO, or ND = 20; PhD = 21; MD–PhD = 25.

Mann–Whitney U test.

Among the 101 research participants, 95 (including 55 cmTBI and 40 HC) were used to study the sensitivity and specificity of the 3D‐MEGNET DL approach for group classification. The remaining six “out‐of‐sample” participants (four cmTBI and two HC) were used to validate the performance accuracy of the DL rs‐MEG approach.

Neuropsychological exams

The neuropsychological assessment (Table 1) focused on tests of executive functions from the Delis‐Kaplan Executive Function System (D‐KEFS; Delis, Kaplan, & Kramer, 2001), and visuospatial processing and psychomotor speed from the Wechsler Adult Intelligence Scale‐Third Edition (WAIS‐III; Wechsler, 1997; Wechsler, 2008), which are sensitive to cognitive decline in mTBI (Robb Swan et al., 2015). The neurocognitive correlation analyses presented herein focused on the following subset of measures from a larger battery that was administered, namely those showing significant correlations with MEG measures in mTBI and/or control groups in previous studies (Huang et al., 2019; Huang et al., 2020; Robb Swan et al., 2015). The D‐KEFS Number‐Letter Switching subtest of the Trail Making Test measures cognitive flexibility. Subtests of the D‐KEFS Verbal Fluency Test require participants to generate, as quickly as possible, words beginning with specific letters (Letter Fluency subtest) and words in specific semantic categories while shifting between categories (Category Switching subtest). Letter Fluency tests phonemic processing, which is sensitive to language difficulties, whereas Category Switching tests cognitive flexibility during semantic processing. The WAIS‐III Digit Symbol Coding subtest utilizes a number of cognitive processes and is primarily used as a measure of visuospatial processing and psychomotor speed (Wechsler, 2008). The WAIS‐III Processing Speed Index is a composite standard score derived from the Symbol Search and Digit Symbol Coding subscales.

MEG data acquisition and signal pre‐processing

Resting‐state MEG data were collected at the UCSD MEG Center using the VectorView™ whole‐head MEG system (Elekta‐Neuromag, Helsinki, Finland) with 306 MEG channels. Participants sat inside a multi‐layer magnetically shielded room (Cohen, Schlapfer, Ahlfors, Hamalainen, & Halgren, 2002). MEG recordings were divided into two 5‐min blocks where the participant was instructed to keep his/her eyes closed and empty his/her mind. Data were sampled at 1000 Hz and run through a high‐pass filter with a 0.1 Hz cut‐off, and a low‐pass filter with a 330 Hz cut‐off. Micro eye blinks and eye movements were monitored using two pairs of bipolar electrodes, and heart signals were monitored with another pair of bipolar electrodes.

Substantial efforts were taken to help ensure that participants were alert during the rs‐MEG recordings. Prior to MEG sessions, participants completed a questionnaire about the number of hours they slept the previous night, how rested they felt, and if there was any reason that they might not be attentive and perform to the best of their abilities (e.g., headache, pain). Sessions alternated between eyes‐closed and eyes‐open conditions, and eye blinking and movement were monitored. During MEG recording, participants were viewed on camera and technicians continuously monitored alpha band oscillations, which are consistently associated with tonic alertness (Oken, Salinsky, & Elsas, 2006).

The Supporting Information describes the pre‐processing steps for rs‐MEG and structural MRI data. None of the cmTBI or HC participants had observable deficits on structural MRI as assessed by a board‐certified neuroradiologist (Dr. R. R. Lee).

MEG source magnitude imaging using Fast‐VESTAL

Voxel‐wise MEG source magnitude images were obtained using our high‐resolution Fast‐VESTAL MEG source magnitude imaging method (Huang et al., 2020; Huang, Huang, et al., 2014). This approach requires the sensor waveform covariance matrix. The artifact‐free, eyes‐closed, resting‐state MEG sensor‐waveform datasets were divided into 2.5 s sections. The data in each section were first DC‐corrected and then run through one of the following band‐pass filters for delta‐theta (1–7 Hz), alpha (8–12 Hz), beta (15–30 Hz), and gamma (30–80 Hz) bands. Notch filter at 60 Hz was applied to remove the power line signals.

Waveforms from all 306 sensors were used in the analysis. Sensor‐waveform covariance matrices were calculated in the time domain for individual sections after the band‐pass filtering, in the same way as in our previous studies to preserve the phase information (Huang et al., 2017; Huang, Huang, et al., 2014). Then the final sensor‐waveform covariance matrix was obtained by averaging the covariance matrices across individual sections for the concatenated 10‐min resting‐state data that combined the two 5‐min blocks. From the covariance matrix, whole‐brain MEG source magnitude images for each frequency band were obtained for each participant using the Fast‐VESTAL procedure (Huang et al., 2016; Huang et al., 2020; Huang, Huang, et al., 2014).

3D‐MEGNET DL neural network for classification

Figure 1 diagrams the DL network based on 3D‐MEGNET. In 95 participants and for each frequency‐band model, voxel‐wise whole brain rs‐MEG source images obtained from Fast‐VESTAL were first spatially registered to the MNI‐152 (Grabner et al., 2006) brain‐atlas template using FLIRT, an affine transformation program from FSL software (www.fmrib.ox.ac.uk/fsl/; Smith et al., 2004; Woolrich et al., 2009).

FIGURE 1

3D‐MEGNET deep‐learning diagram. The input images include rs‐MEG source imaging volumes in standard MNI‐152 space across different frequency bands and participants. The rs‐MEG source images are first convolved with a Gaussian kernel. Next, the maximum values from individual functional ROIs are pooled to form the Pooling Layer. The elements in the Pooling Later are subject to Recursive Feature Elimination and then reshaped into the Flatten Layer. During the classification section, two fully connected Dense Layers with ReLU activation function are added. One of these layers is then fully connected with the Dense Output Layer in which SoftMax activation function is used to classify the individuals into either Control or mTBI groups

Spatial convolution and Max‐pooling using functional regions of interest (ROI): Once in MNI‐152 space, the rs‐MEG source imaging data from all frequency bands were first run through a spatial convolution layer. This convolution operation convolved the imaging data with a 3D Gaussian kernel with 5.0 mm full width half maximum (FWHM) to further reduce the inter‐subject variability in anatomy. Then, a Max‐pooling procedure was applied to the spatially convolved rs‐MEG source imaging data. In this procedure, the MEG source imaging voxels were grouped into 184 GM functional ROI variables using the FCONN parcellation with similar sizes (Shen, Tokoglu, Papademetris, & Constable, 2013). In each functional ROI, the maximum activity was obtained across all voxels within such an ROI. After regressing out age and education variables, 184 ROI features for each frequency band were obtained. Additional details are provided in the Supporting Information.

The spatial convolution and Max‐pooling procedure based on functional ROI played important roles in feature selection and dimension reduction, which is similar to the steps in 3D Convolutional Neural Networks (3D‐CNN) for imaging processing (e.g., [Ji, Xu, Yang, & Yu, 2013]), except for two main differences. The first main difference is spatial convolution uses only the Gaussian kernel in 3D‐MEGNET, without any spatial‐gradient based convolutional filters. The second main difference is that the Max‐pooling is based on functional ROI in 3D‐MEGNET, rather than on Euclidean distance in 3D‐CNN. These differences are discussed below.

Only one convolutional filter with Gaussian smoothing kernel was used in 3D‐MEGNET. This is one main difference between 3D‐MEGNET and the typical CNN approaches that use a large number of spatial gradient‐based convolutional filters for processing structural images (Eo et al., 2018; Hammernik et al., 2018; Hyun et al., 2018; Jun et al., 2019; Kwon et al., 2017; Yang et al., 2018). Specifically, the built‐in CNN filters from TensorFlow were not used in the present study. The gradient based convolutional filters are good for detecting edges and shapes in images, and they usually require a large number (e.g., thousands) of samples to train. However, rs‐MEG source magnitude images only contain “hot spots” without edge, shape, and other features. The single layer of convolutional filter with a Gaussian smoothing kernel in 3D‐MEGNET, in combination with our Max‐pooling layer (see below) using 184 functional ROIs, effectively achieved the goal of dimension reduction using a smaller number of samples to train.

The Max‐pooling procedure based on functional ROIs (Shen et al., 2013) was used in 3D‐MEGNET, instead of the geometric Max‐pooling procedure based on Euclidean distance in typical CNN (Eo et al., 2018; Hammernik et al., 2018; Hyun et al., 2018; Jun et al., 2019; Kwon et al., 2017; Yang et al., 2018). It is common that functionally distinct regions are next to each other with small Euclidean distance. The functional ROI based Max‐pooling has the advantage of taking into consideration the boundaries of functionally distinct regions.

Recursive feature elimination (RFE) based on SVM (Lu et al., 2016; Muller & Guido, 2016; Wang, Xiao, & Wu, 2019) was used to obtain the optimal ROI features from the Max‐pooling ROI variables across all frequency bands. The source activity from the subset of ROI variables after the RFE were reshaped into a 1‐dimensial Flatten Layer for classifications (Figure 1). All elements in the Flatten Layer were fully connected to two hidden Dense ReLU Layers. Each layer contained artificial neurons (i.e., nodes) that were connected to all the nodes in the preceding layer, with the rectified linear unit (ReLU) activation function. These two “hidden” ReLU layers combined all the features learned by the previous layers across the rs‐MEG ROI source images to optimize the weights of the features for classification. The third fully connected dense layer containing nodes with Softmax activation function led to the final classification the MEG source images. The output classes as shown in Figure 1 were the cmTBI and HC groups. More details of the 3D‐MEGNET setting were provided in Supporting Information.

The performances of 3D‐MEGNET were examined for each individual frequency band, by changing the rs‐MEG source imaging data in the Input Layer from the combined all frequency bands data to that from each frequency band. The rest of the network design remained the same.

3D‐MEGNET procedure and parameters setting

In the present study, 85 (89.5%) individuals were randomly selected for a training data set, and the remaining 10 individuals (10.5%) formed the testing set. Such setting is somewhat different from the typical DL settings of 80–20% splitting ratio between the training and testing data sets (Muller & Guido, 2016). The main reason of our choice of splitting ratio was the moderate number of individuals in the present study (see Section 4).

RFE based on SVM (Lu et al., 2016; Muller & Guido, 2016; Wang et al., 2019) was used to obtain the optimal ROI features from the Max‐pooling ROI variables across all frequency bands. Initially, all K variables present in our ROI dataset were included in training the network model. Then, the performance of the RFE‐SVM model was calculated. Next, we computed the performance of the model after eliminating each variable (K times). Specifically, we dropped one ROI variable every time and trained the model on the remaining K‐1 variables. The ROI variable whose removal had produced the smallest (or no) change in the performance of the model was dropped. This process was repeated until no ROI variable could be eliminated. The source activity from the subset of ROI variables after the RFE were reshaped into a 1‐dimensial Flatten Layer for classifications (Figure 1).

All elements in the Flatten Layer were fully connected to a Dense Layer with 2,048 artificial neurons (i.e., nodes) for the all‐band model or 4,096 nodes for individual frequency band models. A fully connected dense layer contained nodes that were connected to all the neurons in the preceding layer. Batch normalization was used in this layer. Using the rectified linear unit (ReLU) activation function, artificial neurons in the Dense ReLU Layer nonlinearly transformed the activations to speed up the network training and reduce the sensitivity to network initialization. A dropout rate of 0.15 was used in this hidden layer to prevent overfitting. Then, the outputs were fed to the second ReLU Dense Layer with 64 nodes. These two “hidden” ReLU layers combined all the features learned by the previous layers across the rs‐MEG ROI source images to optimize the weights of the features for classification.

The third fully connected dense layer led to the final classification the MEG source images. This Dense Softmax Output Layer used softmax activation function to normalize the output of the second fully connected Dense ReLU Layer. The output of the Softmax Layer consisted of positive numbers that sum to one, which was be used as classification probabilities to assign the individuals to one of the mutually exclusive classes and compute the loss. The output classes in Figure 1 are the cmTBI and HC groups.

Sensitivity and specificity of 3D‐MEGNET for classification

The pipeline for studying training and testing accuracies of the 3D‐MEGNET was as follows. First, the rs‐MEG data from 95 participants in the two groups (55 cmTBI and 40 HC) were randomly split into a training set of 85 individuals (89.5%) and a testing set of 10 individuals (10.5%). The training set was then used to train the 3D‐MEGNET model, and the training accuracy was examined. To control for label imbalance, the number of subjects in any group (bmTBI or HC) was at least 30% in each testing set. After the network was trained, the testing data set, which did not overlap with the training set, was used to independently assess the testing accuracy of the network. The testing accuracies were divided into: (a) Sensitivity = TP / P, where TP is the true positive (i.e., correctly identified mTBI subjects) and P is number of real mTBI subjects in the data; and (b) Specificity = TN/N, where TN is the true negative (i.e., correctly identified HC subjects) and N is the number of real HC subjects in the data. To examine the robustness of the performances of 3D‐MEGNET, the procedures in Figure 1 were repeated for 1,000 different combinations of the 85–10 splits.

The 95 participants were randomly split to 85 in‐training datasets for training 3D‐MEGNET, and then the trained network was applied to the remaining 10 subjects' in‐testing datasets to assess classification accuracies in each iteration. This process was repeated 1,000 times, with each iteration using a different combination of the 85–10 split, such that each participant had the chance of being randomly partitioned into the testing dataset ~100 times (i.e., 10% × 1,000), thereby forming a pool of ~100 samples of training results per participant. To assess the sensitivity, Monte‐Carlo analyses were performed for the 55 cmTBI participants by randomly sampling each individual's pool of ~100 training results, which contained categorical classifications into a group and probability estimates for group membership. The Monte‐Carlo procedure was repeated 1,000 times to reliably calculate classification accuracy or sensitivity for the cmTBI group. The same process was then carried out for the 40 HC participants to calculate specificity.

Training 3D‐MEGNET

The network was trained with Adam (Kingma & Ba, 2017) as the optimizer with a fixed learning rate of α = 0.001, β = 0.9, β = 0.999. Adam is computationally efficient with less memory demand compared to other optimizers (e.g., stochastic gradient optimization). The number of training iterations, in each of which all training data goes through once, is called an epoch. The 3D‐MEGNET was trained for 100 epochs with a batch size the same as the size of the training data set (i.e., 85). Categorical cross‐entropy was used as the loss function. The final model was trained when the epoch number reached 100 or when there was no improvement in training loss for five consecutive epochs, whichever was less. The hardware and software settings are listed in Supporting Information.

ROIs that contributed to the 3D‐MEGNET classifications

Additional analyses were conducted for features that were chosen by the RFE‐SVM process in the Flatten Layer of 3D‐MEGNET. Post hoc ANOVAs tested for statistically significant group differences in the rs‐MEG magnitude of each of these ROIs for a specific frequency band. Pearson correlations analyses also tested for relationships between rs‐MEG activity in each ROI and neuropsychological test performances separately for each group of participants.

3D‐MEGNET versus SVM, Fast‐VESTAL versus Beamformer

To evaluate the performance of DL based 3D‐MEGNET over more a traditional machine learning approach, we examined the classification accuracy of the SVM approach. In the SVM approach, the two hidden layers in the diagrams of Figure 1 were replaced with the SVM, and the rest of the diagrams remained the same. Here, the comparisons were between two different machine learning approaches (i.e., SVM vs. 3D‐MEGNET), based on the same Fast‐VESTAL source magnitude data input from the 55 cmTBI and 40 HC participants.

In addition, we examined the performances of 3D‐MEGNET using the source imaging data from the Beamformer source analysis as the input data, rather than using Fast‐VESTAL source magnitude images as the input. Here, the machine learning structure did not change (i.e., 3D‐MEGNET), and the comparisons were between two different source imaging approaches for the input layer (i.e., Beamformer vs. Fast‐VESTAL) in the same 55 cmTBI and 40 HC participants. Specifically, the Fast‐VESTAL input data were replaced by those from the Beamformer source analysis approach, and the rest of the diagrams in Figure 1 remained essentially the same (see Section S.1.5 in Supporting Information for details).

Further validation of using Out‐of‐Sample data sets

The performance of 3D‐MEGNET and SVM algorithms was further validated using rs‐MEG data from six “out‐of‐sample” participants (four cmTBI and two HC). In these analyses, the rs‐MEG data from the 95 participants were used to train the machine‐learning algorithms using Fast‐VESTAL source magnitude images as the input. Then the output from this trained algorithm was tested on the six out‐of‐sample participants for classification accuracy.

RESULTS

Classification accuracies and ROC analyses of 3D‐MEGNET

The classification accuracies of 3D‐MEGNET in the testing sets are shown in upper panel of Table 2 and Figure 2a,b, for each model. Fast‐VESTAL source magnitude imaging data were used as the input layer. Using the diagnostic classification outcome measure (cmTBI vs. HC), the all‐band model had excellent classification accuracies for both groups for sensitivity (99.92%) and specificity (98.95%). Classification accuracies for the all‐band model were markedly better than for models using individual frequency bands (Mann–Whitney U test: p < .0001). Figure 2a graphs the classification accuracies of the individual frequency band models for Fast‐VESTAL input, showing that the rank order of performance from high‐to‐low in both cmTBI and HC groups was gamma > delta‐theta > alpha ~ beta.

TABLE 2

Upper panel: In 3D‐MEGNET analysis of Fast‐VESTAL data input, the categorical classification accuracies (sensitivity and specificity) and the Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) curve in testing data sets for different frequency‐band models. Lower Panel: The corresponding values in SVM analysis of Fast‐VESTAL data input

DL Fast‐VESTAL	All bands (%)	δ‐θ band (%)	α band (%)	β band (%)	γ band (%)
55 mTBI (sensitivity)	99.92 ± 0.38	89.99 ± 2.34	85.56 ± 2.77	85.48 ± 2.42	93.17 ± 2.28
40 HC (specificity)	98.95 ± 1.54	84.00 ± 3.62	79.39 ± 3.17	78.34 ± 3.63	84.40 ± 3.54
AUC of ROC (%)	99.99	93.54	88.77	87.26	94.35
95% CI for AUC	99.90, 100.00	91.77, 95.18	86.59, 90.95	85.00, 89.59	92.64, 96.14

Note: The All‐band Model (bold values) was the main focus of this study.

FIGURE 2

Upper panel: Classification results from the five 3D MEGNET models with Fast‐VESTAL input. (a): 3D‐MEGNET's Percent accuracy in testing data sets in mTBI and control groups, plotted for all frequency bands combined (i.e., all‐band model, green bars), and for individual frequency bands separately. The inverted “U” shapes indicate the all‐band model was statistically more accurate than each of the individual frequency band models at p < .0001 (Mann–Whitney U test). (b): Operating characteristic (ROC) curves for all frequency bands combined (i.e., all‐band model, green curve), and for individual frequency bands separately. The True Positive Rate (i.e., Sensitivity on y‐axis) was plotted as the function of False Positive Rate (i.e., 1—Specificity on x‐axis). The dashed line represents a naive / non‐discretionary classifier. Probabilistic Classification data were used to calculate the ROC curves. Lower panel: Corresponding classification results for SVM models with Fast‐VESTAL input. (c) Percent accuracy in testing data sets for the SVM approach. (d) ROC curves for the SVM approach

The receiver operating curve (ROC) analyses evaluated the overall accuracy of models by analyzing the area under the curve (AUC) for the sensitivity and specificity distributions. Figure 2b shows the ROC curves for all Fast‐VESTAL input models. All curves markedly outperformed the naive / nondiscretionary classifier (diagonal dashed line). However, the all‐band model (green curve) showed nearly perfect performance with both high sensitivity and high specificity. At the level of 0.8 True Positive Rate (sensitivity), the performance order for the False Positive Rate from low to high (lower is better) was all‐band < gamma < delta‐theta < alpha ~ beta. The AUC (Table 2) for all models was greater than 80%, with near perfect performance for the all‐band model at 99.99%, followed by gamma > delta‐theta > alpha ~ beta. Table S.2 details the number of PCS symptoms for each individual participant and the probability of correctly classifying individual participants into their diagnostic group for each of the rs‐MEG models.

By comparison, when Beamformer input was used as the source imaging approach, classification accuracies for each frequency band (i.e., sensitivity and specificity) were significantly lower than those corresponding to models using Fast‐VESTAL input (Section S.2 of the Supporting Information; Figure S1 and Table S5).

Classification accuracies and ROC analyses of SVM for Fast‐VESTAL input

Next we compared the performance of DL‐based 3D‐MEGNET over the more traditional SVM learning approach. The lower panel of Table 2 and Figure 2c,d show the classification accuracies of SVM models in the testing sets. The SVM approach used the same Fast‐VESTAL source magnitude data input from the same 55 cmTBI and 40 HC participants. As in the DL‐based 3D‐MEGNET, classification accuracies with SVM for the all‐band model (i.e., sensitivity at 90.04%, specificity at 89.88%, Figure 2c) were significantly higher than for models using individual frequency bands (Mann–Whitney U test: p < .0001). However, for each frequency band, the classification accuracies (i.e., sensitivity and specificity values) with SVM were significantly lower than those corresponding ones with 3D‐MEGNET (Mann–Whitney U test: p < .0001).

Figure 2d shows the ROC curves for all models using SVM with Fast‐VESTAL input. All curves markedly outperformed the naive / nondiscretionary classifier (diagonal dashed line). However, each of these ROC curves underperformed the corresponding one for 3D‐MEGNET (Figure 2b).

ROIs contributing to All‐Band classifications

The spatial distribution of brain regions that contributed to the diagnostic classification in the all‐band model is displayed in Figure 3 (upper panel) for each frequency band. The numbers of contributing classifiers were 28 (delta‐theta), 29 (alpha), 27 (beta), and 23 (gamma), respectively. The total number of contributing classifiers was 107 from 83 independent ROIs, with some regions identified by multiple bands. Post hoc ANOVAs were performed to explore if 3D‐MEGNET mainly included ROIs that showed significant group differences based on univariate analyses. As these analyses were exploratory, p values were uncorrected for multiple comparisons. The results revealed regions for which activity significantly increased (red), decreased (blue), or showed no difference (green) in the cmTBI relative to the HC group. Some ROIs that contributed to the all‐band classification model showed significantly increased delta‐theta (7.5%) and gamma band (13%) activity in the cmTBI group. Several regions also showed significantly increased or decreased alpha (12%) and beta (6.5%) activity in the cmTBI group. However, across frequency bands, activity in the majority of ROI classifiers (61%) did not differ between the groups. Table S.1 lists anatomic regions of the ROIs that contributed to the all‐band model.

FIGURE 3

Classifiers for the all‐band model with Fast‐VESTAL input. For the all‐band model, the ROIs that contributed to the mTBI‐HC classification. Red: the ROIs contributing to the classification that also showed significant increases (p < .05) in rs‐MEG activity in mTBI than HC; Blue: the ROIs contributing to the classification that also showed significant decreases in rs‐MEG activity in mTBI than HC; Green: the ROIs contributing to the classification but did not show significant group differences in rs‐MEG activity. (a)–(l): locations of representative ROIs are indicated as yellow arrows. The scatter plots of these ROIs are presented in Figure 4

FIGURE 4

Correlations between 3D‐MEGNET classifiers and neuropsychological test performances with Fast‐VESTAL input. First two columns: scatter plots from representative ROIs with significant group differences in rs‐MEG activity that also showed significant correlations with neuropsychological exams in the mTBI (red stars), but not in HC (blue circles) groups. Third column: ROIs without significant group differences in rs‐MEG activity that showed significant correlations with neuropsychological exams in mTBI. In (f) and (l), but not (c) or (i), rs‐MEG in HC also showed significant correlations with neuropsychological exams. The first, second, third, and fourth rows show associations between delta‐theta, alpha, beta, and gamma activity, respectively, and cognitive performances. In (b), (c), and (e), Spearman correlations were used. The locations of the ROIs are indicated by yellow arrows in Figure 3

Correlations between rs‐MEG activity and cognition

Figure 4 displays scatter plots of representative correlations between classifiers and performances on neuropsychological tests of executive functioning, separately for each group. In Figure 3, arrows highlight the anatomy of each classifier shown on scatter plots in Figure 4. The normality of the neuropsychological variables was checked using the Shapiro–Wilk and Shapiro–Francia normality test. All neuropsychological variables passed the normality test, except Number‐Letter Sequencing. Thus, non‐parametric Spearman correlations were used in the analyses that involved Number‐Letter Sequencing. For all the other analyses, Pearson correlations were used. Correlations were unadjusted for familywise error rates, as they were exploratory. As such, they should be interpreted cautiously. We first investigated neurocognitive associations for classifiers that showed significant group differences in rs‐MEG activity (Figure 4, first two columns). Then neurocognitive correlations were evaluated for classifiers that did not show significant group differences in rs‐MEG activity (Figure 4, third column). Significant correlations between other classifiers and executive functioning are detailed in Table S.3. With only a few exceptions, most neurocognitive associations were significant in the cmTBI group, but not in the HC group. Classifiers that correlated with executive functioning were largely in the frontal cortex (Figure 4; Table S.3).

For regions showing abnormal delta‐theta hypoactivity in the cmTBI group (Figure 4, first row), right dorsolateral prefrontal cortex (DLPFC, Figure 3a) and left ventral medial precentral cortex (vmPrecentral, Figure 3b) activity respectively correlated negatively with Category Fluency and Number‐Letter Sequencing scores only in the cmTBI group (Figure 4a,b). Although delta‐theta activity in the right fusiform gyrus did not differ between the groups (Figure 3c), it negatively correlated with Number‐Letter Sequencing scores in only the cmTBI group (Figure 4c). Increased delta‐theta band activity in the above regions predicted worse cognitive functioning.

For regions showing abnormal alpha hypoactivity in cmTBI (Figure 4, second row), right dorsal anterior cingulate cortex (dACC, Figure 3d) and right ventromedial prefrontal cortex (VMPFC, Figure 3e) activity respectively correlated positively with Category Switching and Number‐Letter Sequencing scores only in the cmTBI group (Figure 4d,e). Although there were no group differences in left putamen alpha activity (Figure 3f), it positively correlated with Digit Symbol Coding scores in both groups (Figure 4f). Alpha band decreases in the above regions predicted worse cognitive functioning.

As for regions showing abnormal beta hypoactivity in the cmTBI group (Figure 4, third row), left VMPFC (Figure 3g) and right rostral anterior cingulate (rACC, Figure 3h) beta activity respectively correlated positively with Letter Fluency and Processing Speed scores only in the cmTBI group (Figure 4g,h). Despite an absence of group differences in left precuneus beta activity (Figure 3i), it positively correlated with Category Fluency scores in the cmTBI, but not in the HC group (Figure 4i). Beta band decreases in the above regions predicted worse cognitive functioning.

For regions showing abnormal gamma hyperactivity in the cmTBI group (Figure 4, fourth row), left frontal pole (FP, Figure 3j) and left ventrolateral prefrontal cortex (VLPFC, Figure 3k) activity respectively correlated negatively with Letter Fluency Digit Symbol Coding scores only in the cmTBI group (Figure 4j,k). Although the groups did not differ in right posterior superior frontal gyrus (pSFG, Figure 3l) gamma activity, it negatively correlated with Category Switching scores in both groups (Figure 4l). Increased gamma band activity in these regions was associated with worse cognitive functioning.

Validation results on the “Out‐of‐Sample” data sets

To further validate the performance of 3D‐MEGNET and SVM algorithms, rs‐MEG Fast‐VESTAL source magnitude imaging data from six out‐of‐sample participants (four cmTBI and two HC) were used as testing data, after the machine‐learning algorithms were trained using the 95 participant cohort. The results show that 3D‐MEGNET (all‐band model) accurately classified all 6 “out‐of‐sample” participants, whereas SVM correctly classified five out of the six participants with one cmTBI participant misclassified as a HC.

DISCUSSION

The DL 3D‐MEGNET algorithm demonstrated excellent sensitivity (>99%), specificity (>98%), and overall accuracy (>99%) in identifying optimal rs‐MEG source imaging features that distinguished Veterans and active‐duty service members with diagnosed cmTBI from those without based on Veterans Affairs and Department of Defense diagnostic criteria. The DL 3D‐MEGNET algorithm also outperformed the more traditional SVM machine‐learning approach. Furthermore, within 3D‐MEGNET, the input layer using the Fast‐VESTAL source magnitude imaging approach outperformed the Beamformer source analysis approach. Modeling of all regional frequency bands together also rendered significantly superior classifications relative to models that included only a single frequency band. This finding underscores the diversity in abnormal regional oscillatory frequencies in cmTBI. We also discovered that optimal classifiers in the all‐band model were from brain regions that did and did not show group differences in activity. The latter finding likely relates to the diversity of brain injuries that individuals sustain, some of which do not characterize the cmTBI participants as a group, but bolster the discrimination of individual participants from HCs. For both types of classifiers, regional activity was cognitively relevant, with hyperactivity from delta‐theta and gamma bands and hypoactivity from alpha and beta bands correlating with poorer executive functioning and visuospatial processing/psychomotor speed, typically only in cmTBI individuals. However, these neurocognitive associations require validation in future studies of larger samples of cmTBI cohorts.

Source amplitude imaging for diagnostic classification

3D‐MEGNET modeling of source amplitude images directly revealed specific brain regions that contribute to mTBI classification, identifying features for which activity did and did not differ significantly between groups in univariate analyses. This finding highlights the power of 3D‐MEGNET, particularly for diagnostic classification of conditions in which there are considerable individual differences in underlying pathology. Our approach substantially differs from machine learning methods based on rs‐MEG functional connectivity, which have been used with very good success in classifying mTBI (Dimitriadis et al., 2015; Vakorin et al., 2016) and mild cognitive impairment (Dimitriadis et al., 2018; Yang, Bornot, Wong‐Lin, & Prasad, 2019). Functional connectivity is a powerful tool and well suited for characterizing pathology at the level of coordinated networks. However, the interpretation of abnormal functional connectivity can be challenging in terms of pinpointing the source(s) of injury, which may lie in only one of two identified regions, both regions, and/or direct and indirect pathway(s) connecting the regions. Future DL studies that combine abnormalities from source magnitude and functional connectivity measures may advance a more complete understanding of regional and network features that optimally characterize neuronal pathology in cmTBI.

Source imaging based classification also differs from the sensor based approach (Dimitriadis et al., 2015), wherein the sensor waveform cannot be directly linked with specific brain regions, since one MEG sensor can receive signals from multiple brain regions and signals from any brain region can be detected by multiple MEG sensors. Although our study and Vakorin and colleagues (28) both used MEG source imaging, 3D‐MEGNET modeling in our study used high resolution Fast‐VESTAL (Huang, Huang, et al., 2014), which is highly accurate in recovering source time‐course activity (Huang, Huang, et al., 2014). This differs from Beamformer source imaging (Vakorin et al., 2016), which can lead to signal leakage due to its assumption that brain activity is uncorrelated (Huang, Huang, et al., 2014). This is likely the reason that the Beamformer models in the present study underperformed those for Fast‐VESTAL (Table 2 and Figure 2), since we believe the correlated source activity may also contribute to the group classifications between cmTBI and HC cohorts. Despite these differences, it is clear from our study and others that rs‐MEG exhibits high sensitivity and specificity in detecting mild brain injuries (Dimitriadis et al., 2015; Huang et al., 2012; Huang, Nichols, et al., 2014; Vakorin et al., 2016). Moving forward, it will be critical to determine the degree to which classifiers from different DL algorithms can reliably distinguish different cohorts of patients with mTBI or cmTBI from HC.

Challenges of combat‐related mTBI

Compared with the classification of non‐combat civilian mTBI and HC subjects (Dimitriadis et al., 2015; Vakorin et al., 2016), differentiating cmTBI from HC individuals with similar combat experiences is more challenging. This is because cmTBI and HC combatants have both experienced exposure to high levels of stress and potential exposure to blasts. Deployment‐related head injury exposures and symptoms are typically available only through self‐reports, which are prone to the vagaries of memory (Alosco et al., 2016). Moreover, PCS symptoms are non‐specific and overlap with diagnostic criteria of other medical conditions, increasing the challenge in making accurate, universally agreed upon cmTBI diagnoses (Radigan, McGlinchey, Milberg, & Fortier, 2018). Active duty and Veteran participants referred to this study were assigned to a diagnostic category, either cmTBI or HC, based on the current gold standard, Veterans Affairs and Department of Defense diagnostic criteria defined as head injury exposure accompanied by altered mental state, post‐trauma amnesia or loss of consciousness, and the reporting of at least three current PCS (The Management of Concussion/mTBI Working Group, 2009). Despite the challenges inherent in cmTBI diagnosis, we showed for the first time that 3D‐MEGNET DL applied to rs‐MEG source magnitude imaging can distinguish cmTBI individuals from combat‐deployed HCs with excellent sensitivity and specificity.

Diffusion axonal injury in gray matter, mechanisms of Delta‐theta and gamma activity

DAI plays a major role in mTBI. Traditionally and intuitively, it is assumed that white‐matter tracts and brain regions sensitive to tissue shearing are primarily vulnerable to DAI, which produces cortical network disconnection (see reviews in [Asken et al., 2018; Hannawi & Stevens, 2016]). Yet even sophisticated diffusion‐based MRI techniques for detecting white‐matter abnormalities in mTBI are not sufficiently sensitive for meaningful clinical applications (Asken et al., 2018; Douglas et al., 2015). Moreover, recent animal studies challenge this view by showing that GM is at least and potentially more vulnerable to DAI (see references in (Vascak et al., 2018) and discussion below). As such, rs‐MEG may play a particularly important role in revealing regional gray‐matter vulnerabilities in cmTBI, beyond the traditional focuses on white‐matter tracts and regions sensitive to shearing injury.

Interestingly, gamma and delta‐theta band activity were better classifiers than alpha or beta bands, when modeled individually. Abnormal increases in gamma and delta‐theta activity are consistent with previous rs‐MEG studies of mTBI. For example, aberrantly increased rs‐MEG slow‐wave markers are sensitive in distinguishing mTBI patients with persistent PCS from neurologically intact individuals (Huang et al., 2012; Huang, Nichols, et al., 2014; Lewine et al., 1999; Lewine et al., 2007; Robb Swan et al., 2015). Neurophysiological studies in animals have established a solid connection between pathological delta‐wave (1–4 Hz) generation in GM and injuries. Polymorphic delta‐waves produced by physical lesions in cats were localized to the GM of cortex overlying the lesion (Ball, Gloor, & Schaul, 1977; Gloor, Ball, & Schaul, 1977). Abnormal delta‐waves can also be induced by the administration of atropine (Schaul, Gloor, Ball, & Gotman, 1978), a competitive antagonist of acetylcholine (ACh) receptors that can block and/or limit ACh. These experiments concluded that cortical deafferentation was a key factor in abnormal delta‐wave production, owing to blockages/limitations in the cholinergic pathway (Schaul, 1998). Delta wave generation may be a protective mechanism that isolates injured regions from the rest of the brain for neural healing (Naviaux, 2019). Delta wave activity is a normal feature of deep (“slow‐wave”) sleep in which parasympathetic autonomic tone is high and associated with metabolic patterns that underpin healing and repair in both the CNS and periphery (Javaheri & Redline, 2012; Naji, Krishnan, McDevitt, Bazhenov, & Mednick, 2019; Xie et al., 2013; Yüzgeç, Prsa, Zimmermann, & Huber, 2018). In the periphery, slow‐wave sleep is required for normal protective and restorative decreases in blood pressure and the sleep‐associated increase in healthy heart rate variability (Boudreau, Yeh, Dumont, & Boivin, 2013; Javaheri & Redline, 2012). In genetic forms of mitochondrial disease that preferentially affect neuronal function like Alpers syndrome, CNS healing efforts are also associated with increased delta wave activity that persists during wakefulness (Naviaux et al., 1999).

Increased gamma band activity in the cmTBI cohort is consistent with our recent report of striking rs‐MEG gamma‐band hyperactivity in mTBI (Huang et al., 2020), which we suggested may be a consequence of injury to GABA‐ergic parvalbumin‐positive (PV+) interneurons. Specifically, DAI plays a major role in brain dysfunction in mTBI because it produces an imbalance in excitatory/inhibitory neural activity. In GM, DAI is also directly associated with injuries to GABA‐ergic inhibitory interneurons, specifically the parvalbumin‐positive (PV+) interneurons (Vascak et al., 2018). DAI in GM is perisomatic, near the soma of the PV+ interneurons (Vascak et al., 2018), or with a degradation of the perineuronal net (PNN), a specialized extracellular structure enwrapping cortical PV+ inhibitory interneurons (Hsieh et al., 2017). Fast‐spiking (FS) PV+ inhibitory interneurons are the most common type of GABA‐ergic cells that express the calcium‐binding protein PV+ and receive N‐methyl‐D‐aspartate (NMDA)‐dependent excitatory input from pyramidal cells (Carlén et al., 2012). FS‐PV+ interneurons regulate the activity of neural networks through GABA‐ergic inhibition of local excitatory neurons, and synchronous activity of FS‐PV+ interneurons generates gamma oscillations (30–80 Hz; Carlén et al., 2012). Dysfunction or injury to PV+ interneurons in animals causes disinhibition in the neural network, which upregulates spontaneous gamma activity (and maybe beta), owing to a lack of inhibition to pyramidal and other excitatory neurons (Carlén et al., 2012; Cho et al., 2015; Kalemaki, Konstantoudaki, Tivodar, Sidiropoulou, & Karagogeos, 2018). Our results are compatible with GABA‐ergic inhibitory mechanisms of brain injury, since MEG is primarily sensitive to neuronal signals from GM.

Classifiers correlate with cognitive functioning

Our results further demonstrated the cognitive relevance of some diagnostic features of cmTBI. Although mean test performance was within the normal range for both groups, the striking associations with activity within different frequency bands illustrates the power of MEG to elucidate relatively subtle changes in cognitive functioning. Delta‐theta and gamma band hyperactivity in cmTBI group was associated with worse cognition, whereas alpha and beta band hypoactivity was associated with worse cognition. For classifiers that did not show group differences in rs‐MEG activity, similar neurocognitive associations were also found across the frequency bands, mainly within the cmTBI group, with a few exceptions. This finding emphasizes the cognitive relevance of these classifiers, which may be sensitive to brain injuries in cmTBI participants with injuries that are less characteristic of the entire cmTBI cohort. Neurocognitive associations were overwhelmingly found only in the cmTBI, which suggests that they were injury related, rather than characteristic of normal cognitive activity. In addition, cognitively relevant classifiers were mostly found for frontal cortical regions (Figure 4; Table S.3), consistent with the emphasis of our neuropsychological test battery on executive functions.

Increased delta‐theta and gamma band activity was associated with poorer cognitive flexibility (Number‐Letter Sequencing), executive control (Category and Letter Fluency), and visuospatial processing/psychomotor speed (Digit Symbol Coding). Classifiers in these frequency bands were dominant in the lateral frontal cortex (e.g., DLPFC, FP, VLPFC, pSFG), but also found in medial prefrontal cortex (vmPrecentral), areas that govern higher‐order cognitive and motor functions (Bludau et al., 2014; Fuster, 2015; Koechlin, 2011; Owen, 2000; Owen, McMillan, Laird, & Bullmore, 2005). Both delta‐theta and gamma hyperactivity in the frontal cortex may signify a stronger engagement of cognitive‐control processes (Gillis & Hampstead, 2015; McAllister et al., 1999; Medaglia et al., 2012; Perlstein et al., 2004; Phillips, Parry, Mandalis, & Lah, 2017; Scheibel et al., 2003; Turner, McIntosh, & Levine, 2011), especially in individuals who exhibit signs of executive dysfunction. An exception was the association between poorer Number‐Letter Sequencing and elevated fusiform delta‐theta activity, which is compatible with its role in processing higher‐order visual information (Weiner & Zilles, 2016).

In contrast, decreased alpha and beta band activity was associated with poorer cognitive flexibility, executive control, and visuospatial processing/psychomotor speed. Classifiers in these frequency bands were notable in medial frontal areas (dACC, rACC, VMPFC), which modulate many cognitive functions. For example, dACC and rACC are routinely engaged during error detection, conflict monitoring, attention, and emotional processing (Bush, Luu, & Posner, 2000; Nieuwenhuis, Ridderinkhof, Blom, Band, & Kok, 2001; Pardo, Pardo, Janer, & Raichle, 1990; Weissman, Gopalakrishnan, Hazlett, & Woldorff, 2005), although rostral and dorsal divisions may be more or less involved depending on cognitive (Bush et al., 2002; Polli et al., 2005; Taylor et al., 2006) or emotional (Bush et al., 2000) processing demands. The VMPFC also controls higher‐level processes involved in decision making, emotion regulation, and cognitive flexibility (Kim, Johnson, Cilles, & Gold, 2011) (see review in [Hiser & Koenigs, 2018]). Hypo‐activity in this region is also common in individuals who suffer post‐traumatic stress disorder (PTSD) (Huang et al., 2014; Hughes & Shin, 2011; Rauch, Shin, & Phelps, 2006; Rauch, Shin, Whalen, & Pitman, 1998), and is a key component in an influential neurocircuitry model of PTSD (Hughes & Shin, 2011; Rauch et al., 2006). Altogether, across frequency bands, cognitively relevant classifiers were primarily located in the frontal cortex, in alignment with their relationships to executive functions and visuospatial processing and psychomotor speed.

Sample size, and training‐to‐test ratio

In the present study, the training‐to‐test ratio in our data sets was 85–10 (or 89.5–10.5%). Such setting is somewhat different from the typical DL settings of 80–20% splitting ratio between the training and testing data sets (Muller & Guido, 2016). The rationale for our splitting ratio was the moderate number of individuals in the present study cohort. Although a sample of 95 individuals is considered “large” from neuroimaging point of view, it is much smaller than the number of samples used in typical DL classification approaches. The typical 80% ‐ 20% splitting ratios, sometimes referred as the Pareto rule, is usually associated with thousands of samples in DL. In preliminary analyses, we had found the DL network coefficients became unstable, with reduced classification accuracy, when the number of individuals in the training data set was less than 80, and this phenomenon was independent of the splitting ratio. In contrast, the 85–10 (actual) or (89.5–10.5% proportion) splitting ratio provided a good balance of a stable training data set and accuracy prediction in the testing data set.

Limitations and conclusions

To our knowledge this is the first study to demonstrate that a DL neural network application (3D‐MEGNET) to rs‐MEG source magnitude imaging data produces excellent classification of Veterans and active‐duty service members with and without cmTBI. The all‐frequency band model in 3D‐MEGNET using Fast‐VESTAL source analysis significantly outperformed the individual frequency band models, underscoring the multifaceted nature of combat‐related neuronal injuries. Thus, mTBI is better characterized by considering all frequencies, rather than limiting approaches to single bands. In addition, the DL‐based 3D‐MEGNET outperformed the more traditional SVM machine‐learning approach, and 3D‐MEGNET models using Fast‐VESTAL source analysis also significantly outperformed those using Beamformer source analysis. Importantly, across all frequency bands some classifiers were cognitively relevant, particularly in the cmTBI group, suggesting that they were sensitive to declines in executive functioning and visuospatial processing/psychomotor speed. One limitation was that an intelligence quotient was not available for many participants, which could confound the interpretation of some neurocognitive associations. Future studies in civilian mTBI may include HC data from open access databases like OMEGA (https://www.mcgill.ca/bic/resources/omega) which is not applicable to the present study due to lack of active‐duty service members or Veterans with matched combat experience to our cmTBI participants in these databases. In addition, with the increased interest in conducting multisite clinical trials, future studies are needed that systematically investigate the effects of differences in sensor and noise signatures across different MEG sites, systems, and manufactures on machine‐learning based classification algorithms. These signatures are due to substantial differences in sensor configurations, noise signatures, and noise‐reduction software packages that are used across different MEG sites. Insight into the potential effects of these nuisance variables on classification could aid in refining imaging software/protocols and machine‐learning algorithms. It will also be essential to evaluate the generalizability of classification algorithm accuracy to different mTBI cohorts, who may or may not differ in demographic or clinical features (e.g., etiology, co‐morbidities). In this regard, studies are also required to evaluate the diagnostic accuracy of the obtained classifiers in distinguishing different cmTBI and HC cohorts. This may entail refinement of DL algorithms along with testing different features of rs‐MEG, to identify the most stable classifiers of mTBI or cmTBI conditions. If successful, the 3D MEGNET approach could also be extended to the prediction of recovery from PCS using longitudinal study designs or to classification of other clinical conditions that can be challenging to diagnose such as PTSD, which was recently explored using rs‐MEG functional connectivity measures and an SVM machine‐learning approach (Zhang, Richardson, & Dunkley, 2020).

Appendix S1. Supporting Information.

Click here for additional data file.