Three-way parallel group independent component analysis: Fusion of spatial and spatiotemporal magnetic resonance imaging data.

Advances in imaging acquisition techniques allow multiple imaging modalities to be collected from the same subject. Each individual modality offers limited yet unique views of the functional, structural, or dynamic temporal features of the brain. Multimodal fusion provides effective ways to leverage these complementary perspectives from multiple modalities. However, the majority of current multimodal fusion approaches involving functional magnetic resonance imaging (fMRI) are limited to 3D feature summaries that do not incorporate its rich temporal information. Thus, we propose a novel three-way parallel group independent component analysis (pGICA) fusion method that incorporates the first-level 4D fMRI data (temporal information included) by parallelizing group ICA into parallel ICA via a unified optimization framework. A new variability matrix was defined to capture subject-wise functional variability and then link it to the mixing matrices of the other two modalities. Simulation results show that the three-way pGICA provides highly accurate cross-modality linkage estimation under both weakly and strongly correlated conditions, as well as comparable source estimation under different noise levels. Results using real brain imaging data identified one linked functional-structural-diffusion component associated to differences between schizophrenia and controls. This was replicated in an independent cohort, and the identified components were also correlated with major cognitive domains. Functional network connectivity revealed visual-subcortical and default mode-cerebellum pairs that discriminate between schizophrenia and controls. Overall, both simulation and real data results support the use of three-way pGICA to identify multimodal spatiotemporal links and to pursue the study of brain disorders under a single unifying multimodal framework.

INTRODUCTION

Advances in brain imaging acquisition technologies have enabled the collection of more than one imaging modality from each subject, for example, magnetoencephalography (MEG), positron emission tomography, and magnetic resonance imaging (MRI), which includes structural MRI (sMRI), functional MRI (fMRI), and diffusion MRI (dMRI). Each of them provides limited yet unique views of the functional, structural or dynamic temporal features of the brain. For example, sMRI provides higher spatial resolution on specific brain tissue types (e.g., gray matter, white matter, cerebrospinal fluid) but lacks temporal information, while MEG measures brain electrical activity with high temporal resolution, fair spatial resolution, but no tissue specificity. The limitations of each individual modality and their evident mutual complementarities suggest that advances in multimodal data fusion could be key to maximizing clinical and scientific benefit from large multimodal imaging studies (Baltrusaitis, Ahuja, & Morency, 2019; X. Zhang, Pan, et al., 2020).

Several multimodal fusion methods have been proposed with the purpose of taking full advantage of complementary brain imaging modalities (Sui, Adali, Yu, Chen, & Calhoun, 2012). A representative two‐modality (two‐way) fusion approach named parallel independent component analysis (pICA) (Liu et al., 2009) maximizes both intramodality independence and intermodality correlation to identify linked independent component (IC) pairs between two modalities. An important feature of pICA is that each modality is characterized by its own mixing matrix. For three‐way fusion, joint ICA (jICA) (Calhoun et al., 2006) maximizes the independence among the concatenated multimodal features but generates the same mixing matrix for all modalities. Multiset canonical correlation analysis (mCCA) (Li Yiou, Wang, & Calhoun, 2009) maximizes intersubject covariation across any number of modalities but without an independence constraint, and limited to orthogonal features, which is extended to grouped sparse CCA (X. Zhang, Pan, et al., 2020). The mCCA + jICA (Sui et al., 2011) method maximizes intramodality independence and intersubject covariation by combining mCCA and jICA using a two‐step process, where mCCA acts as preprocessing for jICA. Another two‐stage fusion is multiple regression (De Martino et al., 2010), in which source separation was performed first, then followed by a cross‐modality regression. Linked ICA (Groves, Beckmann, Smith, & Woolrich, 2011) is based on Bayesian framework, which is extend to identify phenotypes from large population sample (W. K. Gong, Beckmann, & Smith, 2021) and to extract spatially and temporally shared components of complex‐valued fMRI by shift‐invariant canonical polyadic decomposition (X. F. Gong, Lin, Cong, & De Lathauwer, 2018; Kuang et al., 2020). Meanwhile, the three‐way pICA (Vergara et al., 2014) method maximizes intramodality independence and intermodality correlation under a unified optimization framework. Another multiway fusion method, multidataset independent subspace analysis (Silva, Plis, Adali, Pattichis, & Calhoun, 2020) accomplishes the fusion aim from a subspace independence perspective.

While the aforementioned methods use an unsupervised fusion paradigm, supervised models also exist, offering alternative models for hypothesis testing. The pICA with reference (Chen et al., 2013) imposes an additional constraint upon spatial maps to minimize the distance between an IC and the reference. The reference‐based mCCAR + jICA (Qi, Calhoun, et al., 2018) method adds a correlation constraint on the mixing matrix that can identify underlying co‐varying patterns that correlate with the reference (Qi, Yang, et al., 2018; Qi et al., 2020; Qi, Schumann, et al., 2021). Supervised big FMRIB's linked ICA can use multiple nonimage derived phenotypes to supervise multimodal MRI fusion (Gong Weikang, Ying‐Qiu, Smith Stephen, & Beckmann Christian, 2021). Importantly, all the fusion methods mentioned so far are restricted to second‐level fMRI features (e.g., 3D regional homogeneity or connectivity maps, subjects × imaging variables), rather than first‐level 4D fMRI (subjects × voxels × time points) (Plis et al., 2018). Under this 3D framework, the valuable temporal dynamic information in fMRI data cannot be fully utilized.

The temporal variation conveys important information in fMRI signal (Yan et al., 2019). However, most multimodal methods rely on preprocessing fMRI with dimensionality reduction approaches that collapse the temporal information and yield 3D images. Reducing 4D fMRI to 3D spatial data in such manner before data fusion does not guarantee a full, efficient utilization of the temporal information during the fusion step. Although there are methods incorporating the first level 4D fMRI in single subject case which assume concurrent EEG/fMRI data (Martinez‐Montes, Valdes‐Sosa, Miwakeichi, Goldman, & Cohen, 2004; Van Eyndhoven et al., 2021) and multisubject case (Chatzichristos, Davies, Escudero, Kofidis, & Theodoridis, 2018; Christos Chatzichristos, De Lathauwer, Theodoridis, & Van Huel, 2020; Jonmohamadi et al., 2020), most of them are limited to fusion of EEG and fMRI. We previously proposed a two‐way parallel GICA + ICA (Qi et al., 2019) to incorporate first‐level 4D fMRI data by parallelizing group ICA (GICA) and ICA, aiming to enable direct fusion of first‐level fMRI features, but it was limited to two modalities. In this study, we extended parallel GICA + ICA to enable three‐way fusion by integrating GICA into a three‐way pICA framework.

In addition to denoising the BOLD fMRI signal in preprocessing (Churchill, Raamana, Spring, & Strother, 2017), ICA can further isolate artifact components from fMRI (Du et al., 2016; X. F. Gong, Wang, & Lin, 2015; Niu et al., 2021). This is particularly noticeable in GICA (Calhoun, Adali, Pearlson, & Pekar, 2001), which deal with 4D fMRI from multiple subjects and not only extracts subject‐ and group‐level ICs but also their time courses. In contrast, three‐way pICA identifies linked component triangles among all three modalities by simultaneously maximizing both intermodality association (at the subject expression level) and intramodality independence (at the group component level). Building on this three‐way framework, we exploit GICA's success with temporal information by letting it take the place of the customary ICA of 3D fMRI data from traditional three‐way pICA. Fittingly, our approach enables analysis of temporal information included 4D fMRI together with the other two 3D modality features, in a single three‐way pICA setting, thus keeping the original properties of both GICA and three‐way pICA.

In order to combine GICA and pICA within a single model, a variability matrix that calculate the L2‐distance between subject‐specific and group maps was defined. We then maximize the correlation among variability matrices estimated from each of the three modalities in a three‐way parallel ICA decomposition. The summarized variability matrices of each modality portray how much the subject‐specific components different from group‐common maps. In turn, direct correlations can be measured between every modality pair to achieve data fusion of 4D fMRI (subjects × voxels × time points), plus two other 3D modality features, which in this work are gray matter volume (GMV) from sMRI and fractional anisotropy (FA) from dMRI (both subjects × voxels). Based on this intuition, a novel three‐way parallel GICA (three‐way pGICA) method that leverages the temporal information from 4D fMRI ensues, enabling the identification of a linked fMRI–sMRI–dMRI component triangle was proposed.

Notably, the fMRI decomposition obtained with the proposed three‐way pGICA is imbued with knowledge from the other two modalities and, naturally, differs from that obtained with simple GICA. Therefore, we can also evaluate its functional network connectivity (FNC) (Damaraju et al., 2014; Gonzalez‐Castillo & Bandettini, 2018; Karahanoglu & Van De Ville, 2015; Preti, Bolton, & Van De Ville, 2017; Y. Zhang, Zhang, et al., 2020) to identify abnormal cross‐network FNCs that are also related with diffusion and structural covariations. In this study, we verified the effectiveness of three‐way pGICA in both simulation and human brain data. The Function Biomedical Informatics Research Network (fBIRN) dataset (N = 254) (Keator et al., 2016) was used as a discovery and the Center for Biomedical Research Excellence (Jorge Nocedal, 1999) (COBRE, N = 89) was an independent replication. Results show that three‐way pGICA can generate a linked fMRI–sMRI–dMRI component triangle in both simulation and human brain data that can be validated in an independent cohort. These results support the use of three‐way pGICA for stable identification of reliable spatiotemporal intermodality linkage among three modalities to study brain disorders under a single unifying multimodal framework.

METHODS AND MATERIALS

GICA of fMRI

represents concatenated 4D fMRI over subjects () in the temporal direction (), where is the data matrix of subject , represents the number of voxels, and has dimensions of . The mean was removed from each subject data before principal component analysis (PCA). PCA is used to reduce subject‐ and group‐level dimensions for fMRI. Let be the reduced data matrix () of subject , where is whitening matrix () obtained by PCA decomposition at subject‐level (the notation indicates that is the pseudoinverse of the rank subject dewhitening matrix , based on the top eigenvalues of the eigenvalue decomposition of the subject‐specific covariance matrix ), and represents the rank from PCA decomposition. Then, the PCA‐reduced data is concatenated over subjects along the temporal dimension. After that, a group‐level PCA is performed to further reduce the concatenated group data to the number of components (), as in (1): where is the group‐level whitening matrix generated from group‐level PCA ( is the pseudoinverse of the rank group dewhitening matrix ), and is the reduced fMRI matrix.

Following ICA decomposition of (details provided further below), we get , where is the estimated aggregate group‐level fMRI IC matrix (one component per row) and is the estimated mixing matrix. Substituting into (1) and multiplying (group‐level dewhitening matrix) on both sides, we obtain the rank approximation: The group dewhitening matrix can be partitioned by subject, as in the following equation: Taking each partition , we can write the following least‐squares approximation for subject‐specific ICs: where consists of spatial fMRI ICs specific to subject and indicates the pseudoinverse. Likewise, multiplying two sides of (4) by (subject‐level dewhitening matrix) yields the final data reconstruction approximation for subject by least squares: which reveals the final ICA decomposition of by means of its least squares approximation (). The matrix contains subject‐specific ICs, while the mixing matrix contains the corresponding subject‐specific time courses. These are the classic linear back‐reconstruction equations from GICA1 (Calhoun et al., 2001).

Three‐way parallel GICA

In order to take advantage of the first‐level temporal fMRI dynamics in a combined GICA and pICA framework, we define a variability matrix () for three‐way parallel fusion that captures component‐wise subject‐level functional variability to link with the other two modalities (Qi et al., 2019; Qi, Plis, et al., 2021). Specifically, each element is defined for the ith subject as the L2‐distance between the jth group‐level IC () and the jth subject‐level IC (), where , . Naturally, the dimensionality of is equal to and (subject expression profiles) from sMRI and dMRI, respectively. Following the three‐way pICA framework (8), we maximize the independence within each modality while also adding constraint terms that maximize the modality linkage among functional ( and structural (, ) variabilities (Figure 1c).

FIGURE 1

Framework of three‐way parallel group independent component analysis (GICA) fusion method. (a) Preprocessed first‐level spatiotemporal 4D functional magnetic resonance imaging (fMRI) (), and second‐level voxelwise gray matter volume (GMV) from 3D structural MRI (sMRI) () and voxelwise FA from 3D diffusion MRI (dMRI) (). (b) Feature matrix for each modality. (c) Three‐way parallel GICA (pGICA) that maximize the independence for all the three modalities based on independent component analysis (ICA) and GICA portions individually, and maximize the correlation between subject expression profiles of sMRI and dMRI from ICA and the variability matrix of fMRI from GICA. (d,e) Group‐level and subject‐level variability matrix and spatial components generated from GICA and ICA portions. (f) Functional network connectivity (FNC) analysis for time courses of fMRI modality. GPCA: group‐level PCA; GICA: group ICA; : fMRI dewhitening matrix from GPCA; : fMRI dewhitening matrix generated from subject‐level PCA (sPCA); and : sMRI and dMRI dewhitening matrix from PCA, respectively; : jth group component; : jth component of subject ; : variability matrix for fMRI; and : mixing matrix for sMRI and dMRI. , , and represent the Infomax gradient. , , and are the between‐modality linkage‐regularized gradient

The intersubject functional variability matrix is estimated as where and is the jth row of . is normalized as to avoid the component scaling.

Suppose that the matrix and the matrix represent the sMRI and dMRI data, respectively ( and represent the number of voxels in each modality). Akin to fMRI, the mean was also removed from sMRI () and dMRI () modality. A PCA reduction step ( for sMRI, for dMRI) followed by ICA decomposition yields the following for the remaining modalities (assuming the subject‐level data accounts for only one row of or ): where and are the group‐level dewhitening matrices, and are the and reduced matrices for sMRI and dMRI, respectively, and likewise for the corresponding and mixing matrices. Here, () and () are the modality‐specific IC matrices, each with components.

We then formulate the cost function of the proposed three‐way parallel GICA method as: where and

with computed element‐wise (), and are in dimension of . is the differential entropy, is the expected value, and is the Pearson correlation coefficient. are the IC indices for the components selected (the top correlated component pair were selected) in each optimization iteration. is the probability density function of and represents the bias vector for each modality.

Note that the objective function (9) looks like the three‐way pICA (Vergara et al., 2014); however, represent GICA portion on 4D fMRI, but not the ICA on 3D fMRI. Moreover, we redefined (variability matrix as in Equation (7)) for fMRI to link with sMRI and dMRI. This estimates how much the subject component different from the group‐common component. Thus, whether this difference is associated with the other two modalities can be investigated. While in traditional three‐way pICA, the three modalities are linked by maximizing the correlation among mixing matrices generated from each ICA portion. The first three terms in (9) are solved parallel by the Infomax (Makeig, Bell, Jung, & Sejnowski, 1996). The last three regularization terms (correlation among selected columns of the variability matrices, based on selected ICs) are optimized using the steepest ascent method. Finally, based on the definitions given above, we obtain the following updating rules.

For the first three terms (major updates for , , and , using the relative gradient (Amari, 1998; Zarzoso & Hyvärinen, 2010)): and correspondingly for the bias weights: where is the learning rate, which is annealed periodically. For the last three terms (minor updates for , , and based on , , and ): where and represents the mean of the mth column of (likewise for below), where

and where

Following the updates on and , and are updated as , and . Thus, the above updated on and eventually becomes an update on and , due to the relationship between and (, and ). is the learning rate for terms of fMRI, sMRI, dMRI and is the step size of terms. In summary, there are two different procedures to update , , and , one is from the regular Infomax framework and the other is from the intermodality correlation regularizer. We adaptively tune the values of and in Equations ((13), (14), (15), (16), (17)) to control the weight in the cross‐modality correlation regularizer. The criterion is that when the maximum value in is larger than (the predefined maximum weight), then annealing is triggered with the update (likewise for ) to prevent values in from blowing up. All in all, our approach adds a further small adjustment on , , and by maximizing the cross‐modality correlations. Hence, the proposed three‐way pGICA is using an “alternating” approach to optimize (9), that is, it iteratively alternates between Infomax and Correlation optimization. Note that and (corresponding to terms) are updated in mini‐batches. Following, at the end of each epoch, , , and (corresponding to terms) are evaluated once in a single full batch update. The number of mini‐batches is the same for all modalities and the stopping criteria is whether all modalities converge or meets the maximum ().

Simulation

In this simulation, a comparison among the proposed three‐way pGICA, separate GICA/ICA, and separate GICA/pICA was conducted to estimate the ability of source separation accuracy and cross‐modality linkage detection. Eight nonoverlapping spatial brain maps were generated by the simTB toolbox (Erhardt, Allen, Wei, Eichele, & Calhoun, 2012) (https://trendscenter.org/software/SimTB) for fMRI and sMRI. Eight fiber bundles from the Johns Hopkins University white matter atlas were selected to simulate dMRI. Thus, the true sources for fMRI in dimension of , for sMRI in dimension of , and for dMRI in dimension of were generated. The number of time points was set to 100. time course matrices were generated from simTB. Mixing matrices (sMRI) and (dMRI) were normally distributed elements in a size of . FMRI variability matrix is generated by the distance between subject‐ and group‐level brain maps. The simulated connections among the three modalities were carefully designed by selecting one column (the third, fifth, and eighth columns for fMRI, sMRI, and dMRI, respectively) from each of , , and to be correlated under certain noise conditions. Therefore, the simulated data contains only one linkage triangle, and all other components are not linked.

Here, we used two ways to design the simulation. One is fixing noise: under peak signal‐to‐noise ratio (PSNR) = 5, the linkage between modalities were varied from weak (r = .2) to strong (r = .9) for both fMRI with sMRI and fMRI with dMRI pairs, with a fixed strong correlation between sMRI and dMRI (r = .8). The other one is fixing modality linkages: r = .7 between fMRI and sMRI, r = .18 (not statistically significant at uncorrected p > .05) between sMRI and dMRI, and r = .6 between fMRI and dMRI, and changing noise levels to PSNR values from −10 to 17 dB. The final observation data matrices result from the linear combination of sources and mixing matrices plus a noise term: . For fMRI, each is generated by . This resulted in 100 simulated subjects, with 10,000; 22,500; and 40,000 voxels for fMRI, sMRI, and dMRI, respectively.

Real brain data

In real data application, two independent cohorts were used. For the discovery dataset we used fBIRN (Keator et al., 2016), including 123 SZ patients (38.8 ± 11.7, 30F/93M) and 131 gender‐age matched healthy controls (HCs) (36.8 ± 10.9, 39F/92M). FBIRN demographic and cognitive scores are available in Supplementary Table 1. All participants are adults between 18 and 60 years. DSM‐IV (SCID) (M. B. First, Spitzer, Gibbon, & Williams, 2002) was used to diagnose SZ patients. Cognition was measured by the Computerized Multiphasic Interactive Neurocognitive System (CMINDS) (van Erp et al., 2015). There is no gender (p = .30) or age (p = .17) differences between HC and SZ for the discovery fBIRN.

We also included a validation dataset called COBRE (Aine et al., 2017) consisting of 47 SZs (39.6 ± 13.1, 12F/35M) and 42 gender‐age matched HCs (37.0 ± 11.8, 10F/32M). Cognition was estimated by the Measurement and Treatment Research to Improve Cognition in Schizophrenia Consensus Cognitive Battery (MCCB) (Green, Kern, & Heaton, 2004). COBRE demographic and cognitive scores are shown in Supplementary Table 2. There were no gender (p = .85) or age (p = .32) differences in the COBRE data. Written informed consent was obtained for both fBIRN and COBRE which is approved by the Institutional Review Boards. Resting fMRI, sMRI, and dMRI were collected for fBIRN and COBRE. Preprocessing steps and imaging parameters are in Supplementary “Imaging preprocessing.” sMRI GMV and dMRI FA were used as features for the fusion input.

RESULTS

Results of simulation

The proposed three‐way pGICA (red) was compared with separate GICA/ICA (green), separate GICA/pICA (blue) on simulated data to test the ability in detecting between modality linkages. Figure 2a–c shows the results for detecting intermodality association under fixed noise levels PSNR = 5 for different methods. The black lines in Figure 2a–c represent the true modality linkage (fMRI–sMRI and fMRI–dMRI: [0.2–0.9]; sMRI–dMRI: 0.8). Results show that three‐way pGICA can get relatively high cross‐modality linkage estimation comparing with separate GICA/ICA, and separate GICA/pICA without fusion. Separate GICA/ICA gets the worst linkage estimation (Figure 2c), since there is no linkage regularization between sMRI and dMRI in separate GICA/ICA without fusion. Due to the lack of regularization between fMRI and the other two modalities, the linkage estimation accuracy is lower for fMRI‐sMRI (Figure 2a) and fMRI‐dMRI (Figure 2c) than three‐way pGICA when comparing separate GICA/ICA, separate GICA/pICA. Figure 2d–f shows the source accuracy estimation for fMRI, sMRI, and dMRI under changing modality linkages. The correlation between the true sources and the estimated components was used as source estimation accuracy. It is evident that the proposed three‐way pGICA can achieve comparable source estimation accuracy relative to separate GICA/pICA, and relative to separate GICA/ICA. This means that without losing the source separation ability in ICA portions, three‐way pGICA can get more accurate modality linkage estimation, indicating the benefit of the proposed method. Collectively, these results show that three‐way pGICA provides accurate intermodality linkage detection under both weak and strong correlations with comparable source decomposition.

FIGURE 2

Simulated three‐way data fusion comparison between three‐way parallel group independent component analysis (GICA) (red), separate GICA/independent component analysis (ICA) (green), and separate GICA/parallel ICA (pICA) (blue). (a–c) Modality linkage estimation under fixed noise (peak signal‐to‐noise ratio [PSNR] = 5) with different level of associations (r = .2–.9). (d–f) Source estimation under different modality linkage. The black lines in (a–c) represent the true simulated modality linkage

We also compared three‐way pGICA with its alternatives when fixing the modality linkage but changing noise levels, as shown in Figure 3. It is clear that three‐way pGICA can achieve higher estimation accuracy of between modality linkages (Figure 3a–c), and comparable estimation accuracy for source separation (Figure 3d–f) under different noise levels. As expected, estimation accuracy decreases as noise level increases.

FIGURE 3

Comparison of modality linkage and source separation under different noise levels with fixed intermodality correlation among three‐way parallel group independent component analysis (GICA) (red), separate GICA/independent component analysis (ICA) (green), and separate GICA/parallel ICA (pICA) (blue). The black lines in (a–c) represent the true simulated modality linkage

Results on real data application

Linked fMRI–sMRI–dMRI triangle

In the fBIRN data application, the preprocessed images from resting‐state fMRI, GMV from sMRI and FA from dMRI were used as feature input for the proposed three‐way pGICA to identify three‐way linked fMRI–sMRI–dMRI components. The component number at subject‐level () was set to 80 given the higher dimensionality of the fMRI data, and the group‐level component number () was set to 20 based on the MDL criterion (Li, Adali, & Calhoun, 2007) for the fMRI modality. Following the three‐way pGICA decomposition, two‐sample t tests were applied on the estimated mixing matrix within each modality to compare group differences between patients of SZ and HC within the identified component triangle.

Among the 20 generated ICs, the seventh fMRI IC, the sixth sMRI IC, and second dMRI IC (Figure 4a) were the linked components triangle (r = .25, p = 4.0e‐05* between fMRI_IC7 and GM_IC6; r = .30, p = 9.1e‐07* between fMRI_IC7 and FA _IC2; r = .58, p = 9.1e‐53* between GM _IC6 and FA _IC2, Figure 4b). The symbol * means FDR correction for multiple comparisons. The loadings of the component triangle also show significant group difference (p = 7.9e‐05*, p = 1.3e‐07*, p = 1.6e‐04*, Figure 4c) for fMRI_IC7, GM _IC6, and FA_IC2, respectively. For the aggregate fMRI component fMRI_IC7, the red regions indicate a higher variability in HCs than SZs and the blue regions are opposite. For sMRI and dMRI, the red brain regions indicate a higher contribution weight in HCs than SZs. In SZ, higher functional activity in thalamus, para‐hippocampus, with lower activity in cerebellum, are accompanied with lower GMV in anterior cingulate cortex (ACC), insula, prefrontal and para‐hippocampus, plus lower FA in forceps major and forceps minor. Detailed anatomical information of the identified multimodal brain regions is summarized in Supplementary Table 3.

FIGURE 4

Linked independent component (IC) triangle that also keep group difference in functional magnetic resonance imaging (fMRI), structural MRI (sMRI), and diffusion MRI (dMRI). (a) The Z‐scored brain maps (|Z| > 2). (b) Correlation between loadings among different modality features. (c) Group difference between SZ and healthy control (HC) on the loadings for the linked IC triangle

Modality specific component

Apart from the linked three‐way fMRI–sMRI–dMRI components, we also identified one fMRI component (fMRI_IC18, Figure 5), containing the left attention network, whose time courses were anticorrelated with attention scores (r = −.2, p = .004). This result highlights the ability of the proposed method in detecting both linked multimodal components as well as modality specific component within the joint model.

FIGURE 5

Modality specific left attention network in fMRI_IC18, whose time courses correlated with attention scores

Correlation with cognition

We found that the loadings of the identified linked GM_IC6 and FA_IC2 are also correlated with major cognitive domains, including speed of processing (r = .40, p = 2.2e‐10* for GM; r = .30, p = 2.4e‐06* for FA), working memory (r = .33, p = 3.2e‐07* for GM), verbal learning (r = .35, p = 5.0e‐08* for GM; r = .30, p = 2.1e‐06* for FA), visual (VIS) learning (r = .36, p = 1.7e‐08* for GM; r = .31, p = 1.0e‐06* for FA) and composite scores (r = .39, p = 1.1e‐09* for GM; r = .34, p = 1.4e‐07* for FA), as displayed in Figure 6.

FIGURE 6

Association with Computerized Multiphasic Interactive Neurocognitive System (CMINDS) cognitive domains of linked component pairs

FNC analysis

One of the advantages of the proposed three‐way pGICA compared with traditional three‐way fusion methods is that, we can also calculate FNC by subject‐specific time courses of the identified fMRI components after the fusion analysis. These fMRI components were then manually classified into seven categories: cerebellar, auditory, VIS, sensorimotor, default‐mode network (DMN), subcortical (SC), and cognitive control (CC) networks based on the brain regions in the components. The mean FNC matrices of SZ and HC are displayed in Figure 7a,b. The group difference FNC matrix between SZ and HC is displayed in Figure 7c. Values in Figure 7c are presented as . Results show that the VIS–SC (fMRI_IC7, Figure 4c) and DMN–CC are group discriminating FNC pairs (, FDR corrected).

FIGURE 7

Functional magnetic resonance imaging (fMRI) functional network connectivity (FNC) matrix based on results of three‐way parallel group independent component analysis (pGICA). (a) SZ mean FNC matrix. (b) Healthy control (HC) mean FNC matrix. (c) Group discriminating FNC matrix. The notation “xxx” denotes elements where (FDR corrected). The components' indexes are displayed in the diagonal cells of the FNC matrix

Linkage replication in COBRE dataset

Here, we further tested the stability of our proposed method in detecting modality linkage, that is, whether the linkage among fMRI_IC7, sMRI_IC6, and dMRI_IC2 extracted in fBIRN can be replicated in the independent COBRE dataset, by performing cross‐site linear projection analysis. Starting with sMRI, in the fBIRN cohort. For the validation cohort, we obtained the mixing matrix by linear projection as: . The same approach was used for dMRI projection from fBIRN to COBRE.

For fMRI, we start with representing the fMRI data for fBIRN cohort. After subject‐level and group‐level PCA, we have where () is the reduced data matrix, () and () are the group‐level and subject‐level whitening matrices, respectively. Here, the same L, which means the same mask, for COBRE and fBIRN. After ICA decomposition, we can calculate , where is the mixing matrix for fMRI and is the group‐level matrix of ICs. Thus, we obtained a group component for COBRE by linear projection as: The subject‐specific components for COBRE were obtained as: where is the ith partition of the dewhitening matrix (, where is number of subjects).

Finally, we can calculate the between‐subject functional variability matrix from Equation (7). Correlation analyses were further performed on the loadings of the projected target component between , , and . Results (Figure 8) show that the triangle linkage detected in fBIRN can be replicated in COBRE cohort (r = .29, p = .001* between fMRI_IC7 and GM_IC6; r = .23, p = .02 between fMRI_IC7 and FA_IC2; r = .44, p = 1.6e‐05* between GM_IC6 and FA_IC2), suggesting that the proposed three‐way pGICA provides a reliable association detection among fMRI–sMRI–dMRI, that is, this association can be replicated in an independent cohort. Furthermore, this replicable fMRI–sMRI–dMRI component triangle also supports group differences between SZ and HC (two‐sample t tests: p = .003* for fMRI, p = .001* for sMRI, and p = .006 for dMRI).

FIGURE 8

Linkage replication by projecting components from Function Biomedical Informatics Research Network (fBIRN) data to Center for Biomedical Research Excellence (COBRE) data. The intermodality linkage detected in fBIRN cohort can be replicated in COBRE cohort, which also supports group differences

In addition to linear projection, we also performed the proposed method on COBRE cohort separately to test the stability in both modality linkage detection and pattern replication. The similarity of component between fBIRN and COBRE was calculated. We calculated the spatial correlation of the identified linked component between two cohorts using only voxels masked at |Z| > 2. First, the spatial maps were transformed into Z scores and masked at |Z| > 2. Then, we obtained two masks from FBIRN (mask_FBIRN) and COBRE (mask_COBRE) respectively, which were used to perform the voxel selection. Only voxels that fell in the union of the masks () were used to calculate the cross‐cohort spatial similarity. Thus, the total number of voxels for calculating the spatial correlation is greatly reduced, for example, from n = 153,594 (whole brain voxels) to m = 2,135 (for the fMRI component). Spatial similarity was finally performed on these commonly identified voxels between two cohorts. In this way, the spatial similarity between fBIRN and COBRE was measured (r = .56, r = .50, and r = .76 for fMRI, GMV, and FA components, respectively), and the modality linkage between the identified component for COBRE was also replicated (Figure 9). It is clear that the proposed three‐way GICA can extract a stable linked component triangle across cohorts. More importantly, PCC and caudate in fMRI, bilateral insula and ACC in GMV, and forceps major and forceps minor in FA were replicated individually in COBRE. This means that not only the modality linkage but also the linked three‐way spatial pattern can both be validated.

FIGURE 9

The proposed three‐way parallel group independent component analysis (GICA) was performed on Center for Biomedical Research Excellence (COBRE) cohort separately to test whether both the pattern and the modality linkage identified in Function Biomedical Informatics Research Network (fBIRN) can be replicated

DISCUSSION AND CONCLUSION

In this study, we proposed a novel approach to combine spatial and spatiotemporal MRI data via a parallel fusion, called three‐way pGICA. We incorporate a regularizer that maximizes the correlation among mixing matrices for fMRI, sMRI, and dMRI, aiming to extract three‐way linked components. Compared with existing three‐way multimodal fusion methods (mCCA, jICA, three‐way pICA), a key benefit of three‐way pGICA is that it can work directly with first‐level 4D fMRI to fuse it with 3D modalities such as sMRI and dMRI. Another benefit is that fMRI FNC matrix can be calculated based on the results of three‐way pGICA. To the best of our knowledge, this is the first proposed method that can incorporate 4D fMRI in a three‐way fusion framework. This method provides a new way to interpret linked patterns among 4D fMRI and two other 3D modality features.

We note that several conditions were explored when designing the simulations for comparing three‐way pGICA with separate GICA/pICA and separate GICA/ICA without fusion, in identifying linked fMRI–sMRI–dMRI component triangle. Simulation results indicate that the proposed three‐way pGICA provides relatively high intermodality linkage estimation accuracy under both strong and weak modality correlation, as well as different noise levels. This means that the proposed method does not inflate the link artificially while achieving comparable accuracy on source map estimation. Due to the lack of correlation regularization between fMRI and the other two modalities, the modality linkage estimation of separate GICA/pICA and separate GICA/ICA decreases significantly when the real association is weak, or the noise level is high, as compared with three‐way pGICA. This demonstrates the advantages of three‐way pGICA in detection cross‐modality association with weak linkage and high noise levels, which are always the case in real human brain data.

In real brain imaging application, we combined data from 4D fMRI, GMV and FA from SZ patients and HCs. One linked component triangle (fMRI_IC7‐GM_IC6‐FA_IC2) was identified. This triangle correlated between each modality pair and also presented significant group difference between SZ and HC. In SZ patients, lower functional activity in subcortical brain areas including thalamus and para‐hippocampus, and higher activity in the cerebellum are identified in the linked fMRI_IC7, consistent with widely reported thalamus‐related network dysfunction and cerebral dysconnections (Pinault, 2011). GM_IC6 and FA_IC2 are correlated with several main cognitive domains, which is consistent with the salience network (including ACC, insula, and prefrontal cortex) detected in sMRI. The covaried forceps major and forceps minor identified in dMRI are involved in multiple high‐order cognitive functions (Kochunov et al., 2017), including attention, working memory (Hegde et al., 2020), and learning (Davidson, 2019). Forceps major has also been reported in ENIGMA large‐scale coordinated study of white matter microstructural differences in schizophrenia (Kelly et al., 2017). Both modality linked components (Figure 4) and modality specific components (Figure 5) can be detected based on the proposed three‐way pGICA. The left attention network identified in fMRI_IC18 is consistent with the association observed between time courses of fMRI_IC18 and attention scores. After multimodal fusion analysis, we also performed post FNC analysis of fMRI modality. Two abnormal FNC pairs, VIS‐SC and DMN‐CC were identified discriminating between SZ and HC groups, which is consistent with the VIS distortions impairment (Butler et al., 2007) and self‐related processing in SZ (Potvin, Gamache, & Lungu, 2019). More importantly, the extracted linked fMRI–sMRI–dMRI triangle can be replicated in the independent COBRE cohort by both linear projection and separately performing three‐way GICA, suggesting that the proposed three‐way pGICA provides a reliable association detection and pattern extraction among the three modalities. The reproducibility of neuroimaging analysis methods has become a point of critical concern in brain imaging studies (Botvinik‐Nezer et al., 2020; Poldrack, Gorgolewski, & Varoquaux, 2019). Here, we further demonstrate that the proposed three‐way pGICA method offers reproducible and reliable linked three‐way component triangle detection among the three MRI modalities, signaling its great potential for extracting linked multimodal imaging biomarkers in brain disorders.

Apart from MRI modalities, single nucleotide polymorphisms, methylation data, and behavioral assessments could also be used as one of the modality fused along with the temporal information of fMRI via the proposed three‐way pGICA. This highlights the flexibility of our method for general multimodal fusion. Moreover, in addition to static FNC analysis, as shown in this work, dynamic FNC can also be evaluated based on three‐way pGICA results. A possible limitation is that the linear back reconstruction for the fMRI modality has be computed in each Infomax iteration to get the group‐ and subject‐specific fMRI components that generate the functional variability matrix C1 linking to the other two modality features, which is time consuming. However, with current advanced computing server clusters, the linked fMRI–sMRI–dMRI triangle can be extract within a few hours.

In sum, we proposed a new temporal information aware three‐way multimodal fusion method called three‐way pGICA, and verified its effectiveness in simulation and real brain imaging data. To the best of our knowledge, this is the first proposed three‐way fusion method that can directly link 4D fMRI with the other two 3D MRI, seeking for potential covaried multimodal biomarkers for brain disorders. Based on the proposed three‐way pGICA, we identified one linked fMRI–sMRI–dMRI triangle that was associated with SZ deficits in major cognitive domains in a discovery dataset and was replicated in an independent dataset, highlighting the promise of the proposed method to detect joint multiway multimodal biomarkers and capture novel information that may be useful to characterize and predict brain disorders. Future work needs to expand the flexibility of the detected association structure to retrieve multiple triangles, as well as potentially missing and higher dimensional linkages.

CONFLICT OF INTEREST

The authors declare no potential conflict of interest.

Appendix S1. Supporting Information

Click here for additional data file.