Explainable fuzzy clustering framework reveals divergent default mode network connectivity dynamics in schizophrenia

Introduction Dynamic functional network connectivity (dFNC) analysis of resting state functional magnetic resonance imaging data has yielded insights into many neurological and neuropsychiatric disorders. A common dFNC analysis approach uses hard clustering methods like k-means clustering to assign samples to states that summarize network dynamics. However, hard clustering methods obscure network dynamics by assuming (1) that all samples within a cluster are equally like their assigned centroids and (2) that samples closer to one another in the data space than to their centroids are well-represented by their centroids. In addition, it can be hard to compare subjects, as in some cases an individual may not manifest a state strongly enough to enter a hard cluster. Approaches that allow a dimensional approach to connectivity patterns (e.g., fuzzy clustering) can mitigate these issues. In this study, we present an explainable fuzzy clustering framework by combining fuzzy c-means clustering with several explainability metrics and novel summary features. Methods We apply our framework for schizophrenia (SZ) default mode network analysis. Namely, we extract dFNC from individuals with SZ and controls, identify 5 dFNC states, and characterize the dFNC features most crucial to those states with a new perturbation-based clustering explainability approach. We then extract several features typically used in hard clustering and further present a variety of unique features specially designed for use with fuzzy clustering to quantify state dynamics. We examine differences in those features between individuals with SZ and controls and further search for relationships between those features and SZ symptom severity. Results Importantly, we find that individuals with SZ spend more time in states of moderate anticorrelation between the anterior and posterior cingulate cortices and strong anticorrelation between the precuneus and anterior cingulate cortex. We further find that individuals with SZ tend to transition more rapidly than controls between low-magnitude and high-magnitude dFNC states. Conclusion We present a novel dFNC analysis framework and use it to identify effects of SZ upon network dynamics. Given the ease of implementing our framework and its enhanced insight into network dynamics, it has great potential for use in future dFNC studies.

Introduction: Dynamic functional network connectivity (dFNC) analysis of resting state functional magnetic resonance imaging data has yielded insights into many neurological and neuropsychiatric disorders.A common dFNC analysis approach uses hard clustering methods like k-means clustering to assign samples to states that summarize network dynamics.However, hard clustering methods obscure network dynamics by assuming (1) that all samples within a cluster are equally like their assigned centroids and (2) that samples closer to one another in the data space than to their centroids are wellrepresented by their centroids.In addition, it can be hard to compare subjects, as in some cases an individual may not manifest a state strongly enough to enter a hard cluster.Approaches that allow a dimensional approach to connectivity patterns (e.g., fuzzy clustering) can mitigate these issues.In this study, we present an explainable fuzzy clustering framework by combining fuzzy c-means clustering with several explainability metrics and novel summary features.
Methods: We apply our framework for schizophrenia (SZ) default mode network analysis.Namely, we extract dFNC from individuals with SZ and controls, identify 5 dFNC states, and characterize the dFNC features most crucial to those states with a new perturbation-based clustering explainability approach.We then extract several features typically used in hard clustering and further present a variety of unique features specially designed for use with fuzzy clustering to quantify state dynamics.We examine differences in those features between individuals with SZ and controls and further search for relationships between those features and SZ symptom severity.
Results: Importantly, we find that individuals with SZ spend more time in states of moderate anticorrelation between the anterior and posterior cingulate cortices and strong anticorrelation between the precuneus and anterior cingulate cortex.We further find that individuals with SZ tend to transition more rapidly than controls between low-magnitude and high-magnitude dFNC states.

Introduction
Resting state functional magnetic resonance imaging (rs-fMRI) dynamic functional network connectivity (dFNC) data has historically been used to give insight into a variety of neurological (1) and neuropsychiatric disorders (2)(3)(4)(5)(6) and cognitive functions (7,8).A common dFNC analysis approach involves applying a hard clustering approach (e.g., k-means clustering) to assign dFNC samples to a set of dynamical states that are supposedly representative of the overall time series (2,(9)(10)(11)(12)(13)(14)(15).Features can then be extracted based on the identified time-resolved states that can give insight into various aspects of the state dynamics.This method has been widely applied but makes a critical assumption that could obscure useful disease-related dynamics.Namely, it assumes that all samples assigned to a state equally belong to a state.This is problematic given that samples very near to one another in the data space may be assigned to two distant cluster centroids.A few studies have presented fuzzy clustering approaches that indicate the degree to which samples belong to different states, but this area is still highly understudied (11,15).Although explainability methods have been developed specifically to help characterize identified states, the use of hard clustering methods can make explainability difficult.In this study, we present a novel explainable fuzzy clustering framework for fMRI dFNC that identifies fuzzy states and assigns samples a probability of belonging to each state.We further present novel dynamical features that use the output probabilities and demonstrate their utility by applying them to identify differences in the dynamics of default mode network (DMN) activity between individuals with schizophrenia (SZs) and healthy controls (HCs).
Several modalities have been used for insight into the effects of neurological and neuropsychiatric disorders upon brain dynamics.These include electroencephalography (EEG), magnetoencephalography (MEG), and fMRI.All three modalities -EEG (16-18), MEG (19), and fMRI (13,(20)(21)(22)(23)(24) -have been used extensively for SZ analysis.EEG and MEG capture much higher resolution temporal information.However, localizing the region of the brain associated with MEG and EEG signals can be challenging.In contrast, fMRI has much higher spatial resolution at a lower temporal resolution relative to EEG and fMRI.A select group of studies have also combined multimodal EEG, MEG, or fMRI data for insight into disorders (9).However, in general, fMRI is more often used in SZ analysis than the other modalities.Within fMRI analysis, both task (25,26) and resting state (10,12,13,(27)(28)(29)(30)(31)(32) data are frequently used.However, resting state data offers several advantages.Specifically, the majority of brain activity is spontaneous (i.e., better captured by resting state), so resting state analysis provides an avenue to understand how the brain operates under most circumstances (32).Additionally, task performance in healthy controls relative to individuals with schizophrenia often varies greatly, introducing a potential confounder into any neuroimaging analyses (32).
Many studies have analyzed resting state fMRI data within the context of SZ and other neuropsychiatric disorders.For example, studies have used independent components (ICs) (24), spectral features (9,33), and functional network connectivity.Functional network connectivity offers a unique benefit in that it provides insights into the interaction of different brain regions and networks.Early in the use of rs-fMRI functional network connectivity, it was more common to analyze the correlation between brain regions across a whole recording.This approach is referred to as static functional network connectivity (sFNC) (34).Although sFNC had widespread use, a number of studies found that dFNC (i.e., functional connectivity captured in windows over time) offered insights into brain interactions that would otherwise be obscured by sFNC analysis (32,34).Both sFNC and dFNC have been used to gain insight into a variety of neurological and neuropsychiatric disorders and cognitive functions, including Alzheimer's disease (1,35,36), major depressive disorder (3-6), schizophrenia (2, 37), cognition (7), and spatial orientation (8).However, dFNC offers greater opportunities to learn about the brain than sFNC.As discussed in (32,34), early studies extracted functional network connectivity from time-series extracted using one seed per brain region, multiple seeds from a given region of interest (ROI), multiple seeds from subregions of multiple ROIs, and seeds from whole-brain ROIs.However, in recent years, the fully automated, group independent component analysis (ICA)-based Neuromark pipeline has been developed as a approach for extracting time-series used to calculate functional network connectivity (38).It yields components that are reproducible across datasets and studies and contains components from a variety of brain networks and subregions.Furthermore, it has been used in many studies (1,3,24,(35)(36)(37)(38)(39)(40)(41)(42).
Multiple approaches have been used to analyze dFNC extracted using the NeuroMark pipeline or other approaches.Many studies have used classification approaches (21,22,43) or a combination of clustering and classification (39).However, a described in (32), a many studies have used clustering approaches (9)(10)(11)(12)(13)(14)(15).These clustering approaches involve assigning dFNC samples to states that summarize the dFNC time-series.Features can then be extracted to quantify various aspects of the state time-series.By far the most common clustering approach used for identifying states of dFNC activity is k-means clustering (44).K-means clustering involves randomly initializing cluster centroids, calculating the average of the samples nearest each centroid, updating the cluster centroids to be equal to the average of the nearest samples, and iteratively repeating the process up to the point that the cluster centroid stops moving at each iteration above a pre-defined threshold.The main advantage of k-means clustering for use with dFNC is that it is easy to implement using existing libraries (45).However, while k-means clustering is widely used, it does have some disadvantages.First, it can yield low quality clusters.Previous studies have proposed approaches to try to address this problem (14).Additionally, it does not assign each sample a probability of belonging to each cluster centroid or resting-state, which involves an implicit assumption in subsequent analysis that the identified states are able to adequately summarize the dFNC time-series (46).This is a big assumption given that samples very near one another in the data space can be assigned to completely different states.In addition, in the hard clustering approach it can be hard to compare individuals as it is likely that some individuals may never enter a given state strongly enough to enter the cluster.
As such, a simple fuzzy clustering-based approach could go a long way towards helping future studies in the field of dFNC clustering to uncover new aspects of disorder-related brain dynamics by better summarizing brain activity.
As the field of rs-fMRI FNC clustering has developed, a few studies have introduced explainability approaches with the goal of helping characterize the differences between identified states (7,23,47).In contrast to earlier efforts involving statistical hypothesis testing (38), these methods quantify the importance of FNC features in a manner that acknowledges the multi-dimensional nature of the underlying clustering.However, these methods have an important shortcoming that arises from their use with hard clustering methods.Specifically, they perturb FNC features and examine the sensitivity of the underlying clusters to perturbation.They quantify the effect of this perturbation upon the clustering by calculating the percentage of samples that completely switch clusters.Given that only a small minority of samples switch clusters following perturbation, it is safer to say that the methods quantify the effects of perturbation upon that subset of samples rather than upon the overall clustering.Similar to how a fuzzy clustering approach could help future studies better summarize brain activity and uncover new aspects of disorders, fuzzy clustering could also contribute to the development of explainability approaches that are better able to quantify the effects of perturbation.
In this study, we present fuzzy c-means clustering (48) as an approach for the identification of fMRI dFNC fuzzy states.Fuzzy cmeans has been used widely across application areas like emotion recognition from speech data (49), customer segmentation (50), and  4) We applied several explainability approaches for insight into the dFNC features characterizing the states.(5,6) We extracted dynamical and stability features.(7) We performed t-tests and trained interpretable machine learning models for insight into the features that differed between HCs and SZs.(8) We performed linear regression analyses controlling for age and gender to identify relationships between symptom severity and dynamical and stability features.brain image segmentation (51).Furthermore, fuzzy c-means is simple to use given its inclusion in several publicly available Python packages (52) and MATLAB (53).As such, it has the potential to be widely used within the domain of FNC clustering.
Our approach outputs a probability that each sample belongs to each fuzzy state.After performing clustering, we present a pair of new explainability approaches to characterize the identified states, comparing them to an existing approach (47).We then present a variety of novel dynamical and stability features that highlight different aspects of the fuzzy state time-series while also showing that our approach is compatible with dynamical features that have been used in previous hard clustering studies.We further show how the features can be used to differentiate between healthy controls and individuals with schizophrenia in the default mode network (DMN), a network that has previously been associated with SZ (10,12,13,(27)(28)(29)(30)(31)(32), and identify relationships between the features and SZ and SZ symptom severity.

Methods
In this section, we describe and discuss our approach for the study.Figure 1 presents an overview of our approach (1).We used a pre-existing rs-fMRI dataset composed of individuals with schizophrenia (SZs) and healthy controls (HCs).( 2) We preprocessed the data and extracted dFNC.(3) We applied fuzzy c-means clustering to identify 5 fuzzy states.( 4) We sought to identify the important dFNC features for each state by applying and comparing the results for two global clustering explainability methods.(5) We next extracted a number of pre-existing and novel dynamical features to summarize different aspects of the identified states.(6) We applied a novel local cluster explainability approach for insight into the stability of study participants to the perturbation of specific dFNC features.(7) To identify SZ-related differences in the extracted dynamical and stability features, we conducted statistical analyses and trained a series of logistic regression classifiers with elastic net regularization (LR-ENR).( 8) Lastly, we performed statistical analyses to determine whether the dynamical and stability features were related to symptom severity.For the sake of reproducibility, our code can be found at: https:// github.com/cae67/Fuzzy_Clustering/.

Description of dataset
In this study, we used the Functional Imaging Biomedical Informatics Research Network (FBIRN) dataset, consisting of rs-fMRI recordings from 151 SZs and 160 HCs (54).Participant demographics are shown in Table 1.The dataset has been used in many studies, both related to fMRI dFNC clustering and classification (21, 24,43,47).In addition to neuroimaging data, positive and negative symptom severity scores from the Positive and Negative Syndrome Scale (PANSS) (55).Positive symptoms of SZ include hallucinations, delusions, and bizarre behavior.Negative symptoms include alogia, apathy, affective flattening, and asociality (28).The dataset was collected at 7 sites: the University of California at Los Angeles, the University of California at San Francisco, the University of California at Irvine, the University of Iowa, the University of Minnesota, the University of New Mexico, and Duke University/the University of North Carolina at Chapel Hill.Data collection procedures were approved by the institutional review boards of each center, and all study participants gave written informed consent.One site used a 3T GE MR750 scanner, and 6 sites used 3T Siemens TIM Trio Scanners.Data was collected using a T2*weighted AC-PC aligned EPI sequence (TE = 30ms, TR = 2s, slice gap = 1mm, flip angle = 77°, voxel size = 3.4 x 3.4 x 3.4 mm 3 , acquisition time = 5 min and 38s, and number of frames = 162).

Description of data preprocessing
Prior to preprocessing, the first 5 mock scans were removed.Statistical parametric mapping (SPM12, https://www.fil.ion.ucl.ac.uk/ spm/) was used for preprocessing, and head motion was corrected using rigid body motion correction.The data was spatially normalized to an echo-planar imaging template in the standard Montreal Neurological Institute (MNI) space.Following resampling to 3x3x3 mm 3 , a Gaussian kernel with a 6mm full width at half maximum was used to smooth the data.The fully automated Neuromark Pipeline of the Group ICA of fMRI Toolbox (GIFT, http://trendscenter.org/software/gift) involved the use of spatially constrained ICA to extract corresponding components while adaptive to individual datasets.In our case, we used the neuromark_fMRI_1.0template to extract 53 independent components (ICs) with peak activations in the gray matter of various brain networks.Seven of the ICs were associated with the DMN, which has been associated with SZ in multiple studies (10,13,28), and we used those components in this study.The 7 ICs included 3 precuneus (PCN) (56), 2 anterior cingulate cortex (ACC) (57), and 2 posterior cingulate cortex (PCC) (58).After IC extraction, dFNC was estimated using Pearson's correlation with a sliding tapered window.The window consisted of a rectangle with a 40-second step size convolved with a Gaussian (s=3).The window size parameter plays an important role in study outcomes.A shorter window size increases dynamics but can increase noise and susceptibility to artifacts, while longer window sizes decrease dynamics and are more robust to artifacts (59).While a number of window sizes have been utilized in previous studies, a 40-second size is highly common (13,23,37,60).Additionally, a Gaussian with s=3 is also highly common (3, 13,36,37).The dFNC extraction resulted in a dataset composed of 21 dFNC features and 124 time steps per participant.
For the remainder of this paper, correlation between two brain regions is abbreviated using the pattern of IC1/IC2 (e.g., PCC1/ ACC2), where IC1/IC2 is equivalent to IC2/IC1.

Description of clustering approach
After extracting dFNC, we concatenated all dFNC time steps across samples and applied fuzzy c-means clustering (48).Fuzzy cmeans clustering is comparable to k-means clustering, but unlike kmeans clustering, fuzzy c-means assigns each sample a probability of belonging to a given cluster.We used the Python package scikit-fuzzy in our implementation (52) and optimized the random seed used for initialization and the fuzziness parameter, m, based on the fuzzy partition coefficient (61).It is relatively common to set the number of clusters to a number identified in previous studies (62).As such, we used 5 clusters to more easily relate our findings to previous studies (13,62) of SZ dFNC data.Further details on our cluster parameter optimization can be found in the Supplementary Materials.

Description of explainability approaches
We used three global explainability approaches for insight into the relative importance of each dFNC feature to the identified fuzzy states.

Global permutation percent change feature importance
We applied Global Permutation Percent Change (G2PC) feature importance (47), which has been used in several neuroimaging studies (7,23).G2PC extends permutation feature importance (63, 64) to the clustering domain, wherein individual features are permuted and the features that cause the greatest percentage of samples that switch clusters are considered most important.Further details on G2PC can be found in the Supplementary Materials.

Permutation-based distribution divergence and G2PC comparison
We used two variations of a novel approach called, "Permutation-based Distribution Divergence (P2D)".P2D has global (GP2D) and local (LP2D) variations (65).P2D is similar to G2PC.However, rather than computing the percent of samples that switch clusters after permutation, P2D involves calculating the Kullback-Leibler divergence (KLD) between the probabilities of a sample belonging to each cluster before versus after permutation.For GP2D, the median and total KLD across samples were calculated, and the rankings of GP2D feature importance were compared with the G2PC feature importance rankings (66)(67)(68).For our analysis with LP2D, we (1) compared the percentage of samples with non-zero KLD values with the G2PC results to see if P2D was more sensitive than G2PC and (2) compared HC and SZ KLD values to identify differences in cluster stability.Further details on our P2D and comparison analyses can be found in the Supplementary Materials.

Description of dynamical feature extraction
After clustering all of the samples, we extracted features to quantify different aspects of the dynamics of the transition to and from the identified fuzzy states.We used two types of features that have been frequently used in hard cluster-based studies (1,5,36,37) and developed a number of new features uniquely suited for use with fuzzy clustering that give new insight into state dynamics.For traditional hard cluster-based features, we assigned each sample to the cluster for which it had the highest probability and extracted the occupancy rate (OCR, i.e., the percent of time steps spent in a state by a participant) and number of state transitions (NST).For novel fuzzy cluster-based features, we extracted the (1) Kullback-Leibler divergence (KLD) across states and calculated a variety of descriptive statistics summarizing the KLD (5 features total).( 2) We also calculated Shannon entropy over time for each fuzzy state (1 feature per state, so 5 features total).(3) We calculated 3 descriptive statistics for probabilities for each state (15 features total).We calculated (4) the cumulative change between consecutive time points for each state (5 features total) and (5) a measure of the uniformity of probabilities within each state (5 features total).We sought to use these features to understand the effects of SZ upon the DMN dynamics, and statistical analyses using the features will be described in subsequent sections.Further details on each feature can be found in the Supplementary Materials.

Description of class-related statistical analysis
We wanted to determine whether the fuzzy states that we identified were related to SZ dynamics.As such, we performed a series of two-tailed t-tests comparing the dynamical and LP2D stability features belonging to SZs to those belonging to HCs.After performing the t-tests, we performed separate FDR corrections on the OCR features, correlation-based features, cumulative change features, and LP2D-based features to reduce the likelihood of false positive test results.

Description of LR-ENR classificationbased explainability analysis
Our statistical analysis indicated whether there were significant differences in the extracted dynamical features between HCs and SZs and gave insight into the effects of SZ upon DMN dynamics.
However, we also wanted to determine how helpful the extracted features could be for discriminating between HCs and SZs and gain insight into the relative importance of each feature for the classification.As such, we classified SZs and HCs using LR-ENR.LR-ENR is a frequently used (13) interpretable machine learning model.Its coefficients can be visualized for insight into the relative importance of each feature included in the classification.In our Scikit-Learn implementation (45), we feature-wise z-scored each feature and then trained separate LR-ENR classifiers for the traditional, non-uniformity + KLD + entropy, variance, mean, range, correlation, cumulative difference, and LP2D features.We used 10-fold nested cross-validation to optimize the ratio of L1 to L2 normalization and the regularization strength.We lastly calculated the mean and standard deviation of the area-undercurve (AUC) of the receiver-operating-curve, sensitivity (SENS, i.e., true positive rate), and specificity (SPEC, i.e., true negative rate) for the 10 folds corresponding to each set of dynamical features.Details on our parameter optimization approach can be found in the Supplementary Materials.

Description of symptom severity analysis
Lastly, while our earlier analyses gave insight into differences in DMN dynamics between SZs and HCs, we also wanted to determine whether the fuzzy states we uncovered could be used for insight into SZ symptom severity.As such, we performed ordinary least squares regression using age, gender, one-hot encoded site data and the negative PANSS score or age, gender, one-hot encoded site data, and the positive PANSS score as independent variables and the dynamical features as independent variables.We lastly applied FDR correction for the OCR features, entropy features, correlation-based features, cumulative change features, and LP2D-based features separately.

Results
In this section, we describe the 5 fuzzy states that we identified with our clustering approach.We then show the important dFNC features identified using our global explainability analyses.After characterizing the states and the dFNC features that differentiate them, we examine whether there are differences in the dynamical features and LP2D-based stability features between HCs and SZs and whether the features are related to symptom severity.

Identifying 5 fuzzy states of dFNC activity
As shown in Figure 2, we identified 5 distinct fuzzy states of dFNC activity.All states had highly positive intra-ACC dFNC.State 0 had overall highly positive dFNC values.It had moderately highly levels of positive ACC/PCN and ACC2/PCC dFNC.State 1 was characterized as having highly negative ACC/PCN and ACC/PCC dFNC, highly positive intra-network dFNC, and highly positive PCC/PCN dFNC.Relative to the other states, state 1 had more higher magnitude patterns of positive and negative dFNC.State 2 was the most poorly connected state, with relatively low values of positive dFNC.However, PCN1/PCN3, PCN1/PCC1, and intra-ACC had moderate positive values.Similar to state 1, state 3 had highly predominantly negative, though lower magnitude, ACC/ PCN dFNC.State 3 also had higher magnitude positive intra-PCN activity relative to states 2 and 4. It had very low intra-PCC and PCC/ACC dFNC comparable to state 2. Additionally, PCC/PCN and particularly PCC1/PCN were lower magnitude than every state but state 2. Lastly, state 4 had lower magnitude intra-PCN and ACC/PCN dFNC comparable to state 2. It had levels of intra-ACC and intra-PCC dFNC comparable to states 0 and 1.It had moderately positive PCC/ACC dFNC comparable to state 0 and positive PCC/PCN dFNC comparable to states 0 and 1.

Identifying key dFNC features differentiating each state and comparing explainability approaches
Although we visually compared the centroids of the clusters that we identified, the visual comparison did not address the relative importance of each dFNC feature to the identified states in a quantitative manner and did not consider the underlying distribution of the samples in relation to the identified fuzzy states.Figure 3 shows the results for our global explainability analyses.G2PC, total GP2D, and median GP2D all identified PCC2/PCN1 and PCC/PCN3 to be among the most important dFNC features with the other features being of varying importance.G2PC affected around 10% of samples.In contrast, mean LP2D captured the effects of perturbation on 100% of samples.Additionally, based on Kendall's rank correlation between the global feature importance of each method, there was significant high-level agreement in the relative importance (i.e., importance rank) of each dFNC feature across global approaches (e.g., the most important features for one method tended to be among the most important for the other methods).G2PC and total GP2D (p < 0.001) had a correlation coefficient of 0.65.G2PC and median GP2D (p < 0.001) had a correlation coefficient of 0.53, and total GP2D and median GP2D (p < 0.05) had a correlation coefficient of 0.39.Additionally, while G2PC and total GP2D tended to distribute importance widely, median GP2D provided a much sparser importance estimate.

Identifying disorder-related differences in dynamical and stability features
Figure 4 shows the results for our t-tests examining differences in dynamical and LP2D stability features between HCs and SZs, and Supplementary Figures 1-4 show boxplots of the dynamical features for each class.Many features, though not the KLD features, had statistically significant differences between groups.State 3 was of reoccurring importance, as HCs had significantly higher entropy, OCR, average probabilities, range in probabilities, and cumulative difference values in state 3. Interestingly, SZs had greater nonuniformity in their distribution of probabilities across states and greater state 0 and state 1 OCRs and average probabilities.SZs also had greater state 1 range in probabilities and state 0 cumulative differences.Additionally, many state correlations displayed significant differences between groups.HCs had higher state 0/1, 0/2, 1/2, and 1/4 correlations, and SZs had higher state 0/3 and 2/3 correlations.Lastly, SZs had higher levels of sensitivity to ACC1/ PCN1, ACC1/PCN3, and ACC2/PCN3 perturbation than HCs and lower levels of sensitivity to PCC1/PCN2 perturbation than HCs.
Table 2 shows the LR-ENR performance results, and Figure 5 shows the LR-ENR explainability results.As might be expected, given the t-test results, the Correlation, Average, Variance, NST + OCR, and Range feature models had the highest AUC values.The LP2D and Non-Uniformity + KLD + Entropy feature models had near chance-level performance.The Correlation model had the highest overall AUC.For SENS, the NST + OCR and Average

B A
Group-Level Statistical Comparison.(A, B) show the t-test results for the dynamical and LP2D stability features, respectively.In (A), the groups of features we extracted are arranged from left to right, and dark black lines separate each group of features.In (B), results are arranged in the form of a standard connectivity matrix for easier interpretation.Dashed lines separate domain pairs.Both panels are heatmaps showing the t-statistics for each feature.The color bar for (A) is above the heatmap, and the color bar for (B) is to the right of the heatmap.Black and grey asterisks indicate features with significant differences with and without FDR correction, respectively It should be noted that the t-tests were performed as HCs minus SZs.As such, a negative t-statistic indicates that SZs had higher values for a particular feature than HCs.Explainability Results.From left to right, the panels show the G2PC, LP2D, total GP2D, and median GP2D results.The color bar corresponding to each panel is positioned to its right, and the ICs associated with each dFNC feature are arranged on the x-and y-axes and are grouped based upon brain region (i.e., PCN, ACC, PCN).Total GP2D and median GP2D are scaled by their maximal value for easier interpretation.
feature models performed highest, with Variance feature model performing slightly lower.For SPEC, the Range, and Correlation, and LP2D models performed highest.The model coefficients also provided results highly consistent with the t-test results, exhibiting similar directionality and relative magnitude of differences between HCs and SZs.

Identifying relationship between symptom severity and dynamical and stability features
Interestingly, no stability features had significant relationships with or without FDR correction with symptom severity when accounting for age, gender, and study site, so we did not include a figure detailing the results for that analysis.However, accounting for age, gender, and study site, state 2 cumulative difference was positively correlated with the negative PANSS score (p < 0.05) prior to FDR correction.

Discussion
The goals of this study were (1) to present a framework for dFNC state analysis that accounts for the inherent variability of data within states and (2) provide additional insights into the effects of SZ upon DMN dynamics.We identified 5 fuzzy states of DMN activity, characterized the states using a novel explainability approach, identified effects of SZ upon DMN dynamics and state stability, and identified relationships between DMN dynamics and symptom severity.
Interestingly, the centroids for the 5 fuzzy states that we identified diverge greatly from those identified in previous studies (13).Nevertheless, there were some overall similarities in dFNC (13).Multiple states in (13) had highly positive intra-PCN dFNC and highly positive PCN/PCC dFNC.In contrast to the clusters that we identified, most of the centroids in (13) had negative intra-ACC and intra-PCC dFNC, while our clusters had positive values.Additionally, dFNC values seemed to be much more uniform in our clusters, whereas those in (13) had much more variability of dFNC values between independent components in each state.
While not all of our dynamic and stability features were able to obtain LR-ENR performance, comparable to the features extracted in (13), mean AUC for the NST + OCR, correlation, average probability, and variance features were comparable, which supports the utility of our pipeline for identifying key class-specific differences in brain dynamics.Importantly, HCs spent much more time in and had higher average probabilities for a state characterized by moderately negative ACC/PCN, moderately positive intra-PCN, and low magnitude intra-PCC and ACC/PCC dFNC (state 3).As might be expected given that they spent more time within that state, HCs also had much higher variability within that state (i.e., higher levels of variance, entropy, range in probabilities, and cumulative difference) than SZs in the state.In contrast, SZs spent more time in a state with highly negative ACC/PCN, moderately negative ACC/PCC dFNC, highly positive intra-network dFNC, and highly positive PCC/PCN dFNC (state 1).ACC/PCC activity was highly important to SZ dynamics, as they lacked stability to perturbation of ACC/PCC activity.Our findings of more negative or less positive ACC/PCC dFNC in SZ fit with those of multiple studies.Many studies have identified reduced or disrupted ACC/PCC dFNC in SZs (13,30,69,70), and these differences could be related to reduced ACC and PCC grey matter volume in SZs (71) or the disrupted anatomical distance function found in SZ (72).State 1 also had higher intra-PCN activity than state 3, which corroborates the findings of previous studies that identified higher intra-PCN activity in SZs than HCs (31,73).Previous studies have also suggested that this overactivation of the PCN could be a compensatory mechanism to reduce language comprehension deficits found in SZ (73) or that the PCN in SZ could be related to differentiating fantasy from reality (74) and that overactivation tends to be related to increased symptom severity.SZs also tended to be less stable to ACC/PCN perturbation, whereas HCs tended to be less stable than SZs to PCC/PCN perturbation.These findings paired with a visual comparison of state centroids indicates that strongly negative ACC/PCN activity and strongly positive PCC/ PCN activity were important for identifying SZs.Previous studies have found that the PCN strongly inhibits the ACC, which corroborates our finding of strongly negative ACC/PCN activity (75).Although many DMN studies agree with our finding that SZs spend more time in a state with more positive connectivity (31, 76) (i.e., PCN/PCN, PCC/PCN, and PCC/PCC in state 1 versus state 3), many previous studies have found SZ inter-domain connectivity to be less than HCs (10).This difference in results is likely attributable to the use of different regions within the DMN (10,77), the use of whole-brain analyses (77) that have been shown to obscure the activity of the DMN (13,23), or the use of multiple disorders that can disguise the effects of SZ (12).Though there are key differences, the insights that can be obtained from our state correlation features are somewhat related to the hidden Markov model state transition probabilities from (13).While the features in (13) indicate the probability of transitioning between a given state, our correlations features indicate how the similarity of brain dynamics of study participants to each fuzzy state vary over time and relate to multiple states simultaneously.SZs spent more time in states 0 and 1, which had high-magnitude dFNC with the exception of negative ACC/PCC and ACC/PCN in state 1.Some findings related to our correlation features provide helpful insight into SZ, though others are likely explained by aspects of our results that were previously explained.Importantly, similarity to state 0 and state 1 activity tended to be more strongly anticorrelated in SZs, and given that SZs spend a significant amount of time in each of those states, that would indicate that there was a large amount of alternation between the states (i.e., changes in ACC/PCC and ACC/PCN activity).Given that State 2 (low magnitude dFNC) was the dominant state, stronger anticorrelations between states 0/2 and states 1/2 in SZs could be explained by their stronger similarity to states 0 and 1 (i.e., being closer to a state of high magnitude dFNC is being farther from a state of low magnitude dFNC) or by increased transitions between states of high magnitude dFNC (i.e., states 0 and 1) and a state of low magnitude dFNC (state 2).Additionally, states 1 and 4 in SZs had stronger anticorrelation than in HCs, which could be explained by the greater state 1 occupancy of SZs, and SZs had more positive correlation between states 2 and 3, which could be explained by their relatively low similarity to each state 3. Our findings of stronger anticorrelation between states 0 and 3 and states 2 and 3 in HCs could also be explained to indicate that HCs transitioned more between states of moderate and high magnitude dFNC or states of low and moderate dFNC.Together, these findings could suggest that SZs transition more rapidly from low-magnitude dFNC states to high-magnitude states, while HCs transition more gradually from low-magnitude to moderatemagnitude to high-magnitude states.This could corroborate findings of previous studies that have identified the effects of SZ to involve more temporally localized changes in activity (22)(23)(24).
It is unfortunate that our KLD analysis did not uncover any significant differences between SZs and HCs.Previous studies have uncovered significant effects of SZ upon dynamics.However, our lack of findings related to KLD features does not indicate that those features would not provide useful insights within the context of other applications.While we extracted a number of metrics to quantify different aspects of the KLD, it is also possible that other metrics might be applied to obtain other insights into state dynamics.
Several studies have previously identified relationships between DMN dFNC features and symptom severity (12,13).We found that our use of the cumulative difference of fuzzy clustering probabilities uncovered a key relationship with negative symptom severity.Accounting for age, gender, and study site, the state 2 cumulative difference was related to negative symptom severity.This indicates that the more variation SZs had in their similarity to a state of lowmagnitude dFNC, the worse their symptoms may have been.Given that previous studies have related DMN deactivation with SZ symptom severity (78-80), it is necessary to be careful with this conclusion as the result was only found prior to FDR correction.
While our findings hold great significance for the domain of SZ analysis, our proposed explainable fuzzy clustering framework has broader implications.Our use of fuzzy clustering enables insight not only into inter-state dynamics but also into intra-state dynamics.Furthermore, our approach accounts for the inherent variation in how similar samples are to their assigned centroid and to samples assigned to other centroids.These challenges will be present in any dFNC clustering analysis.Our dynamical and stability features also represent key advances, providing new insights into disorderrelated activity.Additionally, as we showed by thresholding the state probabilities and calculating OCR and NST values, traditional k-means-based features can still be extracted with fuzzy c-means.That paired with the ease of implementation of fuzzy c-means using existing Python packages (52) and MATLAB (53) means that our approach could easily replace existing dFNC k-means clustering approaches in future studies.The explainability methods that we present also represent advances for the field, providing a quantitative estimate of the importance of each dFNC feature to the clustering.Relative to G2PC, our approach accounts for the effects of perturbation upon 100% of samples.Additionally, different metrics can be applied to the resulting KLD values to gain different insights into the effects of perturbation or to produce more or less sparse explanations (e.g., in total versus median KLD).

Limitations and future opportunities
Our study and approach have several limitations.Namely, we used a 40-second window when calculating dFNC, and while studies have shown that to be a reasonable window size (60), the window size can affect the dynamics and findings.Additionally, it is possible that our use of a feature extraction approach (clustering followed by extracting dynamical and stability features) applied to the whole dataset may have biased our LR-ENR performance results (81).However, our methods did enable us to directly compare the utility of our novel extracted features to those developed in (13).In future studies, we intend to perform a more robust series of analyses on how to best optimize the number of clusters for our application area.However, in this study, we used 5 clusters for easier comparison to (13).Fuzzy c-means has a couple limitations (1).Similar to k-means clustering, it is affected by outliers.However, there are variations of fuzzy c-means that are capable of addressing this problem (82).Additionally, fuzzy c-means can be more computationally intensive than k-means clustering (83).However, an approach similar to iSparse k-means could be easily adapted to work with fuzzy cmeans (84).Lastly, future research directions might include analyzing the reproducibility of our results across datasets similar to (13) or analyzing multiple disorders in a single analysis like (4).

Conclusion
The analysis of rs-fMRI dFNC data using hard clustering methods to identify states that summarize brain dynamics is a common analysis approach that has provided insights into many neurological and neuropsychiatric disorders.However, the use of hard clustering approaches (e.g., k-means clustering) can obscure key information related to how similar samples are to their respective centroids or to samples assigned to other cluster centroids.As such, the use of hard clustering could obscure disorder-relevant dynamics.In this study, we present a novel explainable fuzzy clustering framework.We present 7 new types of dynamical features and sample stability-based features that provide unique insights into brain dynamics while also demonstrating how traditional dynamical features used in hard clustering analyses can also be included in our analysis.Lastly, we present two novel explainability approaches that help characterize the fuzzy states identified using our clustering approach.We demonstrate our framework within the context of SZ DMN analysis, identifying aberrant dynamics in SZs and uncovering relationships between SZ symptom severity and a state of reduced DMN correlation.Our framework provides greater insight into disease dynamics than traditional hard clustering approaches.Furthermore, it can be implemented with an ease comparable to the standard k-means clustering approach using existing code packages.As such, it represents an ideal method for future widespread use in rs-fMRI dFNC analysis and could lead to an improved understanding of the effects of many neurological and neuropsychological disorders upon brain dynamics.

FIGURE 1
FIGURE 1 Overview of Methods.Red dots indicate each step of the methods.(1) We recorded rs-fMRI data from HCs and SZs.(2) We extracted dFNC data.(3) We performed fuzzy c-means clustering, identifying 5 states.(4) We applied several explainability approaches for insight into the dFNC features characterizing the states.(5,6)We extracted dynamical and stability features.(7)We performed t-tests and trained interpretable machine learning models for insight into the features that differed between HCs and SZs.(8)We performed linear regression analyses controlling for age and gender to identify relationships between symptom severity and dynamical and stability features.

Fuzzy
States and Example State Trajectories.(A) shows the centroids for each of the 5 fuzzy states that we identified.Each subpanel shares the same color bar to the right of the panel for Fuzzy State 4. The ICs associated with each dFNC feature are arranged on the x-and y-axes and are grouped based upon brain region (i.e., PCN, ACC, PCN).(B) shows example state trajectories for an SZ participant (left) and HC participant (right).The y-axis indicates the probability of belonging to each state, and the different colors of lines correspond to the state numbers shown in the legends.Note the variation in probability of belonging to each cluster.

FIGURE 5 LR
FIGURE 5 LR-ENR Explainability Results.Each panels shows the LR-ENR coefficient values for the groups of features indicated in their respective titles.Specific feature names are indicated on the x-axis of each panel.Coefficient values are shown on the y-axis.Higher magnitude coefficients indicate greater importance to their respective model.Because a label of 0 was used for SZs and a label of 1 was used for HCs, a negative coefficient value indicates that an increase in that feature corresponded to an increase in likelihood of belonging to the SZ class.

TABLE 1
Description of Study Participants.