We use data from a large, multi-site eyes-open resting state fMRI study with approximately equal numbers of SZs and HCs (n = 311, nSZ = 150). Imaging data for six of the seven sites were collected on a 3T Siemens Tim Trio System and on a 3T General Electric Discovery MR750 scanner at one site. Resting state fMRI scans were acquired using a standard gradient-echo echo planar imaging paradigm: FOV of 220 × 220 mm (64 × 64 matrix), TR = 2 s, TE = 30 ms, FA = 770, 162 volumes, and 32 sequential ascending axial slices of 4 mm thickness and 1 mm skip. Subjects had their eyes closed during the resting state scan. The data were preprocessed with a standard pipeline and decomposed with group-independent component analysis (GICA) into 100 group-level functional network spatial maps (SMs) with corresponding subject-specific TCs. Through a combination of automated and manual pruning, N = 47 functionally identifiable resting state networks (RSNs) were retained (Figure 1). The remaining network TCs were detrended, despiked, and orthogonalized with respect to estimated subject motion parameters. Subject-specific SMs and TCs were obtained from the group level SMs via spatiotemporal regression. The TCs were detrended, despiked, and subjected to some additional post-processing steps [details in Damaraju et al. (2014)]. The retained RSNs fall into seven functional domains: subcortical (SC) (five networks), auditory (AUD) (two networks), visual (VIS) (10 networks), sensorimotor (SM) (seven networks), cognitive control (CC) (13 networks), default mode network (DMN) (eight networks), and cerebellum (CB) (two networks). In all figures showing functional network connectivity matrices in heatmap form, the functional domains appear in the indicated order along both axes. All subjects in the study signed informed consent forms. Symptom scores for patients were determined using the Positive and Negative Syndrome Scale (PANSS) (Kay et al., 1987). The focus in this paper will be on the six symptoms (delusions, grandiosity, hallucinations, suspiciousness/persecution, preoccupation, and unusual thought) that load on the “positive factor” in a heavily used factor analysis study of PANSS symptoms (Lindenmayer et al., 1995).

Dynamic functional connectivity (dFNC) between RSN TCs was estimated using sliding window correlations. Following the protocols from published studies on dynamic connectivity (Damaraju et al., 2014), we employed a tapered rectangular window length of 22 TRs (44 s), advanced 1 TR at each step, and computed pairwise correlations between RSN TCs within these windows. After dropping the first three and final TRs, this procedure yields a 47(47−1)/2 = 1,081-dimensional dFNC measure on each of 136 windows of length 22TRs for each subject. Clustering this collection of time-resolved connectivity observations using MATLAB’s implementation of k-means clustering (Euclidean distance, 2,000 iterations, 250 repetitions, five clusters chosen according the elbow criterion) produces five non-varying cluster centroids, often referred to as “dynamic states” or dFNC states (Figure 2) reported in previous studies (Damaraju et al., 2014; Rashid et al., 2014, 2019; Yu et al., 2015; Miller et al., 2016b,c; Yaesoubi et al., 2017; Espinoza et al., 2019b; Mennigen et al., 2019; Weber et al., 2020; Faghiri et al., 2021; Fu et al., 2021). We mention them because they are referenced later in the results section of this paper. To distinguish time-blind non-varying states from EVOdFNCs, we will refer to them as non-varying “snapshot” dFNC states (SNAPdFNC). There are two modularized SNAPdFNC states dominated by strong positive auditory-visual-sensorimotor (AVSM) connectivity: one (State 1) has strong negative connectivity between the default mode network (DMN) and the AVSM networks; the other (State 2) presents weak DMN-to-AVSM connectivity. HCs spend more time in both of these AVSM-dominant modularized states than the patients do. A third modularized state (State 3) is defined by a contrast between strong negative connectivity between the DMN and the rest of the brain, with diffuse positive connectivity in all other blocks of the connectome. We will use the shorthand, DMNneg, in future references to this particular modularized SNAPdFNC state, which is more occupied by patients than controls. The other two SNAPdFNC states (States 4 and 5) show much less modularity: one (State 4) is diffusely hyperconnected and more occupied by controls; the other (State 5) is diffusely disconnected and more occupied by patients (Damaraju et al., 2014).

We apply a MATLAB implementation (McInnes et al., 2018; Meehan et al., 2020) of UMAP to embed all 1, 081-dimensional dFNCs into the plane. Although both UMAP and t-distributed stochastic neighbor embedding (t-SNE) preserve high-dimensional local structure in the lower-dimensional embedding, UMAP holds onto more global structure than t-SNE, exhibits greater stability across runs, and is also significantly more efficient (see Supplementary Files for experiments clarifying the differences between UMAP and t-SNE and the role of parameters in stabilizing the embedding). Key user–chosen parameters in UMAP are the number of nearest-neighbors (n_neighbors) over which the high-dimensional uniformized neighborhoods are defined, and the minimum distance (min_dist) which parameterizes proximity in the low-dimensional embedding. The stability and “reasonableness” of UMAP embeddings, especially with respect to preservation of global structure, can be sensitive to initialization (Kobak and Linderman, 2021), although our specific data were relatively stable under re-initialization of the algorithm and we further stabilized the final embedding by averaging over multiple runs (Supplementary Figure 4). Because dFNCs are computed on sliding windows that advance one TR at a time, they exhibit considerable temporal smoothness in their native high-dimensional space (mean elementwise squared distance between successive dFNC measures is less than 0.0006). The high-dimensional trajectories are also somewhat “densely packed,” i.e., for any fixed dFNC observation, the average mean squared elementwise distance between that observation and its 10,000 nearest neighbors from different subjects is less than 0.01. The ideal embedding for our purposes would preserve intra-subject trajectory continuity and inter-subject trajectory proximity. This within-subject smoothness is an intrinsic feature of the actual dynamics, which we experimentally optimized UMAP’s two main parameters (n_neighbors = 25; mindist = 0.75) toward conserving in the planar embedding (Supplementary Figure 3). It is worth quickly noting that linear dimension-reduction methods such as PCA or ICA produced diffuse two-dimensional (2D) clouds lacking intra-subject temporal continuity (Figure 3), whereas t-SNE performed more than 200-fold more slowly than UMAP on these data and, thus, was not practical for achieving a planar embedding of the entire dataset.

Because the UMAP embeddings are not fully deterministic, we run the algorithm R = 25 times on the input dFNC data 𝒟 and use the average ℰ^(𝒟)={(1R∑r=1Rxr(t),1R∑r=1Ryr(t)),s=1,2,..,n,t=1,2,…,T}, of the R individual 2D embeddings ℰr(𝒟)={(xrs(t),yrs(t)),s=1,2,…,n,t=1,2,…,T}. This averaging procedure stabilizes the embedding but also de-densifies the group level geometry (Supplementary Figure 4), reducing its qualitative fidelity to the native high-dimensional setting.

Because of our choice of UMAP parameters within intervals that preserve intra-subject trajectory continuity, the majority of each subject’s high-dimensional dFNC trajectory Γ(t) = (v₁(t), v₂(t), … , v₁₀₈₁ (t)) embeds, under the r^th run of UMAP, into an approximately continuous segment γ_r(t) = (x_r (t), y_r (t)) in 2D (Figure 4). Continuity is preserved under summation, so the average γ^(t)=(x^(t),y^(t))=(1R∑r=1Rxr(t),1R∑r=1Ryr(t))∈ℰ^(𝒟) will also be continuous. We collect all continuous trajectory sub-segments (CTSs) α^j(γ^)={γ^(t),t∈[tj,tj+τ]},j=1,2,…,T-τ along γ^(t) of temporal duration τ = 44 (twice the window length used to estimate the dFNC observations in 𝒟). Each α^j(γ^) can be approximated by a line L_j, yielding a reduced characterization (x^j,y^j,m^j,ℓ^j) of αj^(γ^) in terms of its length ℓ^j=maxt,t′∈[tj,tj+τ]||γ^(t)-γ^(t′)||2 , its geometric midpoint (x^j,y^j), and the slope m^j of its linearization L_j (Figure 5). The linearized trajectory segment (LTS) triples: (x^j,y^j,m^j) are then clustered (Figure 10) with k-means [MATLAB implementation, squared Euclidean distance, 2,000 iterations, 250 repetitions, slope bounded between positive and negative max ({|x|, |y|}) to keep arbitrarily large near-vertical slopes from dominating the spatial (x, y) components of the clustering] where the number of clusters, K = 10, was chosen according to the elbow criterion. From each of the K = 10 LTS cluster centroids (x^j,y^j,m^j), we induce a 2D line segment of length d_i (equal to the mean length of all LTSs in that cluster) and slope m^i centered at (x^i,y^i) (Figures 4, 5). These segments are, roughly speaking, gradients of the CTSs, which are, in turn, averages of continuous embedded segments of the high-dimensional dFNC dynamics. The collection of segments induced by LTS cluster centroids will be called linear trajectory exemplars or just exemplars. Although UMAP is a non-linear dimension reduction algorithm, implying, among other things, that line segments in the planar embedding are not necessarily linear segments in the source domain – the line segments moving along dominant embedded curve directions, capturing local spatial trends from t to t + 1, will have non-linear inverse images that similarly reflect local trends in the high-dimensional dynamics as locality is highly preserved under UMAP.

Because UMAP is not straightforwardly invertible, the 2D linear trajectory exemplars above cannot easily be mapped back directly into the native dimension of the dFNC space. To obtain the high-dimensional dFNC trajectory segment ϒ of integer temporal duration τ corresponding to a given 2D linear trajectory exemplar ϒ, (i.e., the “data driven inverse” of υ), we average the high-dimensional observations corresponding to the n = 25 nearest 2D neighbors of each of τ evenly spaced points along υ (Figure 7). The points utilized are the 25 spatial nearest embedded neighbors to each of the τ evenly spaced points along v. For the purpose of inverting planar points not associated with samples, the slope of the linear trajectory segments centered at these points are not taken into account because the embedding itself is working in spatial terms exclusively and its “data-driven” inverse should work in the same terms as the embedding. The number of 2D neighbors used for this inversion of the UMAP embedding, n = 25, was chosen to match the number of nearest neighbors parameter that we employed for running the UMAP algorithm. The operation thus described effectively “lifts” a localized 2D linear trajectory exemplar back into high-dimensional dFNC space. The 2D linear trajectory exemplars are, by construction, concentrated in more densely occupied parts of the plane, and the continuity preserving parameterization of UMAP encourages the high-dimensional data-driven inverse of each 2D linear trajectory exemplar to exhibit naturalistic smoothness (Figure 11). These high-dimensional inverse images of the 2D linear trajectory exemplars are each multiframe representations of evolving functional network connectivity (EVOdFNCs).

The high-dimensional EVOdFNCs obtained from 2D exemplar linear trajectories are more properly viewed as representations of dominant directional trends than as a hard segmentation of the observed dynamics. We, therefore, employed them as basis objects through which to parameterize observed high-dimensional evolving connectivity dynamics. This is done by correlating successive time points (i.e., t = 1, 2, … , τ)of each length-τ window w_i in a subject’s observed dFNC sequence with the corresponding set of t^th time points from the K length-τ EVOdFNCs (Figure 8). For each w_i and each t ∈ {1, 2, … , τ}, we identify the EVOdFNC k ∈ {1, 2, … , K} whose t^th time point has the highest correlation with the t^th time point of observed dFNC window w_i, yielding a length-τ vector ρ = (ϱ₁, ϱ₂, … , ϱ_τ), ϱ_j ∈ {1, 2, … , K}. The weight λ_k(w_i) ∈ [0, 1] on each k ∈ {1, 2, … , K} is the proportion of time points ⌊{j∈{1,2,…,τ}:ϱj=k}⌋τ in w_i that most resemble (are most correlated with) the same time point in EVOdFNC k. This approach both normalizes the weights in [0, 1] and accommodates the possibility that within any given time window of length τ, a subject’s dFNC sequence might contain subintervals that are best represented by the corresponding subinterval of different EVOdFNCs. From this, we obtain a K-variate time series of EVOdFNC weights Λ={λk(wi)}k,i=1K,(T-τ) for each subject (Figure 8, bottom left), capturing the relative RI of each EVOdFNCs to successive length-τ intervals of the subject’s observed high-dimensional dFNC trajectory. To be clear, the weight attached to a given EVOdFNC k, for a given length-τ window w of observed dFNCs, will be λ ∈ [0, 1], denoting the proportion of the τ dFNC observations in w that are more correlated to the corresponding frame in EVOdFNC k than they are to the corresponding frames in other EVOdFNCs. This assigns weight to each EVOdFNC according to the level at which its ordered sequences of frames represent (correlationally) represent ordered subsequences within the window of dFNC observations.

Because EVOdFNCs are inverse images of 2D segments in a planar embedding of very high-dimensional data, using weighted combinations of the EVOdFNCs (as opposed, e.g., to a binary “occupancy” approach) is arguably a better strategy for capturing the evolving variability of high-dimensional dFNC observations. Toward this end, we cluster (across all windows and subjects) the time-indexed weight vectors from the multivariate EVOdFNC time series into M = 10 clusters (elbow criterion, MATLAB’s k-means implementation, squared Euclidean metric, 2,000 iterations, and 500 repetitions). We then induce M meta-EVOdFNCs as centroid-weighted sums of the basis EVOdFNCs (Figure 9). Subject-specific occupancy rates of each meta-EVOdFNC are computed as the fraction of the subject’s time-indexed weight vectors that belong to the cluster whose centroid defines the meta-EVOdFNC.

All reported SZ effects are from a multiple regression on gender, age, head motion (mean frame displacement), and SZ. Reported positive PANSS symptom score effects are from a multiple regression on the six positive PANSS symptoms (delusions, grandiosity, hallucinations, suspiciousness/persecution, preoccupation, and unusual thought) from the Lindenmayer five-factor PANSS model (Lindenmayer et al., 1995), in addition to age, gender, and head motion. The purpose of this model is to identify the effect of each positive symptom corrected for the contributions of the others, in a manner indifferent to, i.e., summing over, the various profiles of negative and general symptoms that subjects might manifest. Whereas negative symptoms in the SZ subjects analyzed here are highly intercorrelated [mean correlation between negative symptom pairs = 0.31 (std. dev = 0.12)], positive symptoms are not [mean correlation between positive symptom pairs = 0.12 (std. dev = 0.13)]. Schizophrenia effects are only reported when significant at the p < 0.05 level after correction for multiple comparisons. All displayed PANSS symptom effects are significant at the p < 0.05 level and remain significant at this level after correction for multiple comparisons where specifically indicted.

We found widespread schizophrenia effects on the RI of each of the 10 EVOdFNCs in subject data (Figure 11). The evolving motifs fluidly move through transient states of connectivity that resemble familiar formations obtained from the basic time-blind clustering into single transiently realized connectivity patterns (Figure 2). Consistent with published results (Damaraju et al., 2014; Miller et al., 2016b; Yaesoubi et al., 2017; Espinoza et al., 2019b; Rashid et al., 2019) on occupancy rates of time-blind SNAPdFNC states, we found the following: strongly modularized and hyperconnected patterns feature more prominently in EVOdFNCs with greater RI in controls (1, 5, 7, and 8) and in EVOdFNCs whose RI in controls is not statistically distinguishable from that in patients (4, 6, and 9); weak connectivity and modularized negative DMN-to-other (DMNneg) patterns (2 and 3) feature more prominently in EVOdFNCs with significantly higher RI in patients. A novel modularized pattern of functional organization, not seen in time-blind SNAPdFNC states, appears in EVOdFNC 10, which features a persistent stretch of strong modularized negative SM-VIS/CC/DMN connectivity. EVOdFNC 10 has significantly greater RI in patients, and its distinctive modularity eventually dissolves into unstructured weak connectivity. EVOdFNC 7, with greater RI in controls, transitions from a weakly connected negative DMN-to-other structure (more characteristic of patients in time-blind analysis, e.g., SNAPdFNC State 3) into another novel pattern that is strongly modularized with positive DMN-VIS/SM connectivity. EVOdFNC 2 (+SZ, RI) and EVOdFNC 5 (+HC, RI) show roles for disconnected periods leading into and out of modularized structures that are more characteristic of SZ and HC, respectively. EVOdFNC 2 (+SZ, RI) also exhibits a transition from something more like SNAPdFNC State 1 (+HC, OCR) to SNAPdFNC State 3 (+SZ, OCR), presenting a path from strong AVSM-focused modularity (+HC, SNAPdFNC OCR) to the DMNneg configuration (+SZ, SNAPdFNC OCR) to weak connectivity (+SZ, SNAPdFNC OCR).

EVOdFNCs 1 and 8 (both +HC, RI) show the role of disconnected periods in transitioning between AVSM-centered modularity and hyperconnectivity in healthy people: in EVOdFNC 1, we see modularity dissolving into diffuse dysconnectivity, which then uniformly inflates toward diffuse hyperconnectivity, whereas, in EVOdFNC 8, diffuse hyperconnectivity gradually erodes toward greater dysconnectivity away from the AVSM core and then starts manifesting negative DMN/CC-AVSM connectivity from a substrate of dysconnectivity. Although EVOdFNCs 2 and 3 both have greater RI in patients and both contain patterns similar to SNAPdFNC State 3 (+SZ, OCR), it is the weaker, less dynamically varying EVOdFNC 3 whose RI is significantly further elevated in more delusional patients (Figure 12). The RI of EVOdFNCs 4 and 9 is not significantly affected by subject diagnosis, but they reveal patterns of functional organization that do not appear in time-blind dFNC clustering. Specifically, EVOdFNC 4 evolves from strong AVSM-centered modularity (+HC, SNAPdFNC OCR) through a DMNneg configuration (+SZ, SNAPdFNC OCR) connectivity into a modularized pattern, in which intra-AVSM and intra-CC connectivity is at approximate parity. EVOdFNC 9 passes through a stage in which intra-CC connectivity is even stronger than intra-AVSM connectivity. Transient periods of whole-brain connectivity where connectivity strength within the CC domain equals or exceeds that within AVSM domain block do not appear in time-blind dFNC clustering of these data. EVOdFNCs 4 and 9 both have strong population-wide RI (Supplementary Figure 2) indicating that these patterns play a meaningful role in dynamic connectivity despite not appearing in SNAPdFNC clusters centroids.

Three of the 10 EVOdFNCs start with intervals characterized by modularized DMNneg patterning, i.e., strong negative connectivity between DMN and all other brain areas with positive connectivity everywhere else. The SNAPdFNC state with this pattern is more occupied by SZs than HCs. Two of the EVOdFNCs (2 and 3) have higher RI in SZs, one (EVOdFNC 7) in HCs. Two have significant positive PANSS symptom effects (Figure 12): One of these (EVOdFNC 3) has higher RI in SZs, and the other (EVOdFNC 7) has higher RI in HCs. Although all three EVOdFNCs pass through a DMNneg type of pattern, each embeds this pattern in a different dynamic context (Figure 13). EVOdFNC 2 presents a high-contrast DMNneg pattern dynamically dissolving into diffuse dysconnectivity. Its temporal mean (bottom left) resembles that of EVOdFNC 3, which features persistent, low-contrast, unvarying DMNneg. Both EVOdFNCs 2 and 3 have higher RI in SZs, but only the lower contrast, less dynamic version, EVOdFNC 3, has significant further relationships to positive symptoms within the patient population (Figure 12). The temporal averages of EVOdFNCs 2 and 3 (Figure 13, bottom row) are very similar, suggesting that the differential sensitivity of EVOdFNCs 2 and 3 to positive symptoms resides in their temporal patterning. Finally, although, in time-blind SNAPdFNC analysis, the DMNneg pattern is significantly more occupied by SZs, we found EVOdFNC 7 (Figure 13, top right) featuring DMNneg, leading into AVSM-dominant, lightly modularized hyperconnectivity to have significantly higher RI in HCs and also to be sensitive to positive symptoms in SZs (Figure 12).

The 2D linear exemplars that lift to high-dimensional EVOdFNCs are relatively sparse in the planar embedding. Many local linearizations of embedded subject dFNC trajectories are not geometrically proximal to an exemplar or are not directionally aligned with the nearest exemplar. Moreover, the local linearizations are only approximations to the embedded curves, which, in turn, are imperfect representations of the high-dimensional ground truth. This suggests that a multivariate characterization of the dynamics in terms of replicable patterns, i.e., clusters, of concurrent multivariate EVOdFNC RI would induce meta-EVOdFNCs (Figure 9) that capture more of the high-dimensional data variability than the RI of individual EVOdFNCs. In practice, for the data evaluated in this study, the resulting weight-vector centroids (Figure 14, center right) were each highly concentrated on one EVOdFNC, yielding meta-EVOdFNCs that strongly resemble the dominantly weighted EVOdFNC in the corresponding weight-vector cluster centroid (Figures 11, 14). However, although focused on one dominant EVOdFNC, the meta-dFNCs include content from all EVOdFNCs, which allows for more flexible representation of the data when warranted. Moreover, the occupancy rate forces a harder segmentation than RI, yielding slightly different results. Delusions, for example, show positive effects on the RI of EVOdFNC 3 and on the occupancy rate of meta-EVOdFNC 2 (Figure 15), which is heavily weighted toward EVOdFNC 3. Delusions, however, also exhibit negative effects on the occupancy rate of meta-EVOdFNC 4, which is heavily weighted on EVOdFNC 7 whose RI is not significantly affected by Delusions.

It is also possible to explore diagnosis and symptom effects directly in the 2D embedding, considering the linear exemplar cluster membership of each local linearization to an embedded subject trajectory (Figure 10). This approach “trusts” the embedding to distill some important relationships in the data while unavoidably distorting others. Surprisingly, there are no significant SZ effects on exemplar cluster membership of 2D local linearizations: nothing significant at the p < 0.05 level after correction for multiple comparisons, nor any raw uncorrected p-values less than 0.05. However, in the case of positive PANSS symptoms, we found a number of significant relationships between positive PANSS scores and linear exemplar cluster membership (Figure 16), one of which, the negative effect of Hallucinations on the linear exemplar that lifts to EVOdFNC 8, is significant at the p < 0.05 level after correction for multiple comparisons.

High-dimensional EVOdFNCs are inverses or lifts of 2D linear exemplars that summarize the planar embedding of group level high-dimensional dynamics. The exemplar membership of successive linear approximations to the embedded curves offers a “first pass” view of the alignment of subject connectivity dynamics with directional trends in the embedded data. As an embedding-intrinsic measure, exemplar cluster membership of 2D linear approximations is less concretely empirical than high-dimensional EVOdFNC RI. For positive SZ symptoms, this embedding-intrinsic point of view strengthens the uncorrected significance of, e.g., the Hallucinations effects on EVOdFNCs 3 and 8 while identifying several new relationships between positive symptoms and dimensionally reduced dFNC dynamics (Figure 16).