Unraveling Diagnostic Biomarkers of Schizophrenia Through Structure-Revealing Fusion of Multi-Modal Neuroimaging Data

Fusing complementary information from different modalities can lead to the discovery of more accurate diagnostic biomarkers for psychiatric disorders. However, biomarker discovery through data fusion is challenging since it requires extracting interpretable and reproducible patterns from data sets, consisting of shared/unshared patterns and of different orders. For example, multi-channel electroencephalography (EEG) signals from multiple subjects can be represented as a third-order tensor with modes: subject, time, and channel, while functional magnetic resonance imaging (fMRI) data may be in the form of subject by voxel matrices. Traditional data fusion methods rearrange higher-order tensors, such as EEG, as matrices to use matrix factorization-based approaches. In contrast, fusion methods based on coupled matrix and tensor factorizations (CMTF) exploit the potential multi-way structure of higher-order tensors. The CMTF approach has been shown to capture underlying patterns more accurately without imposing strong constraints on the latent neural patterns, i.e., biomarkers. In this paper, EEG, fMRI, and structural MRI (sMRI) data collected during an auditory oddball task (AOD) from a group of subjects consisting of patients with schizophrenia and healthy controls, are arranged as matrices and higher-order tensors coupled along the subject mode, and jointly analyzed using structure-revealing CMTF methods [also known as advanced CMTF (ACMTF)] focusing on unique identification of underlying patterns in the presence of shared/unshared patterns. We demonstrate that joint analysis of the EEG tensor and fMRI matrix using ACMTF reveals significant and biologically meaningful components in terms of differentiating between patients with schizophrenia and healthy controls while also providing spatial patterns with high resolution and improving the clustering performance compared to the analysis of only the EEG tensor. We also show that these patterns are reproducible, and study reproducibility for different model parameters. In comparison to the joint independent component analysis (jICA) data fusion approach, ACMTF provides easier interpretation of EEG data by revealing a single summary map of the topography for each component. Furthermore, fusion of sMRI data with EEG and fMRI through an ACMTF model provides structural patterns; however, we also show that when fusing data sets from multiple modalities, hence of very different nature, preprocessing plays a crucial role.

1. Introduction

Multiple neuroimaging techniques provide complementary views of neural structure and function. For instance, one of the most commonly used neuroimaging methods, electroencephalography (EEG), measures the electrical activity with high temporal but low spatial resolution, while functional magnetic resonance imaging (fMRI) records the changes in the blood flow with high spatial but low temporal resolution (Bunge and Kahn, 2009; Uludag and Roebroeck, 2014). Therefore, joint analysis of signals from multiple neuroimaging modalities is of interest in order to better understand neural activity and to discover reliable diagnostic biomarkers for psychiatric disorders, such as schizophrenia (Daunizeau et al., 2009; Sui et al., 2012, 2018; Dahne et al., 2015; Liu et al., 2015).

With the advances in technology, vast amounts of neuroimaging data has been generated; however, data mining or signal processing methods so far have limited success in terms of finding reliable diagnostic imaging biomarkers for many psychiatric disorders (Sui et al., 2012; Takahashi and Suzuki, 2018). One of the reasons for this limited success has been the fact that data fusion is a particularly challenging task when the goal is to extract reproducible and interpretable patterns. Data from different sources consists of both shared (or common) and unshared (or distinct) underlying patterns (Alter et al., 2003; Daunizeau et al., 2009; Lock et al., 2013; Uludag and Roebroeck, 2014), and even the definition of “sharedness” is a topic of current research (Farias et al., 2016; Smilde et al., 2017). Furthermore, data sets from different modalities may be of different orders, such as multi-channel EEG signals from multiple subjects can be represented in the form of a third-order tensor with modes: subject, time, and channel, while fMRI data is often represented as a subject by voxel matrix (Figure 1). Similar challenges have been observed in other disciplines targeting biomarker discovery as well, e.g., in joint analysis of omics data (Acar et al., 2015), where the ultimate goal is to discover significant metabolites, genes, etc. as potential biomarkers.

FIGURE 1

Figure 1. A third-order tensor representing multi-channel EEG signals is coupled with fMRI and sMRI data in the form of matrices in the subject mode.

The common approaches for fusion of multi-modal neuroimaging data are based on matrix factorizations, such as joint independent component analysis (jICA) (Calhoun et al., 2006b), parallel ICA (Calhoun et al., 2009) and independent vector analysis (IVA)-based techniques (Adali et al., 2015a), where signals from multiple modalities are represented as matrices, e.g., fMRI data in the form of a subject by voxel matrix, and EEG signals as a subject by time matrix (Adali et al., 2015a). Matrix factorization-based fusion methods require additional constraints to recover patterns uniquely (Alter et al., 2003; Klami et al., 2013; Lock et al., 2013; Adali et al., 2015a) and a common practice in neuroscience is to assume that extracted patterns (i.e., biomarkers, or spatial/temporal patterns) are statistically independent. Drawbacks of the traditional methods are 2-fold: (i) in the presence of multi-channel EEG signals, which can naturally be represented as a third-order tensor, data is either matricized in the form of a subject by time-channel matrix (Swinnen et al., 2014) or only the signal from a single channel is analyzed (Adali et al., 2015a), ignoring the potential multilinear structure of multi-channel EEG signals, (ii) statistical independence might be a too strong constraint to impose on the patterns; therefore, methods may fail to capture the true patterns (Acar et al., 2015).

In contrast, coupled matrix and tensor factorizations (CMTF), introduced more recently, have proven useful in terms of addressing the drawbacks of matrix factorization-based fusion methods by jointly analyzing data sets in the form of matrices and higher-order tensors without imposing constraints on the components when the higher-order tensors have a defined multilinear structure (Acar et al., 2015). CMTF-based approaches factorize higher-order tensors using a tensor factorization model while simultaneously factorizing the data sets in the form of matrices, and enable the exploration of the potential multilinear structure inherent to, for instance, multi-channel EEG signals. Previously, analyzing multi-channel EEG signals using tensor factorizations has shown promising performance in terms of capturing spatial, spectral and temporal signatures of epileptic seizures (Acar et al., 2007; De Vos et al., 2007) as well as providing better understanding of brain activity patterns (Möcks, 1988; Miwakeichi et al., 2004; Mørup et al., 2007), see also (Cong et al., 2015) for a review. Another review (Hunyadi et al., 2017) discusses the studies making use of tensor factorizations to analyze EEG and fMRI signals, with a particular focus on epilepsy. Therefore, recent studies analyzing neuroimaging signals from multiple modalities have arranged neuroimaging data sets as higher-order tensors, and used CMTF-type methods to jointly analyze, e.g., EEG and magnetoencephalography (Becker et al., 2012; Naskovska et al., 2017) (both arranged as higher-order tensors), EEG and electro-ocular artifacts (Rivet et al., 2015) (both arranged as higher-order tensors), and EEG and fMRI (Karahan et al., 2015; Hunyadi et al., 2016; Eyndhoven et al., 2017) in the form of a matrix coupled with a third-order tensor, or as coupled higher-order tensors as in Chatzichristos et al. (2018), or arranged multiple diffusion tensor imaging modalities as a third-order tensor and coupled that with gray matter maps (Groves et al., 2011). However, when jointly analyzing a matrix coupled with a higher-order tensor, CMTF-based methods assuming that coupled data sets have only shared components, may fail to capture the underlying patterns in the presence of both shared and unshared components (Acar et al., 2014; Lathauwer and Kofidis, 2017); therefore, they are not ideal for biomarker discovery.

In this paper, we use a CMTF-based approach to jointly analyze neuroimaging signals from multiple modalities, more specifically, fMRI, sMRI and EEG data, collected during an auditory oddball (AOD) task from a group of subjects consisting of patients with schizophrenia and healthy controls with the goal of unraveling potential diagnostic biomarkers for schizophrenia. To the best of our knowledge, this is the first comprehensive study of a CMTF-based method for biomarker discovery for a psychiatric disorder discussing both strengths and limitations of the proposed framework, building onto our preliminary results in Acar et al. (2017a,b). Furthermore, due to the reproducibility and uniqueness requirements of such an application, we use a structure-revealing CMTF model, known as the advanced CMTF (ACMTF) model (Acar et al., 2014), to estimate weights of the components in each modality in order to identify shared/unshared components and quantify the contribution from each modality. Our preliminary studies have shown the promise of the ACMTF model in terms of capturing neural patterns that can differentiate between patients with schizophrenia and healthy controls by jointly analyzing EEG-fMRI signals (Acar et al., 2017b) and EEG-fMRI-sMRI data (Acar et al., 2017a); however, those two studies used only a subset of electrodes, making it difficult to evaluate the added value of each modality in terms of biomarker discovery. Also, in this paper, we include an additional metric to study the additive value of each modality, and evaluate the performance of the models in terms of clustering subjects from different groups, whereas the previous studies only used the interpretation and statistical significance of extracted patterns in terms of differentiating between groups. Clustering results complement univariate statistical significance tests and show whether combinations of potential biomarkers provide meaningful clusters. We show that EEG analysis using a CP [also known as Canonical Decomposition (CANDECOMP) and Parallel Factor Analysis (PARAFAC)] tensor model and joint analysis of EEG, fMRI as well as EEG, fMRI and sMRI reveal statistically significant and biologically meaningful components in terms of differentiating between patients with schizophrenia and healthy controls. In comparison to the results when only the EEG data is analyzed, the incorporation of fMRI signals results in clearer spatial maps and better clustering performance. With the incorporation of sMRI, we obtain structural patterns in addition to temporal and spatial patterns of functional activity without degrading the clustering performance. ACMTF models with different parameter settings have been compared, and based on detailed experiments, we observe that ACMTF consistently reveals similar significant patterns, which provide a concise summary of the topography, while being sensitive to certain parameters for uniqueness.

2. Materials and Methods

2.1. Background

In this section, we briefly discuss the CP tensor model as well as structure-revealing CMTF and jICA models. Let the third-order tensor $\text{X}$ ∈ ℝ^I×J×K with modes subject, time, and electrode, and matrices Y ∈ ℝ^I×M (subject by voxel) and Z ∈ ℝ^I×L (subject by voxel), represent multi-channel EEG, fMRI, and sMRI data, respectively (as in Figure 1).

2.1.1. CANDECOMP/PARAFAC (CP)

The CP model (Carroll and Chang, 1970; Harshman, 1970), also referred to as the canonical polyadic decomposition (Hitchcock, 1927), is one of the most popular tensor factorization models. It is considered as one of the extensions of the matrix singular value decomposition (SVD) to higher-order tensors (N ≥ 3) and represents the tensor as a sum of rank-one tensors, i.e., rank-one components. For a third-order tensor $\text{X}$ ∈ ℝ^I×J×K, the R-component CP model is defined as

$\begin{array}{l}\text{X}\approx \widehat{\text{X}}=\left[\left[\lambda ;\text{A},\text{B},\text{C}\right]\right]={\displaystyle \sum _{r=1}^R}{\lambda }_r{\text{a}}_r◦{\text{b}}_r◦{\text{c}}_r,\end{array}$

where ◦ indicates the vector outer product. The vectors from the rank-one components are collected in the factor matrices $\text{A}\in {ℝ}^{I\times R}=\left[{\text{a}}_1\dots {\text{a}}_R\right],\text{B}\in {ℝ}^{J\times R}=\left[{\text{b}}_1\dots {\text{b}}_R\right]$ and $\text{C}\in {ℝ}^{K\times R}=\left[{\text{c}}_1\dots {\text{c}}_R\right]$. In this definition, columns of all factor matrices are assumed to be normalized to unit 2-norm and the norms are absorbed in the vector λ ∈ ℝ^R×1. For the third-order tensor $\text{X}$ consisting of the EEG data, the factor matrices A, B and C correspond to the extracted factor vectors in the subject, time and electrode modes, respectively. By modeling $\text{X}$ using a CP model, we assume that component r models a brain activity with temporal and spatial patterns represented by b_r and c_r. Multi-channel EEG signals from each subject are a linear mixture of these R brain activities mixed using subject-specific weights. The CP model is also known as a topographic components model (Möcks, 1988). Note that the terms factor and component are used interchangeably throughout the paper, and refer to the rank-one matrices and/or higher-order rank-one tensors.

In contrast to matrix factorizations, the CP model of higher order tensors is unique up to scaling and permutation under mild conditions (Kruskal, 1977; Sidiropoulos and Bro, 2000), without the need to impose additional constraints.

2.1.2. Structure-Revealing Coupled Matrix and Tensor Factorizations

Given the third-order tensor $\text{X}$ coupled with matrices Y and Z in the subject mode we can jointly factorize them using a structure-revealing CMTF model (a.k.a. ACMTF) (Acar et al., 2014) that fits a CP model to tensor $\text{X}$ and factorizes matrices Y and Z in such a way that the factor matrix extracted from the common mode, i.e., subject, is the same in the factorizations of all data sets. An R-component ACMTF model minimizes the following objective function:

$\begin{array}{ll}\begin{array}{l}f(\lambda ,\Sigma ,\Gamma ,\text{A},\text{B},\text{C},\text{D},\text{E})=\parallel \text{X}-\left[\left[\lambda ;\text{A},\text{B},\text{C}\right]\right]{\parallel }^2+\parallel \text{Y}-\text{A}\Sigma {\text{D}}^{\text{T}}{\parallel }^2\\ +\parallel \text{Z}-\text{A}\Gamma {\text{E}}^{\text{T}}{\parallel }^2+\beta \parallel \lambda {\parallel }_1+\beta \parallel \sigma {\parallel }_1+\beta \parallel \gamma {\parallel }_1,\end{array}& (1)\end{array}$

where the columns of factor matrices have unit norm, i.e., ∥a_r∥ = ∥b_r∥ = ∥c_r∥ = ∥d_r∥ = ∥e_r∥ = 1 for r = 1, …, R. λ, σ, γ ∈ ℝ^R×1 are the weights of rank-one terms in $\text{X}$, Y, and Z, respectively. Σ, Γ ∈ ℝ^R×R are diagonal matrices with entries of σ and γ on the diagonal. D ∈ ℝ^M×R and E ∈ ℝ^L×R correspond to factor matrices in the voxel mode in fMRI and sMRI. ∥.∥ denotes the Frobenius norm for matrices/higher-order tensors, and the 2-norm for vectors. ∥.∥₁ denotes the 1-norm of a vector, i.e., $\parallel X{\parallel }_1=\sum _{r=1}^R\left|x_r\right|$ and β > 0 is a penalty parameter. Imposing penalties on the weights in Equation (1) sparsifies the weights so that unshared components have weights close to 0 in some data sets. The model is illustrated in Figure 2. By jointly analyzing neuroimaging data using an ACMTF model, we assume that each component extracted from $\text{X}$ models a brain activity with certain temporal (b_r) and spatial (c_r) signatures, and the corresponding component in Y models related brain activity with higher spatial resolution using d_r while the component in Z provides information about the tissue type at a very high spatial resolution using e_r. Since the same factor matrix A is extracted from the subject mode from all data sets, subject-specific coefficients in all modalities are assumed to be the same. The ACMTF model inherits uniqueness from CP (Sorensen and De Lathauwer, 2015), as long as all factors are shared, and provides reproducible and interpretable factors. Note that, in the presence of both shared/unshared components, 1-norm penalties on the weights help to obtain unique solutions, which has been so far only shown experimentally (Acar et al., 2014).

FIGURE 2

Figure 2. Modeling of tensor $\text{X}$ coupled with matrices Y and Z in the subject mode using a structure-revealing CMTF model.

2.1.3. Joint ICA

An alternative approach to jointly analyze $\text{X}$, Y and Z is to use a matrix factorization-based fusion approach by matricizing $\text{X}$ in the subject mode as a subject by time–electrode matrix denoted as X₍₁₎. Joint ICA (Calhoun et al., 2006b) concatenates matrices representing the data from different modalities and models the constructed matrix using an ICA model as follows:

$\begin{array}{ll}\left[{\text{X}}_{(1)}\text{Y}\text{Z}\right]=\text{A}\text{S}& (2)\end{array}$

where, for an R-component ICA model, A ∈ ℝ^I×R corresponds to the mixing matrix, similar to the factor matrix in Equation (1), and S ∈ ℝ^R×(JK+M+L) represents the source signals. Note that subject covariations across all data sets, i.e., modalities, are assumed to be the same in jICA as in ACMTF, since the same mixing matrix is shared across the data sets. However, in this case the model does not include an adaptive estimation of contributions from each modality as in ACMTF, and though this can be captured to a degree within the weights from the estimated components from each modality, it represents a more constrained model. The rows of S corresponding to patterns of brain activity are assumed to be statistically independent.

2.2. Experiments

We make use of EEG, fMRI, and sMRI data collected from patients with schizophrenia and healthy controls to show the use of ACMTF models to discover potential diagnostic biomarkers for schizophrenia. Our experiments focus on joint analysis of EEG and fMRI data, and discuss the effect of different modeling choices, i.e., number of components (R), the penalty parameter (β), preprocessing, and use of subsets of electrodes. We also discuss the performance of ACMTF in comparison with jICA. Furthermore, the analysis of only EEG signals and joint analysis of EEG, fMRI, and sMRI have been studied to show the information gain with each modality and potential issues due to the use of additional modalities.

2.2.1. Data

The participants in this study were recruited through advertisements, presentations at local universities, and word-of-mouth. The 32 participants provided written and informed IRB-approved consent at Hartford Hospital. The participants were compensated for their participation. Patients met the criteria for schizophrenia as defined in the DSM-IV on the basis of a SCID diagnosis as well as a review of their case file (First et al., 2002b). The healthy participants were screened to ensure that they were free from DSM-IV Axis I or Axis II psychopathology, as assessed using the SCID (First et al., 2002a) prior to their inclusion in the study. Exclusion criteria for the healthy participants also included a history of psychosis in any first-degree relatives. This was assessed through an interview. The participants in both groups had self-reported normal hearing and were able to perform the AOD task successfully in the trial sessions prior to scanning.

The EEG, fMRI and sMRI data were separately collected from 21 healthy controls and 11 patients with schizophrenia during an auditory oddball task (AOD), where subjects pressed a button when they detected an infrequent target sound within a series of auditory stimuli. The details of the AOD task used in this work have been published previously (Kiehl et al., 2005). The task consisted of two runs of 244 pseudo-randomly distributed stimuli each, following an initial practice block of 10 trails to ensure that the participants understood the instructions. The stimuli were three different types of tones: a standard stimulus (1,000 Hz tone, probability = 0.8), a target stimulus (1,500 Hz tone, probability = 0.1), and novel stimuli (non-repeating random digital noises, such as tone sweeps or whistles, probability = 0.1). Each auditory stimulus was presented for 200 ms with an inter-stimulus interval of 2,000 ms. The intervals between the non-standard stimuli consisted of 3–5 standard stimuli (approximately 8–12 s) (Kiehl et al., 2005). Two scans were performed, each lasting approximately 9 min. Thus, the total duration of the scan lasted approximately 18 min.

2.2.2. Data Preprocessing

For the fMRI data, task-related spatial activity maps for each subject were computed by the general linear model-based regression approach using the statistical parametric mapping toolbox (SPM2)¹. The fMRI preprocessing pipeline used for this data has been described previously (Calhoun et al., 2006a). More precisely, for each subject, each image in the scan was realigned using the INRIalign algorithm (Freire and Mangin, 2001). The aligned data was then normalized into Montreal Neurological Institute (MNI) space and spatially smoothed with a 10 × 10 × 12 mm³ full width at half maximum (FWHM) Gaussian kernel. The normalized data was then subsampled to 3 × 3 × 3 mm, thus resulting in brain volumes of dimensions 53 × 63 × 46 voxels. The data for each subject was then analyzed through a multivariate regression where the regressors are formed from the auditory stimuli, their temporal derivatives, and an intercept. The regressors for the auditory stimuli are formed by modeling the stimuli as delta functions convolved with the default SPM2 hemodynamic response function. We use the regression coefficients (beta-values) corresponding to the target tone as the feature for these analyses (Calhoun et al., 2006a). By making use of these features and removing voxels not corresponding to brain regions, we constructed a matrix of 32 subjects by 60,186 voxels representing the fMRI signals.

The details of the preprocessing pipeline used for the EEG data have been previously published in Calhoun et al. (2006b). Briefly, the EEG channels were amplified (20,000 gain) with a passband between 0.01 and 100 Hz. The signals were digitized at a rate of 500 samples per second. The EEG data was then preprocessed using ICA to remove any ocular artifacts (Jung et al., 2000). The data were then filtered with a 20 Hz low-pass filter in order to reduce electromyographic activity. EEG epochs from -115 to 788 ms (451 time points) around stimulus onsets were used for the event-related potentials (ERP). ERP were formed only from those trials in which the subjects correctly identified the target stimulus (Calhoun et al., 2006b). Out of 64 electrodes in total, we used 62 electrodes by excluding the two corresponding to vertical and horizontal electrooculography (EOG) electrodes. Multi-channel EEG signals were then arranged as a third-order tensor: 32 subjects by 451 time samples by 62 electrodes. In order to assess the modeling assumptions, we also used a subset of electrodes, i.e., AF3, AF4, Fz, T7, C3, Cz, C4, T8, Pz, PO3, PO, and in that case, formed a third-order tensor with 11 electrodes as in (Acar et al., 2017a,b). This subset of electrodes was determined with the goal of selecting a small set of electrodes that would record the relevant functional activity related to the AOD task. The selected subset includes electrodes from the frontal, motor and parietal areas that are expected to be involved in the motor and auditory responses as well as the planning stage for the given task.

For the sMRI data, we computed probabilistically segmented gray matter images for each subject and by using these features formed a matrix of 32 subjects by 306,640 voxels. The details of the sMRI preprocessing pipeline have also been presented previously in Calhoun et al. (2006a). The image files were first normalized using a 12 parameter affine model to the 152 average T₁ MNI template. The images were next segmented into gray matter, white matter, and cerebrospinal fluid. Then the segmentation was smoothed with an 8-mm FWHM Gaussian kernel and the segments were then averaged to create the gray matter, white matter, and cerebrospinal fluid templates. Each subject's segmented T₁ images were then normalized to the customized gray matter templates. Then, for each subject, the warped T₁ images were segmented into gray matter, white matter, and cerebrospinal fluid maps using a model clustering algorithm resulting in probabilistic segmentations of gray matter (Ashburner and Friston, 2000). Finally, the resulting gray matter images were smoothed with a 12-mm FWHM Gaussian kernel (Good et al., 2001).

2.2.3. Experimental Setting

Before the analysis, we centered the third-order EEG tensor across the time mode, and scaled within the subject mode by dividing each horizontal slice by its standard deviation (see Bro and Smilde, 2003 for further details on preprocessing of higher-order tensors). The fMRI and sMRI data were also preprocessed by centering each row (subject-wise) and dividing each row by its standard deviation. When fitting the ACMTF model, each data set was also divided by its Frobenius norm to give equal importance to the approximation of each data set.

In order to demonstrate the information gained by the addition of each modality and sensitivity of the fusion approach to various modeling choices, the following experiments are carried out:

• Individual analysis of the EEG tensor using a CP model,

• Joint analysis of the EEG tensor coupled with fMRI using an ACMTF model (i) by leaving out signals from one subject at a time, (ii) for 11-electrode vs. 62-electrode case, (iii) in comparison with jICA, (iv) with different number of components, R, (v) with different sparsity penalty parameters, β, (vi) with/without additional centering across the subject mode.

• Joint analysis of the EEG tensor coupled with fMRI and sMRI using an ACMTF model.

The CP model is fit using CP-OPT (Acar et al., 2011) from the Tensor Toolbox version 2.5² using the non-linear conjugate gradient algorithm (NCG). For the ACMTF model, we use ACMTF-OPT (Acar et al., 2014) from the CMTF Toolbox version 1.1³, also using NCG to fit the model. Multiple random initializations are used to fit the models, and the solution corresponding to the minimum function value is reported. Furthermore, the ACMTF model is experimentally validated to be unique by obtaining the same minimum function value⁴ a number of times and checking the uniqueness of model parameters, i.e., factor matrices and weights of the components (up to permutation)⁵. For jICA, we unfold the EEG tensor in the subject mode forming a matrix of 32 subjects × 27,962 (time-electrodes), and concatenate the resulting matrix with the fMRI matrix. The concatenated matrix is decomposed using an ICA algorithm based on entropy bound minimization (ICA-EBM)⁶ (Li and Adali, 2010), which makes use of a flexible density model that is a better fit to data formed by concatenating signals from different modalities (Adali et al., 2015a). We fit the model using multiple random initializations and report the most stable run determined by a minimum spanning tree-based approach (Du et al., 2016).

2.3. Performance Evaluation

The performance is assessed both qualitatively and quantitatively. The qualitative assessment relies on the interpretation of the extracted temporal and spatial patterns as well as comparisons with the previous findings in the literature on schizophrenia. For the quantitative assessment, we perform the following:

• Two-sample t-test: Since the ultimate goal of any factorization of this data is the discovery of latent factors that can differentiate between patients with schizophrenia and healthy controls, we can quantify the performance of a method based on its ability to produce factors that can provide such a differentiation. With the assumption of unequal variance for the healthy control and patient groups, a two-sample t-test is applied on each column of the factor matrix extracted from the subject mode, which is of size 32 by R. Out of R columns, those that have subject mode vectors that are statistically significant, i.e., with p-values < 0.05, are identified and corresponding temporal and spatial patterns are reported as potential biomarkers.

• Clustering: Subjects are clustered into two groups based on the factor matrix corresponding to the subject mode using k-means clustering, where k-means is performed 100 times with different initializations and the most consistent cluster assignments are used. Unlike the t-test based approach that is performed on each column individually, clustering is performed on all possible combinations of the columns of the factor matrix and the performance of the best combination is reported. Therefore, this approach provides a more global view of the discriminatory power of the resulting factorization than the t-test based approach. The clustering performance is assessed in terms of accuracy and F₁-score, where F₁-score $=\frac{2\mathrm{\text{Precision}}\times \mathrm{\text{Recall}}}{\mathrm{\text{Precision}}+\mathrm{\text{Recall}}}$. Precision and Recall are defined as $\mathrm{\text{Precision}}=\frac{\mathrm{\text{\# of true positives}}}{\mathrm{\text{\# of true positives}}+\mathrm{\text{\# of false positives}}}$ and $\mathrm{\text{Recall}}=\frac{\mathrm{\text{\# of true positives}}}{\mathrm{\text{\# of true positives}}+\mathrm{\text{\# of false negatives}}}$, and a patient being clustered as a patient is considered a true-positive.

When assessing different modeling choices, we also report the model fit defined as:

$\begin{array}{l}\text{Fit}=100\times \left(1-\frac{{\parallel \text{X}-\hat{\text{X}}\parallel }^2}{\parallel \text{X}{\parallel }^2}\right),\end{array}$

where X stands for the raw data (e.g., EEG tensor or fMRI/sMRI matrix), and $\widehat{\mathbf{\text{X}}}$ denotes the model. A fit of 100% means that the data is fully explained by the model. The fit shows whether the model explains the data well and indicates the unexplained part left in the residuals. Also, the change in model fit for different number of components shows whether there is a significant gain, in terms of explaining the remaining part in the residuals, by adding more components.

Finally, we compare the similarity of the significant components extracted by different models using a similarity score called the factor match score (FMS). The FMS of component k from two models ${\widehat{\text{X}}}^{(1)}$ and ${\widehat{\text{X}}}^{(2)}$ of the EEG tensor $\text{X}$ is defined as

$\begin{array}{l}\text{FMS}{}_k({\widehat{\text{X}}}^{(1)},{\widehat{\text{X}}}^{(2)})=\frac{{|{\text{a}}_k^{(1)}}^{\text{T}}{\text{a}}_k^{(2)}|}{\parallel {\text{a}}_k^{(1)}\parallel \parallel {\text{a}}_k^{(2)}\parallel }\times \frac{\left|{{\text{b}}_k^{(1)}}^{\text{T}}{\text{b}}_k^{(2)}\right|}{\parallel {\text{b}}_k^{(1)}\parallel \parallel {\text{b}}_k^{(2)}\parallel }\times \frac{\left|{{\text{c}}_k^{(1)}}^{\text{T}}{\text{c}}_k^{(2)}\right|}{\parallel {\text{c}}_k^{(1)}\parallel \parallel {\text{c}}_k^{(2)}\parallel },\end{array}$

and from two models ${\widehat{\mathbf{\text{Y}}}}^{(1)},{\widehat{\mathbf{\text{Y}}}}^{(2)}$ of the fMRI matrix Y as

$\begin{array}{l}\text{FMS}{}_k({\hat{\text{Y}}}^{(1)},{\hat{\text{Y}}}^{(2)})=\frac{\left|{{\text{a}}_k^{(1)}}^{\text{T}}{\text{a}}_k^{(2)}\right|}{\parallel {\text{a}}_k^{(1)}\parallel \parallel {\text{a}}_k^{(2)}\parallel }\times \frac{\left|{{\text{d}}_k^{(1)}}^{\text{T}}{\text{d}}_k^{(2)}\right|}{\parallel {\text{d}}_k^{(1)}\parallel \parallel {\text{d}}_k^{(2)}\parallel },\end{array}$

where ${\text{a}}_k^{(i)},{\text{b}}_k^{(i)},{\text{c}}_k^{(i)},{\text{d}}_k^{(i)}$ correspond to the kth column of the factor matrix corresponding to subject, time, electrode, and voxel mode of the ith model, respectively, after finding the best matching factors for the two models. When components are compared for the models with mismatching dimensions, such as number of subjects or number of electrodes, the mismatching dimension is omitted in the product. An FMS close to one implies similarity of the compared components, while very different components will have an FMS close to zero. FMS is used to quantify the reproducibility of the extracted patterns in addition to qualitative interpretations based on the plots.

For visualization of the extracted components, patterns from fMRI and sMRI voxel modes are plotted as z-maps, thresholded at |z| ≥ 2.7, where red indicates an increase in controls over patients and blue indicates an increase in patients over controls. This threshold value was selected after considering several alternatives in order to aid in the visual interpretation of the components. Patterns extracted from the electrode mode of the EEG tensor are plotted using the topoplot function from the EEGLAB v13.6.5b (Delorme and Makeig, 2004).

3. Results

3.1. Individual Analysis of the EEG Tensor Using a CP Model

As shown in Figure 3, the CP model of the EEG data in the form of a subject by time by electrode tensor constructed using 62 electrodes has captured significant components in terms of differentiating healthy controls and patients with schizophrenia. The model is fit using R = 3 components, and reveals factor matrices A ∈ ℝ^{32 × 3}, B ∈ ℝ^{451 × 3}, and C ∈ ℝ^{62 × 3} corresponding to subject, time, and electrode modes, respectively. The fit of the model is 62% indicating that using only three components a major part of the data can be explained. The t-tests performed on the columns of A indicate that all three components are significant. The first component, whose factor vectors in time and electrode modes are represented in Figure 3A, corresponds to the third positive peak (P3) and is heavily weighted by central electrodes. The second component with the factor vectors in time and electrode modes represented in Figure 3B, refers to the N1-P2 as well as the N2-P3 transitions and is heavily weighted by central and parietal electrodes. The third component with the factor vectors in time and electrode modes represented in Figure 3C, refers to the N2 as well as a negative peak after P3 and is heavily weighted by frontal and central electrodes. CP models with different number of components have been fitted to the data as well but those either revealed fewer components with less significance or are degenerate, i.e., a CP model with that many components is not an appropriate model for the data (see Kolda and Bader, 2009 for more on degeneracy). Table 1 shows that subjects can be clustered into two groups with 81% accuracy using the factor matrix A. Note that clustering performance is similar to the one achieved using the CP model of a tensor constructed using only a subset of electrodes indicating that the assumption of same subject coefficients and temporal patterns across all electrodes is not decreasing the performance. This also may indicate that the additional electrodes are not providing much added information beyond that which is contained by a subset of the electrodes.

FIGURE 3

Figure 3. Temporal and spatial patterns from the statistically significant components of the CP model. (A) Component 1 corresponds to the P3 peak mainly represented by central electrodes, (B) Component 2 refers to the N1-P2 as well as N2-P3 transitions, with high contributions from central and parietal electrodes, (C) Component 3 refers to the N2 as well as a negative peak after P3, heavily weighted by frontal and central electrodes. The corresponding p-values are 2.1 × 10⁻³, 1.6 × 10⁻², 1.4 × 10⁻⁴, respectively. Columns of the factor matrix in the time mode are in red while green plots show signals from individual electrodes averaged across all subjects.

TABLE 1

Table 1. Performance in terms of clustering for different modeling values as well as the factor match scores in comparison to the 10-component ACMTF model (no centering, 62 electrodes).

3.2. Joint Analysis of EEG and fMRI

Shown in Figure 4, the joint analysis of the EEG tensor and fMRI matrix using an ACMTF model has revealed significant components in terms of differentiating between healthy controls and patients while also providing spatial patterns in much higher resolution and improving the clustering performance compared with the CP model of the EEG tensor. The 10-component ACMTF model extracts factor matrices A ∈ ℝ^{32 × 10}, B ∈ ℝ^{451 × 10}, C ∈ ℝ^{62 × 10}, and D ∈ ℝ^{60186 × 10} corresponding to subject, time, electrode, and voxel modes, respectively, as well as weights of the components in EEG (λ ∈ ℝ^{10 × 1}) and fMRI (σ ∈ ℝ^{10 × 1}). The sparsity penalty parameter is set to β = 10⁻³. The fit is 79% and 65% for EEG and fMRI, respectively, indicating that the extracted factors, which have high weights in both data sets indicating shared factors, account for a large part of both data sets. The t-test on the columns of A reveals that out of ten components, only two, components 1 and 9, are statistically significant. Figures 4A,B illustrate the factor vectors of the two significant factors in time (b_r), electrode (c_r) and voxel (d_r) modes. Figure 4C shows component weights in each data set. From Figure 4, we see that though both significant components have a contribution from both EEG and fMRI, the contribution from EEG to each of these components is greater. This means that the discriminatory information plays a larger part in the EEG data than it does in fMRI (see Acar et al., 2014 for further information on how to interpret the weights of the components). Note that the statistical significance of the components and weights of components are different concepts. Statistical significance is determined based on a two-sample t-test on the columns of the factor matrix in the subject mode and indicates whether a component can be considered to be significant in terms of separating patients and controls. On the other hand, the weight of a component indicates how much that component contributes to the model fit. Therefore, a component with a high weight may not be statistically significant in terms of distinguishing between patients and controls.

FIGURE 4

Figure 4. Temporal and spatial patterns from the statistically significant components of the ACMTF model of the EEG tensor with 62 electrodes and fMRI data, with R = 10. (A) Component 1 refers to the P2-N2 transition as well as the P3 peak, heavily weighted by the frontal and central electrodes in the EEG, and the fMRI shows increased activation of controls over patients in the sensorimotor cortex, (B) component 9 is related to the P3 peak, heavily weighted by the frontal electrodes in the EEG and the fMRI indicates a decrease in activation of controls vs. patients in some regions of the sensorimotor cortex and parietal lobe, (C) weights of the components in each data set. The corresponding p-values are 6.2 × 10⁻³, 1.9 × 10⁻², respectively. Columns of the factor matrix in the time mode are in red while green plots show signals from individual electrodes averaged across all subjects.

The first component, whose factor vectors in the time, electrode, and voxel modes are shown in Figure 4A, is similar to the component shown in Figure 3C, and refers to the P2-N2 transition as well as the P3 peak and is heavily weighted by the frontal and central electrodes. A similar activation pattern is seen in the positive activations in the fMRI, though with greater spatial resolution. The second component, with the factor vectors in time, electrode, and voxel modes shown in Figure 4B, shares some similarity with the component shown in Figure 3A, since both are related to the P3 peak but this component is more heavily weighted by the frontal electrodes. The fMRI part in Figure 4B indicates a decrease in activation of controls vs. patients in parts of the anterior sensorimotor cortex. We should note that there are some similarities between the areas highlighted in the topographic maps and the regions highlighted in the fMRI. The areas of increased activation of controls over patients in the fMRI, namely frontal and sensorimotor, generally correspond to the greatest weights in the topographic maps. Similar components have been found previously in other analyses of similar data (Acar et al., 2017b), thus increasing our confidence in the results. In comparison to the individual analysis of the EEG tensor, the clustering performance is also much higher, i.e., 88% accuracy and F₁ score of 0.82, as shown in Table 1. This indicates that the ACMTF reveals more discriminatory factors through the inclusion of complementary information from the fMRI.

3.2.1. Sensitivity Analyses

3.2.1.1. Leave-one-out

The patterns captured in different modes using an ACMTF model are reproducible in case of changes in data sets. In order to evaluate the consistency of the results to changes in the original data, we leave out one subject at a time and fit the ACMTF model using the same parameters (i.e., R = 10, β = 10⁻³). In other words, we construct 32 different EEG-fMRI data set pairs (with 31 subjects) and compare the significant factors extracted using the ACMTF model of each pair with the model derived using 32 subjects. Table 2 shows that average FMS for component 1 (Figure 4A), which is the most significant factor, is 0.98 for EEG and 0.95 for fMRI indicating close to exact recovery of the same patterns. Average FMS for the less significant component, i.e., component 9 (Figure 4B), is around 0.90 indicating similar patterns. Furthermore, the average clustering performance is the same as the performance of the original model estimated using data from 32 subjects.

TABLE 2

Table 2. Leave-one-out sensitivity analysis: average values (standard deviation) of FMS, clustering performance and fit of the models built on data sets with 31 subjects.

3.2.1.2. 11 electrodes vs. 62 electrodes

When jointly analyzing EEG and fMRI, ACMTF achieves a slightly better performance using a subset of electrodes from certain regions of interest during the construction of the EEG tensor than the case where all 62 electrodes are used to construct the tensor. In our preliminary studies (Acar et al., 2017b), we observed similar components when comparing the 11-electrode case with 62-electrode case while achieving higher statistical significance and better interpretability using 11 electrodes. These observations are also supported by our findings in this study on a slightly different set of subjects (38 subjects in Acar et al., 2017b vs. 32 subjects in this paper). Figure 5 illustrates the factor vectors, in the time and voxel modes, of the two most significant components, captured by a 10-component ACMTF model of the EEG tensor with 11 electrodes and fMRI data⁷. The fit is 80% and 65% for EEG and fMRI, respectively. As previously observed, both components have higher significance thus indicating that the additional electrodes are not contributing much additional discriminatory information compared with the original 11 electrodes. This may also be related to the fact that the most contributing electrodes to the components in Figure 4 are the electrodes that are part of the set of 11 electrodes. The first component, shown in Figure 5A, is similar to the component shown in Figure 4A and refers to P2-N2 transition as well as the P3 peak though with more parietal activation in the fMRI. The second component, with the factor vectors in time and voxel modes shown in Figure 5B, is similar to the component shown in Figure 4B and refers to the P3 peak though has more parietal activation in the fMRI, similar to the default mode network. When significant components from 11- and 62-electrode cases are compared, FMS is 0.82 for EEG (excluding the electrode mode) and 0.73 for fMRI for the component given in Figure 4A, and 0.80 for EEG and 0.72 for fMRI for Figure 4B. These scores indicate that components are similar to some extent but are not identical. Table 1 indicates slightly higher clustering performance for the 11-electrode case. To summarize, considering that there is minimal difference in performance beyond a slight increase in significance, and that the factors are similar, using all electrodes is preferable to choosing a subset of electrodes, as the latter requires prior knowledge about the functionally relevant electrodes to select and may also introduce a bias by targeting specific regions.

FIGURE 5

Figure 5. Temporal and spatial patterns from the statistically significant components of the ACMTF model of the EEG tensor with 11 electrodes and fMRI data, with R = 10. (A) Component 10 corresponds to the P2-N2 transition as well as the P3 peak in the EEG with an increase in sensorimotor and parietal activation of controls over patients in the fMRI, (B) component 5 refers to the P3 peak with a decrease in default mode activity of controls vs. patients in the fMRI, (C) weights of the components in each data set. The corresponding p-values are 4.3 × 10⁻³, 8.0 × 10⁻³, respectively. Columns of the factor matrix in the time mode are in red while green plots show signals from individual electrodes averaged across all subjects.

3.2.1.3. ACMTF vs. JICA

The traditional fusion approach jICA can also capture components that can differentiate between healthy controls and patients; however, jICA provides less interpretable patterns. For jICA, the EEG tensor unfolded in the subject mode is concatenated with the fMRI data resulting in a 32 (subject) by 88,148 (time × electrode–voxel) matrix. When this matrix is modeled using jICA with R = 10, 15, 20 components, the 10-component model reveals a single component that may be considered statistically significant but the p-value is 0.05. The 15-component model reveals a more significant component as illustrated in Figure 6. JICA captures neither a single temporal pattern for all electrodes nor a spatial pattern for all time points, making the interpretation of the components more difficult. In order to provide insight into the topology, spatial patterns as in ACMTF can be computed post hoc from the analysis, e.g., by using peak value of each channel to construct a spatial map for each component (Liu et al., 2009); however, that comes with additional assumptions and does not reveal the underlying patterns as naturally as a tensor factorization-based approach. The component, whose source signals corresponding to the time and voxel parts are shown in Figure 6, is related to the P2-N2 and N2-P3 transition and the fMRI map includes the expected temporal lobe and default mode regions. We should note that using a 20-component model, a component similar to the one in Figure 6 is captured, indicating that jICA has some stability in regards to the value of R. Table 1 shows that while the clustering performance of jICA is lower compared with ACMTF models of EEG and fMRI data sets for R = 10 and R = 15, it is similar for R = 20, indicating that methods with different assumptions may perform the best with different number of components.

FIGURE 6

Figure 6. The statistically significant component captured by the jICA model fitted to the concatenation of the unfolded EEG tensor (with 62 electrodes) and fMRI, with R = 15. The corresponding p-value is 5.2 × 10⁻³. Parts of s_r, i.e., rth column of the source matrix S in Equation (2), corresponding to the time samples for each electrode in EEG and voxels in fMRI are plotted. The EEG part is related to the P2-N2 and N2-P3 transitions. The fMRI indicates some increased activation in the temporal lobe of controls vs. patients as well as some posterior cingulate representing the default mode network. In the EEG plot, green dashed plots show signals from individual electrodes averaged across all subjects.

3.2.1.4. Parameter selection

The ACMTF model is sensitive to two parameters, i.e., the number of components, R, and the sparsity penalty parameter, β. So far, R and β are set to R = 10 and β = 10⁻³. In order to probe the effect of the model order, R, we have increased the number of components until the model fails to give a unique solution. We find that as we increase the number of components, the most significant component (i.e., Figure 4A) is still consistently captured; however, with a decreasing level of similarity. With both R = 11 and R = 12, the ACMTF model is still unique and reveals significant components that can differentiate between patients and healthy controls. The fit is 81% for EEG and 68% for fMRI with R = 11 while 82% for EEG and 70% for fMRI with R = 12, indicating that additional components do not explain much additional information. Table 1 shows that a component with a FMS score around 0.90 when compared with the most significant component in a 10-component model (Figure 4A) has been revealed by both models. The less significant component (Figure 4B), on the other hand, has limited similarity of around FMS 0.60 with the captured components. The clustering performance for both R = 11 and R = 12 is similar to that of R = 10. Figure S1 illustrates the factor vectors, in the time, electrode, and voxel modes, of the significant components captured by the ACMTF model with R = 12 components. When R is increased any further, we cannot obtain a unique solution.

The sensitivity of an ACMTF model to different values of β has been studied in Acar et al. (2014) using simulated data sets, and it has been shown that in the presence of both shared and unshared components, small values, such as β = 10⁻³ or β = 10⁻⁴ are effective in terms of uniquely recovering the underlying patterns. For β = 0, which corresponds to a CMTF model, the model fails to give a unique solution in the presence of unshared components (Acar et al., 2014). For larger values, such as β = 10⁻², it is still possible to find the true solution but the algorithm is very sensitive to the initialization. When EEG and fMRI data sets are jointly analyzed using a 10-component model, weights of the components shown in Figure 4C indicate that all components are shared. However, even in the presence of only shared components, β = 10⁻⁴ fails to give a unique solution, which is due to the fact that without an effective sparsity penalty term, the model can become degenerate. For β = 10⁻², it is not possible to reach to the same function values even with many random initializations due to the sensitivity to initialization.

3.2.1.5. Preprocessing

In addition to the preprocessing steps described in section 2, in data fusion studies, further preprocessing may be needed, in particular, when the average behavior across subjects accounts for a large variation in one of the data sets vs. the other. SMRI is such a data set and we perform additional centering when we include sMRI in the analysis. Here, when we only consider the joint analysis of EEG and fMRI, an additional centering step across the subjects mode does not affect the clustering performance of the ACMTF model and the significant component in Figure 4A has also been captured with FMS between 0.80 and 0.85 (for different number of components R = 9, 10, 11) as shown in Table 1. FMS drops for R = 12. The less significant component (Figure 4B), on the other hand, is also estimated, but with FMS within the range 0.55–0.69. Figure S2 illustrates the significant components captured by an ACMTF model with R = 12 components. It is important to note that in this case, despite the low FMS values, temporal and spatial patterns are similar to the ones observed in Figure 4 and the interpretation of these two components is the same.

3.3. Joint Analysis of EEG, fMRI, and sMRI

Inclusion of the sMRI data introduces several issues highlighting challenges in data fusion, in particular, preprocessing. If a joint analysis of EEG, fMRI, and sMRI data is carried out using an ACMTF model after the preprocessing steps described in section 2.2.3, the model has two statistically significant components in terms of differentiating healthy controls and patients. Figure 7 illustrates the temporal patterns as well as functional/structural spatial patterns revealed by the significant components. Both components are similar to the components shown in Figure 4 with FMS values between 0.84 and 0.93 as given in Table 1 as well as components shown in Figure 3. However, now information from the two functional modalities, EEG and fMRI, has been combined with information from the structural modality, sMRI. In Figure 7A, the sMRI portion of the component shows an increase in concentration of gray matter in controls over patients in sections of the parietal lobe and cerebellum. In Figure 7B, the sMRI portion of the component shows increases in gray matter for controls over patients in multiple portions of the frontal and parietal lobes. Overall these components are similar to components found previously in other analyses of data from the same subjects but with 11 electrodes (Acar et al., 2017a).

FIGURE 7

Figure 7. Temporal and spatial patterns from the statistically significant components of the ACMTF model of the EEG tensor with 62 electrodes, fMRI and sMRI data, with R = 10. (A) Component 10 refers to the P2-N2 transition as well as the P3 peak, heavily weighted by the frontal and central electrodes in the EEG, an increase in sensorimotor and parietal activation of controls over patients in the fMRI, and an increase in concentration of gray matter in controls over patients in sections of the parietal lobe and cerebellum in the sMRI, (B) component 6 corresponds to the P3 peak heavily weighted the frontal electrodes in the EEG, the fMRI indicates a decrease in activation of controls vs. patients in some regions of the sensorimotor cortex and parietal lobe, and increases in gray matter for controls over patients in multiple portions of the frontal and parietal lobes in the sMRI. The corresponding p-values are 8.0 × 10⁻³, 4.0 × 10⁻³, respectively. Columns of the factor matrix in the time mode are in red while green plots show signals from individual electrodes averaged across all subjects.

Weights of the components, shown in Figure 8, indicate that sMRI does not contribute much to the significant components. In most components the weights of the components in the sMRI are low, indicating potentially unshared factors. However, a closer look at the model reveals that the model fit is 97% for sMRI, while it is 79% for EEG and 64% for fMRI (see Figure S3 for the singular value spectrum of each data set), and components with high weights in sMRI are mainly modeling the average structure across subjects with highly correlated components, i.e., the correlation of component vector in the voxel mode is 0.95 for components 1 and 5. Therefore, data sets must be centered across the subject mode to incorporate information other than the mean from sMRI into the analysis. Furthermore, as a result of including components from sMRI, which do not differentiate between the groups, we also observe a drop in the clustering performance in Table 1. Note that the clustering performance is based on all possible combinations of the columns of the factor matrix in the subject mode. However, we do note that the significance of both components have increased, thus indicating that there is additional discriminatory information that the sMRI is providing.

FIGURE 8

Figure 8. Weights of the components in EEG, fMRI, and sMRI extracted by the ACMTF model of the EEG tensor with 62 electrodes, fMRI and sMRI data, with R = 10 and no additional centering across the subject mode.

When all data sets are centered across the subject mode, the ACMTF model has three statistically significant components, which are illustrated in Figure S4. The two most significant ones, with the factor vectors in time, electrode, fMRI voxel, sMRI voxel modes shown in Figures S4A,B, are similar to the components in Figure 4, also indicated by the FMS values in Table 1. The factor vector of the third component in the time mode represents the N2 peak as well as the P3 peak in the EEG. The topographic map indicates activation in the parietal lobe, while the fMRI part shows increased activation of controls over patients in the sensorimotor cortex and a decrease in activation of controls vs. patients in the occipital lobe. The sMRI portion of the component indicates changes to gray matter concentration throughout the frontal and parietal lobes. Note that the clustering performance of the overall model has improved significantly by modeling more relevant structure in sMRI compared to the case where there is no centering. The model fit is 69, 47, and 70% for EEG, fMRI, and sMRI, respectively. In order to increase the model fits, in particular for fMRI, when we increase the number of components to R = 15, a unique model can still be obtained with model fits 78, 64, and 79% for EEG, fMRI, and sMRI, respectively. In that case, however, only a single statistically significant component (p-value=1.4 × 10⁻⁴) that is similar to Figure 4B, is captured. The clustering performance of the 15-component ACMTF model is similar to the 10-component case. These observations indicate that the model is consistent to some degree across models with different numbers of components, by still revealing one of the significant components. The additional components explain some of the remaining information in the data sets but do not reveal any additional significant components in terms of distinguishing between the groups. The results show that they still help with increasing the significance of the relevant component by modeling the structure not contributing to the discrimination of patients and controls.

4. Discussion

In this paper, we have jointly analyzed multi-modal neuroimaging signals, namely, EEG, fMRI, and sMRI, collected from patients with schizophrenia and healthy controls, using a structure-revealing CMTF model. The model captures temporal as well as functional/structural spatial patterns that can differentiate between patients and healthy controls. Compared to traditional fusion approaches, such as jICA, the structure-revealing CMTF model enables us to exploit the multilinear structure of multi-channel EEG signals providing both interpretable patterns and improved uniqueness properties without imposing additional constraints on the extracted patterns. Through joint analysis of EEG, and fMRI, the following temporal and spatial patterns are identified as potential biomarkers:

• Pattern 1: The temporal part referring to the P2-N2 transition as well as the P3 peak, and the functional spatial part showing increased activation of controls over patients in the sensorimotor cortex.

• Pattern 2: The temporal part referring to the P3 peak, and the functional spatial part indicating a decrease in activation of controls vs. patients in some regions of the sensorimotor cortex and parietal lobe.

The biomarkers that are extracted using the ACMTF model correspond to signals observed in previous investigations of the structural and functional impacts of schizophrenia. The EEG signals are similar to those observed in previous schizophrenia research (Calhoun et al., 2010; Kayser et al., 2010). Additionally, the regions in spatial patterns have also been shown to be affected in patients with schizophrenia previously (Schroder et al., 1999; Wolf et al., 2008; Minzenberg et al., 2009). Through the incorporation of the sMRI data, these patterns have been complemented with the following structural spatial parts: (i) Pattern 1, the structural spatial part indicating an increase in concentration of gray matter in controls over patients in sections of the parietal lobe and cerebellum. (ii) Pattern 2, the structural spatial part showing increases in gray matter for controls over patients in multiple portions of the frontal and parietal lobes. All three regions have been shown to be impacted in patients with schizophrenia (Shenton et al., 2010; Olabi et al., 2011; Brent et al., 2013; Lungu et al., 2013). These patterns are reproducible and have been revealed even in the case of changes in data sets, as we have illustrated by leaving out data from one subject at a time and in our preliminary studies on a slightly different set of subjects (Acar et al., 2017b) and using a subset of electrodes (Acar et al., 2017a).

Any method targeting biomarker discovery must capture the underlying patterns corresponding to the potential biomarkers uniquely; therefore, in this paper, we have used the structure-revealing CMTF model that focuses on unique identification of underlying patterns when jointly analyzing multi-modal data sets with shared and unshared factors, rather than other CMTF methods that have proved useful in missing data estimation applications (where uniqueness of underlying patterns is not of interest) (Zheng et al., 2010; Ermis et al., 2015). In addition to the patterns interpreted as potential biomarkers, the structure-revealing CMTF model also reveals weights of the patterns that can be used to identify shared/unshared patterns in each data set and quantify the contribution from each data set. In joint analysis of EEG and fMRI, all components, including the statistically significant components differentiating between patients and controls, correspond to shared components. Similarly, in joint analysis of EEG, fMRI, and sMRI, as long as necessary preprocessing steps, such as centering are carried out as discussed in section 3.3, all components are shared among three modalities as shown in Figure S5.

The structure-revealing CMTF model relies on the hard coupling assumption, i.e., the factor vectors in the coupled (e.g., subject) mode are exactly the same in different modalities for the shared factors. That is a strong assumption, which may be relaxed using soft coupling approaches (Seichepine et al., 2013; Farias et al., 2016), e.g., having different factor matrices in the coupled mode and penalizing the distance between the factor matrices based on various norms, as in Rivet et al. (2015) when jointly analyzing EEG and electro-ocular artifacts, or as in Chatzichristos et al. (2018) when jointly analyzing EEG and raw fMRI signals. However, a hard coupling-based approach may be preferable in noisy cases (Adali et al., 2015b). Therefore, whether or not it is a limitation is an open research question.

One potential drawback of the structure-revealing CMTF model is sensitivity to its parameters, i.e., the number of components, R, and the sparsity penalty parameter, β. Despite the sensitivity, we have consistently observed similar statistically significant patterns for different number of components, as shown in Table 1 in terms of FMS values. Note that FMS takes into account every entry in the factor vectors (i.e., many voxels in the fMRI) and it is rather a strict measure. Therefore, we have observed that even for lower FMS values, interpretations of the captured patterns are the same visually (e.g., Figure 4 vs. Figure S2). Another challenge as a result of sensitivity to model parameters is that the model must be experimentally validated to be unique. An important future research direction is to study the landscape of the optimization problem and develop ways to make the problem less sensitive to parameters as well as the initialization. A theoretical understanding of the uniqueness properties of the structure-revealing CMTF model would also significantly mitigate this challenge. Recent work on uniqueness (Lathauwer and Kofidis, 2017) of coupled matrix and tensor factorizations provides a step in this direction studying the indeterminacies in CMTF in the presence of unshared factors. Furthermore, while we have used the same sparsity penalty parameter for all data sets in this paper, in some applications, this parameter may need to be data-specific depending on the structure of each data set.

While data fusion methods are of interest in many disciplines, preprocessing steps have not been well-studied within the framework of data fusion. In this paper, we have shown that while centering across the subject mode does not affect the joint analysis of EEG and fMRI data, it has a dramatic effect when the sMRI data is incorporated, and the interpretation of the component weights changes significantly. In addition to such preprocessing steps, there are further steps that should be carefully incorporated to data fusion methods, such as outlier removal, residual analysis, which may also enable better tools for selecting the number of components.

This paper is a systematic study of the structure-revealing CMTF model for biomarker discovery but with limited number of subjects. In order to see the real promise of the method as a biomarker discovery approach and assess the validity of the potential biomarkers for schizophrenia, joint analysis of EEG, fMRI, and sMRI signals must be carried out on a larger set of subjects and also including patients with different neurological disorders. There are increasingly more studies exploiting the multilinear structure of different neuroimaging signals (Andersen and Rayens, 2004; Davidson et al., 2013; He et al., 2014; Ferdowsi et al., 2015; Madsen et al., 2017; Chatzichristos et al., 2018, 2019; Williams et al., 2018), and understanding the advantages and limitations of CMTF-based approaches, that can exploit those multilinear structures, is crucial as the fusion of multimodal neuroimaging data holds promise to enhance biomarker discovery.

Ethics Statement

The protocol was reviewed and approved by the IRB at Hartford Hospital and all participants provided written, informed, consent.

Author Contributions

EA and TA conceived the project. EA, YL-S, and TA designed the experiments. EA and CS performed the data analysis. YL-S, TA, and VC interpreted the extracted patterns. EA, CS, YL-S, and TA were involved in the preparation of the manuscript. All authors have given approval to the final version of the manuscript.

Funding

This work has been funded in part by the following grants: NSF-CCF 1618551, NSF-NCS 1631838, NIH-NIBIB R01EB005846, and NIH-T32 HD049311.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2019.00416/full#supplementary-material

Footnotes

1. ^https://www.fil.ion.ucl.ac.uk/spm/

2. ^http://www.sandia.gov/tgkolda/TensorToolbox/

3. ^http://www.models.life.ku.dk/joda/CMTF_Toolbox

4. ^Up to the sixth decimal place.

5. ^Depending on the difficulty of the problem, 160-200 random initializations are used to check for uniqueness.

6. ^The ICA-EBM code is available at http://mlsp.umbc.edu/ica_ebm.html.

7. ^T-tests reveal, in total, four statistically significant components for this model. However, the other two components, not illustrated in Figure 5, are less significant and have lower FMS values.

References

Acar, E., Bingol, C. A., Bingol, H., Bro, R., and Yener, B. (2007). Multiway analysis of epilepsy tensors. Bioinformatics. 23, i10–i18. doi: 10.1093/bioinformatics/btm210

PubMed Abstract | CrossRef Full Text | Google Scholar

Acar, E., Bro, R., and Smilde, A. K. (2015). Data fusion in metabolomics using coupled matrix and tensor factorizations. Proc. IEEE. 103, 1602–1620. doi: 10.1109/JPROC.2015.2438719

CrossRef Full Text | Google Scholar

Acar, E., Dunlavy, D. M., and Kolda, T. G. (2011). A scalable optimization approach for fitting canonical tensor decompositions. J. Chemometr. 25, 67–86. doi: 10.1002/cem.1335

CrossRef Full Text | Google Scholar

Acar, E., Levin-Schwartz, Y., Calhoun, V., and Adali, T. (2017a). “ACMTF for fusion of multi-modal neuroimaging data and identification of biomarkers,” in Eusipco'17: Proceedings of the 25th European Signal Processing Conference (Kos), 673–677. doi: 10.23919/EUSIPCO.2017.8081286

CrossRef Full Text | Google Scholar

Acar, E., Levin-Schwartz, Y., Calhoun, V. D., and Adali, T. (2017b). “Tensor-based fusion of EEG and FMRI to understand neurological changes in schizophrenia,” in ISCAS '17: Proceedings of IEEE International Symposium on Circuits and Systems (Baltimore, MD), 314–317. doi: 10.1109/ISCAS.2017.8050303

CrossRef Full Text | Google Scholar

Acar, E., Papalexakis, E. E., Gurdeniz, G., Rasmussen, M. A., Lawaetz, A. J., Nilsson, M., et al. (2014). Structure-revealing data fusion. BMC Bioinform. 15:239. doi: 10.1186/1471-2105-15-239

PubMed Abstract | CrossRef Full Text | Google Scholar

Adali, T., Schwartz, Y. L., and Calhoun, V. D. (2015a). Multimodal data fusion using source separation: application to medical imaging. Proc. IEEE. 103, 1494–1506. doi: 10.1109/JPROC.2015.2461601

CrossRef Full Text | Google Scholar

Adali, T., Schwartz, Y. L., and Calhoun, V. D. (2015b). Multimodal data fusion using source separation: two effective models based on ica and iva and their properties. Proc. IEEE. 103, 1478–1493. doi: 10.1109/JPROC.2015.2461624

CrossRef Full Text | Google Scholar

Alter, O., Brown, P. O., and Botstein, D. (2003). Generalized singular value decomposition for comparative analysis of genome-scale expression data sets of two different organisms. PNAS. 100, 3351–3356. doi: 10.1073/pnas.0530258100

PubMed Abstract | CrossRef Full Text | Google Scholar

Andersen, A. H., and Rayens, W. S. (2004). Structure-seeking multilinear methods for the analysis of fMRI data. Neuroimage. 22, 728–739. doi: 10.1016/j.neuroimage.2004.02.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashburner, J., and Friston, K. J. (2000). Voxel-based morphometry–the methods. Neuroimage. 11, 805 – 821. doi: 10.1006/nimg.2000.0582

PubMed Abstract | CrossRef Full Text | Google Scholar

Becker, H., Comon, P., and Albera, L. (2012). “Tensor-based processing of combined EEG/MEG data,” in Eusipco'12: Proceedings of the 20th European Signal Processing Conference (Bucharest), 275–279.

Google Scholar

Brent, B. K., Thermenos, H. W., Keshavan, M. S., and Seidman, L. J. (2013). Gray matter alterations in schizophrenia high-risk youth and early-onset schizophrenia. Child Adolesc. Psychiatr. Clin. N. Am. 22, 689–714. doi: 10.1016/j.chc.2013.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bro, R., and Smilde, A. K. (2003). Centering and scaling in component analysis. J. Chemometr. 17, 16–33. doi: 10.1002/cem.773

CrossRef Full Text | Google Scholar

Bunge, S. A., and Kahn, I. (2009). Cognition: an overview of neuroimaging techniques. Encyc. Neurosci. 2, 1063–1067. doi: 10.1016/B978-008045046-9.00298-9

CrossRef Full Text | Google Scholar

Calhoun, V., Adali, T., Giuliani, N., Pekar, J., Kiehl, K., and Pearlson, G. (2006a). Method for multimodal analysis of independent source differences in schizophrenia: combining gray matter structural and auditory oddball functional data. Hum. Brain Mapp. 27, 47–62. doi: 10.1002/hbm.20166

PubMed Abstract | CrossRef Full Text | Google Scholar

Calhoun, V., Adali, T., Pearlson, G., and Kiehl, K. (2006b). Neuronal chronometry of target detection: fusion of hemodynamic and event-related potential data. Neuroimage. 30, 544–553. doi: 10.1016/j.neuroimage.2005.08.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Calhoun, V., Wu, L., Kiehl, K., Eichele, T., and Pearlson, G. (2010). Aberrant processing of deviant stimuli in schizophrenia revealed by fusion of FMRI and EEG data. Acta Neuropsychiatr. 22, 127–138. doi: 10.1111/j.1601-5215.2010.00467.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Calhoun, V. D., Liu, J., and Adali, T. (2009). A review of group ICA for fMRI data and ICA for joint inference of imaging, genetic, and ERP data. Neuroimage. 45, S163–S172. doi: 10.1016/j.neuroimage.2008.10.057

PubMed Abstract | CrossRef Full Text | Google Scholar

Carroll, J. D., and Chang, J. J. (1970). Analysis of individual differences in multidimensional scaling via an N-way generalization of “Eckart-Young” decomposition. Psychometrika. 35, 283–319.

Google Scholar

Chatzichristos, C., Davies, M., Escudero, J., Kofidis, E., and Theodoridis, S. (2018). “Fusion of EEG and fMRI via soft coupled tensor decompositions,” in Eusipco'18: Proceedings of the 26th European Signal Processing Conference (Rome). doi: 10.23919/EUSIPCO.2018.8553077

CrossRef Full Text | Google Scholar

Chatzichristos, C., Kofidis, E., Morante, M., and Theodoridis, S. (2019). Blind fMRI source unmixing via higher-order tensor decompositions. J. Neurosci. Methods. 315, 17–47. doi: 10.1016/j.jneumeth.2018.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Cong, F., Lin, Q.-H., Kuang, L.-D., Gong, X.-F., Astikainen, P., and Ristaniemi, T. (2015). Tensor decomposition of EEG signals: a brief review. J. Neurosci. Methods. 248, 59–69. doi: 10.1016/j.jneumeth.2015.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Dahne, S., Biebmann, F., Samek, W., Haufe, S., Goltz, D., Gundlach, C., et al. (2015). Multivariate machine learning methods for fusing multimodal functional neuroimaging data. Proc. IEEE. 103, 1507–1530. doi: 10.1109/JPROC.2015.2425807

CrossRef Full Text | Google Scholar

Daunizeau, J., Laufs, H., and Friston, K. J. (2009). “EEG–fMRI information fusion: Biophysics and data analysis,” in EEG–fMRI, eds C. Mulert and L. Lemieux (Berlin; Heidelberg: Springer), 511–526.

Google Scholar

Davidson, I., Gilpin, S., Carmichael, O., and Walker, P. (2013). “Network discovery via constrained tensor analysis of fMRI data,” in KDD '13: Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (Chicago, IL), 194–202. doi: 10.1145/2487575.2487619

CrossRef Full Text | Google Scholar

De Vos, M., Vergult, A., De Lathauwer, L., De Clercq, W., Van Huffel, S., Dupont, P., et al. (2007). Canonical decomposition of ictal scalp EEG reliably detects the seizure onset zone. Neuroimage. 37, 844–854. doi: 10.1016/j.neuroimage.2007.04.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Delorme, A., and Makeig, S. (2004). EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics. J. Neurosci. Methods. 134, 9–21. doi: 10.1016/j.jneumeth.2003.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, W., Levin-Schwartz, Y., Fu, G. D., Ma, S., Calhoun, V. D., and Adali, T. (2016). The role of diversity in complex ICA algorithms for fMRI analysis. J. Neurosci. Methods. 264, 129–135. doi: 10.1016/j.jneumeth.2016.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermis, B., Acar, E., and Cemgil, A. T. (2015). Link prediction in heterogeneous data via generalized coupled tensor factorization. Data Mining Knowl. Discov. 29, 203–236. doi: 10.1007/s10618-013-0341-y

CrossRef Full Text | Google Scholar

Eyndhoven, S. V., Hunyadi, B., De Lathauwer, L., and Van Huffel, S. (2017). “Flexible data fusion of EEG-fMRI: revealing neural-hemodynamic coupling through structured matrix-tensor factorization,” in Eusipco'17: Proceedings of the 25th European Signal Processing Conference (Kos), 26–30. doi: 10.23919/EUSIPCO.2017.8081162

CrossRef Full Text | Google Scholar

Farias, R. C., Cohen, J. E., and Comon, P. (2016). Exploring multimodal data fusion through joint decompositions with flexible couplings. IEEE Trans. Signal Process. 64, 4830–4844. doi: 10.1109/TSP.2016.2576425

CrossRef Full Text | Google Scholar

Ferdowsi, S., Abolghasemi, V., and Sanei, S. (2015). A new informed tensor factorization approach to EEG–fMRI fusion. J. Neurosci. Methods. 254, 27–35. doi: 10.1016/j.jneumeth.2015.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

First, M. B., Spitzer, R. L., Miriam, G., and Williams, J. B. (2002a). Structured Clinical Interview for DSM-IV: Non-Patient Edition (SCID-NP). New York, NY: Biometrics Research, New York State Psychiatric Institute.

Google Scholar

First, M. B., Spitzer, R. L., Miriam, G., and Williams, J. B. (2002b). Structured Clinical Interview for DSM-IV-TR Axis I Disorders, Research Version, Patient Edition. (SCID-I/P). New York, NY: Biometrics Research, New York State Psychiatric Institute.

Google Scholar

Freire, L., and Mangin, J.-F. (2001). Motion correction algorithms may create spurious brain activations in the absence of subject motion. Neuroimage. 14, 709–722. doi: 10.1006/nimg.2001.0869

PubMed Abstract | CrossRef Full Text | Google Scholar

Good, C. D., Johnsrude, I. S., Ashburner, J., Henson, R. N., Friston, K. J., and Frackowiak, R. S. (2001). A voxel-based morphometric study of ageing in 465 normal adult human brains. Neuroimage. 14, 21–36. doi: 10.1006/nimg.2001.0786

PubMed Abstract | CrossRef Full Text | Google Scholar

Groves, A. R., Beckmann, C. F., Smith, S. M., and Woolrich, M. W. (2011). Linked independent component analysis for multimodal data fusion. Neuroimage. 54, 2198–2217. doi: 10.1016/j.neuroimage.2010.09.073

PubMed Abstract | CrossRef Full Text | Google Scholar

Harshman, R. A. (1970). Foundations of the PARAFAC procedure: models and conditions for an “explanatory” multi-modal factor analysis. UCLA Working Papers Phonetics. 16, 1–84.

Google Scholar

He, L., Kong, X., Yu, P. S., Ragin, A. B., Hao, Z., and Yang, X. (2014). “Dusk: a dual structure-preserving kernel for supervised tensor learning with applications to neuroimages,” in SDM'14: Proceedings of the SIAM International Conference on Data Mining, 127–135. doi: 10.1137/1.9781611973440.15

PubMed Abstract | CrossRef Full Text | Google Scholar

Hitchcock, F. L. (1927). The expression of a tensor or a polyadic as a sum of products. J. Math. Phys. 6, 164–189.

PubMed Abstract | Google Scholar

Hunyadi, B., Dupont, P., Paesschen, W. V., and Huffel, S. V. (2017). Tensor decompositions and data fusion in epileptic electroencephalography and functional magnetic resonance imaging data. WIREs Data Mining Knowl. Discov. 7:e1197. doi: 10.1002/widm.1197

CrossRef Full Text | Google Scholar

Hunyadi, B., Paesschen, W. V., De Vos, M., and Van Huffel, S. (2016). “Fusion of electroencephalography and functional magnetic resonance imaging to explore epileptic network activity,” in Eusipco' 16: Proceedings of the 24th European Signal Processing Conference (Budapest), 240–244. doi: 10.1109/EUSIPCO.2016.7760246

CrossRef Full Text | Google Scholar

Jung, T.-P., Makeig, S., Humphries, C., Lee, T.-W., McKeown, M. J., Iragui, V., et al. (2000). Removing electroencephalographic artifacts by blind source separation. Psychophysiology. 37, 163–178. doi: 10.1111/1469-8986.3720163

PubMed Abstract | CrossRef Full Text | Google Scholar

Karahan, E., Rojas-Lopez, P. A., Bringas-Vega, M. L., Valdes-Hernandez, P. A., and Valdes-Sosa, P. A. (2015). Tensor analysis and fusion of multimodal brain images. Proc. IEEE. 103, 1531–1559. doi: 10.1109/JPROC.2015.2455028

CrossRef Full Text | Google Scholar

Kayser, J., Tenke, C. E., Gil, R., and Bruder, G. E. (2010). ERP generator patterns in schizophrenia during tonal and phonetic oddball tasks: effects of response hand and silent count. Clin. EEG Neurosci. 41, 184–195. doi: 10.1177/155005941004100405

PubMed Abstract | CrossRef Full Text | Google Scholar

Kiehl, K. A., Stevens, M. C., Laurens, K. R., Pearlson, G., Calhoun, V. D., and Liddle, P. F. (2005). An adaptive reflexive processing model of neurocognitive function: supporting evidence from a large scale (n = 100) fmri study of an auditory oddball task. Neuroimage. 25, 899–915. doi: 10.1016/j.neuroimage.2004.12.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Klami, A., Virtanen, S., and Kaski, S. (2013). Bayesian canonical correlation analysis. J. Mach. Learn. Res. 14, 965–1003. Available online at: http://www.jmlr.org/papers/v14/klami13a.html

Google Scholar

Kolda, T. G., and Bader, B. W. (2009). Tensor decompositions and applications. SIAM Rev. 51, 455–500. doi: 10.1137/07070111X

CrossRef Full Text | Google Scholar

Kruskal, J. B. (1977). Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics. Linear Algebra Appl. 18, 95–138. doi: 10.1016/0024-3795(77)90069-6

CrossRef Full Text | Google Scholar

Lathauwer, L. D., and Kofidis, E. (2017). “Coupled matrix-tensor factorizations–the case of partially shared factors,” in ASILOMAR'17: Proceedings of the Asilomar Conference on Signals, Systems and Computers (Pacific Grove, CA). doi: 10.1109/ACSSC.2017.8335436

CrossRef Full Text | Google Scholar

Li, X.-L., and Adali, T. (2010). Independent component analysis by entropy bound minimization. IEEE Trans. Signal Process. 58, 5151–5164. doi: 10.1109/TSP.2010.2055859

CrossRef Full Text | Google Scholar

Liu, J., Kiehl, K. A., Pearlson, G., Perrone-Bizzozero, N. I., Eichele, T., and Calhoun, V. D. (2009). Genetic determinants of target and novelty-related event-related potentials in the auditory oddball response. Neuroimage. 46, 809–816. doi: 10.1016/j.neuroimage.2009.02.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Liu, S., Cai, W., Che, H., Pujol, S., Kikinis, R., et al. (2015). Multimodal neuroimaging feature learning for multiclass diagnosis of Alzheimer's disease. IEEE Trans. Biomed. Eng. 62, 1132–1140. doi: 10.1109/TBME.2014.2372011

PubMed Abstract | CrossRef Full Text | Google Scholar

Lock, E. F., Hoadley, K. A., Marron, J., and Nobel, A. B. (2013). Joint and individual variation explained (jive) for integrated analysis of multiple data types. Ann. Appl. Stat. 7, 523–542. doi: 10.1214/12-AOAS597

PubMed Abstract | CrossRef Full Text | Google Scholar

Lungu, O., Barakat, M., Laventure, S., Debas, K., Proulx, S., Luck, D., et al. (2013). The incidence and nature of cerebellar findings in schizophrenia: a quantitative review of fMRI literature. Schizophr. Bull. 39, 797–806. doi: 10.1093/schbul/sbr193

PubMed Abstract | CrossRef Full Text | Google Scholar

Madsen, K. H., Churchill, N. W., and Mørup, M. (2017). Quantifying functional connectivity in multi-subject fMRI data using component models. Hum. Brain Mapp. 38, 882–899. doi: 10.1002/hbm.23425

PubMed Abstract | CrossRef Full Text | Google Scholar

Minzenberg, M. J., Laird, A. R., Thelen, S., Carter, C. S., and Glahn, D. C. (2009). Meta-analysis of 41 functional neuroimaging studies of executive function in schizophrenia. Arch. Gen. Psychiatry. 66, 811–822. doi: 10.1001/archgenpsychiatry.2009.91

PubMed Abstract | CrossRef Full Text | Google Scholar

Miwakeichi, F., Martínez-Montes, E., Valdés-Sosa, P. A., Nishiyama, N., Mizuhara, H., and Yamaguchi, Y. (2004). Decomposing EEG data into space-time-frequency components using parallel factor analysis. Neuroimage. 22, 1035–1045. doi: 10.1016/j.neuroimage.2004.03.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Möcks, J. (1988). Topographic components model for event-related potentials and some biophysical considerations. IEEE Trans. Biomed. Eng. 35, 482–484.

PubMed Abstract | Google Scholar

Mørup, M., Hansen, L. K., and Arnfred, S. M. (2007). ERPWAVELAB a toolbox for multi-channel analysis of time - frequency transformed event related potentials. J. Neurosci. Methods. 161, 361–368. doi: 10.1016/j.jneumeth.2006.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Naskovska, K., Korobkov, A. A., Haardt, M., and Haueisen, J. (2017). “Analysis of the photic driving effect via joint EEG and MEG data processing based on the coupled CP decomposition,” in Eusipco'17: Proceedings of the 25th European Signal Processing Conference (Kos), 1325–1329. doi: 10.23919/EUSIPCO.2017.8081415

CrossRef Full Text | Google Scholar

Olabi, B., Ellison-Wright, I., McIntosh, A. M., Wood, S. J., Bullmore, E., and Lawrie, S. M. (2011). Are there progressive brain changes in schizophrenia? A meta-analysis of structural magnetic resonance imaging studies. Biol. Psychiatry. 70, 88–96. doi: 10.1016/j.biopsych.2011.01.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Rivet, B., Duda, M., Guerin-Dugue, A., Jutten, C., and Comon, P. (2015). “Multimodal approach to estimate the ocular movements during EEG recordings: a coupled tensor factorization method,” in EMBC'15: Proceedings of 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (Milan). doi: 10.1109/EMBC.2015.7319999

CrossRef Full Text | Google Scholar

Schroder, J., Essig, M., Baudendistel, K., Jahn, T., Gerdsen, I., Stockert, A., et al. (1999). Motor dysfunction and sensorimotor cortex activation changes in schizophrenia: a study with functional magnetic resonance imaging. Neuroimage. 9, 81–87.

PubMed Abstract | Google Scholar

Seichepine, N., Essid, S., Févotte, C., and Cappé, O. (2013). “Soft nonnegative matrix co-factorization with application to multimodal speaker diarization,” in ICASSP'13: Proceedings of the IEEE Int. Conf. on Acoustics, Speech and Signal (Vancouver, BC). doi: 10.1109/ICASSP.2013.6638316

CrossRef Full Text | Google Scholar

Shenton, M. E., Whitford, T. J., and Kubicki, M. (2010). Structural neuroimaging in schizophrenia: from methods to insights to treatments. Dialog. Clin. Neurosci. 12, 317–332.

PubMed Abstract | Google Scholar

Sidiropoulos, N. D., and Bro, R. (2000). On the uniqueness of multilinear decomposition of N-way arrays. J. Chemometr. 14, 229–239. doi: 10.1002/1099-128X(200005/06)14:3<229::AID-CEM587>3.0.CO;2-N

CrossRef Full Text | Google Scholar

Smilde, A. K., Mage, I., Naes, T., Hankemeier, T., Lips, M. A., Kiers, H. A., et al. (2017). Common and distinct components in data fusion. J. Chemometr. 31:e2900. doi: 10.1002/cem.2900

CrossRef Full Text | Google Scholar

Sorensen, M., and De Lathauwer, L. (2015). Coupled canonical polyadic decompositions and (coupled) decompositions in multilinear rank −(l_{r, n}, l_{r, n}, 1) terms — part i: Uniqueness. SIAM J. Matrix Anal. Appl. 36, 496–522. doi: 10.1137/140956853

CrossRef Full Text | Google Scholar

Sui, J., Qi, S., van Erp, T. G. M., Bustillo, J., Jiang, R., Lin, D., et al. (2018). Multimodal neuromarkers in schizophrenia via cognition-guided MRI fusion. Nat. Commun. 9:3028. doi: 10.1038/s41467-018-05432-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Sui, J., Yu, Q., He, H., Pearlson, G. D., and Calhoun, V. D. (2012). A selective review of multimodal fusion methods in schizophrenia. Front. Hum. Neurosci. 6:27. doi: 10.3389/fnhum.2012.00027

PubMed Abstract | CrossRef Full Text | Google Scholar

Swinnen, W., Hunyadi, B., Acar, E., Van Huffel, S., and Vos, M. D. (2014). “Incorporating higher dimensionality in joint decomposition of EEG and fMRI,” in Eusipco'14: Proceedings of the 22nd European Signal Processing Conference (Lisbon), 121–125.

Google Scholar

Takahashi, T., and Suzuki, M. (2018). Brain morphologic changes in early stages of psychosis: implications for clinical application and early intervention. Psychiatry Clin. Neurosci. 72, 556–571. doi: 10.1111/pcn.12670

CrossRef Full Text | Google Scholar

Uludag, K., and Roebroeck, A. (2014). General overview on the merits of multimodal neuroimaging data fusion. Neuroimage. 102, 3–10. doi: 10.1016/j.neuroimage.2014.05.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, A. H., Kim, T. H., Wang, F., Vyas, S., Ryu, S. I., Shenoy, K. V., et al. (2018). Unsupervised discovery of demixed, low-dimensional neural dynamics across multiple timescales through tensor component analysis. Neuron. 98, 1099–1115. doi: 10.1016/j.neuron.2018.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolf, D. H., Turetsky, B. I., Loughead, J., Elliott, M. A., Pratiwadi, R., Gur, R. E., et al. (2008). Auditory oddball fMRI in schizophrenia: association of negative symptoms with regional hypoactivation to novel distractors. Brain Imaging Behav. 2, 132–145. doi: 10.1007/s11682-008-9022-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, V. W., Cao, B., Zheng, Y., Xie, X., and Yang, Q. (2010). “Collaborative filtering meets mobile recommendation: a user-centered approach,” in AAAI'10: Proceedings of the 24th Conf. on Artificial Intelligence (Atlanta, GA), 236–241.

Google Scholar