2.1.1 ADHD-200.

The ADHD-200 dataset is a multi-site combination of neuroimages taken from 8 sites (of which we use 7). It is partitioned into two disjoint sets (training and holdout), where the prediction model is learned from only the training set, and the learned model’s accuracy is measured on the holdout set. The data from the excluded site (Brown University) are not used because subject labels (ADHD vs. control) are not available for that site and we require labels for our machine learning algorithms. Thus, our training set consists of 776 resting state scans: 491 were taken from healthy controls and 279 from patients. To balance our training set, we used all 279 patients and selected 279 healthy controls evenly taken from all the sites as our training set [11, 12, 16], which means that the baseline classification accuracy for the training set is 0.50. The ADHD-200 competition hold-out dataset (excluding Brown University data) consists of 171 subjects: 94 healthy subjects and 77 ADHD cases (baseline accuracy 0.5497). We only use this set for evaluating the quality of the final model; n.b., it is not used in any way during the training process.

The ADHD-200 dataset also includes other non-imaging features for each subject, including gender, age, handedness, site of the imaging, IQ measure, etc. [17]. However, we only use imaging data for our experiments.

2.1.2 ABIDE.

The ABIDE dataset consists of 1111 scans, taken from 17 sites consisting of 573 healthy controls and 538 patients with autism. To evaluate each learning model, we use 800 subjects (70%) for model training and 311 subjects (30%) for hold-out testing. We use the same case/control ratio (0.5157) for both training and test set. The ABIDE dataset provides an array of non-imaging information which includes age, gender, handedness, various IQ scores, site of the imaging and eyestat (which indicates whether the person kept their eyes open or not during the scan). Again, we only use imaging data for our experiments.

2.4.1 Background: Filters, convolution.

While our application deals with 3D images, this section explains the basic ideas using simpler 2D examples.

In image processing, filters are transformations that accentuate certain features within an image. Filters are generally defined on a neighborhood. For example, the 2D image filter (1) is defined on 3 × 3 neighborhood. Below, we view h as a function—here h_contrast: {−1, 0, 1} × {−1, 0, 1} ↦ ℜ.

We then “convolve” a filter against an image. More precisely, a convolution is a mathematical operation on two functions—an image and a filter (where are the integers from −d through + d, inclusive)—to produce a third function that is typically viewed as a modified version of I(x, y), giving the overlap between the two functions. In case of 2-dimensional convolution, I is a 2D image (where I(x, y) is the intensity at position (x, y)), and h(u, v) is the filter that is convolved with the image—such as the h_contrast from Eq 1. The result of the convolution can be interpreted as the similarity measures between each pixel of the image and the filter. For 2D images, this result, at location (x, y), is defined as where is the neighborhood over which the filter h is defined, and (2) is the sigmoid function.

2.4.2 3D patch sampling.

The 3D patch sampling process converts our MRI image into a set of smaller 3D images. This is done by defining a patch size (in our case 5 × 5 × 5) and moving this cube across the image in all dimensions (step size is 1 voxel), extracting the contents of the cube at each point to create a 3D patch. This produces one such patch for every possible 5 × 5 × 5 cube that can fit within the original image—yielding approximately as many patches as there are voxels. Note that “adjacent” patches will overlap.

2.4.3 Using a sparse autoencoder to generate filters.

We then run a sparse autoencoder over the patches, to produce the filters required for convolution. To generate the filters for the convolutional neural network (CNN), we use a single layer unsupervised sparse autoencoder, that seeks to compress the 125-dimensional patches down to just 3-dimensions; see Fig 3. This produces 3 filters (based on the Encoding Weights); see left portion of Eq 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Autoencoder that tries to reconstruct the 5 × 5 × 5 patches (represented as a 125-D vector).

For each of the k = 3 filters, there is a link (“Encoding Weights” θ₁) from each input node to each hidden node, and then from each hidden node to each Output node (“Decoding Weights”, θ₂); we show only a subset for clarity.

https://doi.org/10.1371/journal.pone.0194856.g003

The sparsity aspect of the autoencoder encourages the network to learn different transformations in each of the nodes [19], which is meant to improve the robustness of the classifiers using the filters’ outputs as features. Formally, we seek the parameters [θ₁, θ₂, b₁, b₂] of the network (Fig 3) that minimizes: (3) where D ∈ ℜ^m×n is the data matrix, where each row is a data point x ∈ ℜⁿ, also h ∈ R^k is the hidden representation of the data (i.e., there are k nodes in the hidden layer), is reconstructed data, L(a,b) = ∑_i(a_i − b_i)² is the squared loss, σ(s) is sigmoid function (Eq 2), ρ is sparsity parameter (e.g., ρ = 0.05), is average activation of the j^th hidden neuron, and is the Kullback-Leibler divergence. We used (internal) 5-fold cross-validation to select the hyperparameters β, λ and ρ, and solve Eq 3 using the L-BFGS package (http://users.iems.northwestern.edu/~nocedal/lbfgsb.html).

2.4.4 Single layer unsupervised convolutional neural net.

After learning the k = 3 filters from the sparse autoencoder, the filters are then convolved throughout each entire MRI 3D image (size 79 x 95 x 68), producing a new 3D image (size 74 x 91 x 64) for each filter. (Fig 4 shows a 2D analogue of this.) As we use sigmoid activation functions σ(⋅), these features are a non-linear combination of voxels within the input image [20].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. (a) Median axial slice from MRI scan and (b) the result after convolution with a filter.

https://doi.org/10.1371/journal.pone.0194856.g004

We then use max-pooling after the convolution step. Here, for each volume (produced by each filter), we divide the entire volume into disjoint and exhaustive 5 × 5 × 5 sub-regions, and then compute the maximum value of each sub-region. This step reduces the sensitivity of the feature map to various distortions, for example, by introducing translational invariance for the features within the max-pooling region. We then collect these values into a single vector—here of size 3 × (15 × 18 × 13) = 10 530. We then describe each subject with a single vector that is the concatenation of all k of these size-10,530 vectors. This is the input to the base learning algorithm; see Fig 1.

2.4.5 Evaluating the model produced by LeFM_s.

To evaluate the LeFM_s algorithm, we used 5-fold CV on the training set with internal 5-fold CV to obtain hyperparameters that maximize CV accuracy, and then (after finding the apparent best parameters and model), applying that learned model to the holdout set. The base learner is an SVM.

We applied the resulting learned model to the holdout set, but found that it did not perform as well as state-of-the-art methods that also only use structural MRI features. This process was still relevant, as it identifies relevant features, which will be used within LeFM_SF; see Section 2.6 below.

2.5.1 Dimensionality reduction using 2-step PCA.

An fMRI scan for a single subject has a large number of voxels, which are captured for each of a number of time points. Because combining all the fMRI scans from all subjects for source separation would be intractable, it is important to reduce the number of features [21, 22] while preserving the important subject-level variations in the data. In order to reduce the computational load on group level analysis, we follow Calhoun et al. [21]’s 2-step principal component analysis (PCA) [21, 22], which applies PCA first to each subject separately and then to the all subjects together. We use the implementation of this procedure in the GIFT software package. This 2-step PCA procedure produces a set of components meant to capture the subject level variations and group level commonalities. Each component yields a “spatial map” [23]. This approach is based on the (reasonable) assumption that there are common spatial sources for a set of fMRI scans in a particular study, where subjects differ based on temporal weights of each source.

2.5.2 Subject level PCA.

As shown in Fig 6, we start with a T × V matrix, matrix describing each of the n subjects. T is the number of time points. V is the number of voxels. PCA [21] is used (within GIFT) to reduce this to a T_red × V matrix by selecting T_red largest eigenvalues, to capture 99% of the variance. (For details, see the GIFT documentation http://mialab.mrn.org/software/gift/documentation.html).

2.5.3 Subject concatenation.

To estimate a set of common components across all subjects, we vertically concatenate the PCA-rotated matrices produced for each subject by the above subject-specific PCA step (as shown in Fig 7). Note that vertical concatenation means we concatenate along the temporal direction (as opposed to spatial concatenation, which is not discussed here). Therefore, we are using source separation to decompose each fMRI scan into a set of spatial maps, where each map is derived from a time course defined by a PCA component. As seen in Fig 7, in this context, the inner product of voxel components and time series components reconstructs an approximation of the original fMRI images, so we can view this as a way of representing each voxel’s time series as a linear combination of more general, shared time series.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Temporal vs. spatial concatenation.

Each ‘fMRI data i” box (for i = 1..n) corresponds to one of the n subjects. This shows two ways we can concatenate these boxes from subjects in Fig 6 to produce X, Here, we use temporal concatenation.

https://doi.org/10.1371/journal.pone.0194856.g007

2.5.4 Group level PCA.

The previous steps produced the Concatenated Matrix Y_[(nxTred)xV] (3^rd column in Fig 6). We next run PCA on this matrix, and again included enough rows K to capture 99% of the variance.

2.5.5 Source separation.

Previous studies [21, 24] have explored brain regions that are strongly temporally coherent (co-activated during rest) using source separation techniques like PCA or independent component analysis (ICA). Here, PCA separates the fMRI brain scan into uncorrelated spatial maps or sources based on variations in the time, whereas ICA decomposes the brain fMRI scans into a common set of spatially independent components (a.k.a. sources or spatial maps) and their corresponding user-specific time courses. In this case, ICA assumes that the spatial maps have constant higher order statistics [21, 22, 25].

This paper compares four source separation methods, which are applied after 2-step PCA—including this simple version of ICA, as well as a more complex Non-stationary Spatial Sources Decomposition (NSD) technique—in terms of the accuracy of the models produced by the features obtained from each algorithm. To the best of our knowledge, this is the first study to fully investigate using NSD to produce predictive models from fMRI data.

Principal Component Analysis (PCA) for fMRI The simplest of our 4 approaches just uses the “Representative Matrix” X_[K×V] as our spatial maps, without any further modification (that is, the 4th step of Fig 6 is just the identity transformation). In this case, we use internal cross-validation to determine the appropriate number of components that optimizes the accuracy of the final model, instead of preserving 99% of the variance.

Kernel Principal Component Analysis (k-PCA) for fMRI Here, we use the Radial Basis Function kernel (RBF), with the associated matrix, (4) where x_i and x_j are i^th and j^th column of X respectively. After computing the kernel similarity matrix, we can obtain the eigenvalue/eigenvector pair of the kernel matrix; the projection onto these eigenvectors produces spatial maps (rows of S). Finally, the inner product of each patient’s fMRI scan on these maps results in the time components (columns of A_i). The number of eigenvectors retained for this final inner product is again selected through internal cross-validation.

Independent Component Analysis (ICA) The goal of the standard ICA algorithm is to find independent sources that maximize the log likelihood of the observed data, in the context of temporally concatenated fMRI data—i.e., given the X appearing in 4th column of Fig 6. We also suggest how one might interpret this process.

The standard max likelihood estimate (MLE) interpretation of ICA starts with an observed data matrix X where each row of X (xⁱ) is a single observation. In particular, this model assumes that each xⁱ is a realization of the random variable xⁱ where xⁱ is a linear combination of independent source random variables, {s^j}_{j = 1..K}. We also need a mixing matrix A that defines the contributions of each s^j in producing each xⁱ. Concretely, if we take any column of X (denoted x_j) we find the relationship (5) where W = A^† is pseudo-inverse of A and s_j is a column of S, the matrix whose rows sⁱ are each a realization of sⁱ. Note that A and S are never observed.

In our fMRI setting, shown in Fig 8, X and S each have V columns and each time point in the concatenated data xⁱ is a linear combination of our unknown independent spatial maps S = [s¹, …, s^K]. To help understand ICA applied to temporally concatenated fMRI data, think of the data as snapshots of a someone performing different gestures. For each snapshot, we observe the brain activation associated with a particular gesture (our xⁱ), which can be viewed as a combination of what would be the brain activations associated with more primitive movements S such as rotating a wrist or lifting an arm. In light of this example, ICA is being used to decompose neurological activity into more fundamental components to identify a disease state such as ADHD or Autism.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Temporal concatenation with ICA elaborated from Fig 7.

This corresponds to the ICA decomposition X ≈ AS, except here, for illustrative purposes, we have replaced X with Y (the concatenated fMRI data from Fig 6 before the Group-level PCA).

https://doi.org/10.1371/journal.pone.0194856.g008

Now that we have defined our setting and provided some intuition, we now describe how we use MLE to obtain the (unobserved) independent components of S. Our ultimate goal is to find an optimal inverse mixing matrix W* that maximizes the likelihood of the observed data. Here, we assume that the distribution of the h^th source s^h is given by a density p_h(s^h), and that the joint distribution of the sources are independent of one another; i.e.,

(6)

The probability of the observed signals for j^th observation (x_j) is given by (7) where is the transpose of the h^th column of W [26]. The W that maximizes the log-likelihood given the data X_[K×V], is (8) Note that this p_h(⋅) can be any non-Gaussian distribution [27]. Here, we use GIFT’s default implementation and parameters, which is (9) where each different source h = 1..K will have a different λ^(h). Note this does not depend on the location of the voxel. Given this W*, and the observed x_j for patient j, we then use Eq 5 to solve for the time course A_j.

Non-Stationary Spatial Sources Decomposition (NSD) As mentioned earlier, the probability density of each source is parameterized by a constant value that does not change with location of voxels. NSD generalizes this, by allowing the probability density function for each component to be parameterized by the location of the voxels [28, 29]. Here, we define the distribution for the h^th component S^h as (10) (Compare to Eq 9)

We think that each of these spatial sources should be non-stationary (i.e., have different λ_ijk values at different (i, j, k) locations) as:

1. One source (which corresponds to one row in S [see Fig 6]) may not have the same magnitude and variation throughout the whole brain due to in-homogeneous magnetic susceptibility that depends on the location of voxels in the scan. This differs from commonly used source separation models, which require that the sources have same variation for the whole brain scan.
2. The strength and variability within a particular source depends on the brain tissue type in a particular brain region. For example, the activation values in grey matter and white matter (which depends on the amount of oxygenated blood flow in that tissue) would be different.

In order to develop the theory for nonstationary source decomposition, consider a representative fMRI scan X_[K×V] ≈ A_[K×K] S_[K×V], where S is an K × V matrix of source components, where each row S_r,: (size 1 × V) specifies the contribution of voxel co-ordinates to r^th source. Each spatial activation refers to one row in S. Column h in matrix A will have the time courses for the corresponding h^th spatial component.

In mathematical terms, suppose we have K independent spatial maps (each corresponding to a row in S) and V observations (each corresponding to a 3D brain scan at a time point). Then, we can formulate the covariance matrix at location r = (i, j, k) as R_x(r) = 〈x(r) x(r)^T〉 = A D_s(r) A^T where R_x amd D_s(r) are of size K × K. Further we let x(r) = X(:, N(r)) where N(r) represents any suitably chosen region around position r for which the signals are assumed to have same higher order statistics. For our experiments, we have chosen N(r) to be a 4 × 4 × 4 bounding box (e.g., if the point r = [200, 100, 50], then this box is defined by corners [199, 99, 49] and [202, 102, 52]). We found, empirically, that increasing the bounding box severely degraded the performance of the model.

Assuming the sources are spatially non-stationary, and following [15, 30],

However, as we do not have a perfect estimate for R_x(r), we estimate the covariance matrix R_x(r) for some spatial interval. We denote the sample estimates as , and define the measurement error as

Suppose we have N samples of for r ∈ {r₁, r₂, .., r_N}, then we can estimate the parameters by and estimate the source components as

Then for each patient, we can calculate the time courses A_i from Eq 11. This method decorrelates the spatial sources at different regions of the brain, which is desirable for the reasons described above.

2.5.6 Complete LeFM_F process.

Fig 6 shows the steps. First, the SubjectLevelPCA reduces each subject’s [T × V] matrix Iⁱ to Yⁱ = U_iIⁱ where U_i is the T_red × T reduction matrix. Note that T_red is determined by retaining the independent components that preserve 99% of the variance (see the GIFT software package documentation for complete details on individual and group level PCA reductions). The second step, Concatenation, concatenates {Yⁱ}_{i = 1..n} over all of the n subjects, to obtain the matrix Y of size nT_red × V The third step, GroupLevelPCA, runs PCA on this Y to produce a new matrix X ≈ FY, where F is a K × nT_red reduction matrix. If we divide the F into n sub-blocks F_i, each of size K × T_red, then F_i corresponds to reduction matrix for subject i. We then consider various final steps: (1) Directly using the final F matrix as features. (2) Running kernel PCA on X to obtain a reduction matrix A. (3) Running (stationary) ICA on X to obtain the mixing matrix A. (4) Running NSD on X to obtain the mixing matrix A.

In these last three variants we are using different blind source separation algorithms to express X ≈ A × S, where S is a [C × V] matrix with C spatial maps covering V voxels. Hence, for subject i, (11) where , which is different for each patient (while all use the same S matrix). Here and are pseudo-inverses of U_i and F_i respectively. Thus, Xⁱ is approximated as A_i S.

3.2.1 5-fold cross validation accuracy.

To assess the predictive ability of a model trained only on fMRI features, we ran 5-fold cross validation (CV) over a varying number independent components. We then ran the version with the optimal number on the holdout set. Table 3 shows the results of this CV selection process for both the ADHD and Autism datasets. As shown, for the ADHD data, the model performed best with 45 independent components, achieving an accuracy of 0.6450 compared to the baseline of 0.50 with significance p = 1.83e-12 (via binomial test). For the Autism data, the model performed best with 45 components with a CV accuracy of 0.6225 relative to the baseline 0.5157 (p = 4.95e-10, via binomial test). Note that the 45 components for the ADHD data are unrelated to the 45 components for the Autism data. Fig 9 plots these values, and shows that NSD is better than the other feature extraction methods (PCA, kPCA and ICA).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. 5-fold CV results for classification accuracy (and standard deviation) using various numbers of independent components.

https://doi.org/10.1371/journal.pone.0194856.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. 5-fold cross validation results of classification.

https://doi.org/10.1371/journal.pone.0194856.g009

3.2.2 Model hold-out test accuracy.

Using the optimal parameters selected from cross validation, the resulting hold-out accuracy for ADHD was 0.6491 (baseline 0.5497) and was 0.6233 for ABIDE (baseline 0.5157). The specificity, sensitivity and Jstat for ADHD (resp., ABIDE) are 0.8191, 0.4416 and 0.2607 (respectively, 0.6768, 0.5533 and 0.2301). To our knowledge, our holdout accuracies achieved on ADHD-200 test data and the ABIDE test data are the best known when only using fMRI scans.

3.2.3 Leave-one-out accuracy comparison.

To compare our results with yet other systems, we also computed the leave-one-out accuracy of our ABIDE model; Table 4 shows that our LeFM_F was superior.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Leave-one-out results for autism classification (ABIDE) using only imaging data.

https://doi.org/10.1371/journal.pone.0194856.t004

Note that Abraham et al. [31] report an “inter-site” cross-validation accuracy of 66.8±5.4% on a version of the dataset containing 871 subjects (patients are ommited because they did not pass a visual quality inspection). Here, Abraham et al. [31] use a 16-fold cross validation, done by training on 15 of 16 sites and testing on the remaining site in each fold. By contrast, our model was trained on a smaller set of only 800 instances—recall we trained on only 70% of the original 1111 instances, and left the remaining 30% (311) for testing. It is likely that their reportedly high accuracy is due in part to having a larger CV dataset where samples were heavily scrutinized so as to exclude those that did not meet their criteria.

3.2.4 Spatial maps.

Spatial maps from ADHD-200 data:.

Figs 10 and 11 show the spatial maps found using multiple de-correlation—showing 9 axial slices for each. The component IC1 appears to be an artifact as it assigns very high values to all the voxels—i.e., as all the voxels in the brain are approximately of equal importance in this component, it is shared by all the voxels. We assume this is probably due to noise like motion, breathing and attention signals, which can modulate voxels throughout the brain [32]. We believe that IC3 also appears to be an artifact as it consists of regions from cerebrospinal fluids which means it might be a result of cardiac pulsatility artifacts [33]. Moreover, we found that removing these two artifacts did not change the performance of the model. The other components consist of some common default mode networks (regions that are shown to be active during rest [34]) and apparently some new resting state networks not previously identified in the literature. Some of the resting state networks found are consistent with [35]. For example, IC5 appears similar to the resting state network for peristriate area, and lateral and superior occipital gyrus, which are areas related to visual cortex and might represent spontaneous brain activities like day-dreaming [36, 37]. It is easy to connect this component to ADHD as ADHD patients are more likely to experience mind-wandering [38]. IC6 captured shared functional properties in the frontal and occipital lobe (responsible for planning and many areas of vision respectively). This area is very important for ADHD as well, as ADHD patients may suffer from lack of effective planning [39]. IC24 consists of pons (responsible for eye movement, sleep, and many other vegetal and automatic functions) regions and temporal lobe (for sensory processing, memory formation and higher order association processing). Their usefulness in prediction suggests that some of the physical signals captured by fMRI may also be indicative of a disease state.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Spatial maps (aka components) 1–25 for ADHD-200 using NSD.

Each 3×3 box shows 9 axial slices of one component. The colorbar is same as Fig 4.

https://doi.org/10.1371/journal.pone.0194856.g010

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Spatial maps or components 26–45 for ADHD-200 using NSD.

Each component is shown in a box and 9 axial slices are shown. The colorbar is same as Fig 4.

https://doi.org/10.1371/journal.pone.0194856.g011

Spatial maps from ABIDE data:.

Figs 12 and 13 present the spatial maps found using multiple de-correlation. This IC1, too, appears to be an artifact as it is shared by almost all the voxels [32]. IC3 also appears to be an artifact as it consists of regions from cerebrospinal fluids and is a result cardiac pulsatility artifacts [33]. A deeper investigation into the components shows multiple overlapping components. The components include visual areas (visual cortex, V1 and V2), partial overlapping with some default mode networks (PCC/precuneus, anterior cingulate cortex and frontal lobe) and motor networks. These components are informative: see the cross-validation, test accuracy on the ABIDE data for autism classification.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Spatial maps or components 1-25 for ABIDE using NSD.

Each component is shown in a box and 9 axial slices are shown. The colorbar is same as Fig 4.

https://doi.org/10.1371/journal.pone.0194856.g012

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 13. Spatial maps or components 26-45 for ABIDE using NSD.

Each component is shown in a box and 9 axial slices are shown. The colorbar is same as Fig 4.

https://doi.org/10.1371/journal.pone.0194856.g013