Building Models of Functional Interactions Among Brain Domains that Encode Varying Information Complexity: A Schizophrenia Case Study

2.1. Datasets and Pre-Processing

This study uses two independent datasets, the Function Biomedical Informatics Research Network (fBIRN) dataset (Keator et al., 2016) and the Center of Biomedical Research Excellence (COBRE) dataset (Aine et al., 2017). Subjects with large head motion (≥ 3deg and ≥ 3mm) during the scan and with functional data leading to bad full brain normalization were excluded. After applying this criteria for exclusion, the fBIRN dataset consisted of 160 healthy controls (HC) with age 19 – 59 years (mean 37.04 ± 10.86), 45/115 females/males and 151 subjects who had schizophrenia (SZ) with age 18–62 years (mean 38.77 ± 11.63), 36/115 females/males, whereas the COBRE dataset consisted of 89 healthy controls (HC) with age 18–65 years (mean 38.09 ± 11.67), 25/64 females/males and 68 SZ with age 19–65 years (mean 37.79 ± 14.45), 11/57 females/males. For avoiding confounding effects, HC and SZ subjects of both datasets were matched by age and gender (age: p = 0.1758 (fBIRN), 0.8874 (COBRE); gender: p = 0.3912 (fBIRN), 0.0794 (COBRE)).

Data collection for fBIRN dataset was done using 3-T Siemens Tim Trio scanners for six out of seven sites and 3-T General Electric Discovery MR750 scanner for one site. The same resting-state parameters were used across the scanners a standard gradient echo-planar imaging (EPI) sequence, repetition time (TR)/echo time (TE) = 2000/30 ms, voxel spacing size = 3.4375 × 3.4375 × 4 mm, slice gap = 1 mm, flip angle (FA) = 77°, field of view (FOV) = 220 × 220 mm, number of excitations (NEX) = 1, and number of volumes = 162. Participants had their eyes closed and were instructed to rest quietly during the scanner. The COBRE data was collected on a single site using a 3-T Siemens Trio scanner. For the COBRE data, a gradient-echo EPI sequence was used to acquire T2-weighted functional images with the following parameters: TE =29 ms, TR = 2000 ms, flip angle (FA) = 75°, slice thickness = 3.5 mm, slice gap = 1.05 mm, field of view 240 mm, matrix size = 64 × 64, voxel size = 3.75 mm × 3.75 mm × 4.55 mm and number of volumes = 149. For the duration of the scan, subjects were instructed to keep their eyes open and passively stare at a central cross.

For both the datasets, the statistical parametric mapping (SPM12, http://www.fil.ion.ucl.ac.uk/spm/) toolbox based on the Matlab 2016 platform was used to preprocess fMRI data. To ensure equilibrium of the signal and adaptation of the subjects to scanner noise, the first five scans were removed. The SPM toolbox was used for performing slice timing correction and rigid body head motion correction. Warping of the fMRI data into the standard Montreal Neurological Institute (MNI) space was done using an echo-planar imaging (EPI) template. The data were slightly resampled to 3 × 3 ×3 mm³ isotropic voxels. For smoothing the data, a Gaussian kernel with a full width at half maximum (FWHM) of 6 mm was used. This was followed by further feature extraction detailed in subsequent sections.

2.2. Connectivity Features and Subdomain Interactions

To estimate intrinsic connectivity networks from the fMRI data, the scICA approach described in (Du and Fan, 2013) with the Neuromark (Du et al., 2019a,b) template as reference maps was used. ICA was run for a model order of 100, out of which 53 consistent and reproducible independent non-artifactual components (ICs) were retained. The functional connectivity matrix was computed based on the static correlation between the time courses of pairs of two ICs. All the 53 ICs have been arranged into 7 distinct functional subdomains (disjoint exhaustive sets of ICs), namely the areas falling under the following subdomains: default mode network (DMN), visual (VIS), auditory (AU), cognitive control (CC), sensorimotor (SM), cerebellar (CB) and sub-cortical (SC). Based on this categorization, 28 subdomain interaction (SDI) features were created where each SDI represents the set of intra-network or inter-network connections between all possible pairs of ICs for a given pair functional subdomains or pair of functional subdomains respectively (See Figure 2). As an example, the features in the SDI named DMN-VIS refer to entries from the functional connectivity matrix (fMRI time-series correlations) corresponding to the connections between the DMN and VIS subdomain scICA components. Such connections are termed as inter-network connections. On the other hand, the intra-network connections comprise connections within the same domain, for example, the SDI named DMN-DMN would correspond to the connections where both of the scICA components belong to the DMN subdomain. Thus, connectivity features for the 28 SDIs consist of 21 sub-matrices for the inter-network connections and 7 sub-matrices for intra-network connections. Hence the SDI features are simply sub-matrices of the 53 × 53 functional connectivity matrix created using 53 scICA components.

2.3. Architecture Search Space

A multilayer perceptron (MLP) architecture with multiple branches was used with input layer consisting of the 28 SDI features from the 7 subdomains (sets of ICs) as described in subsection 2.2. At the input level, the architecture has 28 branches for each SDI i.e., a sub-matrix of the 53×53 static functional connectivity matrix. Each input layer is then fed into a variable number of fully connected hidden layers and a late fusion layer, finally followed by a fully connected layer (Figure 1b). In each branch, a constant factor of 0.1 was used to reduce the number of units in subsequent hidden layers, starting with the input layer. Optimizing the variation in the number of fully connected layers in the branch corresponding to each SDI can give insights into understanding the nature of information that is encoded by the component interactions for the subdomains involved in the SDI. SDIs which consistently require no or a small number of layers in an optimized architecture can be said to contain more linear or direct information for the task at hand (schizophrenia classification), whereas the SDIs which require a higher number of fully connected layers can be said to contain more complex and indirect information. For controlling the model complexity, the number of fully connected hidden layers in each branch were limited to vary between 0 to 2. In the exponential search space thus generated with a total of 3²⁸ possible configurations, a 28 length vector x can be used to represent each point, with each element x_i ∈ {0, 1, 2} denoting the number of fully connected layers in the branch for the i-th SDI (Figure 1b).

2.4. Tree-structured Parzen Estimator (TPE) for Hyper-Parameter Optimization

Optimizing hyper-parameters that span an exponential search space has been a well-studied problem in machine learning. In similar lines, Bergstra et al. (2011) had proposed the TPE algorithm which is used for hyper-parameter optimization in this work. TPE is a sequential model-based optimization (SMBO) framework which essentially estimates the conditional and marginal probabilities p(x|y) and p(y) to perform hyper-parameter optimization on hyper-parameters x and cost function y. After spanning an initial set of points selected randomly from the search space and computing the cost function for them, TPE utilizes this information to create a surrogate cost function. The surrogate function is then used to select the next set of points in the traversal sequence based on the estimates of the cost function. Towards this goal, two groups made up of the upper quartile and the lower quartile are created based on a splitting value y* of the cost function at the randomly selected points. The two probability density functions, g(x) and l(x) for the upper and lower quartiles of the cost function respectively is then defined as in Equation 1:

Having estimated l(x) and g(x), the subsequent iterations of the TPE algorithm involve the optimization of the Expected Improvement function defined in Equation 2:

The expected improvement function EI_y*(x) can be shown to be proportional to l(x)/g(x) (Bergstra et al., 2011). This means that if points are sampled from l(x) with high probability and from g(x) with lower probability, then the expected improvement is maximized, thus minimizing the cost function. In each iteration, the algorithm returns the point x* having the largest value of EI_y*(x) and the distributions l(x), g(x) are updated accordingly. An example showing how TPE works is provided in Figure 1c and indicates convergence towards the true optimal on a quadratic cost function (x − 1)² optimized for a real valued parameter x.

2.5. Analysis of Learned Models

Using features from each SDI as inputs, the MLP architecture with 28 branches as described in subsection 2.3 was optimized using the TPE implementation in the HyperOpt python library (Bergstra et al., 2013a). We used a repeated random sub-sampling cross-validation procedure for running the TPE algorithm for 50 repetitions on the functional connectivity data, resulting in a set of optimized architecture vectors ${\left\{{\mathbf{x}}^{(r)}\right\}}_{r=1}^{50}$. The mean across repetitions of the validation accuracy has been plotted against TPE iterations in Figure 5.

The final prediction model, represented by the vector x*, was created using the most frequent number of hidden layers across repetitions for each SDI in the TPE-optimized architectures. The final architecture vector x* is defined as, ${\left\{x_i^{*}\right\}}_{i=1}^{28}={\left\{argma{x}_{l\in \{0,1,2\}}{\sum }_{r=1}^{50}\mathbb{1}\left\{x_i^{(r)}=l\right\}\right\}}_{i=1}^{28}$, where x^(r) represents the TPE-optimized architecture for the r^th repetition and $\mathbb{1}\{.\}$ is the indicator function. The multi-branched MLP architecture generated using the final architecture vector x* was run on a held-out test data for 50 repetitions, resulting in 50 test accuracy values. It should be noted that while the random splits for training and testing were different for each of the 50 repetitions, the same random splits were used across all methods for a given repetition to avoid any unwanted differences in results due to different training data.

2.6. Performance Comparison with Classification Methods

To study the performance of the TPE-optimized architecture, a comparison was made with baseline methods run for 50 repetitions on held-out test data with 20% of the total number of samples. These methods included standard machine learning models, neural network based models, and multi-branched architectures with as well as without flexible branch-depth. The set of methods and corresponding parameters used are detailed subsequently.

3.1. scICA components

After running the scICA framework on the fMRI data, 53 non-artifactual, reproducible independent components (ICs) were obtained. Table 1 shows the peak coordinates for the components along with the brain region corresponding to the component. Following this, static correlations between the time courses of these 53 components were computed resulting in a 53×53 functional connectivity matrix. The 53 ICs were assigned to 7 functional subdomains or sets of components (Table 1,Figure 4). 28 pair-wise subdomain interactions (SDIs) resulting from the 7 subdomains were constructed. The SDIs are defined by the sets of intra-network or inter-network connections between all possible pairs of subdomains; thus, there were 21 SDIs corresponding to inter-network connections and 7 for the intra-network connections (Figure 2).

3.2. Performance using TPE

Starting with an initial set of randomly selected points in the hyper-parameter search space, the TPE optimization procedure subsequently selects new points to traverse using the expected improvement metric and the surrogate cost function which is computationally faster (subsection 2.4). Figure 5 shows the validation accuracy plotted against time (iterations). The validation accuracy, which is the objective function being optimized, increases with the TPE iterations. The procedure returns an optimal architecture for each of the 50 repetitions. The most frequently occurring number of hidden layers for each SDI across repetitions are used to define the final architecture. The test accuracy values on held-out data were obtained for 50 independent repetitions by using the final architectures constructed after running the TPE optimization algorithm (Figure 6). The final TPE models reported a slightly higher prediction accuracy on held-out test data than the baseline standard machine learning methods (SVM, logistic regression and random forest classifier). Specifically, we observed a mean prediction accuracy of 0.81 for fBIRN and 0.78 for COBRE for the proposed TPE-optimized architecture. For both the datasets, TPE also performs significantly better (p < 0.05) than the non-branched neural network based methods (MLP, FNN, BCNN) as well as multi-branch architectures with uniform depth in each branch (UNIF0, UNIF1 and UNIF2). Interestingly, the TPE optimization procedure also results in significantly better (p < 0.05) test accuracy than other optimization techniques for branched architectures with variable branch-depth (RNDS, GRDS). The results have been summarized in Table 2 as well as Figure 6. The results show that treating features from certain subdomains differently in terms of the complexity of architecture required can result in superior prediction performance. In fact, this improvement over the baseline methods also allows for a meaningful analysis of variability in the nature of information carried by the various subdomain interactions in the data.

3.3. Feature Stability

While the TPE algorithm gives the best performance, the next step would be to analyze the subdomain interactions (SDIs) requiring deeper or shallower hidden layers in the multi-branched MLP architecture. However, before checking for this characteristic variation in the SDIs, it is relevant to check whether the features learned by the optimized architecture have consistent discriminative power for each SDI in terms of their importance towards prediction. For this purpose, impurity-based feature prediction power (Louppe et al., 2013), also known as Mean Decrease Impurity (MDI), was computed by running random forest classifier on the parameters learned in the late fusion layer (Figure 1b) of the optimized TPE architecture. These were compared with the impurity-based prediction power by running random forest classifier on the features in the input layer of the network (i.e., the functional connectivity features). The comparison was made by obtaining the prediction power vectors for both the aforementioned cases and computing the cumulative prediction power for elements belonging to each subdomain interaction (SDI). The cumulative prediction power of each SDI for the prediction task, averaged over 50 repetitions of the algorithm on held-out test data is plotted in Figure 7. The results indicate that the prediction power of SDI parameters learned in the late fusion step of the final architecture is highly correlated (0.9 for fBIRN and 0.81 for COBRE) with the prediction power of the same SDI features in the functional connectivity input features. This observation suggests that the TPE-optimized architecture learned appropriate brain representations as the SDI features retained similar properties in terms of the prediction accuracy on the schizophrenia classification task as well as the importance of each SDI in the prediction. With these properties being similar, the TPE procedure that optimizes based on the depth by each SDI in the multi-branched architecture, can be said to additionally provide insights into the complexity of the information stored in the SDIs. The results from the analysis of this variation in the nature of information across the SDIs (as defined in subsection 2.5) is done in the subsequent subsection.

3.4. Analyzing Optimal Models

After having ensured that the TPE procedure for optimizing the multi-branched MLP architecture is consistent in terms of prediction accuracy on held-out test data and also shows a similar trend in the importance of SDIs, an analysis on the variation in depth required by various SDIs in the final architecture returned by TPE was done as described in subsection 2.5. This was done by considering the most frequently occurring number of fully connected layers across all repetitions of TPE for each SDI. The relationship between 7 subdomains (networks) in terms of the number of fully connected layers needed for optimal decoding of information in each pair of networks (SDI) in the final architecture is shown in Figure 8. A notable observation here is that in the final architecture of both the fBIRN and COBRE datasets, the number of fully connected layers required by each SDI is the same for 19 out of the 28 SDIs (Figure 9), which indicates a very similar common pattern across the two independent datasets in terms of complexity of information in subdomain interactions. Note that the number of such common values across two datasets can be modelled using a random variable which is distributed according to the binomial distribution B(k, n, p), where k is the number of successes (=19), n is the number of trials (=28) and p is the probability of a success, which in this case is the probability of having the same depth in both datasets for a particular SDI, given by $\frac{3}{3^2}=\frac{1}{3}$. A binomial test was done to check whether the observed number of common values (= 19) is significantly different from the expected number (= 28/3) in the case of the binomial distribution. This test indicates that the result is significant with a p-value of 2.09 × 10⁻⁴.

3.5. Biological Interpretation

Using the TPE framework, functional connectivity features between brain regions belonging to various functional subdomains were used to analyze them in terms of the complexity involved for uncovering predictive information for schizophrenia classification. While most deep learning architectures do not allow for analyzing specific region-wise functional connectivity input features due to the mixing up in subsequent hidden layers, the TPE framework mixes the region-level features only within the branch corresponding to the functional subdomains that the regions belong to. Thus, the results, upon optimization for complexity required by various subdomain interactions, can be interpreted in terms of the functional connectivity subdomains. Figure 9 summarizes these findings in terms of the 7 functional connectivity subdomains involving brain regions belonging to default mode network (DMN), visual (VIS), auditory (AU), cognitive control (CC), sensorimotor (SM), cerebellar (CB) and sub-cortical (SC) areas.

With schizophrenia being a cognitive disorder involving audio-visual hallucinations occurring as if being sensed in real even while the subject is not actively engaged in a task, it is interesting to note that the same is reflected by the functional subdomains involved in these functions (CC, DMN, SM, VIS, AUD) requiring higher complexity models. Cognitive control deficits are well known to be prevalent in patients diagnosed with schizophrenia, with cognitive control known to regulate a lot of other cognitive systems including behavioral disorganization (Lesh et al., 2011). According to the cognitive control model, failures to retrieve episodic memories are mainly mediated through frontal areas involved in cognitive control (Ragland et al., 2009). Additionally, previous studies have revealed that impaired connectivity in the DMN is associated with reduction in capacity to effectively retrieve episodic memories and process self referential information (Dunn et al., 2014; Kim, 2010). Thus, the associations between CC and DMN domains are evident from these observations. The connectivity of the DMN regions with the SM and VIS regions is also known to be affected in schizophrenia (Wang et al., 2015). The role of connectivity between regions belonging to the SM, VIS and CC subdomains is also known in the case of schizophrenia (Kaufmann et al., 2015; Javitt and Freedman, 2015). Previous studies have found significant differences in the connectivity within the SM/VIS regions as well as their connectivity with the CC regions (Kaufmann et al., 2015; Javitt and Freedman, 2015; Butler and Javitt, 2005; Gaebler et al., 2015). Many studies have shown that inter-subdomain connectivity is affected in schizophrenia for the CC, VIS, SM, DMN and AUD areas (Javitt, 2009; Kim et al., 2009). This is in line with the observation that the inter-subdomain connections between regions from all these subdomains require more complexity (Figure 9).