Relationship between Dynamic Blood-Oxygen-Level-Dependent Activity and Functional Network Connectivity: Characterization of Schizophrenia Subgroups

Abstract

Aim:

In this work, we propose the novel use of adaptively constrained independent vector analysis (acIVA) to effectively capture the temporal and spatial properties of dynamic blood-oxygen-level-dependent (BOLD) activity (dBA), and we efficiently quantify the spatial property of dBA (sdBA). We also propose to incorporate dBA into the study of brain dynamics to gain insight into activity-connectivity co-evolution patterns.

Introduction:

Studies of the dynamics of the human brain using functional magnetic resonance imaging (fMRI) have enabled the identification of unique functional network connectivity (FNC) states and provided new insights into mental disorders. There is evidence showing that both BOLD activity, which is captured by fMRI, and FNC are related to mental and cognitive processes. However, a few studies have evaluated the inter-relationships of these two domains of function. Moreover, the identification of subgroups of schizophrenia has gained significant clinical importance due to a need to study the heterogeneity of schizophrenia.

Methods:

We design a simulation study to verify the effectiveness of acIVA and apply acIVA to the dynamic study of resting-state fMRI data collected from individuals with schizophrenia and healthy controls (HCs) to investigate the relationship between dBA and dynamic FNC (dFNC).

Results:

The simulation study demonstrates that acIVA accurately captures the spatial variability and provides an efficient quantification of sdBA. The fMRI analysis yields synchronized sdBA-temporal property of dBA (tdBA) patterns and shows that the dBA and dFNC are significantly correlated in the spatial domain. Using these dynamic features, we identify schizophrenia subgroups with significant differences in terms of their clinical symptoms.

Conclusion:

We find that brain function is abnormally organized in schizophrenia compared with HCs since there are less synchronized sdBA-tdBA patterns in schizophrenia and schizophrenia prefers a component that merges multiple brain regions. Identification of schizophrenia subgroups using dynamic features inspires the use of neuroimaging in studying the heterogeneity of disorders.

Impact statement

This work introduces the use of joint blind source separation for the study of brain dynamics to enable efficient quantification of the spatial property of dynamic blood-oxygen-level-dependent (BOLD) activity to provide insight into the relationship of dynamic BOLD activity and dynamic functional network connectivity. The identification of subgroups of schizophrenia using dynamic features allows the study of heterogeneity of schizophrenia, emphasizing the importance of functional magnetic resonance imaging analysis in the study of brain activity and functional connectivity to gain a better understanding of the human brain, especially the brain with a mental disorder.

Introduction

Studies of the dynamics of human brain function using neuroimaging data have achieved great success in identifying unique biomedical patterns that enable a better understanding of the functional differences between the healthy and disordered brain (Assaf et al., 2010; Etkin et al., 2009; Hutchison et al., 2013; Lefebvre et al., 2016; Rashid et al., 2014; Uddin et al., 2011).

Functional magnetic resonance imaging (fMRI) is one of the most commonly used imaging modalities and measures the blood-oxygen-level-dependent (BOLD) activity in human brain during either resting- or task-evoked state (Allen et al., 2014; Kucyi and Davis, 2014; Sakoğlu et al., 2010; Verley et al., 2018; Wang et al., 2010). Most dynamic studies obtain a functional connectivity (FC) or functional network connectivity (FNC) matrix, which measures the association between the BOLD activity of different functional regions or coherent networks, as a function of time. A dynamic FNC (dFNC) analysis allows a systematical study of evolving functional patterns by jointly taking multiple brain networks into consideration (Calhoun et al., 2014; Mennigen et al., 2019).

There is rich work showing that both BOLD activity and F(N)C are related to mental and cognitive processes (Britz et al., 2009; Hutchison and Morton, 2016; Marusak et al., 2017; McIntosh et al., 2008; Pereira et al., 2019). Differences in BOLD activity and FNC have been separately reported in multiple mental disorders, especially in schizophrenia, such as reduced amplitude of low-frequency fluctuation (ALFF) in cuneus (Hoptman et al., 2010; Turner et al., 2013), reduced BOLD activation in anterior cingulate gyrus (Baiano et al., 2007; Schultz et al., 2012), dysconnectivity in default mode network (DMN; Mingoia et al., 2012; Van Den Heuvel and Pol, 2010), and dysconnectivity between thalamus and sensory regions (Calhoun et al., 2009; Guller et al., 2012; Kühn and Gallinat, 2013; Malaspina et al., 2004; Zhou et al., 2007). However, the association between time-varying, that is, dynamic BOLD activity (dBA) and dFNC is not well studied and it is desirable to incorporate dBA to gain insight into the activity-connectivity co-evolution by identifying highly correlated patterns between dBA and dFNC.

In Fu et al. (2018), the authors investigate the associations between dBA and dFNC in temporal domain, that is, they measure dBA and dFNC by using the temporal variabilities of functional networks and observe that dBA and dFNC are significantly correlated in time in some cases and patients with schizophrenia show lower or nonexistent associations between dBA and dFNC compared with the healthy controls (HCs). However, similar to Fu et al. (2018), most previous studies conduct a dFNC analysis by using the time courses of BOLD activity by assuming the spatial domain is static (Jafri et al., 2008; Calhoun et al., 2014; Allen et al., 2014; Weber et al., 2020).

Spatial variation of BOLD activity is observed as a change in the volume of a functional network or variations in the activated regions within a functional network and has started to attract attention, since they enrich the dynamic study of brain function (Iraji et al., 2019a,b, 2020). Previous studies have shown that simultaneously considering temporal and spatial changes yields more distinguishable resting-state networks (RSNs) between subject groups (Jie et al., 2018; Kottaram et al., 2018; Iraji et al., 2020). Studies that compute dFNC by using the spatial maps of brain networks have also emphasized the importance of the assumption of spatial variability in a dynamic study, by detecting significant differences between the patients with mental disorder such as between those with schizophrenia and the HCs (Bhinge et al., 2019a,b; Ma et al., 2014).

However, the spatial activation patterns of dBA themselves are not well explored primarily due to the lack of effective quantification strategies of spatial property of dBA (sdBA). An efficient quantification of sdBA does not only enable the investigation of the spatial activation patterns of dBA but also leads to a study of activity-connectivity co-evolution using sdBA and dFNC that is computed by using the spatial maps (sdFNC) of functional networks.

The extraction and efficient quantification of the dynamic features captured through temporal and spatial variation leads one to analyze their importance in studying a mental disorder such as schizophrenia. There has been significant interest in studying the subtypes of schizophrenia (Dwyer et al., 2018; Geisler et al., 2015; Jablensky, 2006, 2010) to gain a better understanding of the uncertainty in whether a precision medicine is needed (Senn, 2018) during clinical diagnosis and treatment. Subtypes of schizophrenia have been studied by using genetic information (Hallmayer et al., 2005; Morar et al., 2018; Sayın et al., 2013) but not yet using neuroimaging modalities such as fMRI data that have been successfully used in the study of schizophrenia (Calhoun and Adalı, 2009; Geisler et al., 2015; Ma et al., 2012).

Our previous work showed that fMRI alone can identify subgroups of schizophrenia that demonstrated significant differences in terms of clinical symptoms and it helps understand the heterogeneity of schizophrenia (Long et al., 2020). However, these subgroups were identified through static imaging features. This motivates an investigation of the effectiveness of dynamic neuroimaging features for identifying and studying the heterogeneity of mental disorders such as schizophrenia.

Both temporal and spatial variabilities of BOLD activity are important to gain a better understanding of brain dynamics. However, existing methods such as group independent component analysis (ICA; Calhoun et al., 2001) and joint ICA (Calhoun et al., 2006) can only effectively capture either temporal variability or spatial variability by making relatively strong assumptions. In this work, we propose the novel use of a recent method, adaptively constrained independent vector analysis (acIVA; Bhinge et al., 2019b), that enables us to simultaneously capture the dFNC and dBA in the spatial and temporal domain, and acIVA also provides an effective quantification of sdBA, as part of the algorithm.

AcIVA is able to precisely preserve spatial variability by adaptively tuning the constraint parameter that controls the spatial association between reference signals and the source estimates. More importantly, measuring the association between reference signals and the source estimates as part of the acIVA algorithm provides an efficient quantification of sdBA. The acIVA method also reduces the undesirable effects of high dimensionality that frequently arise during a dynamic study. We describe the issue of high dimensionality and its relation to the dynamic study in the Data Formation for Dynamic Analysis section. In addition, acIVA eliminates the tough alignment process across multiple decompositions due to the use of reference signals (Bhinge et al., 2019b).

To study the novel use of acIVA to capture the dynamic features, we propose a simulation study and demonstrate that acIVA yields a desirable performance in capturing spatial dynamics by accurately recovering the association between reference signals and the targeted estimates. The adaptively tuned constraint parameter ρ provides the closest lower bound for the association, hence preserving the spatial variation patterns across datasets.

We apply acIVA to extract dynamic features from a resting-state fMRI dataset acquired from 179 subjects (88 schizophrenia subjects and 91 HCs). We quantify the temporal property of dBA (tdBA) by using fractional ALFFs (fALFFs) of network time courses to conduct a comparison between tdBA and sdBA. To study the activity-connectivity co-evolution of dBA and dFNC, dFNC is measured by using the dependence among spatial maps and the correlation between sdBA and sdFNC is computed.

The analysis results show that tdBA and sdBA pairs of the same functional network have higher correlation than those of different networks, and sdBA yields higher correlation with the network connectivity that is associated with the same network. However, schizophrenia subjects show fewer significantly correlated sdBA and tdBA pairs and favor a complex component that merges multiple brain regions—super parietal, visual, and cerebellum (SP-V-C)—when compared with HCs, suggesting that the brain function is abnormally organized in schizophrenia; hence, more brain regions are simultaneously activated for a certain intrinsic function. Most importantly, we identify subgroups of schizophrenia by using the tdBA-sdBA association and the sdBA-sdFNC co-evolution, and we detect significant differences across subgroups in terms of their clinical symptoms that are measured by the positive and negative syndrome scale (PANSS) scores (Kay et al., 1987). This observation inspires further study of the heterogeneity of mental disorders by using neuroimaging modalities such as fMRI.

The rest of the article is organized as follows. The Materials and Methods section presents the framework of using acIVA for dynamic analysis of fMRI data, including data formation by adopting the sliding window approach, a detailed introduction of acIVA algorithm, the extraction of reference signal, and the quantification of dBA and dFNC. The Results section shows the simulation results and the application to real fMRI data. Finally, we conclude in the Discussion section and point out the limitations of this work and possible interesting future directions.

Materials and Methods

Institutional Review Board approval statement

Institutional Review Board approval was obtained for the study. All participants provided informed consent. This work is not a clinical trial.

Data acquisition and preprocessing

We used resting-state fMRI data from the Center of Biomedical Research Excellence (COBRE) that is available on the collaborative informatics and neuroimaging suite data exchange repository (Scott et al., 2011; Çetin et al., 2014; Aine et al., 2017). The data include 88 schizophrenia subjects (average age: $37 \pm 14$ ) and 91 HCs (average age: $38 \pm 12$ ). All schizophrenia subjects had a negative toxicology screen for drugs of abuse at the start of the study. The total olanzapine equivalent dose of these subjects ranges from 2 to 60.2 mg/day, with an average of 14.9 mg/day. The Clinical Core (COBRE Stability Clinic) affiliated with this project determined retrospective stability from relevant psychiatric records, documenting no change in symptomatology or type/dose of psychotropic medications that occurred during the 3 months before the referral.

For this study, the participants were asked to keep their eyes open during the entire scanning period. All images were collected on a single 3-Tesla Siemens Trio scanner with a 12-channel radio frequency coil by using the following parameters: echo time = 29 msec, repetition time = 2 sec, flip angle = 75°, slice thickness = 3.5 mm, slice gap = 1.05 mm, and voxel size $3.75 \times 3.75 \times 4.55$ mm³. Participants were instructed to keep their eyes open during the scan and stare passively at a central fixation cross. Each resting-state scan consists of 150 volumes. To eliminate the T1-related signal fluctuations (T1 effect; Shin et al., 2013), the first 6 volumes are removed in this study, thus 144 volumes remain for each subject. The fMRI data are realigned with the INRIalign algorithm (Freire et al., 2002) for head motion correction, followed by slice-timing correction to account for timing difference by using the middle slice as the reference frame. Then, the fMRI data are spatially normalized to the standard Montreal Neurologic Institute space (Friston et al., 1994) and resampled to $3 \times 3 \times 3$ mm³, resulting in $53 \times 63 \times 46$ voxels. Afterward, the fMRI data are smoothed by using a Gaussian kernel with a full-width at half-maximum of 5 mm.

Data formation for dynamic analysis

The sliding window approach facilitates the dynamic study of fMRI data that are acquired within a certain duration (Ma et al., 2014). In a sliding window approach, the entire scanning period is divided into overlapping windows of length T_w , yielding L windows for each subject, as shown in Figure 1. In our case, this results in $L \times 179$ datasets. Joint analysis of a large number of datasets typically involves estimation of a large number of parameters with a limited number of samples, which degrades the estimation of the spatiotemporal dynamic features. This motivates the use of an algorithm that effectively captures the variability in both the spatial and temporal domains from a large number of datasets.

FIG. 1.

Dynamic study using acIVA with the sliding window approach on a subset of K ₀ subjects. Applying acIVA to the windowed datasets with the spatial maps of common components as reference signals yields the time courses and spatial maps for each dataset. The constrained SCVs contain corresponding components from different windows and subjects. acIVA, adaptively constrained independent vector analysis; SCV, source component vector. Color images are available online.

Independent vector analysis (IVA) is a data-driven technique that extends ICA to multiple datasets and makes effective use of the dependence across datasets (Adalı et al., 2014). The IVA also relaxes the assumptions made in other joint blind source separation solutions such as group ICA (Calhoun et al., 2001), which assumes a common spatial signal space, and joint ICA (Calhoun et al., 2006), which assumes a common temporal signal space. The IVA simultaneously estimates the demixing matrices of all datasets to obtain dataset-specific time courses and spatial maps, effectively capturing both the temporal and spatial variability of functional networks across datasets (Bhinge et al., 2019a,b; Laney et al., 2014; Ma et al., 2014). However, this flexibility comes at a cost that, with a limited number of samples, the performance of IVA degrades as the number of datasets and/or the model order—the number of sources—increases (Bhinge et al., 2019b; Long et al., 2020). This issue is referred to as a curse of high dimensionality in IVA and is addressed by using a novel algorithm, acIVA, that incorporates reference information into the IVA decomposition.

Adaptively constrained IVA

The acIVA algorithm is a semi-blind source separation technique that incorporates reference information regarding the time courses or spatial maps into the IVA decomposition. This reference information guides the decomposition toward a desirable solution in high-dimensional scenarios, thereby addressing the issue of high dimensionality (Bhinge et al., 2019b). In this section, we introduce the general IVA model followed by a description of the acIVA technique.

Given K datasets each containing V samples, IVA assumes that each dataset is a linear mixture of N independent sources, $x^{[k]} (v) = A^{[k]} s^{[k]} (v), 1 \leq k \leq K, 1 \leq v \leq V,$ (1)

where $X^{[k]} = [x^{[k]} (1), x^{[k]} (2), \dots, x^{[k]} (V)] \in ℛ^{N \times V}$ , $S^{[k]} = [s^{[k]} (1), s^{[k]} (2), \dots, s^{[k]} (V)] \in ℛ^{N \times V}$ and $A^{[k]} \in ℛ^{N \times N}$ denote the observed dataset, the set of independent sources, and the invertible mixing matrix, respectively. In addition to the assumption of independence among sources within a dataset, IVA makes effective use of dependence across multiple datasets by defining a source component vector (SCV) as $s_{n} (v) = {[s_{n}^{[1]} (v), s_{n}^{[2]} (v), \dots, s_{n}^{[K]} (v)]}^{T} \in ℛ^{K \times 1}, 1 \leq n \leq N$ , by collecting corresponding components, where $s_{n}^{[k]} \in ℛ^{V \times 1}$ is the nth source from the kth dataset. The IVA finds K demixing matrices by minimizing the mutual information among the SCVs, which results in the following cost function

such that the estimated sources of each dataset are obtained as $y^{[k]} (v) = W^{[k]} x^{[k]} (v)$ for $k = 1, \dots, K$ , where denotes the demixing matrices, $y_{n}$ denotes the estimated SCV, $ℋ (\cdot)$ denotes the (differential) entropy, and $h_{n}^{[k]}$ is a unit vector resulting from the decoupling process that is perpendicular to all rows of $W^{[K]}$ except $w_{n}^{[k]}$ (Bhinge et al., 2017; Li and Zhang, 2007).

The acIVA algorithm guides the decomposition with prior information such as properly selected reference signals for the source components. The IVA decomposition is, hence, achieved by minimizing the cost function in Equation (2) subject to an inequality constraint $g_{n} (y_{n}^{[k]}, d_{n}) = ρ_{n} - |{(y_{n}^{[k]})}^{T} d_{n}| \leq 0$ , where $0 \leq ρ_{n} \leq 1$ is the constraint parameter that provides the lower bound for the similarity between the estimate $y_{n}^{[k]}$ and the reference signal $d_{n}$ that is measured by using Pearson correlation. The cost function of acIVA is defined by incorporating an inequality constraint in the IVA cost function, yielding

where is the IVA cost function as defined in Equation (2), $M (0 \leq M \leq N)$ is the number of source estimates to be constrained, $μ_{m}^{[k]}$ is the regularization parameter, and $γ_{m} > 0$ is the penalty parameter (Bhinge et al., 2019b).

Through an adaptive parameter-tuning process, acIVA yields different values of the lower bound $\hat{ρ}$ for the similarity, which allows the estimate to vary across datasets. Consequently, acIVA is able to effectively capture the variability across datasets that describes sdBA. Note that we refer to the voxel activation conveyed in the spatial maps of estimated components from acIVA as BOLD activity, which is different from the raw BOLD signal in the original fMRI data. The use of reference signals in acIVA also provides us with a useful metric for quantifying spatial variation. Therefore, in this work, we propose a novel use of acIVA to efficiently quantify the spatial variability by using the similarity between the estimates and corresponding reference signals.

The use of reference signals also effectively reduces the effect of high dimensionality, providing a more robust estimation. In addition, it eliminates the alignment problem across multiple decompositions. In application to the dynamic study of resting-state COBRE fMRI data, this makes it possible to divide the 179 subjects into multiple subsets of K ₀ subjects (yielding $K = L \times K_{0}$ datasets) and perform multiple individual IVA decompositions without considering the challenging issue of component alignment. We use an IVA algorithm that uses a multivariate Laplacian distribution model for the sources, such as IVA-L (Kim et al., 2006), and additionally considers the second-order statistics (IVA-L-SOS; Bhinge et al., 2019b), since the Laplacian distribution is a good match for fMRI data and the adding of SOS enables a full statistical characterization of a Laplacian multivariate random vector.

Reference signal extraction

In our implementation of acIVA for extracting and quantifying sdBA from multi-subject resting-state fMRI data, we use the spatial maps of exemplar RSN components as reference signals. Possible choices for the exemplar RSNs include the pre-defined RSN templates (Allen et al., 2011) or the group-level RSN components extracted from the same dataset by using group decomposition algorithms such as group ICA (Bhinge et al., 2019a).

In this work, we apply acIVA to the resting COBRE data collected from a schizophrenia group and an HC group, seeking to gain a better understanding of schizophrenia through a fair comparative study of dynamic features, that is, the association between sdBA and tdBA and the co-evolution of sdBA and sdFNC, between groups. Therefore, instead of using arbitrarily pre-defined or estimated RSNs, it is desirable to use common static RSNs that are shared across the subjects of two groups as the reference signals. We extract common RSNs from the data of all 179 subjects by using a common subspace analysis method (Long et al., 2020) and use their spatial maps as reference signals. All the spatial maps are normalized to have zero mean and unit variance.

Analysis of dBA and dFNC

Spontaneous slow fluctuations of correlated activity are a fundamental feature of the resting brain and can be captured as the BOLD signal that reflects neural synchrony between brain regions (Bluhm et al., 2007; van de Ven et al., 2004). A related measure is fALFF that quantifies the amplitude of these low-frequency oscillations. Therefore, we measure tdBA by computing the fALFF of time courses as the fraction of the square root of power spectrum integrated in a low-frequency band (0.025–0.15 Hz) that summarizes the most dominant frequencies of BOLD signal during resting state (Zang et al., 2007; Zou et al., 2008).

Before computing fALFF, the time courses were detrended, motion parameters were regressed (including their derivatives, their squares, and derivatives of their squares), and finally despiked, which involved detecting spikes as determined by the 3dDespike algorithm that is originally included in the Analysis of Functional NeuroImages toolbox (https://afni.nimh.nih.gov) and replacing spikes by values obtained from third-order spline fit to neighboring clean portions of the data.

The fALFF feature of each component is a vector of length L with each entry as the fALFF value of the associated time course in a single window. In this study, there are M RSN components involved, hence yielding M fALFF features. Since the variabilities of the spatial maps of RSNs reflect sdBA, we quantify sdBA by using correlation between the spatial maps of the estimates and the reference signals, $\tilde{ρ} = |{(y_{n}^{[k]})}^{T} d_{n}|$ . We, hence, introduce an effective metric to quantify the spatial variability of dBA, which enables an investigation of the relationship between the tdBA and sdBA for identifying synchronized patterns. We compute the cross-correlation between $\tilde{ρ}$ and fALFF to study the association between sdBA and tdBA and denote the $\tilde{ρ}$ -fALFF cross-correlation matrix as $C^{[k]}$ for each subject, as shown in Figure 2a.

FIG. 2.

Cross-correlation between sdBA ( $\tilde{ρ}$ ) and tdBA (fALFF) (a) and activity-connectivity co-evolution computed by using sdBA ( $\tilde{ρ}$ ) and sdFNC (b). dBA, dynamic blood-oxygen-level-dependent (BOLD) activity; fALFF, fractional amplitude of low-frequency fluctuation; RSN, resting-state network; sdBA, spatial property of dBA; sdFNC, dFNC that is computed by using the spatial maps; tdBA, temporal property of dBA. Color images are available online.

As there is evidence that both BOLD activity and network connectivity are related to mental and cognitive processes in the brain, it is desirable to study the co-evolution of BOLD activity and network connectivity. In Fu et al. (2018), the authors studied activity-connectivity co-evolution by using temporal characteristics of the BOLD signal. In this work, we make use of an effective metric $\tilde{ρ}$ for sdBA, enabling the study of activity-connectivity co-evolution in the spatial domain by using the spatial variability of the BOLD signal. We measure sdFNC by using the normalized mutual information among the spatial maps. Pearson's correlation between $\tilde{ρ}$ and sdFNC is calculated and the activity-connectivity co-evolution matrix is denoted as $E_{n}^{[k]}$ for each component and each subject, as shown in Figure 2b.

The matrices $C^{[k]}$ and $E_{n}^{[k]}$ are further analyzed to conduct a comparison between the schizophrenia and HC groups. We perform a one-sample t-test on the cross-correlation matrices $C^{[k]}$ across the subjects in the schizophrenia group and HC group separately, yielding a single matrix of t-statistics that are significant for each group. A one-sample t-test is also performed on the co-evolution matrices $E_{n}^{[k]}$ for two groups separately. Note that a co-evolution matrices $E_{n}^{[k]}$ is associated with a single component of a specific subject. The one-sample t-test on $E_{n}^{[k]}$ yields a matrix of t-statistics that is significant for each component in each group. We report the significant statistics with the false discovery rate (FDR) control at p < 0.05.

Study of schizophrenia heterogeneity

We study the heterogeneity of schizophrenia by identifying subgroups of schizophrenia subjects using k-means clustering to cluster the $\tilde{ρ}$ -fALFF cross-correlation matrices $C^{[k]}$ and the activity-connectivity co-evolution matrices $E_{n}^{[k]}$ separately. The unique entries of each matrix form a feature, with a length of 64 for $C^{[k]}$ and 28 for $E_{n}^{[k]}$ . We use the elbow criterion, which is calculated as the ratio of within-cluster distance to between-cluster distance (Allen et al., 2014; Zhang et al., 2018), to determine the number of subgroups $N_{s}$ .

After identifying subgroups of schizophrenia, we perform a multivariate analysis of variance (MANOVA) on $C^{[k]}$ and $E_{n}^{[k]}$ to detect significant differences across the subgroups. Four different test statistic values—Pillai's, Wilks', Hotelling's, and Roy's (Rencher, 2003)—are reported. A MANOVA provides a comprehensive comparison among multiple variables across multiple subgroups. However, it returns a single test statistic value such that we cannot check the correlation patterns. To visually compare the correlation patterns across subgroups, we perform a one-sample t-test on $C^{[k]}$ for each subgroup and a univariate analysis of variance (ANOVA) on $C^{[k]}$ across subgroups with an FDR control at $p < 0.05$ . We then plot the significant statistic values as a matrix to illustrate the correlation patterns for each subgroup.

Simulation study

In the application of acIVA to real fMRI data, we use the spatial maps of RSNs as reference signals to emphasize the importance of spatial variabilities. Therefore, we first use a simulation study to demonstrate that acIVA is able to accurately recover the spatial variabilities due to the adaptive parameter-tuning process.

The acIVA implementation enables the quantification of sdBA by adaptively tuning the amount of correspondence between the estimated functional networks and the reference signals. We use simulated data to demonstrate the ability of acIVA to accurately capture the underlying spatial variability of dBA. To simulate fMRI-like sources that are super-Gaussian distributed, we generate $N = 10$ SCVs that are multivariate generalized Gaussian distributed with randomly selected shape parameter $β \sim U (0.1, 0.8)$ that forms a super-Gaussian distribution, and different correlation matrices $R_{n} \in ℛ^{16 \times 16}, n = 1, \dots, 10,$ as described later in the four cases and shown in Figure 3.

FIG. 3.

Visualization of the correlation matrices and the decomposition results of acIVA in the four cases of simulation. A total of 100 independent realizations are generated for each case. In the ticklabels of x-axis in the boxplots, “H” refers to the results of components that are highly correlated (c ₁) with the reference signal, “L” refers to the results of components that are not highly correlated (c ₂), and “diff” refers to the cases where we summarize the differences between the tuned values of $\tilde{ρ}$ and the ground truth. The boxplot displays the median, the 25th and 75th percentiles of the values with whiskers extending to the 99.3% confidence interval and some outliers in red asterisks beyond whiskers. Color images are available online.

Out of the 10 SCVs, we constrain 6 SCVs, where the first component in each of the first 6 SCVs is used as the reference signal, and the other 15 components are used to generate $K = 15$ datasets. The sources for the kth dataset are obtained by concatenating the $(k + 1)$ th row from all SCVs. Each dataset, $X$ , is computed by mixing the 10 sources using a mixing matrix, whose elements are randomly drawn from a uniform distribution $U (0, 1)$ . We apply acIVA by using the IVA-L-SOS algorithm on the K datasets, with the first $M = 6$ components constrained by using the reference signals. The value of constraint parameter ρ in acIVA is tuned from a default set

The details of four cases where the SCVs are generated to have different correlation matrices are as follows, and the visualization of correlation matrices is illustrated in Figure 3. If the components within an SCV are highly correlated with the reference signal, it means the spatial activation of the component is relatively stationary and does not vary much across datasets. If the components are barely correlated with the reference signal, it means the component varies a lot across datasets in terms of its spatial activation.

Case 1: SCVs 1–3 are generated to have a uniform correlation structure with correlation value $c = 0.9$ , which means that all the components are highly correlated with the reference signal. SCVs 4–6 are generated to have a subset of components that are highly correlated with the reference signal and have high correlation values $c_{1} = 0.9$ . Those components that are not highly correlated with the reference signal have correlation values that are randomly selected from a uniform distribution $c_{2} \sim U (0.1, 0.2)$ . SCVs 7–10 are generated to have a correlation matrix with random entries, $C = {Q Q}^{T}, Q \sim U (- 0.2, 0.8),$ , that allow the components to have different correlation values with the reference signal.

Case 2: SCVs are generated by using the same parameters as in Case 1, except that the low correlation values in SCVs 4–6 are randomly selected from a uniform distribution $c_{2} \sim U (0.001, 0.2)$ . The default set is updated to be

Case 3: SCVs 1–3 are generated to have correlation values $c \sim U (0.8, 0.9)$ . The high correlation values in SCVs 4–6 are $c_{1} \sim U (0.8, 0.9)$ and the low correlation values are $c_{2} \sim U (0.1, 0.2)$ for SCV 4, $c_{2} \sim U (0.2, 0.3)$ for SCV 5, and $c_{2} \sim U (0.3, 0.4)$ for SCV 6. The correlation values are set to be in a specific range to investigate the sensitivity of the adaptive parameter-tuning process to varying values in acIVA. SCVs 7–10 are generated in the same way as in Cases 1 and 2.

Case 4: SCVs 1 and 2 are generated to have correlation values $c \sim U (0.8, 0.9)$ . SCVs 3 and 4 have correlation values $c \sim U (0.6, 0.7)$ . SCVs 5–10 are generated to have a correlation matrix $C = Q Q^{T}, Q \sim U (- 0.2, 0.8)$ . The components in SCVs 5 and 6 have random correlation with the reference signal.

We quantify sdBA by using correlation between the spatial maps of the estimates and the reference signals. This metric is reflected in the tuned value of $\hat{ρ}$ , which we use for measuring the performance of the algorithm. The ground truth of $ρ$ , the lower bound of the correlation, is 0.1 if $c \in [0.1, 0.2]$ , 0.2 if $c \in [0.2, 0.3]$ , and so on. In most cases, $ρ$ is a single value, hence we summarize the tuned values of $\hat{ρ}$ for a direct comparison between $\hat{ρ}$ and $ρ$ . However, the true $ρ$ for components with small correlation values in SCVs 4–6 of Case 2 and all the components in SCVs 5–6 of Case 4 is not a single value, because their correlation values have a large range. For example, the components in SCVs 5–6 of Case 4 have correlation values that change from 0 to 1, hence the true $ρ$ of some components with $c \in [0.1, 0.2]$ is 0.1 whereas for those with $c \in [0.2, 0.3]$ is 0.2. We then summarize the differences between $\hat{ρ}$ and $ρ$ to assess the performance of acIVA in these scenarios.

Results

We summarize the experimental results of both the simulation study and the application to real fMRI data in this section. With the simulation study, we demonstrated that acIVA accurately preserves the spatial variabilities across datasets. When applied to real fMRI data, we investigated the association between the tdBA and sdBA and study the activity-connectivity co-evolution by using spatial variabilities.

Simulation results

We had 100 independent realizations for each of the four cases introduced in the Simulation study section and summarized the results in Figure 3. We used the box plot for comparing values of the tuned $\hat{ρ}$ and the ground truth $ρ$ as well as the correlation between the estimates and true sources. The simulation results demonstrated that the values of $\hat{ρ}$ were very close to the ground truth values, suggesting that acIVA is able to tune the closest lower bound for the correlation between the estimates and reference signals in different scenarios. With the constraint of the reference signals, the estimates were accurately recovered and hence yielding high correlation values with the true sources. When applied to real fMRI data for a dynamic study where the sources come from a sequence of windows, the correlation between the spatial maps of the estimates and the reference signals can be a good measure for sdBA.

Application to fMRI data

We obtained $L = 16$ datasets for each subject by using a sliding window of length $T_{w} = 24$ with a 66.67% (2/3) overlap, yielding 2864 datasets in total. Instead of performing a single acIVA on the whole data, we divided the subjects into 45 subsets of $K_{0} = 4$ subjects and performed acIVA decomposition on each subset. Each subset has $K = L \times K_{0} = 64$ datasets, a value that is higher than the optimal value of the number of datasets that allows a reliable regular IVA decomposition (Long et al., 2020). It should be noted that the last subset of HC only has three subjects, hence 48 datasets. The model order used for each acIVA decomposition was 20, and the dimension of each dataset was reduced from 24 to 20 by using principal component analysis.

Using the common subspace extraction method, we obtained eight components that were shared across all 179 subjects. The spatial maps of $M = 8$ common components were used as reference signals to constrain 8 out of 20 components that were estimated. The eight common components include medial visual (RSN1), sensorimotor (RSN2), cerebellum (RSN3), DMN (RSN4), SP-V-C (RSN5), frontoparietal (RSN6), supplementary motor area (SMA, RSN7), and frontal (RSN8) components, as shown in Figure 4 (anatomical regions are summarized in Supplementary Table S1). An example of the tuned value of $\hat{ρ}$ and the computed value of $\tilde{ρ}$ across 16 windows for one schizophrenia subject and one HC subject is shown in Figure 5. We can clearly see that the spatial variabilities of dBA in the same functional network vary significantly across different subjects.

FIG. 4.

Spatial maps of the eight common RSN components extracted by using IVA-CS for HC and SZ groups. DMN, default mode network; HC, healthy control; IVA-CS, IVA for common subspace analysis; SMA, supplementary motor area; SP-V-C, super parietal, visual, and cerebellum; SZ, schizophrenia. Color images are available online.

FIG. 5.

The tuned values of $\tilde{ρ}$ in acIVA and the computed correlation $\tilde{ρ}$ between the estimates and the reference signals as a function of windows for an HC subject and an SZ subject (a). The spatial activation of the DMN component varies across windows in HC subject (b) and SZ subject (c). The values of spatial activation are normalized to be z-scores and thresholded by using z_t = 2. Color images are available online.

We also computed the normalized mutual information between the estimates and reference signals and compared their values with the correlation values $\tilde{ρ}$ . We plot normalized mutual information values against $\tilde{ρ}$ in Supplementary Figure S1 and showed that they are coherent. In this work, we used $\tilde{ρ}$ for the quantification of sdBA because these values were part of the acIVA to quantify the association between the estimates and reference signals. We used spatial maps of the DMN component as an example to demonstrate how $\tilde{ρ}$ quantifies spatial variability. High values of $\tilde{ρ}$ mean that the estimates are very similar to the reference signal, whereas low values mean that the estimates deviate from the reference signal, indicating higher variability. It shows that the DMN component varies more across windows for the HC individual than for the schizophrenia individual.

Co-evolution of activity and connectivity

We first investigated the association between sdBA and tdBA (referred to as Case 1 for simplicity in the following context) by computing the cross-correlation between $\tilde{ρ}$ and fALFF. The mean $\tilde{ρ}$ -fALFF cross-correlation matrices $\bar{C}$ that were averaged across the subjects in the schizophrenia group and the HC group and the one-sample t-test results are shown in Figure 6a. From Figure 6a, we see that the significant correlation values are mostly computed between sdBA and tdBA of the same network. The results indicate that the temporal variabilities and spatial variabilities of dBA of the same brain functional network are more likely to be significantly correlated than those of different networks. The schizophrenia group has fewer significantly correlated sdBA and tdBA pairs when compared with the HC group, indicating an abnormally organized brain function in schizophrenia.

FIG. 6.

Cross-correlation matrices of two types of dBA (a) and activity-connectivity co-evolution matrices computed by using dBA and dFNC (b) that are averaged across subjects in the SZ and HC groups. The values of t-statistic with an FDR control are shown. FDR, false discovery rate. Color images are available online.

The investigation of the activity-connectivity co-evolution in the spatial domain (referred to as Case 2 for simplicity in the following context) is shown in Figure 6b by demonstrating the mean co-evolution matrices $\bar{E}$ for two groups and the t-statistics of a one-sample t-test. There is a component (RSN5) yielding significant co-evolution between its sdBA and sdFNC in both schizophrenia and HC groups. This component is a complex component that merges multiple brain regions—SP-V-C—and it is the only component that shows significant correlation between sdBA and sdFNC for schizophrenia. Its sdBA activity significantly correlates to more connectivity patterns in sdFNC for the schizophrenia group compared with the HC group. It should be noted that these connectivity patterns are associated with this complex component. The HC group has an additional component (RSN7), the superior and middle frontal component (SMA), that yields significant co-evolution between its sdBA and sdFNC.

Identification of subgroups of schizophrenia using dynamic features

We identified subgroups of schizophrenia subjects by clustering $C^{[k]}$ and $E_{n}^{[k]}$ separately. To determine the number of subgroups $N_{s}$ , we ran the elbow algorithm 100 times, yielding the mean value and standard deviation of 4.31 ± 0.51 and 4.63 ± 0.69 for the two cases, respectively. We conducted an extensive investigation on the influence of $N_{s}$ by forming three scenarios where four, five, and six subgroups of schizophrenia subjects were identified by clustering the dynamic features into four, five, and six clusters.

The mean cross-correlation matrices ${\bar{C}}_{S G i}, i = 1, \dots, 6$ that were averaged across the subjects in each subgroup for the three scenarios are shown in Figure 7 for Case 1. The MANOVA test results show that there exist significant differences across subgroups in all three scenarios with different values of $N_{s}$ . We also plot the statistic values of the one-sample t-test and ANOVA on $C^{[k]}$ across subgroups to illustrate the correlation patterns. The results from Figure 7 show that each subgroup has a unique cross-correlation patterns and within each subgroup, there are statistically significant correlation values, except for subgroup 5 when $N_{s} = 6$ .

FIG. 7.

Mean cross-correlation matrices of each subgroup and the statistical test results (with an FDR control) for three scenarios with different values of N _s. The consistent subclusters of SZ subjects that are denoted by using numbers k = 1, …, 88 are listed under the correlation matrices using different colors. The colored font in bold refers to subclusters that are consistently detected in all three scenarios, and the regular colored font refers to subclusters that are detected in at least two scenarios. The font in black refers to subjects that do not belong to any consistent subclusters. The numbers in brackets denote the number of subjects in each subgroup. It should be noted that in the fifth subgroup (SG5) in the third scenario, all the correlation values are not significant with FDR control. ANOVA, analysis of variance; MANOVA, multivariate analysis of variance. Color images are available online.

We looked in detail at subjects that were clustered into each subgroup, and we identified some subject clusters that were consistently grouped together, as shown in Figure 7 with numbers denoting subject indices (( $k = 1, \dots, 88$ ) under the correlation matrices. We refer to these consistent subject clusters as subclusters. The numbers in bold font and color refer to subclusters that were detected in all three scenarios, and those in regular font and color refer to clusters that were detected in two arbitrary scenarios. The black numbers refer to subjects that did not belong to any consistent subclusters. The numbers in brackets denote the number of subjects in each subgroup. In Case 2 where the co-evolution matrices $E_{n}^{[k]}$ are used to identify the subgroups of schizophrenia subjects, we also detect significant differences across subgroups by using statistical tests and find clusters of subjects that are consistently grouped together in three scenarios of different values for $N_{s}$ . The results are summarized in Supplementary Figure S2.

We also compared the clinical symptoms measured by PANSS scores (Kay et al., 1987) across the schizophrenia subgroups that were identified by using neuroimaging data. There is a total of 7 positive, 7 negative, and 16 general PANSS scales that were measured for each subject. To compare the scores across the identified subgroups, we performed a MANOVA on five statistics—mean, standard deviation, median, minimum, and maximum—of all the 30 PANSS across subgroups. The statistics of each PANSS were computed by using the values of all subjects in a single subgroup. Therefore, in this MANOVA, five variables—mean, standard deviation, median, minimum, and maximum—with 30 samples were tested across the $N_{s}$ subgroups. The MANOVA detected significant differences across subgroups for the scenario of $N_{s} = 6$ with F-score = 2.1148 (p = 0.0012) when using cross-correlation feature and F-score = 1.7942 (p = 0.0100) when using co-evolution feature.

Positive, negative, and general PANSS are three independent categories that can be investigated separately (Kayahan et al., 2005; Peralta and Cuesta, 1994). Therefore, we also investigated the differences across the identified subgroups by looking at the correlation between their positive, negative, and general PANSS for the case when $N_{s} = 6$ . We found that the subgroups demonstrated different trends in the correlation between the positive and negative PANSS, as shown in Figure 8. In Case 1, subgroups 1 and 4 show positive correlation, subgroups 3 and 5 show negative correlation, and the positive and negative PANSS of subgroups 2 and 6 are almost not correlated. Though all subgroups showed positive correlation between their positive/negative PANSS and general PANSS, only the correlation of the two subgroups—subgroups 1 and 4—that showed positive correlations between their positive and negative PANSS was significant. In Case 2, four subgroups showed positive correlation between their positive and negative PANSS scores and the other two subgroups showed negative correlation. More subgroups showed significant positive correlation between their positive and general PANSS than between their negative and general PANSS.

FIG. 8.

Correlation between the positive, negative, and general positive and negative syndrome scale scores for each subgroup when N _s = 6 for Case 1 (a–c) and Case 2 (d–f). The large chart plots the data of six subgroups together for an easier visual comparison, and the small charts on its right plot the data of individual subgroups with the correlation value and its associated significance level listed. Color images are available online.

Discussion

In this work, we proposed the use of a semi-blind source separation method acIVA to extract both the spatial and temporal variations across multiple datasets. We first demonstrated that acIVA yielded preferable performance in preserving the spatial variabilities across multiple datasets when using the spatial maps as reference signals by conducting a simulation study. The use of acIVA in the analysis has a number of advantages. First, acIVA is a flexible model that makes use of reference signals, hence providing a more reliable estimation by reducing the effect of dimensionality issue in IVA. Second, through the use of reference signals, acIVA eliminates the alignment issue of components across multiple decompositions and enables the study of variabilities of a set of target components.

Although there are a significant number of works studying the functional dynamics in the human brain, most of these have a focus on temporal variabilities (Fu et al., 2018; Jafri et al., 2008; Weber et al., 2020). This is true, primarily because it is straightforward to withdraw temporal variability by applying the window strategy to the time series of raw fMRI or time courses of functional networks that are estimated from the original fMRI data. The extraction of spatial variability requires effective analysis of a large number of windowed datasets, such as the use of acIVA on the windowed datasets in this study.

There is evidence that functional patterns change in both spatial and temporal domains, thus also taking spatial variability into account is expected to enrich the dynamic study of the brain (Iraji et al., 2019a,b, 2020; Jie et al., 2018; Kottaram et al., 2018). Data-driven source separation techniques enable the extraction of dataset-specific variabilities (Bhinge et al., 2019a,b; Laney et al., 2014; Ma et al., 2014). However, a consistent estimation of functional networks from windowed datasets is challenging due to the iterative nature of data-driven algorithms and the alignment issue across multiple decompositions. The AcIVA desirably addresses these issues through the use of a set of reference signals and is able to extract both temporal and spatial variations from a large number of windowed datasets, allowing for a more comprehensive study of brain dynamics.

In the application of acIVA to the dynamic study of fMRI data, we investigated the effect of using different values as the window length and decided to obtain the windowed datasets by using a sliding window of length 24TR with a 66.67% (8TR) overlap. We ensured the window length to be in the suggested range of 30–60 sec for capturing brain dynamics (Allen et al., 2014; Damaraju et al., 2014; de Lacy et al., 2017; Rashid et al., 2014; Shirer et al., 2012). In our recent study using the same dataset (Long et al., 2020), there were 14 common components identified for subjects with schizophrenia and 16 for the HCs. Eight of these common components are aligned across the two groups. We decided to use 20 as the model order when applying acIVA to the windowed datasets to allow the estimation of all of the common components for each group as well as a small number of subject-specific components. Therefore, the smallest window length we can use in this study is 20TR (40 sec).

We designed a hybrid simulation study similar to that in Long et al. (2020) to investigate the effect of changes in window length and compare three scenarios where 20TR, 24TR (48 sec), and 28TR (56 sec) were used as window lengths. We summarized the results in Supplementary Figure S3. The results illustrate that there is no significant difference across the three scenarios. Therefore, we finally chose 24TR as the window length in this work since it is in the middle of the suggested range for studying brain dynamics.

Most previous studies compute dF(N)C matrices and do not require a decomposition of a large number of windowed datasets. They identify discrete dF(N)C states by using these matrices; hence, a sliding step of 1TR avoids missing a single potential state. In this work, we did not identify dFNC states but extracted dynamic features such as dBA and dFNC to study the correlation among the extracted dynamic features. To extract dynamic features, we performed acIVA on the windowed datasets.

We chose a relatively large sliding step of 8TR when dividing the data into overlapping windows to avoid a significant increase in the computational cost. With 24TR as the window length, we compared the $\tilde{ρ}$ -fALFF cross-correlation matrices for the cases when 1TR and 8TR was used as the sliding size separately. The results demonstrated that the cross-correlation values did not yield a significant difference between the two cases. Therefore, we chose 8TR as the sliding size, which helps to significantly decrease the computational cost without affecting the statistical analysis results. A relatively large sliding step was also commonly used in previous studies, especially those that require following processes such as decomposition for feature extraction (Bhinge et al., 2019b; Demirtaş et al., 2016). Nevertheless, it is of interest to conduct future work that uses reduced sliding steps to explore brain dynamics at faster time scales.

When applied to real fMRI data, we demonstrated that acIVA extracted both the temporal and spatial variabilities and provided an efficient measure of the sdBA in the study of brain dynamics. A desirable measure of sdBA enables us to investigate the cross-correlation between the two types of dBA as well as a first attempt in studying the activity-connectivity co-evolution patterns in the spatial domain. Our study of the cross-correlation between sdBA and tdBA illustrated that spatial and temporal variations of BOLD signal synchronously reflected the functional dynamics.

The results showed that the schizophrenia group had fewer significantly correlated sdBA and tdBA pairs when compared with the HC group, which may illustrate the dysconnectivity in the brain of patients with schizophrenia. The investigation of the activity-connectivity co-evolution in the spatial domain yielded a complex component that merges multiple brain regions—SP-V-C. This was the only component that yields significant co-evolution patterns for schizophrenia. Compared with HC, this complex component even yielded more significant co-evolution patterns; hence, it can be interpreted as performing a more centralized role in the brain for schizophrenia. The preference of this complex component suggested the reduced efficiency of brain function in schizophrenia, because more regions were simultaneously activated for a certain intrinsic function in their brain.

The additional SMA component (RSN7) that yielded significant activity-connectivity co-evolution patterns for the HC group may reveal that HCs show better functional organization in terms of the spatial variation patterns in the superior and middle frontal cortex compared with the schizophrenia group. The frontal cortex has been widely studied to better understand schizophrenia, and there is rich work showing the functional deficits in the frontal lobe for schizophrenia (Callicott et al., 2000, 2003; Carter et al., 1998; Hahn et al., 2018; Manoach, 2003; Mubarik and Tohid, 2016; Perlstein et al., 2001; Sheffield and Barch, 2016).

The dynamic features—cross-correlation and co-evolution patterns—were used to study schizophrenia heterogeneity by identifying schizophrenia subgroups. We validated the results by detecting significant differences across the identified subgroups and summarized the results in Figure 7 and Supplementary Figure S1. A study of the relationship between the findings of imaging data and clinical data is of great interest. PANSS scores are important clinical diagnostic data of schizophrenia and have been widely used in studies to gain a better interpretation of schizophrenia (Kayahan et al., 2005; Peralta and Cuesta, 1994).

We performed statistical tests on all 30 PANSS scores and compared the three independent PANSS categories across the schizophrenia subgroups that were identified by using imaging data. As shown in Figure 8a, schizophrenia subgroups 3 and 5 demonstrate a negative correlation between their positive and negative PANSS. These two subgroups yielded two extreme cases regarding the $\tilde{ρ}$ -fALFF cross-correlation matrices that are shown in Figure 7. Subgroup 3 yielded the most significant correlation values, whereas subgroup 5 yielded no significant correlation values. More specifically, the component that yielded the most significant correlation values in subgroup 3 was the DMN component (RSN4), highlighting the potential important role of DMN in illuminating some of the heterogeneity in schizophrenia.

There has been rich work showing that schizophrenia subjects demonstrate dysconnectivity in DMN (Mingoia et al., 2012; Van Den Heuvel and Pol, 2010) and this work provides novel evidence that the dynamics of DMN also help better interpret schizophrenia. More importantly, our work shows that biological data suggest multiple subgroups that also correspond to unique aspects of the clinical symptom scores. Although it is tempting to label these as schizophrenia subtypes, additional work is needed to further confirm this. Replication of these results would be important, and ideally, we would like to also evaluate a broader sampling of the psychosis spectrum for a further validation. However, this preliminary work showcases the potential of an approach that lets the biological data drive the categorizations by using the (often ignored) spatial dynamic information and suggests a way forward for linking self-reported symptoms to biological markers. With the increase in emphasis of personalized therapy for schizophrenia, a joint analysis of multi-domain data that helps gain a better understanding of the heterogeneity may help improve treatment strategies.

Given the interesting results discussed earlier, there are certain limitations of our work and some of the limitations lead to promising future research directions. First, we used dynamic features of eight aligned components that are common across subjects of both the schizophrenia and HC groups to identify the schizophrenia subgroups for the study of heterogeneity in schizophrenia since this work involves a comparative study between the two groups. However, it might be also reasonable to use components that are only common within the schizophrenia group and in this case more components might be involved, hence yielding more comprehensive observations. Second, we studied the association between two types of dBA and the co-evolution of dBA and dFNC by using the resting-state COBRE data and we drew interesting conclusions. The co-evolution of dBA and dFNC is a very promising research topic and should be further investigated in more datasets and other mental disorders.

Conclusion

We investigated the effectiveness of dynamic functional features extracted from fMRI in studying the heterogeneity of schizophrenia by identifying subgroups of schizophrenia subjects. Besides the dFNC information that is widely used in dynamic studies, BOLD activity is also shown to be related to mental and cognitive processes in the human brain, making a strong case for the desirability of using both BOLD activity and FC when studying dynamics. We proposed a novel use of a recent method, acIVA, to effectively capture tdBA and sdBA and to quantify sdBA as part of the algorithm. The efficient quantification of sdBA allows the study of the association between tdBA and sdBA as well as the activity-connectivity co-evolution in the spatial domain by using sdBA and sdFNC.

We conducted a simulation study to demonstrate that the adaptive parameter-tuning process in acIVA helps accurately capture the spatial variabilities and provides a nice metric to efficiently quantify the sdBA. The application of acIVA to the dynamic study of multi-subject resting-state fMRI data shows that dBA demonstrates synchronized patterns in temporal and spatial domains, and dBA and dFNC are significantly correlated in the spatial domain. The results illustrate that the brain function is differently organized in schizophrenia subjects compared with HCs. In addition, we identify unique subgroups of schizophrenia subjects that demonstrate different dynamic patterns by using the cross-correlation and co-evolution matrices separately. More importantly, significant differences are detected across the subgroups in terms of their clinical symptoms that are measured by PANSS. This observation inspires further study of the heterogeneity of mental disorders by using neuroimaging modalities such as fMRI.

Footnotes

Authors' Contributions

Q.L. conducted the project, including software development, method development, formal data analysis, visualization of results, and the writing of the article. S.B. contributed to software development, method development, and the writing of the article. V.D.C. contributed to project conceptualization, funding acquisition, and writing of the article, and provided data and related resources. T.A. conceptualized the project, supervised and administered the project, provided research resources and funding, and contributed to writing of the article.

Acknowledgments

The authors thank the research staff from the Mind Research Network COBRE study who collected, preprocessed, and shared the data. The authors appreciate the valuable feedback provided by the members of Machine Learning for Signal Processing Laboratory at the University of Maryland, Baltimore County. The hardware used in the computational studies is part of the UMBC High Performance Computing Facility (HPCF). The facility is partially supported by the U.S. National Science Foundation through the MRI program (OAC–1726023). See hpcf.umbc.edu for more information on HPCF and the projects using its resources.

Author Disclosure Statement

No competing financial interests exist.

Funding Information

Q.L., S.B., V.D.C., and T.A. were supported by NSF grants CCF 1618551 and NCS 1631838, and NIH grant R01MH 118695. V.D.C. was also supported by NIH grant R01EB 020407.

Supplementary Material

Supplementary Figure S1

Supplementary Figure S2

Supplementary Figure S3

Supplementary Table S1

References

Adali

, Anderson

, Fu

G-S

. 2014. Diversity in independent component and vector analyses: identifiability, algorithms, and applications in medical imaging. IEEE Signal Process Mag, 31:18–33.

Aine

, Bockholt

, Bustillo

, et al. 2017. Multimodal neuroimaging in schizophrenia: description and dissemination. Neuroinformatics, 15:343–364.

Allen

, Damaraju

, Plis

SMet al

. 2014. Tracking whole-brain connectivity dynamics in the resting state. Cerebral Cortex, 24:663–676.

Allen

, Erhardt

, Damaraju

, et al. 2011. A baseline for the multivariate comparison of resting-state networks. Front Syst Neurosci, 5:2.

Assaf

, Jagannathan

, Calhoun

, et al. 2010. Abnormal functional connectivity of default mode sub-networks in autism spectrum disorder patients. NeuroImage, 53:247–256.

Baiano

, David

, Versace

, et al. 2007. Anterior cingulate volumes in schizophrenia: a systematic review and a meta-analysis of MRI studies. Schizophr Res, 93:1–12.

Bhinge

, Long

, Calhoun

, et al. 2019a. Spatial dynamic functional connectivity analysis identifies distinctive biomarkers in schizophrenia. Front Neurosci, 13:1006.

Bhinge

, Long

, Levin-Schwartz

, et al. Non-orthogonal constrained independent vector analysis: application to data fusion. In 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), New Orleans, LA, USA, 2017, pp. 2666–2670.

Bhinge

, Mowakeaa

, Calhoun

, et al. 2019b. Extraction of time-varying spatiotemporal networks using parameter-tuned constrained IVA. IEEE Trans Med Imag, 38:1715–1725.

10.

Bluhm

, Miller

, Lanius

, et al. 2007. Spontaneous low-frequency fluctuations in the BOLD signal in schizophrenic patients: anomalies in the default network. Schizophr Bull, 33:1004–1012.

11.

Britz

, Landis

, Michel

. 2009. Right parietal brain activity precedes perceptual alternation of bistable stimuli. Cerebral Cortex, 19:55–65.

12.

Calhoun

, Adalı

. 2009. Feature-based fusion of medical imaging data. IEEE Trans Inform Technol Biomed, 13:711–720.

13.

Calhoun

, Adalı

, Pearlson

, et al. 2001. A method for making group inferences from functional MRI data using independent component analysis. Hum Brain Mapp, 14:140–151.

14.

Calhoun

, Adalı

, Pearlson

, et al. 2006. Neuronal chronometry of target detection: fusion of hemodynamic and event related potential data. NeuroImage, 30:544–553.

15.

Calhoun

, Eichele

, Pearlson

. 2009. Functional brain networks in schizophrenia: a review. Front Hum Neurosci, 3:17.

16.

Calhoun

, Miller

, Pearlson

, et al. 2014. The chronnectome: time-varying connectivity networks as the next frontier in fMRI data discovery. Neuron, 84:262–274.

17.

Callicott

, Bertolino

, Mattay

, et al. 2000. Physiological dysfunction of the dorsolateral prefrontal cortex in schizophrenia revisited. Cerebral Cortex, 10:1078–1092.

18.

Callicott

, Mattay

, Verchinski

, et al. 2003. Complexity of prefrontal cortical dysfunction in schizophrenia: more than up or down. Am J Psychiatry, 160:2209–2215.

19.

Carter

, Perlstein

, Ganguli

, et al. 1998. Functional hypofrontality and working memory dysfunction in schizophrenia. Am J Psychiatry, 155:1285–1287.

20.

Çetin

, Christensen

, Abbott

, et al. 2014. Thalamus and posterior temporal lobe show greater inter-network connectivity at rest and across sensory paradigms in schizophrenia. NeuroImage, 97:117–126.

21.

Damaraju

, Allen

, Belger

, et al. 2014. Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia. NeuroImage, 5:298–308.

22.

de Lacy

, Doherty

, King

, et al. 2017. Disruption to control network function correlates with altered dynamic connectivity in the wider autism spectrum. NeuroImage, 15:513–524.

23.

Demirtaş

, Tornador

, Falcón

, et al. 2016. Dynamic functional connectivity reveals altered variability in functional connectivity among patients with major depressive disorder. Hum Brain Mapp, 37:2918–2930.

24.

Dwyer

, Cabral

, Kambeitz-Ilankovic

, et al. 2018. Brain subtyping enhances the neuroanatomical discrimination of schizophrenia. Schizophr Bull, 44:1060–1069.

25.

Etkin

, Prater

, Schatzberg

, et al. 2009. Disrupted amygdalar subregion functional connectivity and evidence of a compensatory network in generalized anxiety disorder. Arch Gen Psychiatry, 66:1361–1372.

26.

Freire

, Roche

, Mangin

J-F

. 2002. What is the best similarity measure for motion correction in fMRI time series?. IEEE Trans Med Imag, 21:470–484.

27.

Friston

, Holmes

, Worsley

, et al. 1994. Statistical parametric maps in functional imaging: a general linear approach. Hum Brain Mapp, 2:189–210.

28.

, Tu

, Di

, et al. 2018. Characterizing dynamic amplitude of low-frequency fluctuation and its relationship with dynamic functional connectivity: an application to schizophrenia. NeuroImage, 180:619–631.

29.

Geisler

, Walton

, Naylor

, et al. 2015. Brain structure and function correlates of cognitive subtypes in schizophrenia. Psychiatry Res Neuroimag, 234:74–83.

30.

Guller

, Tononi

, Postle

. 2012. Conserved functional connectivity but impaired effective connectivity of thalamocortical circuitry in schizophrenia. Brain Connect, 2:311–319.

31.

Hahn

, Robinson

, Leonard

, et al. 2018. Posterior parietal cortex dysfunction is central to working memory storage and broad cognitive deficits in schizophrenia. J Neurosci, 38:8378–8387.

32.

Hallmayer

, Kalaydjieva

, Badcock

, et al. 2005. Genetic evidence for a distinct subtype of schizophrenia characterized by pervasive cognitive deficit. Am J Hum Genet, 77:468–476.

33.

Hoptman

, Zuo

X-N

, Butler

, et al. 2010. Amplitude of low-frequency oscillations in schizophrenia: a resting state fMRI study. Schizophr Res, 117:13–20.

34.

Hutchison

, Morton

. 2016. It's a matter of time: reframing the development of cognitive control as a modification of the brain's temporal dynamics. Dev Cogn Neurosci, 18:70–77.

35.

Hutchison

, Womelsdorf

, Gati

, et al. 2013. Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques. Hum Brain Mapp, 34:2154–2177.

36.

Iraji

, Deramus

, Lewis

, et al. 2019a. The spatial chronnectome reveals a dynamic interplay between functional segregation and integration. Hum Brain Mapp, 40:3058–3077.

37.

Iraji

, Fu

, Damaraju

, DeRamus

, et al. 2019b. Spatial dynamics within and between brain functional domains: a hierarchical approach to study time-varying brain function. Hum Brain Mapp, 40:1969–1986.

38.

Iraji

, Miller

, Adali

, et al. 2020. Space: a missing piece of the dynamic puzzle. Trends Cogn Sci, 24:135–149.

39.

Jablensky

2006. Subtyping schizophrenia: implications for genetic research. Mol Psychiatry, 11:815.

40.

Jablensky

2010. The diagnostic concept of schizophrenia: its history, evolution, and future prospects. Dialog Clin Neurosci, 12:271.

41.

Jafri

, Pearlson

, Stevens

, et al. 2008. A method for functional network connectivity among spatially independent resting-state components in schizophrenia. NeuroImage, 39:1666–1681.

42.

Jie

, Liu

, Shen

. 2018. Integration of temporal and spatial properties of dynamic connectivity networks for automatic diagnosis of brain disease. Med Image Analysis, 47:81–94.

43.

Kay

, Fiszbein

, Opler

L. A

. 1987. The positive and negative syndrome scale (PANSS) for schizophrenia. Schizophr Bull, 13:261–276.

44.

Kayahan

, Ozturk

, Veznedaroglu

, et al. 2005. Obsessive-compulsive symptoms in schizophrenia: prevalance and clinical correlates. Psychiatry Clin Neurosci, 59:291–295.

45.

Kim

, Eltoft

, Lee

T-W

. Independent vector analysis: an extension of ICA to multivariate components. In International Conference on Independent Component Analysis and Signal Separation, London, UK, 2006, pp. 165–172.

46.

Kottaram

, Johnston

, Ganella

, et al. 2018. Spatio-temporal dynamics of resting-state brain networks improve single-subject prediction of schizophrenia diagnosis. Hum Brain Mapp, 39:3663–3681.

47.

Kucyi

, Davis

. 2014. Dynamic functional connectivity of the default mode network tracks daydreaming. NeuroImage, 100:471–480.

48.

Kühn

, Gallinat

. 2013. Resting-state brain activity in schizophrenia and major depression: a quantitative meta-analysis. Schizophr Bull, 39:358–365.

49.

Laney

, Westlake

, Ma

, et al. Capturing subject variability in data driven fMRI analysis: a graph theoretical comparison. In 48th Annual Conference on Information Sciences and Systems (CISS), Princeton, NJ, USA, 2014, pp. 1–6.

50.

Lefebvre

, Demeulemeester

, Leroy

, et al. 2016. Network dynamics during the different stages of hallucinations in schizophrenia. Hum Brain Mapp, 37:2571–2586.

51.

X-L

, Zhang

X-D

. 2007. Nonorthogonal joint diagonalization free of degenerate solution. IEEE Trans Signal Process, 55:1803–1814.

52.

Long

, Bhinge

, Calhoun

, et al. 2020. Independent vector analysis for common subspace analysis: application to multi-subject fMRI data yields meaningful subgroups of schizophrenia. NeuroImage, 2016:116872.

53.

, Calhoun

, Eichele

, et al. 2012. Modulations of functional connectivity in the healthy and schizophrenia groups during task and rest. NeuroImage, 62:1694–1704.

54.

, Calhoun

, Phlypo

, et al. 2014. Dynamic changes of spatial functional network connectivity in healthy individuals and schizophrenia patients using independent vector analysis. NeuroImage, 90:196–206.

55.

Malaspina

, Harkavy-Friedman

, Corcoran

, et al. 2004. Resting neural activity distinguishes subgroups of schizophrenia patients. Biol Psychiatry, 56:931–937.

56.

Manoach

DS.

2003. Prefrontal cortex dysfunction during working memory performance in schizophrenia: reconciling discrepant findings. Schizophr Res, 60:285–298.

57.

Marusak

, Calhoun

, Brown

, et al. 2017. Dynamic functional connectivity of neurocognitive networks in children. Hum Brain Mapp, 38:97–108.

58.

McIntosh

, Kovacevic

, Itier

. 2008. Increased brain signal variability accompanies lower behavioral variability in development. PLoS Comput Biol, 4:e1000106.

59.

Mennigen

, Fryer

, Rashid

, et al. 2019. Transient patterns of functional dysconnectivity in clinical high risk and early illness schizophrenia individuals compared with healthy controls. Brain Connect, 9:60–76.

60.

Mingoia

, Wagner

, Langbein

, et al. 2012. Default mode network activity in schizophrenia studied at resting state using probabilistic ICA. Schizophr Res, 138:143–149.

61.

Morar

, Badcock

, Phillips

, et al. 2018. The longevity gene Klotho is differentially associated with cognition in subtypes of schizophrenia. Schizophr Res, 193:348–353.

62.

Mubarik

, Tohid

. 2016. Frontal lobe alterations in schizophrenia: a review. Trends Psychiatry Psychother, 38:198–206.

63.

Peralta

, Cuesta

. 1994. Psychometric properties of the positive and negative syndrome scale (PANSS) in schizophrenia. Psychiatry Res, 53:31–40.

64.

Pereira

, Direito

, Sayal

, et al. 2019. Self-modulation of premotor cortex interhemispheric connectivity in a real-time functional magnetic resonance imaging neurofeedback study using an adaptive approach. Brain Connect, 9:662–672.

65.

Perlstein

, Carter

, Noll

, et al. 2001. Relation of prefrontal cortex dysfunction to working memory and symptoms in schizophrenia. Am J Psychiatry, 158:1105–1113.

66.

Rashid

, Damaraju

, Pearlson

, et al. 2014. Dynamic connectivity states estimated from resting fMRI identify differences among schizophrenia, bipolar disorder, and healthy control subjects. Front Hum Neurosci, 8:897.

67.

Rencher

AC.

2003. Methods of Multivariate Analysis, Vol. 492. New Jersey, USA: John Wiley &. Sons.

68.

Sakoğlu

, Pearlson

, Kiehl

, et al. 2010. A method for evaluating dynamic functional network connectivity and task-modulation: application to schizophrenia. MAGMA, 23:351–366.

69.

Sayın

, Yüksel

, Konac

, et al. 2013. Effects of the adverse life events and disrupted in schizophrenia-1 (DISC1) gene polymorphisms on acute symptoms of schizophrenia. DNA Cell Biol, 32:73–80.

70.

Schultz

, Koch

, Wagner

, et al. 2012. Reduced anterior cingulate cognitive activation is associated with prefrontal-temporal cortical thinning in schizophrenia. Biol Psychiatry, 71:146–153.

71.

Scott

, Courtney

, Wood

, et al. 2011. COINS: an innovative informatics and neuroimaging tool suite built for large heterogeneous datasets. Front Neuroinform, 5:33.

72.

Senn

2018. Statistical pitfalls of personalized medicine. Nature, 563:619–621.

73.

Sheffield

, Barch

. 2016. Cognition and resting-state functional connectivity in schizophrenia. Neurosci Biobehav Rev, 61:108–120.

74.

Shin

, Ahn

, Hu

. 2013. Correction for the T1 effect incorporating flip angle estimated by Kalman filter in cardiac-gated functional MRI. Magn Reson Med, 70:1626–1633.

75.

Shirer

, Ryali

, Rykhlevskaia

, et al. 2012. Decoding subject-driven cognitive states with whole-brain connectivity patterns. Cerebral Cortex, 22:158–165.

76.

Turner

, Damaraju

, Van Erp

, et al. 2013. A multi-site resting state fMRI study on the amplitude of low frequency fluctuations in schizophrenia. Front Neurosci, 7:137.

77.

Uddin

, Supekar

, Ryali

, et al. 2011. Dynamic reconfiguration of structural and functional connectivity across core neurocognitive brain networks with development. J Neurosci, 31:18578–18589.

78.

van de Ven

, Formisano

, Prvulovic

, et al. 2004. Functional connectivity as revealed by spatial independent component analysis of fMRI measurements during rest. Hum Brain Mapp, 22:165–178.

79.

Van Den Heuvel

, Pol

HEH

. 2010. Exploring the brain network: a review on resting-state fMRI functional connectivity. Eur Neuropsychopharmacol, 20:519–534.

80.

Verley

, Torolira

, Pulido

, et al. 2018. Remote changes in cortical excitability after experimental traumatic brain injury and functional reorganization. J Neurotrauma, 35:2448–2461.

81.

Wang

, Yu

, Chen

, et al. 2010. Dynamic functional reorganization of the motor execution network after stroke. Brain, 133:1224–1238.

82.

Weber

, Johnsen

, Kroken

, et al. 2020. Dynamic functional connectivity patterns in schizophrenia and the relationship with hallucinations. Front Psychiatry, 11:227.

83.

Zang

Y-F

, He

, Zhu

C-Z

, et al. 2007. Altered baseline brain activity in children with ADHD revealed by resting-state functional MRI. Brain Dev, 29:83–91.

84.

Zhang

, Li

, Wang

, et al. 2018. Abnormal dynamic functional connectivity between speech and auditory areas in schizophrenia patients with auditory hallucinations. NeuroImage, 19:918–924.

85.

Zhou

, Liang

, Tian

, et al. 2007. Functional disintegration in paranoid schizophrenia using resting-state fMRI. Schizophr Res, 97:194–205.

86.

Zou

Q-H

, Zhu

C-Z

, Yang

, et al. 2008. An improved approach to detection of amplitude of low-frequency fluctuation (ALFF) for resting-state fMRI: fractional ALFF. J Neurosci Methods, 172:137–141.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.20 MB

2.87 MB

0.10 MB

0.09 MB