Enhancing multi-center generalization of machine learning-based depression diagnosis from resting-state fMRI

Subjects

One hundred sixty-three patients with MDD (age 20-75, average 44.1 ± 12.2) were recruited by 5 sites (the Psychiatry Department of Hiroshima University and collaborating medical institutions, Table 1). They were screened using the Mini International Neuropsychiatric Interview (M.I.N.I,) [13, 14], which enables medical doctors to identify psychiatric disorders according to DSM-IV criteria [15]. Patients had an initial MRI scan before or after starting medication within 0-2 weeks.

As a control group, 195 healthy subjects (age 20-75, average 39.1 ± 15.5) with no history of mental or neurological disease were recruited from the local community. All control subjects underwent the same self-assessment and examination administered to the MDD group. Thereafter, for the melancholic MDD classifier, the dataset was limited to have the subtype of melancholia (based on M.I.N.I.). The numbers of patients and healthy controls were set to be equal for each site except for the subjects from the National Institutes for Quantum and Radiological Science and Technology, in order to develop a classifier unbiased toward either group (see Table 1). The subjects from the National Institutes for Quantum and Radiological Science and Technology were used only for the test data (replication cohort).

The study protocol in this study was approved by the Ethics Committee of Hiroshima University and the Radiation Drug Safety Committee and by the institutional review board of the National Institutes for Quantum and Radiological Science and Technology, in accordance with the ethical standards laid down in the 1964 Declaration of Helsinki and its later amendments. All necessary patient/participant consent has been obtained and the appropriate institutional forms have been archived. The detailed properties of the subjects are shown in Table 1.

Acquisition and preprocessing of functional MRI data

fMRI scanners were used to generate magnetic resonance images. Functional data were collected using gradient echo planar imaging (EPI) sequences. High-resolution T1-weighted magnetization-prepared rapid gradient echo images were also acquired before scanning the functional data. In the scan room with dimmed lights, participants were instructed not to think of anything in particular, not to sleep, and keep looking at a cross mark at the center of the monitor screen. The first 4 to 7 images were discarded to allow magnetization to reach equilibrium. All participants underwent an approximately 5 to 10 min resting-state scan. The scanners and imaging parameters are different depending on the site (see Table 1 and 2).

All the resting-state fMRI data were preprocessed using the identical procedures described in [7]. T1-weighted structural image and resting-state functional images were preprocessed using SPM8 (Wellcome Trust Centre for Neuroimaging, University College London, UK) on Matlab R2014a (Mathworks inc., USA). The functional images were preprocessed with slice-timing correction and realignment to the mean image. Thereafter, using the normalization parameters obtained through the segmentation of the structural image aligned with the mean functional image, the fMRI data was normalized and resampled in 2 ⨯ 2 ⨯ 2 mm³ voxels. Finally, the functional images were smoothed with an isotropic 6mm full-width half-maximum Gaussian kernel. After these preprocessing steps, the scrubbing procedure [16] was performed to exclude any volume (i.e., functional image) with excessive head motions, based on the frame-to-frame relative changes in time series data.

Calculation of functional connectivity

For each individual, the time course of fMRI data was extracted for each of 137 regions of interest (ROIs), anatomically defined in the Brainvisa Sulci Atlas (BSA; http://brainvisa.Info) [6, 16] covering the entire cerebral cortex without a cerebellum. After applying a band-pass filter (0.008– 0.1 Hz), the following nine parameters were linearly regressed out: the six head motion parameters from realignment; the temporal fluctuation of the white matter; that of the cerebrospinal fluid; and that of the entire brain. Pair-wise Pearson correlations between 137 ROIs were calculated to obtain a matrix of 9,316 functional connectivities for each participant.

Data standardization, removal of nuisance variables, and feature selections

For functional connectivity data, we applied the Fisher z-transformation to each correlation coefficient. We then regressed them on training subjects’ ages, sex, and dummy variables for each site. The resulting residuals of the regression are used for the subsequent procedures as the features controlling for age, sex, and site effects. We also used the regression coefficients to remove nuisance variables such as age, sex, and site from test subjects. Features like functional connectivities were high-dimensional. We then used a rank-sum test for dimensional reduction; we used the threshold of 0.05 unless noted. The procedures are summarized in Figure 1.

• Download figure
• Open in new tab

Figure 1. The procedure of the removal of site bias and feature selection.

LOOCV

Use one subject as test data and the remaining subjects as training data, and repeat for all subjects LOSOCV

Use subjects from one site as test data and subjects from the remaining 3 sites as training data, and repeat for all subjects

Independent site test

Use subjects from four sites as training data, and use subjects from the independent site as test data. In some cases, a part of data was used for the removal of nuisance variables

Classification algorithms

We then applied machine learning classification using the above features. We used ensemble learning (bagging and AdaBoost), support vector machine (SVM), or sparse logistic regression (SLR) [17] to classify depressed patients and healthy controls, or melancholic patients and healthy controls. Ensemble learning is a method to make a classifier by combining weak classifiers. While bagging makes the decision based on voting by parallel weak classifiers, AdaBoost combines weak classifiers sequentially. The hyperparameters adjusting tree depth (minimum leaf size and the maximum number of splits), and the other hyperparameters such as the number of ensemble learning cycles, and learning rates were optimized using Bayesian optimization described in the following section. Ensemble learning has advantages in dealing with small sample size, and high-dimensionality [18]. SVM is a widely used classifier and performs classification by finding the hyperplane that differentiates the two classes. Types of the kernel and penalty parameter were determined by the optimization in the following section. SLR is a Bayesian extension of logistic regression in which a sparseness prior is imposed on the logistic regression. SLR has the ability to select features objectively that are related to classifying depressive disorders [2, 7, 17].

Evaluation and estimation of the hyperparameters of Classifiers

The performance of the classification was evaluated using Leave-One-Subject-Out cross-validation, Leave-One-Site-Out cross-validation within discovery cohort or independent dataset from a replication cohort.

For Leave-One-Subject-Out cross-validation, one subject was used for the test and the rest subjects were used for the training and this process was repeated until all subjects were used as a test data. The hyperparameters of the AdaBoost, Bagging, and SVM were optimized with Bayesian optimization using Gaussian process [19].

For the optimization, the whole data was divided into 9-folds keeping an equal amount of (diagnosis, gender, and site) combinations per fold. The optimization was performed using the 8-folds with 8-folds cross-validation (inner 7-folds are used for the training and the remaining inner one-fold was used for the validation). The test data was selected in the outer 1-fold. The rest from outer 1-fold and the outer 8-folds were used for the training of the classifier [7].

For Leave-One-Site-Out cross-validation, subjects from three out of the four sites were used as the training data and subjects from the remaining sites were used as test data. In the case that independent dataset was used as a test data, all subjects from the four sites were used for training. For both cases, the optimization of the hyperparameters was done using 8-folds cross-validation within the training data.

Generalization for different sites

There are some studies trying to overcome site bias. Some studies used independent component analysis or sparse canonical component analysis to obtain site-independent brain activity [23]. These methods are time-consuming and it is possible to lose important features that characterize patients with depression. On the other hand, regressing out is a simple way to reduce site bias and retains all features.

In our study, SVM showed good classification performance for independent site data. Depressive disorder is known to be heterogeneous. The melancholic depressive disorder is a typical and severe type of disorder. This would be a reason why the classification performance is better when we focus on melancholic patients. Targeting other types of depressive disorders would be future work. There are some studies about the subtyping the depressive disorder by using biological markers including functional connectivity [1, 6]. If a new subtype of depressive disorder that showed characteristic neural behavior is found, the diagnosis for this subtype using machine learning would be improved.

For clinical diagnosis

When we consider an application for clinical diagnosis, site bias is a large problem. If traveling subjects, the same persons visit multiple sites, are available, the calibration would be more effective and better performance of the diagnosis of patients with major depression would be expected [2]. Even if traveling subjects are not available, our proposed procedure would be helpful. Site bias should be calibrated using data from healthy subjects before clinical diagnosis, using the fMRI scanner used for the clinical diagnosis, which would be a more practical way for an application at a novel site.