Model Order Effects on ICA of Resting-State Complex-Valued fMRI Data: Application to Schizophrenia

Background Component splitting at higher model orders is a widely accepted finding for independent component analysis (ICA) of functional magnetic resonance imaging (fMRI) data. However, our recent study found that intact components occurred with subcomponents at higher model orders.

To date, most studies on model order effects on ICA have considered magnitude-only fMRI data.fMRI data initially are acquired as a bivariate complex-valued form containing magnitude and phase image pairs (Rowe, 2005;Rowe et al., 2007;Adalı and Calhoun, 2007;Calhoun and Adalı, 2012;Hagberg and Tuzzi, 2014).Although phase data tend to be noisy, phase fMRI data may contribute unique biological A C C E P T E D M A N U S C R I P T information.For example, phase data previously have been used to reduce unwanted draining vein (Nencka and Rowe, 2007) and macrovascular (Menon, 2002;Tomasi and Caparelli, 2007) effects on blood oxygenation level-dependent signals.In addition, phase data show task-related changes (Hoogenraad et al., 1998(Hoogenraad et al., , 2001;;Arja et al., 2010).The phase change in activated brain areas can be explained via the sphere of Lorentz effect (Zhao et al., 2007;Feng et al., 2009).Multiple studies also have shown that meaningful brain activations can be extracted from phase fMRI data (Calhoun et al., 2002;Castro et al., 2014;Yu et al., 2015;Chen and Calhoun, 2016;Chen et al., 2017b).fMRI analyses based on evaluation of magnitude-only data can benefit from including information conveyed by phase data (Rodriguez et al., 2011(Rodriguez et al., , 2012;;Castro et al., 2014;Hagberg and Tuzzi, 2014;Yu et al., 2015).For instance, Chen and Calhoun (2016) extracted a task-evoked brain functional map from phase fMRI data using an ICA-based brain functional magnetic susceptibility mapping method.Castro et al. (2014) achieved a 5% improvement in classification accuracy rate for HCs and patients with schizophrenia (SZs) by incorporating features, which were extracted by ICA from phase data, into a multiple kernel learning algorithm.Rodriguez et al. (2011Rodriguez et al. ( , 2012) ) achieved better sensitivity and specificity than the magnitude-only method by proposing quality map phase denoising in the complex-valued method.Yu et al. (2015) detected more (139% for task-related and 331% for default mode) expected activations in a complex-valued method (as opposed to a magnitude-only method) by using a spatial phase mask for post-ICA denoising.Chen et al. (2017b) analyzed functional connectivity in phase data using group ICA and obtained different results from magnitude-only data, which provided a new arena for more comprehensive brain function analysis.
A fixed model order typically is used when performing ICA of complex-valued fMRI data.The model order is larger than that used for magnitude-only data because additional phase information is incorporated (Adalı and Calhoun, 2007;Li et al., 2007;Xiong et al., 2012;Yu et al., 2015).Motivated by model order effects seen in ICA of magnitude-only fMRI data, we were interested in how model order would affect ICA of complex-valued fMRI data.In a preliminary study, we investigated model order A C C E P T E D M A N U S C R I P T effects on independent vector analysis using 16 task-related complex-valued fMRI data sets (Kuang et al., 2017a).Results indicated that complex-valued ICA analysis also detected component splitting at higher model orders but differed from magnitude-only analysis in that an intact component and its subcomponents existed simultaneously.
In this study, we examined model order effects on ICA using resting-state complex-valued fMRI data from 82 subjects (40 HCs and 42 SZs).In addition, we investigated underlying causes of component-splitting differences between complex-valued data and magnitude-only data.We chose DMN-, visual-, and sensorimotor-related networks as components of interest because of clear component-splitting differences at higher model orders in magnitude-only analysis (Abou-Elseoud et al., 2010).We repeated ICA 10 times at varying model orders from 10 to 140 to gather complex-valued data, magnitude-only data, and phase data from each subject.We selected a best run of ICA results to provide DMN-, visual-, and sensorimotor-related components for further analyses.To ensure consistency of ICA results across runs, Du et al. (2014Du et al. ( , 2016) ) proposed a best run selection method using the assignment problem and a minimum spanning tree (MST) to sort components.We identified the best run by the highest average correlation between all estimated components and their corresponding t-maps, which we obtained by performing a one-sample t-test on each component estimate from all subjects.Furthermore, we have proposed an improved method for selecting the best ICA run.The proposed method jointly performs averaging of the spatial maps of a component across all runs and statistical significance tests for each voxel in the spatial maps, which provides a more concise spatial reference.We constructed and used three spatial references, which corresponded to intact DMN-, visual-, and sensorimotor-related components, to select the best run in terms of the highest sum of correlation coefficients between the three reliable spatial references and their corresponding spatial components.We conducted experimental tests to evaluate the efficacy of the proposed best run selection method.
To evaluate model order effects, we obtained the group components by averaging a component of interest in the best run across 82 subjects and then we used these components for assessment.The group components in complex-valued analysis were further denoised using group phase masks (Yu et al., 2015;Kuang et al., 2017b), and the group components in phase analysis were further spatially smoothed for denoising.We examined the model order effects on ICA of complex-valued data by comparing these effects with ICA of magnitude-only data and phase data.We made these comparisons in terms of qualitative observations of DMN-, visual-, and sensorimotor-related group components and quantitative indexes, which included the number of activated voxels and cluster quality.In addition, we demonstrated the advantage of intact networks obtained from complex-valued fMRI data in distinguishing HCs and SZs.
We selected the intact DMN as an example and compared HCs and SZs in terms of the number of activated voxels within several ROIs in group DMN components as well as the significance of differences in voxel number, voxel-wise t-maps (obtained by a two-sample t-test), and voxel-wise variance maps across subjects.

Selection of the best run
Principal component analysis (PCA)-reduced fMRI data of a single-subject are described as  ∈ ℂ × , where N is the model order and V is the number of in-brain voxels.The spatial ICA model of fMRI data is  =  , where  = [ 1 , ⋯ ,   ] ∈ ℂ × includes N spatial maps. = [ 1 , ⋯ ,   ]  ∈ ℂ × is the mixing matrix containing N time course information, and T denotes transpose.The model order N could vary from 1 to 146, which is the number of scans.By learning a demixing matrix  ∈ ℂ × under the assumption of statistical independence of spatial components (McKeown et al., 1998;Calhoun et al., 2001b) where _(•) denotes the p-value of a one-sample t-test on   1 (), ⋯ ,    (), and  ℎ is a p-value threshold (we used  ℎ = 0.01 ).Note that we removed the voxels with negative values from magnitude-only spatial references, because the magnitude maps in complex-valued analysis have only voxels with positive values.Removing the negative values did not substantially affect the results.In equation (1), the spatial reference  ̅  takes advantage of denoising by averaging ICA estimates over all runs and of the consistency of activations using statistical significance tests.

ICA analyses with varying model order
We changed the model orders from 10 to 140 (in increments of 10).We used the entropy bound minimization (EBM) algorithm (Li and Adalı, 2010) to perform complex-valued ICA in complex-valued fMRI analysis, and we used the Infomax algorithm (Bell and Sejnowski, 1995) to perform real-valued ICA in magnitude-only and phase fMRI analyses.The parameters used for EBM and Infomax were the default parameters used in the GIFT toolbox (Medical Image Analysis Lab, University of New Mexico, Albuquerque, NM, US; http://mialab.mrn.org/software/gift/).We repeated ICA 10 times for each subject and selected the best run using our proposed method in Section 2.2.1.
Given the best run of ICA for each subject, we first selected the DMN-, visual-, and sensorimotor-related components by using their spatial reference networks.Because some DMN-, visual-, and sensorimotor-related subcomponents emerged at higher model orders because of component splitting, we also constructed subcomponent references by dividing the reference networks (see Section 2.2.1) into several subnetworks.As a result, we obtained: (1) six total reference networks for DMN-related components, which included the intact DMN as well as subcomponents DMNA, DMNP, IPL, right IPL (RIPL), and left IPL (LIPL), (2) six visual reference networks, which included the intact visual network as well as subnetworks MVS, OPVS, LVS, right LVS (RLVS), and left LVS (LLVS), and (3) four sensorimotor reference networks, which included the intact sensorimotor network as well as SMA, RPMA, and LPMA.We selected these DMN-, visual-, and sensorimotor-related components from the best run of ICA via maximal Pearson correlation coefficients between spatial map estimates and corresponding references as well as visual inspection.
Next, we performed post-ICA processing to obtain the DMN-, visual-, and sensorimotor-related spatial group components for subsequent evaluation.We built group components in complex-valued, magnitude-only, and phase analyses via the following three steps: 1. Averaging each component over all 82 subjects.For complex-valued analysis, we denoised the averaged components using the group phase denoising method (Yu et al., 2015;Kuang et al., 2017b).
For phase analysis, we also removed negative voxels and then spatially smoothed the group components using an 8 × 8 × 8 mm³ FWHM Gaussian kernel.
3. Removing voxels with relatively smaller magnitudes by thresholding at Z > 0.5.
Previous studies have demonstrated that the stability of magnitude-only analysis decreases as model order increases.In this study, we examined the stability of complex-valued analysis, magnitude-only analysis, and phase analysis under increasing model order.We used the cluster quality index in ICASSO, which is software for investigating reliability of ICA estimates by clustering and visualization (Himberg et al., 2003(Himberg et al., , 2004)), to evaluate consistency of a spatial component estimated from different runs of ICA.More specifically, we used the mean cluster quality index  ̅ averaged first over selected components and then over all subjects for evaluation (Kuang et al., 2017a).The definitions of  ̅ and the standard deviation  ̅ averaged across 82 subjects are provided in Appendix B.

Complex-valued DMN component differences between HCs and SZs
To investigate the advantage of intact components obtained in complex-valued analysis at higher model orders, we used the intact group DMN as an example to show its capacity for distinguishing HCs and SZs.
We selected a higher model order N = 120 for evaluation.When compared with magnitude-only analysis, we analyzed a magnitude subcomponent containing dominant posterior areas and weak IPL areas since this subcomponent was the most similar to the intact DMN component.In addition, we compared the intact DMN components obtained at lower model orders using both complex-valued and magnitude-only analyses.Specifically, we selected N = 20 for magnitude-only analysis and N = 40 for complex-valued analysis.We used a larger model order in complex-valued analysis because of the additional phase information that had been incorporated compared with magnitude-only data (Adalı and Calhoun, 2007;Li et al., 2007;Xiong et al., 2012;Yu et al., 2015).Based on the DMN estimates at selected model orders for each subject, we calculated the group DMN components for 40 HCs and 42 SZs using the method given in Section 2.2.2.
We counted the number of activated voxels within several ROIs of group DMN components and calculated the difference between the HC and SZ groups.We selected the ROIs according to Brodmann areas (BA), which included the precuneus (BA 7), medial frontal gyrus (BAs 8, 9), ventromedial prefrontal cortex (BA 10), orbitofrontal areas (BA 11), lingual gyrus (BAs 17,18,19), and anterior cingulate cortex (BA 32).To investigate group differences, we determined the number of activated voxels within ROIs for each subject and performed a two-sample t-test for 40 HCs and 42 SZs.
Next, we aligned single-subject DMN estimates by masking them with a group DMN mask for each model order.In the group DMN mask, a voxel = 0 if its value in the group component over 82 subjects was smaller than 0.5; otherwise, a voxel = 1.Based on these aligned single-subject DMN estimates, we computed the difference t-map and difference in variance maps between HCs and SZs to assess activation differences and subject variability.We formed the difference t-map by performing voxel-wise two-sample

Best run selection
We first present an example of the spatial DMN references generated using the proposed method, which combines across-subject averaging and a one-sample t-test (Z > 0.5; t > 3.250; p < 0.01; df = 9) with a comparison of references constructed by an across-subject one-sample t-test (i.e., t-map; t > 3.250; p < 0.01; df = 9) and averaging (Z > 0.5), respectively.Fig. 1 shows results for ICA of complex-valued and We further compared our proposed method, MST-based method, and ICASSO for magnitude-only analysis and complex-valued analysis with varying model orders (from 10 to 140 with 10 intervals).The quantitative indexes were means and standard deviations of Pearson correlation coefficients between the estimated DMN-, visual-, and sensorimotor-related components selected from the best run and their corresponding spatial references across all 82 subjects.Fig. 2 shows results of magnitude-only analysis and complex-valued analysis.We observe from Fig. 2 that the proposed method achieved higher means and lower standard deviations of correlation coefficients for all DMN-, visual-, and sensorimotor-related components as compared with the MST-based method and ICASSO.This result suggests that our best run selection method benefited from both subject averaging and the one-sample t-test.In addition, we found that component splitting resulted in decreasing means of correlation coefficients in magnitude-only analysis when we increased model orders from 60 to 140 (see magenta lines in Fig. 2).Conversely, means of correlation coefficients in complex-valued analysis showed increasing trends (see green lines in Fig. 2), which suggests intact DMN-, visual-, and sensorimotor-related components (i.e., close to their corresponding reference networks) at higher model orders in complex-valued analysis.

ICA analyses with varying model orders
We obtained spatial group components using the method described in Section 2.2.2 for complex-valued fMRI data, magnitude-only fMRI data, and phase fMRI data.Our first observation from Figs. 3-5 is that complex-valued analysis required a higher model order to extract meaningful components than magnitude-only and phase analyses.Specifically, magnitude-only analysis separated DMN and visual components at model order N = 10 and sensorimotor components at N = 20, whereas complex-valued analysis could not extract DMN and visual components until N ≥ 20 and could not extract sensorimotor components until N ≥ 30.This result is consistent with previous findings that the model order of complex-valued fMRI data is higher than that of magnitude-only fMRI data (Li et al., 2007;Xiong et al., 2012;Kuang et al., 2017a) due to including additional phase fMRI data.
Our second observation is that, with increasing model order, magnitude-only analysis showed component splitting, whereas phase analysis and complex-valued analysis showed component integration.LIPL at N = 50, and RIPL and IPL at N = 60.The model orders for detecting these subcomponents in complex-valued analysis were either identical to or higher than those of magnitude-only or phase analyses.
For example, the model order was the same for detecting DMNP1 (N = 20) in all three analyses, but it was higher for extracting DMNA than in phase analysis (complex: N = 30; phase: N = 10).At each model order (e.g., N = 140), essential subcomponents to form the intact complex-valued DMN existed in the magnitude-only (e.g., DMNP2) and phase (e.g., "DMNP2+DMNA" and IPL) analyses.
We also observed splitting of magnitude-only components and integration of phase components for visual networks as shown in Fig. 4 and sensorimotor components as shown in    Conversely, as shown in Fig. 6(1), in complex-valued analysis, the number of activated voxels in the three intact group components showed an initial increase and then maintained a nearly constant value as model order increased.This quantified the existence of the intact components in complex-valued analysis at higher model orders.A similar tendency was observed in the phase analysis shown in Fig. 6(3).The representative intact components showed a roughly constant voxel number, which indicated that phase analysis also detected intact components.Regarding the other subcomponents in complex-valued and phase analyses, the number of voxels showed a slight increase as model order increased (except for DMNA and LPMA in complex-valued analysis).

Stability of ICA analysis
Fig. 7 shows ICA stability for complex-valued, magnitude-only, and phase analyses with increasing model order.We compared the mean cluster quality index  ̅ and the standard deviation  ̅ averaged across 82 subjects (threshold at  ℎ = 0.7).We found that  ̅ values kept decreasing in magnitude-only analysis as the full range of model orders increased from 10 to 140.In contrast,  ̅ values remained stable at higher model orders of N ≥ 60 after initially decreasing at lower model orders of N < 60 in complex-valued and phase analyses.Note that  ̅ values of complex-valued analysis were lower than  -13 ) analyses.This lower value may have been because stability analysis used the unmixing matrix of ICA, which corresponded to the initial component estimates without any post-ICA processing (such as post-ICA denoising for complex-valued fMRI analysis).
Regarding the standard deviation  ̅, the three analyses first demonstrated similar tendencies in that  ̅ values were ascending at lower model orders of N < 60.Magnitude-only analysis then exhibited slightly increasing  ̅ values, whereas complex-valued and phase analyses showed steady  ̅ values at higher model orders of N ≥ 60.This was consistent with previous findings that magnitude-only analysis demonstrated decreasing stability with increasing model order (Abou-Elseoud et al., 2010).In contrast, complex-valued analysis yielded a relatively stable separation at higher model orders.

Number of activated voxels within ROIs
As stated in Section 2.2.3, we used the intact group DMN to show its capacity for distinguishing HCs and SZs.We calculated the number of activated voxels within several ROIs of group DMN and tested for group-level HC-SZ differences using a two-sample t-test for 40 HCs and 42 SZs (Table 1).In general, compared with SZs, HCs had a higher number of activated voxels in the medial frontal gyrus (BAs 8, 9), ventromedial prefrontal cortex (BA 10), orbitofrontal areas (BA 11), and anterior cingulate cortex (BA 32).In contrast, compared with HCs, SZs had a higher number of activated voxels in the precuneus (BA 7) and lingual gyrus (BAs 17,18,19).
For a higher model order of N = 120, at which complex-valued analysis had an intact DMN but magnitude-only analysis had only a DMNP2 subcomponent, the complex-valued intact DMN not only

Table 1
Figure 7 A showed much larger HC and SZ differences in activated voxel numbers within all ROIs, but it also showed significant subject-level differences in that all p-values were less than 0.05.In contrast, the magnitude-only DMNP2 component showed much smaller voxel number differences for HCs and SZs (BAs 17-19: a larger voxel number difference but insignificant subject-level difference).
For lower model orders (complex: N = 40; magnitude: N = 20), it was evident that the complex-valued intact group DMN was comparable to the magnitude-only intact group DMN in finding differences between HCs and SZs and presented larger differences in the number of activated voxels and extremely significant subject-level differences for the ventromedial prefrontal cortex (BA 10; HC-SZ = 94; p < 1.6 × 10 -7 ) and orbitofrontal areas (BA 11; HC-SZ = 136; p < 5.2 × 10 -13 ).These two ROIs demonstrated similar performance at a higher model order of N = 120 for complex-valued analysis.Specifically, the ventromedial prefrontal cortex (BA 10) had HC-SZ = 212 (p < 2.9 × 10 -9 ), and the orbitofrontal areas (BA 11) had HC-SZ = 167 (p < 2.7 × 10 -21 ).In fact, ROIs BAs 8-11 and BA 32 were concentrated in anterior DMN areas.These results are supported by previous findings.For example, HCs had stronger effects in the bilateral medial frontal cortex (BAs 8, 9) and left anterior cingulate cortex (BA 32) compared with SZs (Mingoia et al., 2012), and SZs exhibited moderate increases in the posterior cingulate and parietal lobe (Calhoun et al., 2008).By using complex-valued fMRI data, BA 10 and BA 11 exhibited highly significant subject-level differences for HCs and SZs.

T-maps of spatial group differences
The difference DMN t-map was formed by performing voxel-wise two-sample t-tests (p < 0.05) on 40 aligned HC components and 42 aligned SZ components (see Section 2.2.3), and results are shown in Fig. 8(3).The group DMN spatial components of HCs and SZs are shown in Fig. 8(1)-8(2).We first compared the difference DMN t-maps for complex-valued and magnitude-only analyses at a higher model order of N = 120, as shown in Fig. 8A(3)-8B(3).It is evident that the difference DMN t-map in complex-valued analysis included similar changes as those observed in the magnitude-only difference t-map.That is, the posterior areas in SZs were significantly higher than those areas in HCs.Furthermore, the difference DMN t-map in complex-valued analysis demonstrated unique and additional differences in anterior areas where HCs were significantly higher than SZs.This can be readily seen in the group DMN components of HCs and SZs at N = 120, as shown in Fig. 8A(1)-8A(2) and Fig. 8B(1)-8B(2).The group DMN of HCs showed much larger activations in anterior areas when compared with SZs in complex-valued analysis, whereas anterior areas showed subtle activations and differences in the group DMNs of HCs and SZs in magnitude-only analysis.
Next, we compared the difference DMN t-maps for complex-valued and magnitude-only analyses at lower model orders (complex: N = 40; magnitude: N = 20).We observed that the difference t-map in the complex-valued analysis showed more differences in the anterior areas than the magnitude-only analysis, which was consistent with the results at higher model orders.The difference found in posterior areas (SZs were significantly higher than HCs) at higher model orders was not observed at lower model orders, as shown in Fig. 8C(3)-8D(3).The results at lower model orders are reasonable in that ICA can extract more activations from complex-valued fMRI data than magnitude-only fMRI data because of the inclusion of phase fMRI data (Rodriguez et al., 2011(Rodriguez et al., , 2012;;Yu et al., 2015).

Difference variance maps
Fig. 9 shows variance maps of HCs, variance maps of SZs, and difference variance maps for the four cases.As shown in Fig. 9(3), we computed the difference in the DMN variance maps by calculating voxel-wise differences between the variance maps of HCs and SZs.We generated the variance maps of HCs and SZs shown in Fig. 9(1)-9(2) by computing voxel-wise variance values across DMN components of HCs and SZs.Briefly, the difference variance maps were similar for complex-valued and magnitude-only analyses at lower model orders.For the higher model order (N = 120), however, complex-valued analysis achieved better results in that SZs showed higher variance differences in posterior areas of the DMN than HCs, and HCs showed higher variance differences in anterior areas of the DMN than SZs.This finding is consistent with previous results of the difference DMN t-map.In sum, the intact DMNs obtained at the higher model orders in complex-valued analysis appeared to be useful for differentiating HCs and SZs.

Discussion
This study examined model order effects on ICA of complex-valued fMRI data based on resting-state complex-valued fMRI data sets from 82 subjects (40 HCs and 42 SZs).We analyzed DMN-, visual-, and sensorimotor-related networks as components of interest.We proposed a best run selection that combines across-subject averaging and a one-sample t-test to select components of interest for varying model orders from 10 to 140.Magnitude-only analysis exhibited component splitting at higher model orders (N ≥ 70), whereas complex-valued and phase analyses showed component integration with increasing model order.
Complex-valued analysis detected an intact component at all model orders and extracted some subcomponents at a higher model order compared with magnitude-only or phase analyses.The complex-valued intact DMN component obtained at higher model orders exhibited highly significant group differences for HCs and SZs.
We could infer the difference between model order effects on complex-valued fMRI data and magnitude-only fMRI data from incorporation of phase data.Results of phase analysis together with magnitude-only analysis provided quantitative supports for splitting differences between complex-valued and magnitude-only data at higher model orders.The white boxes at the bottom of Figs.3-5 show one example of the complementary role of magnitude-only and phase data in generating an intact component.
In fact, this can be found at any model order where all essential subcomponents for forming the complex-valued intact component can be found in either magnitude-only or phase analysis.For example, when DMNA was invisible in DMNP2 at N ≥ 50 in magnitude-only analysis, a DMNA component was extracted (from N = 10) in phase analysis (Fig. 3).When OPVS was split out of the intact visual component in magnitude-only analysis at N = 30, phase analysis detected OPVS at N = 30 (Fig. 4).When RPMA and LPMA were invisible in SMA2 in magnitude-only analysis at N ≥ 60, phase analysis extracted RPMA and LPMA from N = 20 (Fig. 5).As a result, the existence of an intact component is reasonable at higher model orders.Along this line, the coexistence of IPL, PIPL, and LIPL shown in Fig.
3A was also acceptable.Considering that we always detected an intact component at varying model order, the subcomponents in the complex-valued analysis may have been generated by magnitude-only or phase analysis instead of by component splitting.The results given in Section 3.2 support the finding that the model order for detecting subcomponents in complex-valued analysis was either identical to or higher than magnitude-only or phase analysis.
From a different point of view, we can say that complex-valued analysis can detect more voxels than magnitude-only analysis, and these additional voxels are meaningful as they have increasing correlations with the reference networks in Smith et al. (2009).This can be derived from the results shown in Figs.6and 2. In comparison with magnitude-only analysis showing component splitting with increasing model order, complex-valued analysis detected an increasing number of voxels for the intact DNM, visual, and sensorimotor components (Fig. 6).The correlation coefficients between these intact complex-valued components and their reference networks also increased with increasing model order, similar to those shown in Fig. 2. Therefore, compared with the magnitude components, the extra voxels in the intact complex-valued components have increasing correlations with the reference networks as model order increases.These additional and meaningful voxels detected by complex-valued analysis at higher model orders can provide additional brain information beyond magnitude-only analysis, which may be useful for applications such as differentiating HCs and SZs, as demonstrated in Section 3.5.
The uniqueness of phase information has been verified in multiple publications (Hoogenraad et al., 1998(Hoogenraad et al., , 2001;;Calhoun et al., 2002;Zhao et al., 2007;Feng et al., 2009;Arja et al., 2010;Castro et al., 2014;Yu et al., 2015;Chen and Calhoun, 2016;Chen et al., 2017b).Nevertheless, this study elucidated additional details about model order effects on ICA of phase data.First, component integration in phase analysis is the opposite of component splitting in magnitude-only analysis.We compared eigenvalues (sorted in a descending order) of the covariance matrix of phase data with those of magnitude-only data.Fig. 10 shows the normalized and averaged results over all 82 subjects.It is evident that the eigenvalues of the phase data decreased much more slowly than those of magnitude-only data for all subjects.Note the eigenvalues of phase data were also much smaller than those of magnitude-only data (maxima: 4.29 vs. 3.94 × 10 9 ).Therefore, it is reasonable to infer that, with increasing model order, increasing phase information is retained, which compensates for the loss of magnitude information caused by splitting and enables generation of integrated components.The relationship of the number of activated voxels between magnitude-only analysis and phase analysis, as shown in Fig. 6, supports this point.Second, we separated different components from phase fMRI data compared with those from magnitude-only data (e.g., "DMNP2+DMNA" in Fig. 3C, visual component in Fig. 4C, and sensorimotor component in Fig. 5C).
These components in phase analysis were nearly intact at varying model orders (see the first column in Figs.3C-5C).Moreover, we extracted the same components for magnitude-only and phase data at different model orders (e.g., DMNA emerged at N = 30 in magnitude-only analysis, but at N = 10 in phase analysis).In this case, it is safe to conclude that phase data provide unique and complementary brain activations in complex-valued fMRI analysis.Thus, complex-valued fMRI data require higher model orders than magnitude-only fMRI data for ICA to extract meaningful components.
The DMN component has been widely implicated in studies of many diseases including Alzheimer's disease, schizophrenia, and epilepsy (Broyd et al., 2009).Multiple studies have found significant abnormalities in activity and deactivation within the DMN (Garrity et al., 2007;Calhoun et al., 2008Calhoun et al., , 2009;;Pomarol-Clotet et al., 2008;Whitfield-Gabrieli et al., 2009;Jeong et al., 2010;Karbasforoushan and Woodward, 2012;Mingoia et al., 2012;Manoliu et al., 2014).Spatial differences in the DMN have been widely found in schizophrenia based on both task-related and resting-state fMRI data (Garrity et al., 2007;Calhoun et al., 2008Calhoun et al., , 2009;;Sambataro et al., 2010;Karbasforoushan and Woodward, 2012;Mingoia et al., 2012;Manoliu et al., 2014;Gopal et al., 2016;Hu et al., 2017).In this study, we examined HC-SZ group differences in terms of the number of activated voxels, spatial t-maps, and variance maps based on the complex-valued intact DMN.As for the number of activated voxels within ROIs of the DMN, we obtained results similar to those of Mingoia et al. (2012), who showed that HCs had stronger effects in the bilateral medial frontal cortex (BAs 8-9) and left anterior cingulate cortex (BA 32) compared with SZs.Our difference DMN t-maps in magnitude-only analysis were consistent with previous studies by Manoliu et al. (2014) and Calhoun et al. (2008).Manoliu et al. reported that the DMN at a higher model order showed significantly decreased activity in the bilateral precuneus of DMNP for SZs.Calhoun et al. detected more activations in anterior areas of HCs when compared with SZs, whereas SZs presented more activity in the posterior cingulate and parietal regions in a group average DMN component at a lower model order during an auditory oddball task (Calhoun et al., 2008).Other studies also found decreased task-related suppression in the medial and dorsolateral prefrontal cortexes (Pomarol-Clotet et al., 2008;Whitfield-Gabrieli et al., 2009;Jeong et al., 2010) in schizophrenia.Note that Garrity et al. (2007) and Harrison et al. (2007) both detected abnormally increased deactivation in the prefrontal cortex and anterior cingulate cortex for SZs.For difference variance maps, our results are supported by those reported in Gopal et al. (2016).Gopal et al. captured greater variance for HCs than for SZs in a few anterior areas of the DMN at a high model order (N = 75), although they detected less significant activity in the DMN.Moreover, our difference variance maps in complex-valued analysis additionally captured clearly higher variance for SZs than for HCs in posterior areas.Compared with magnitude-only analysis, the ventromedial prefrontal cortex (BA 10) and orbitofrontal areas (BA 11) exhibited highly significant subject-level differences for HCs and SZs.SZs showed higher subject-level spatial variability in posterior areas of the DMN, whereas HCs showed higher spatial variability in anterior areas of the DMN.In sum, our complex-valued results at a higher model order (N = 120) not only support previous results, but provide additional information that might serve as biomarkers for distinguishing SZs from HCs.This suggests that additional unique biological information can be extracted from complex-valued fMRI data to understand better the disease and distinguish HCs from SZs.
We note that some subcomponents for the visual component (Gopal et al., 2016;Van de Ven et al., 2017) and temporal lobes (Calhoun et al., 2004(Calhoun et al., , 2008;;Gopal et al., 2016) also show spatial differences between HCs and SZs.We found that the intact temporal lobe component had similar benefits for distinguishing SZs from HCs, but the intact visual and sensorimotor components did not show significant differences.
Additional quantitative indices may highlight more differences in their complex-valued form and this is a topic of future work.Future study of functional connectivity among brain networks at higher model orders for complex-valued fMRI data is also promising.Finally, we used across-subject averaging to generate group components.Group ICA and independent vector analysis (Garrity et al., 2007;Calhoun et al., 2008Calhoun et al., , 2009;;Calhoun and Adalı, 2012;Gopal et al., 2016;Chen et al., 2017a;Kuang et al., 2017c) are good candidates for providing group components with new characteristics to investigate model order effects.

Phase correction and denoising methods
Given a complex-valued spatial component    ∈ ℂ  ( = 1,2,3;  = 1, ⋯ , ) , we assume its corresponding time course is  ̃  = { ̃  ()} ∈ ℂ  ( = 1, ⋯ , , and T is the number of time points), which is obtained by decompressing the vector    ∈ ℂ  of the estimated mixing matrix .Phase correction and denoising are done as follows (Yu et al., 2015): Step 1: Construct  ̃  ∈ ℝ 2× using the real part  ̃, Step 2: Calculate the covariance of  ̃  : Step 3: Perform eigenvalue singular decomposition on cov( ̃  ) =   , where  is the eigenvalues matrix, and  is the eigenvector matrix: Calculate the rotation angle    from the elements of .
Step 4: Rotate    as follows: =        . (A.4) Step 5: Remove sign ambiguity of    according to the correlation between    and its reference  ,ref :   () denotes the set of indices belonging to the nth cluster;  − () is the set of indices that does not belong to the nth cluster, and 〈•〉 is size.The cluster quality index is  , ∈ [0,1].Thus, cluster quality improves as  , increases.
Let  , ′ ∈ { , >  ℎ } denote the threshold at  , ;  ℎ is the threshold, and ′ is the cluster number for  , ′ .The mean and standard deviation of  , ′ for subject k are calculated as follows (Kuang et al., 2017a): Finally,  ̅  and σ  are each averaged across 82 subjects to assess overall stability of each ICA analysis (Kuang et al., 2017a): Tables Table 1.HC and SZ differences in activated voxel number within several ROIs of group DMN at lower model orders (complex: N = 40; magnitude: N = 20) and a higher model order (N = 120).At N = 120, the subcomponent DMNP2 in magnitude-only analysis was analyzed.2) averaging (Z > 0.5), and (3) the proposed method (Z > 0.5; t > 3.250; p < 0.01) for ICA of complex-valued (A) and magnitude-only (B) fMRI data from subject 2 at model order N = 60.The magnitude parts of complex-valued spatial maps are shown.
Figure 2. Comparison of our proposed method, the MST-based method, and ICASSO for magnitude-only analysis (magenta) and complex-valued analysis (green) with varying model order (from 10 to 140 with 10 intervals) in terms of correlation coefficients (CC) between component estimates selected from the best run and their corresponding reference networks in Smith et al. (2009).We calculated means and standard deviations of correlation coefficients across all 82 subjects.We analyzed DMN (A), visual (B), and sensorimotor (C) components.Sensorimotor was an integrated component containing RPMA and LPMA regions.This integration is illustrated in a blue box.In complex-valued analysis, an intact sensorimotor component was detected at N = 30.LPMA was detected at N = 40, and SMA1 appeared at N = 50.RPMA emerged at N = 70.Finally, the bottom white box highlights one example of a complementary form of magnitude-only and phase data in complex-valued analysis.The intact sensorimotor component was integrated by magnitude SMA2 and phase sensorimotor.and SZs (2) are shown in addition to the difference DMN t-maps (3) (two-sample t-test, p < 0.05).A higher model order of N = 120 was used for both complex-valued analysis (A) and magnitude-only analysis (B).A lower model order of N = 40 was used for complex-valued analysis (C), and N = 20 was used for magnitude-only analysis (D).The magnitude parts of complex-valued group components and difference t-maps are displayed.


Phase data has a complementary role in preserving the integrity of brain networks. Intact group DMN at higher orders showed significant difference between HCs and SZs.We proposed an improved best run selection to select a run of ICA results.
-state complex-valued fMRI data were collected from 40 HCs and 42 SZs with written subject consent overseen by the University of New Mexico Institutional Review Board.During the scan, all participants were instructed to rest quietly in the scanner, keep their eyes open without sleeping, and not think of anything in particular.fMRI scans were acquired using a 3.0 Tesla Siemens Allegra scanner (Siemens Medical Solutions USA, Inc., Malvern, PA, US) equipped with 40 mT/m gradients and a standard quadrature head coil.The functional scan was acquired using gradient-echo echo-planar imaging with the following parameters: TR = 2 s, TE = 29 ms, field of view = 24 cm, acquisition matrix = 64 × 64, flip angle = 75°, slice thickness = 3.5 mm, slice gap = 1 mm.Data preprocessing was performed using the Statistical Parametric Mapping software package (Wellcome Trust Center for Neuroimaging, UniversityCollege of London, London, England; http://www.fil.ion.ucl.ac.uk/spm).Five scans were excluded because of steady-state magnetization effects.Functional images were motion corrected and then spatially standard Montreal Neurological Institute space.Following spatial normalization, the data were slightly subsampled to 3 × 3 × 3 mm³, which resulted in 53 × 63 × 46 voxels.Both magnitude and phase images were spatially smoothed with an 8 × 8 × 8 mm³ full-width half-maximum (FWHM) Gaussian kernel.Phase images were motion corrected using the transformations computed from magnitude-only data.Complex division of phase data by the first time point reduced the need for phase unwrapping, and spatial normalization of phase images used the warp parameters computed from magnitude-only data.A total of 146 scans were entered into the ICA.
t-tests (p < 0.05) on 40 aligned HC components and 42 aligned SZ components.The difference variance map was computed by calculating the voxel-wise difference values between the variance maps of HCs and SZs (HC-SZ).We generated the variance map of HCs (or SZs) by computing the variance values at each voxel across 40 aligned HC components (or 42 aligned SZ components).
fMRI data from subject 2 at model order N = 60.The spatial references generated by the proposed method [Fig.1A(3) and 1B(3)] included fewer voxels than those generated by the one-sample t-test [Fig.1A(1) and 1B(1)] and averaging [Fig.1A(2) and 1B(2)].We retained voxels with high mean Z-values and high t-values in the proposed spatial DMN reference.

Fig. 5 .
Three splitting cases are marked in Fig. 4B.The intact visual component was split into OPVS and an "MVS1+LVS" component, which included medial and lateral areas at N = 30."MVS1+LVS" further branched into MVS1 (a medial visual component) and LVS at N = 60, and LVS segmented into LLVS at N = 100 and RLVS at N = 110.One integrating case marked in the phase analysis demonstrated that an integrated LVS component containing both RLVS and LLVS areas was detected at N = 110 when LLVS emerged at N = 110 (RLVS appeared at N = 70), as shown in Fig. 4C.In the complex-valued analysis (Fig. 4A), we detected an intact visual component at N ≥ 20 and subcomponents OPVS and MVS2 (another medial visual component) at N = 30 as well as LVS at N = 70.We also detected these subcomponents at the same or lower model order in magnitude-only analysis (OPVS at N = 30; LVS at N = 60) and phase analysis (OPVS at N = 30; MVS2 at N = 10).At N = 140, the intact visual component demonstrated an example integration by magnitude MVS1, phase visual, phase LVS, and phase OPVS. Figure5

Figure 3 .
Figure3.DMN-related group components showing the effect of increasing model order for complex-valued analysis (A), magnitude-only analysis (B), and phase analysis (C).Magnitude-only analysis shows component splitting, whereas phase and complex-valued analyses show component integration.In magnitude-only analysis, an intact DMN was detected at N = 10, and DMNP1 appeared at N = 20.The intact DMN first branched into DMNP2, DMNA, and IPL at N = 30 (DMNA region was visible in DMNP2) and then branched into DMNP2, DMNA, RIPL, and LIPL at N = 50 (DMNA region was invisible in DMNP2), whereas IPL split into RIPL and LIPL at N = 50.Regions of bilateral IPL within DMNP2 gradually became invisible from N = 30 to N = 140.Splitting is illustrated in orange boxes.For phase analysis, "DMNP2+DMNA" and DMNA were detected at N = 10.DMNP1 appeared at N = 20, and LIPL emerged at N = 40.RIPL emerged at N = 70, and an integrated component IPL containing both LIPL and RIPL areas was detected.This integration is shown in the blue box.In complex-valued analysis, intact DMN and DMNP1 were detected at N = 20.DMNA was extracted at N = 30, and LIPL appeared at N = 50.RIPL emerged at N = 60, and an integrated component IPL containing LIPL and RIPL areas was detected.This integration also is shown in the blue box.Finally, the bottom white box highlights one example of a complementary form of magnitude-only and phase data in complex-valued analysis.The intact DMN was integrated by magnitude DMNP2, phase "DMNP2+DMNA", and phase IPL.

Figure 4 .
Figure 4. Visual-related group components showing the effect of increasing model order for complex-valued analysis (A), magnitude-only analysis (B), and phase analysis (C).Magnitude-only analysis shows component splitting, whereas phase analysis shows component integration.In magnitude-only analysis, an intact visual network was detected at N = 10.At N = 30, the intact visual network branched into "MVS1+LVS" and OPVS.MVS2 emerged at N = 40."MVS1+LVS" branched into MVS1 and LVS at N = 60.LVS further branched into LLVS at N = 100 and RLVS at N = 110.The three orange boxes show this splitting.For phase analysis, visual and MVS2 were detected at N = 10.OPVS appeared at N = 30, and RLVS emerged at N = 70.LLVS emerged at N = 110, and an integrated component LVS containing both RLVS and LLVS areas was detected.This integration is shown in the blue box.In complex-valued analysis, an intact visual network was detected at N = 20.OPVS and MVS2 A C C E P T E D M A N U S C R I P T

Figure 5 .
Figure5.Sensorimotor-related group components showing the effect of increasing model order for complex-valued analysis (A), magnitude-only analysis (B), and phase analysis (C).Magnitude-only analysis shows component splitting, whereas phase analysis shows component integration.In magnitude-only analysis, an intact sensorimotor component was detected at N = 20.The intact sensorimotor component first branched into SMA2, LPMA, and SMA1 at N = 40 and then branched into SMA2, RPMA, LPMA, and SMA1 at N = 50.Splitting is shown in an orange box.For phase analysis, RPMA and LPMA were detected at N = 20, and sensorimotor and SMA1 emerged at N = 50.Sensorimotor was an integrated component containing RPMA and LPMA regions.This integration is illustrated in a blue box.In complex-valued analysis, an intact sensorimotor component was detected at N = 30.LPMA was detected at N = 40, and SMA1 appeared at N = 50.RPMA emerged at N = 70.Finally, the bottom white box highlights one example of a complementary form of magnitude-only and phase data in complex-valued analysis.The intact sensorimotor component was integrated by magnitude SMA2 and phase sensorimotor.

Figure 7 .
Figure 7.Comparison of ICA stability for complex-valued, magnitude-only, and phase analyses with increasing model order from 10 to 140 in terms of the mean cluster quality index  ̅ and the standard deviation σ ̅ averaged over 82 subjects.

Figure 8 .
Figure8.Difference DMN t-maps between HCs and SZs.Group DMN spatial components of HCs (1) and SZs (2) are shown in addition to the difference DMN t-maps (3) (two-sample t-test, p < 0.05).A higher model order of N = 120 was used for both complex-valued analysis (A) and magnitude-only analysis (B).A lower model order of N = 40 was used for complex-valued analysis (C), and N = 20 was used for magnitude-only analysis (D).The magnitude parts of complex-valued group components and difference t-maps are displayed.

Figure 9 .
Figure 9. Results of difference variance maps between HCs and SZs.Variance maps of HCs, variance maps of SZs, and difference variance maps of HCs and SZs (HC-SZ) (threshold at 0.5) are shown in subfigures (1)-(3), respectively.Complex-valued analysis at N = 120 (A) and N = 40 (C) as well as magnitude-only analysis at N = 120 (B) and N = 20 (D) are presented.The magnitude parts of complex-valued variance maps and difference variance maps are displayed.

Figure 10 .Figure 1 Figure 1 .Figure 2 Figure 2 .Figure 3 Figure 3 .
Figure 10.Normalized and averaged eigenvalues (in descending order) of the covariance matrix of magnitude-only data and phase data over all 82 subjects.A C C E P T E D M A N U S C R I P T