We obtained neuroimaging data of ADHD patients from the ADHD-200 MRI dataset (Biswal et al., 2010), which is publicly available at the 1,000 Functional Connectomes Project website http://fcon_1000.projects.nitrc.org/indi/adhd200. These T1-weighted structural and resting state functional MRI scans were acquired at six different institutes. The ages of the participants ranged from 7 to 14 years. The sites were Brown University (BU), New York University Child Study Center (NYU), Beijing Normal University (BNU), Kennedy Krieger Institute (KKI), Oregon Health, and Science University (OHSU), and Washington University in St. Louis (WU). All participants have been scanned under 3.0-Tesla scanners. The other technical details about the scanner parameters from each participating site are available at the above-mentioned URL of the ADHD-200 global competition. In addition, all the sites contributing to ADHD-200 had been approved by their local institutional review board (IRB) and complied with local IRB protocols (Qureshi et al., 2016).

Cortical reconstruction and volumetric segmentation were performed with the FreeSurfer v5.3.0 image analysis suite, which is documented and freely available for download online at http://surfer.nmr.mgh.harvard.edu/. The technical details of these procedures have been described elsewhere (Dale and Sereno, 1993; Dale et al., 1999; Fischl et al., 1999a,b, 2001, 2002, 2004a,b; Fischl and Dale, 2000; Ségonne et al., 2004; Han et al., 2006; Jovicich et al., 2006; Reuter et al., 2010, 2012). Briefly, the preprocessing procedure includes motion correction of volumetric T1-weighted images, removal of non-brain tissue using a hybrid watershed/surface deformation procedure (Ségonne et al., 2004), and automated Talairach transformation. Moreover, the procedure includes the segmentation of the subcortical white matter and deep gray matter volumetric structures including hippocampus, amygdala, caudate, putamen, ventricles (Fischl et al., 2002, 2004a), intensity normalization (Sled et al., 1998), tessellation of the gray matter-white matter boundary, and automated topology correction (Fischl et al., 2001; Ségonne et al., 2007). In addition, surface deformation following intensity gradients was performed to optimally place the gray/white and gray/cerebrospinal fluid borders at the location where the greatest shift in intensity defines the transition to the other tissue class (Dale and Sereno, 1993; Dale et al., 1999; Fischl and Dale, 2000). Once the cortical models are complete, a number of deformable procedures can be performed for further data processing and analysis. These procedures include surface inflation (Fischl et al., 1999a), registration to a spherical atlas that is based on individual cortical folding patterns to match the cortical geometry across subjects (Fischl et al., 1999b). Moreover, parcellation of the cerebral cortex into units with respect to gyral and sulcal structure (Fischl et al., 2004b; Desikan et al., 2006) and creation of a variety of surface based data including maps of curvature and sulcal depth are part of the procedure. This method uses both intensity and continuity information from the entire three-dimensional MR volume in segmentation and deformation procedures to produce representations of cortical thickness, calculated as the closest distance from the gray/white boundary to the gray/CSF boundary at each vertex on the tessellated surface (Fischl and Dale, 2000). The maps are created using spatial intensity gradients across tissue classes rather than simply relying on absolute signal intensity. The maps produced are not restricted to the voxel resolution of the original data, and are thus capable of detecting submillimeter differences between groups. Procedures for the measurement of cortical thickness have been validated against histological analysis (Rosas et al., 2002) and manual measurements (Kuperberg et al., 2003; Salat et al., 2004). FreeSurfer morphometric procedures have been demonstrated to show good test-retest reliability across scanner manufacturers and across field strengths (Han et al., 2006; Reuter et al., 2012). A cortical surface-based Desikan-Killiany-Tourville (DKT) atlas Klein and Tourville, 2012) was mapped to a sphere aligning the cortical folding patterns, which provided accurate matching of the morphologically homologous cortical locations across subjects. For each of the DKT31 protocol-based segments, FreeSurfer calculated the nine different measures, including number of vertices, surface area, gray matter volume, average cortical thickness, cortical thickness standard deviation, cortical mean curvature, cortical Gaussian curvature, cortical folding index, and cortical curvature indices (Colby et al., 2012). For the subcortical regions, FreeSurfer calculated the area and volume of the whole segment, white matter volume, intensity and overall volume of the whole brain divisions including cerebrospinal fluid (CSF), intracranial volume (ICV), gray matter (GM), and white matter (WM). Two of the selected measures are the most common features in the structural studies (Arbabshirani et al., 2017). The surface area was calculated by computing the area of every triangle in a standardized spherical surface tessellation. The local curvature was computed using the registration surface based on the folding patterns (Qureshi et al., 2016).

Preprocessing of functional MRI data was based on Analysis of Functional Neuroimages AFNI software; http://afni.nimh.nih.gov/afni/ (Cox, 1996). Every single echo planar image (EPI) volume was co-registered to the corresponding anatomical image of the subject and mapped to Talairach coordinates space with the TT_N27+tlrc template. We excluded the first six images from each EPI volume to achieve the MR steady state. In addition, slice-timing correction was performed. We censored and cut out time points based on their number of outliers and head motion magnitude. Same number of slices were cut out for all subjects. Slice alignment was applied by using the local Pearson's correlation (LPC) cost function. The correction of head motion along with averaging the EPI volumes was performed to obtain a mean functional image. The dataset was already reoriented to right posterior inferior (RPI) in the ADHD-200 dataset repository. Each EPI volume underwent the linear multiple regression to regress the motion derivatives and effects of the white matter and cerebrospinal fluids. Spatial smoothing was performed by using a Gaussian kernel with a blur size of 6-mm full width at half maximum (FWHM). A polynomial detrending was applied. AFNI program 3dmaskave was used to calculate the beta values for the correlation analysis. A run of 10,000 Monte Carlo simulations was conducted with AlphaSim program. The cluster size of 10 voxels was determined at a family-wise error rate corrected with p < 0.000001 to avoid the problem of multiple comparisons. As a result, the clusters consisting of lower than 10 voxels were discarded from the analysis.

After preprocessing, the structural data were converted into surface maps for ANOVA analysis. Heat kernel smoothing was applied on the dataset after conversion. An AFNI program 3dANOVA2 (Cox, 1996) was used for this analysis.

First, to examine the differences among the three groups, the multiple comparisons test was performed using the AFNI program 3drefit with the FDR correction on both the structural and the functional data. In addition, for the structural data, the AFNI program SurfClust was used to find regions that have significant differences between groups with an area of more than 50 mm². The p-value threshold was set to be 0.001 (Qureshi et al., 2016). This was followed by a post-hoc t-test to evaluate the actual differences among the three groups. For the functional data, we used AFNI program 3dclust to generate the mask of significant ROIs of ANOVA analysis. In addition, 3dROIstats was used with this mask to acquire the average connectivity value from each subject. Since these data were very significant, we set the p-value at a very low level of 0.000001. Finally, we extracted 12 significant regions from the structural data and nine significant regions from the functional data that served as the ANOVA-based features for classification. It is worth noting that the ANOVA analysis was performed only on the training dataset to avoid the potential bias (Arbabshirani et al., 2017) of the overall accuracy of the classifier.

The Pearson correlation coefficients were calculated between structural and functional ANOVA-based significant ROIs and age separately. For the structural data, we extracted the average cortical thickness of each ROI for each individual subject using the AFNI program 3dmaskave. The p-value of the correlation coefficients was calculated based on our sample size (Rosenhead, 1963). For the functional data, we calculated the average beta value of each ANOVA-based ROI of each individual subject using AFNI program 3dmaskave.

In neuroimaging, pattern recognition is sometimes referred to as multi-voxel pattern analysis (MVPA), as voxels are often used directly as features. The features can be categorized as voxel, region, or network-derived set of values. In a voxel-based feature set, features are extracted on the voxel level. In a region-based feature set, features are derived by parcellation of brain images from predefined regions based on anatomical or functional brain atlas. In a network-based feature set, features are extracted by combining voxels across networks, like the functional connectivity based features (Wolfers et al., 2015). We used the last two approaches (region and network driven features) in this study in addition to an uncommon approach of statistical measurement (ANOVA) based feature extraction. The region-based feature selection is the most common approach used in the fMRI studies (Arbabshirani et al., 2017).

In this study, the multi-measure features were acquired from the significant regions obtained by ANOVA analysis and atlas-based ROIs separately for both structural and functional MRI datasets. Finally, we gathered a big pool of features as shown in Table 2. Therefore, optimal feature selection was necessary to get a high classification accuracy.

Two classifiers were used in this study, namely, linear extreme learning machine (ELM) and support vector machine (SVM) with linear and radial basis function (RBF) kernels. The classifiers were used in both binary and multi-class settings. In the case of multi-class, we used the one versus all (OVA) approach. In addition, we also mentioned the ELM-based classification results without applying the feature selection to validate the significance of the proposed feature selection framework. A brief description of both classifiers used in this study is as following.

The classifier performance was measured in terms of classification accuracy and Kappa score for multi-class classification. Classification accuracy was determined as

Where TP is true positive, TN is true negative, FP is false positive and FN is false negative. These values were obtained by computing the confusion matrix (Qureshi et al., 2016). Kappa score is considered as a robust statistical measure for the assessment of multi-class classifiers. It is defined by Cohen (1960) as

Where Pr(a) is the relative observer agreement between raters and Pr(e) is hypothetical probability of chance agreement. The value of k ranges between +1 and −1, where k = 1 represents perfect classification, k = 0 represents chance level, and k = −1 represents completely erroneous classification (Wang et al., 2012; Qureshi et al., 2016).

A Permutation test was used to assess the statistical significance of the ELM classifier performance (Golland and Fischl, 2003). Briefly, it works as follows. First, we choose the actual test accuracy as the test statistic of the classifier, the class labels for testing dataset permuted randomly and given to the classifier and check for the cross validation. Generally, the lower p-value of the permuted prediction rate against the prediction rate with the original data labels indicates the higher significance of the classifier performance. There is no fixed rule for setting the number of permutations. We have permuted the data 10,000 times in the current study (Qureshi et al., 2016).

Figure 2 depicts the overview of the classification framework of this study.

ANOVA analysis was applied to the cortical thickness measures acquired by FreeSurfer and further smoothed with Heat kernel by AFNI preprocessing pipelines. We used only cortical thickness measures for the ANOVA analysis as it is the most important structural imaging biomarker (Shaw et al., 2007; Narr et al., 2009). A basic description of this results was presented briefly in a recent previous work (Qureshi et al., 2016). In the present study, we are using the ANOVA results for extracting the features for classification. The ANOVA tests showed (see Figure 3) 10 areas of significant thickness change in the left hemisphere and two areas in the right hemisphere (corrected at p < 0.001, cluster size 50 mm²). AlphaSim was applied on the result for correction. Table 3 summarizes our findings.

The ANOVA tests of resting state fMRI data showed nine areas of significant change (corrected at p < 0.000001, cluster size 11 mm²) (see Figure 4) from four significant functional networks of the human brain including the ventral attention network, sensorimotor network, visual network, and cerebellar network. Table 4 summarizes the details of the significant regions.

Our classification results confirm that these regions can be used as potential biomarkers/features for the multi-class classification of ADHD s/fMRI dataset. The ANOVA-based features gave high classification accuracy both as standalone modality and in combination with ROI-based features.

Our results suggest that the optimized features are robust and yield higher classification accuracies both in binary and multi-class classification settings. We achieved a maximum of 76.190% accuracy (p < 0.0001) on the test data by using the proposed hierarchical feature extraction framework with (ELM) in multi-class settings. Table 5 summarizes the overall multi-class (one vs. all settings) accuracies, kappa score, and p-values of structural, functional, multi-modal, and multi-measure data.

Binary classification results are presented in the current study to demonstrate the fact that, ideally, multi-class classification results show equaling and generally lower accuracy than results from the binary classification of the same dataset (Qureshi et al., 2016).

The major result of the current study is the highest prediction accuracy of the ADHD dataset using all the measures and modalities together as features; therefore, we performed the binary classification on the combined measures only. In binary classification settings, we achieved a maximum of 92.857% (p < 0.0001) prediction accuracy. It is interesting to note that the proposed feature selection framework selected only structural features from multi-modal input to acquire this accuracy. Most significant feature measures were mean cortical thickness, volume, area, and standard deviation of cortical thickness respectively. Even though the functional features were not selected in the final selection but they assist the selection algorithm to extract the most discriminative structural features (both ROI and ANOVA-based) from the pool. When we used only the structural features for binary classification we could not reach the same accuracy as we did while using the multi-modal features. In addition, we were not able to achieve this high accuracy in the previous work (Qureshi et al., 2016) because all of the cortical feature measure types utilized in the present study were not used in it. We also learned from the previous study (Qureshi et al., 2016) that many of the feature measures do not have any predictive power including folding and curvature indices, therefore we did not include those measures in this study. Briefly, we achieved the following binary accuracies with each individual modality. With only structural ANOVA features 89.29%; with only functional ANOVA features 82.14%; with only functional ROIs 85.17%, and with only structural ROIs 85.30%.

Table 6 summarizes the classification scores of all the disease groups along-with sensitivity, specificity, F1-score, recall, and precision measures for all the multi-modal and multi-measure features.

It is important to note that two of the three ELM classifier-based binary classification scores (ADHDC-TDC and ADHDI-TDC) without feature selection overlap exactly with the accuracy of the linear SVM with the most discriminative features. In other words, we can say that the lowest accuracy of the ELM classifier overlaps with the highest accuracy of the SVM classifier.

After applying the Bonferroni correction method with p < 0.000463, we found two regions that were negatively correlated with age in TDC group (Table 7). These regions were left middle temporal gyrus and left superior frontal gyrus. On the other hand, some individual IQ and symptoms severity score data were not included in the original database. Therefore, we only use the subjects with available IQ and symptoms severity score information to perform the calculation of the correlation coefficients. Those subjects with missing IQ information included nine from the TDC group and 13 from the ADHDC group. The subjects with missing symptoms severity score were included 29 from ADHDI group, 11 from ADHDC group and 36 from TDC group. In contrast, there was no missing age information.

For the functional data, the beta value for each of the ANOVA-based significant cluster was calculated by using AFNI program 3dmaskave. The Pearson correlation coefficients were calculated between the beta values and age, IQ, and symptoms severity scores.

After applying the Bonferroni correction method with p < 0.000617, no significant region was found.

Highest subgroup classification accuracies were obtained by using the significant features as shown in Table S2 (Supplementary Material).

In recent years, many studies about classification-based diagnosis of ADHD have been reported. Most of these studies were binary and the classification scores were ranging from about 72 to 95% (Igual et al., 2012; Lim et al., 2013; Peng et al., 2013; Johnston et al., 2014; Deshpande et al., 2015; Hammer et al., 2015; Iannaccone et al., 2015; Wolfers et al., 2015; Qureshi and Lee, 2016; Qureshi et al., 2016; Xiao et al., 2016). Except for one study, all used feature selection methods prior to classification (Wolfers et al., 2015). In Xiao et al. (2016), a different ADHD dataset was used with only 47 subjects, compared to the current study. Although the theme of this study by Xiao et al. (2016) is similar to our present work, the results are not comparable due to dataset differences and their use of only cortical thickness ROIs. In addition, our main focus was the multi-class, multi-measure, and multi-modal classification of the ADHD-200 dataset. Previously, a very few multi-class (Fair et al., 2012; Qureshi et al., 2016) and multi-modal (Colby et al., 2012; Dai et al., 2012) classification studies have been published with ADHD-200 dataset. However, to the best of our knowledge, the present study is the first multi-modal, multi-measure, and multi-class classification study of the ADHD-200 dataset. In addition, the multi-class classification accuracy of 76.190% (p < 0.0001) is the highest reported accuracy on the ADHD-200 dataset in multi-class settings. However, all the aforementioned studies, except a recent study (Qureshi et al., 2016), used different subject selection compared to us. Therefore, the present results are not comparable to the previous studies. The multi-class classification score of the structural part 73.810% (p < 0.0001) of the present study is much higher than the results (60.78%) in Qureshi et al. (2016), which utilized the exact same subject selection for classification. The present study is also comparable to Qureshi et al. (2016) in that it also utilized sparse hierarchical feature extraction framework with basic ELM, while Qureshi et al. (2016) used H-ELM classification framework, which is a similar approach in principal. However, the method of the present study outperformed all the results in Qureshi et al. (2016) because the proposed hierarchical feature extraction framework can extract features with the higher predictive power as compared to the sparse autoencoder. In addition, in the present study we utilized many additional feature measures that were not used in Qureshi et al. (2016). We believe that these measures played a vital role in increasing the classification accuracy. However, the exact causes of the boosted accuracy are not yet determined.

We two significant negative Pearson correlation coefficients between the cortical thickness of the regions and age. During adolescence, normal children undergo the cortical thinning phase (Sowell et al., 2004). Consistent with these findings, the two cortical regions in our results showed a decreasing tendency of cortical thickness with age in the TDC group. However, this pattern seemed to be weaker in the ADHDC and ADHDI groups: no regions in each had negative correlation coefficients. In other words, the patients exhibited some deficits in cortical thinning in the middle temporal gyrus and superior frontal gyrus. While considering the middle temporal gyrus; it can be included in ventral attention networks and superior frontal gyrus is related to the dorsal attention networks, therefore, the structural thinning deficits might be used as important markers to discriminate patients and controls (Aboitiz et al., 2014). Thus, we postulate that normal cortical thinning development was interrupted in the ADHD groups. Based on previous findings suggesting that ADHD is related to cortical thinning, the patients might experience abrupt thinning of some brain regions related to attention and motor functions (Makris et al., 2007; Narr et al., 2009).

Moreover, it is evident from our classification results that those significant statistical differences may underlie the high classification accuracy. The results in Vaidya (2011), McLeod et al. (2014), Tegelbeckers et al. (2015) further endorse our results of the correlation analysis of ANOVA-based significant regions with age.

The significant regions from both structural and functional ANOVA results play an important role in acquiring highly accurate differential diagnosis when used in combination with the ROI based significant features that did not form part of the outcome of the ANOVA analysis. The following description of the regions obtained from our ANOVA analysis may have clinical significance:

Most of the neuroimaging studies addressing ADHD emphasize that the selected regions in the present study has a strong connection with the neurodevelopmental basis of this disorder. Specifically, the cortical thinning problem and the attention and motor network hypoconnectivity have been suggested to be involved in ADHD (McCarthy et al., 2014; McLeod et al., 2014). These regions showed significant differences between the ADHD and control groups (Qureshi et al., 2016). The functional ANOVA results again endorse the significance of our noteworthy results. A few of the significant regions from our functional ANOVA results overlap with those of structural ANOVA results. They include the left precentral gyrus, the right precentral gyrus, the left middle frontal gyrus, and the right superior temporal gyrus, which are all related to the motor or attention networks.

Our significant feature results showed that structural information is more crucial for the classification. It means that combined gross classification by sMRI and supporting information from fMRI can increase accuracy. In our structural ANOVA results, multiple areas of the frontal cortex, including the middle frontal gyrus (MFG), the precentral gyrus, and the postcentral gyrus in the parietal cortex, showed significantly higher cortical thickness in the TDC group. Our findings are in accordance with previous results (Durston, 2003; Shaw et al., 2007; Valera et al., 2007), and this might contribute to the high classification accuracy of our model. On the other hand, although not much is known about the temporal lobe (Lorberboym et al., 2004), it is evident that the cortical thickness differs significantly in ADHD patients as compared to the normal control group (Fernández-Jaén et al., 2014). Cortical thickness changes over time, therefore we focused on the dynamic changes of thickness according to the ages. Generally, ADHD has cortical thinning compared to controls after development. In Narr et al. (2009) the mean age was of 11.7 years. However, ADHD has cortical thinning problems during the development which are related to the delayed cortical maturation. The mean age of peak brain thickness is 10.5 years for ADHD and 7.5 years for controls. It means after 10.5 years of age, ADHD experiences fast decreases of cortical thickness while controls have enough time to maturate the brain. In line with these evidences, we can say that ADHD has delayed cortical thinning, but finally they can get thinner cortex than controls. In addition, many brain areas can maturate with different speeds, therefore that the results were confusing for us. Our results from the functional ANOVA reconfirm this finding. The middle temporal gyrus is associated with language abilities, visual perception, multi-modal sensory integration, and semantic memory processing (Qureshi et al., 2016). Our findings from the functional ANOVA result are consistent with these results. The bilateral postcentral gyrus (Martinussen et al., 2005), right declive (Ortiz et al., 2015), and Left middle occipital gyrus were reported to be affected in ADHD. Visual perception ability is also related to the visuospatial working memory which is the common symptom of ADHD (Valera et al., 2010; Vaidya, 2011; McLeod et al., 2014). ADHD is related to the semantic memory and language problems (Felton et al., 1987). Our functional ANOVA results endorse these findings, thus making our observations highly beneficial to clinicians.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.