All the participants were scanned with a Siemens TIM 3.0-Tesla scanner. A multi-echo MPRAGE (MEMPR) sequence was used with the following parameters: TR/TE/TI = 2530/[1.64, 3.5, 5.36, 7.22, 9.08]/900 ms, flip angle = 7°, FOV = 256 × 256 mm, Slab thickness = 176 mm, Matrix = 256 × 256 × 176, Voxel size = 1 × 1 × 1 mm, Number of echoes = 5, Pixel bandwidth = 650 Hz, Total scan time = 6 min. Using 5 echoes, the TR, TI, and time to encode partitions for the MEMPR are similar to that of a conventional MPRAGE, resulting in similar GM/WM/CSF contrast.

Cortical reconstruction and volumetric segmentation were performed with the FreeSurfer v. 5.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 previously (Dale and Sereno, 1993; Dale, 1999; Fischl et al., 1999a,b, 2001, 2002, 2004; Fischl and Dale, 2000; Ségonne et al., 2004; Han et al., 2006; Reuter et al., 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 subcortical white matter and deep gray matter volumetric structures including the hippocampus, amygdala, caudate nucleus, putamen, ventricles (Fischl et al., 2002, 2004); intensity normalization, tessellation of the gray-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 (CSF) borders at the location where the greatest shift in intensity defines the transition to the other tissue class (Dale and Sereno, 1993; Dale, 1999; Fischl and Dale, 2000). Once the cortical models are complete, several deformable procedures can be performed for further data processing and analysis. These procedures include surface inflation (Fischl et al., 1999a)and 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., 2004; 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 shown 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 9 different measures, including the number of vertices, surface area, gray matter volume, mean cortical thickness, cortical thickness standard deviation, mean cortical curvature, Gaussian cortical 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, and intensity and overall volume of the whole brain divisions, including of the cerebrospinal fluid (CSF), intracranial volume (ICV), gray matter (GM), and white matter (WM). Two of the selected measures are the most common features in 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, 2017b).

We used measures including the mean cortical thickness, cortical thickness standard deviation, surface area, volume, mean curvature, white matter volume, subcortical segment volume, subcortical intensity, and overall brain volume and intensity as the structural features. All features were acquired as morphological statistics for each subject during the FreeSurfer based preprocessing of the structural MRI data. In addition, after the preprocessing, FreeSurfer's QA Tools were used for the detection and removal of outliers and negative features (Qureshi et al., 2017b).

Resting state functional MRI data were collected with single-shot full k-space echo-planar imaging (EPI) with ramp sampling correction using the intercomissural line (AC-PC) as a reference (TR: 2 s, TE: 29 ms, matrix size: 64 × 64, 32 slices, voxel size: 3 × 3 × 4 mm³). During the acquisition process, all subjects were instructed to keep their eyes open and stare at the fixation cross.

Preprocessing of functional MRI data was based on Analysis of Functional Neuroimages AFNI software; http://afni.nimh.nih.gov/afni/ (Cox, 1996). Every single EPI volume was coregistered to the corresponding anatomical image of the subject and mapped to the 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 excluded time-points based on the number of outliers and head motion magnitude. The same number of slices was excluded for all subjects. Slice alignment was applied using the local Pearson's correlation (LPC) cost function. Correction of head motion along with averaging of EPI volumes was performed to obtain a mean functional image. Each EPI volume underwent linear multiple regression to regress the motion derivatives and effects of the white matter and cerebrospinal fluid. Spatial smoothing was performed using a Gaussian kernel with a blur size of 6-mm full width at half maximum (FWHM). A polynomial detrending was applied (Qureshi et al., 2017b).

A global connectivity measure was used to calculate the average brain-wise correlation coefficients (GCOR) of all the possible combinations of voxel time series. The GCOR estimation of cortical regions is a computationally expensive process. It involves the calculation of M (M−1)/2 correlation estimates for an M voxel volume (Saad et al., 2013). AFNI simplifies the process by taking reduced time series of each voxel and scaling it by its Euclidean norm. In addition, it averages the scaled time series over the whole brain mask and finally the length l₂ − norm of this averaged series represents the GCOR. We used the AFNI program 3dTcorrMap to get the global functional connectivity maps. This involved the calculation of correlation (Pearson's r) of the residual time series in each voxel with every other voxel in their brain mask and recording the mean correlation back in the voxel. These connectedness values were further transformed with Fisher's z-transformation to yield normally distributed measures (Gotts et al., 2012). We used the DKD_ Desai_PM Atlas (Desikan et al., 2006) of the AFNI package to acquire the ROI measures from 102 cortical and subcortical regions across the whole brain. Figure 1 depicts a typical global functional connectivity map for one of the control subjects from the COBRE dataset.

Resting state functional MRI data were preprocessed again with FSL (FMRIB Software Library, www.fmrib.ox.ac.uk/fsl) 6.0 for acquiring the group independent component analysis (gICA) based connectivity measures. FSL Multivariate Exploratory Linear Optimized Decomposition into Independent Components (MELODIC) version 3.14 was utilized to perform a single-session ICA. The number of independent components was set as 30, as a model with an order higher than 20 components is usually required for the detection of some components (e.g., S1, S2, striatum), whereas a higher order model shows a decrease in ICA repeatability (Abou-Elseoud et al., 2010). We used variance normalization and thresholded the independent component maps with an alternative hypothesis test that was based on the fitting of a Gaussian/gamma mixture model to the distributions of the voxel intensities within the spatial maps and controlling the local false-discovery rate at p < 0.5 (Smith et al., 2004; Beckmann et al., 2005). Among the 30 independent components, 11 were discarded as noise and/or artifacts upon visual inspection by an experienced clinical psychiatrist. Finally, 19 components were identified as functional networks or subnetworks. Similar networks were used in our previous study (Oh et al., 2017). Figure 2 depicts some of the 19 selected components; they represent the well-known resting state functional networks: the default mode network, sensorimotor network, visual network, and frontal network.

Using these components, we applied the FSL dual regression technique (Filippini et al., 2009; Littow et al., 2010; Veer et al., 2010), to acquire subject-specific spatial maps. During this process, we utilized the group spatial maps in a linear model that was fit against the separate fMRI datasets and then acquired the time-course matrices depicting the temporal dynamics for each subject and component. After that, we could acquire the subject-specific spatial maps with the time-course matrices. Next, we calculated the correlation coefficient between the selected 19 components. It resulted in a 19 × 19 correlation based connectivity matrix. We isolated the upper diagonal elements of the connectivity matrix and after vectorizing them, we used them as the gICA features for each subject. Figure 3 depicts a typical ICA connectivity matrix for one of the control subjects from the COBRE dataset.

Finally, we gathered four types of functional connectivity based features, namely: cortical ROI-based global connectivity, subcortical ROI-based global connectivity, whole brain-based global connectivity, and gICA connectivity features. The first three types of features were based on the DKD_ Desai_PM Atlas ROI's (Desikan et al., 2006) and all of them were acquired by considering the average measure of global connectivity of the atlas ROIs.

Table 2 presents a brief overview of the structural and functional measures and number of features used in this study.

The extreme learning machine originally proposed by Huang et al. (2006) has been adopted in many previous neuroimaging studies (Termenon et al., 2013; Zhang and Zhang, 2015; Qureshi et al., 2016), and (Qureshi et al., 2017b) in the binary and multiclass settings. The ELM randomly assigns the weights and bias to the input data to compute the output weight matrix. This random assignment of weights makes the ELM algorithm very fast compared to other gradient-based classifiers such as the SVM. A more detailed discussion of the classifier can be found elsewhere (Huang et al., 2006; Qureshi et al., 2016). The ELM classifier requires the “number of nodes” as the only hyper-parameter that has to be tuned for achieving maximum performance in terms of accuracy. In this study, we used the Matlab implementation of the ELM. A greedy search method was used to tune this parameter for achieving maximum test accuracy. In this study, the search scale for selecting this parameter was set to N = [90, 91, …, 220]. A careful selection of the number of nodes is necessary to prevent the classifier from overfitting. All the features were normalized and scaled to values between –1 and 1 for improving the performance of the ELM classifier (Qureshi et al., 2016, 2017b).

The SVM classifier, which was originally proposed by Cortes and Vapnik (1995), has been one of the most popular machine-learning tools in the neuroscience domain in the last decade. It is a supervised classification algorithm. It maps features in higher dimensional space using linear and non-linear functions known as kernels. In this study, we used both linear and non-linear (RBF) SVM. In case of SVM-RBF, value of scaling parameter “σ” was tuned in the range σ = [0.000001, 0.00001, …, 100000, 1000000] and box-constraint parameter “C” was tuned in range of C = [0.005, 0.01, 0.1, …, 100000, 1000000] using logarithmic grid search method.

The LDA classifier, originally proposed by Duda and Hart (1973), is a linear classifier that utilizes hyperplanes to discriminate the data. These hyperplanes maximize the interclass mean while keeping the interclass variance minimized.

Random forest trees, originally proposed by Breiman (1996), are a method of building a forest of uncorrelated trees with randomized node optimization and bagging. Out-of-bag errors are used as an estimate of the generalization error. The random forest (RF) measures variable importance through permutation (Liaw and Wiener, 2002). The general bootstrap aggregation algorithm was used for training. In RF implementation, with decision trees as the learner type, only the number of trees and the number of splits for prediction need to be defined. In this study, we used 30 trees and 20 splits.