Aberrant Recovery of Timescale-Aligned Amplitude Balance Links to Symptoms and Cognition in Schizophrenia

9 Understanding how metabolic energy is distributed across brain networks is essential for 10 elucidating healthy brain function and neurological disorders. Research has established 11 the link between blood flow changes and glucose metabolic processes that fuel neural 12 activity. Here, we introduce a novel framework based on the normalized dynamic time 13 warping algorithm robust to neural temporal variability, enabling reliable insights into 14 metabolic energy demands using functional magnetic resonance imaging data. Our 15 findings indicate that healthy brains maintain balanced energy distribution, whereas 16 imbalances are more pronounced in schizophrenia with links to both positive and negative 17 symptoms, particularly during rapid neural processes. Additionally, we identified a 18 dynamic state that supports the brain criticality theory and is associated with higher-order 19 cognitive abilities, demonstrating our framework’s functional and clinical relevance. By 20 linking metabolic energy distribution to neural dynamics, this framework provides a novel 21 way to estimate and quantify the brain’s maintenance of functional balance in a broadly 22 applicable manner for studying brain health and disorders. 23

brain energetics, dynamic time warping , high -frequency BOLD, brain 24 criticality, schizophrenia 25 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 2 1. INTRODUCTION. 26 The human brain, despite constituting about 2% of the body’s mass, consumes roughly 27 20% of its total energy 1—a cost driven by intensive neuronal processing and synaptic 28 activity essential for complex cognitive functions 2. Functional magnetic resonance 29 imaging (fMRI) leverages this metabolic demand by measuring the blood -oxygenation-30 level-dependent (BOLD) signal, which indirectly reflects neural activity through changes 31 in cerebral blood flow and oxygenation2,3. The detectability of the BOLD signal arises from 32 the retention of oxygen in the increased blood flow, underscoring glucose as the primary 33 energy source for neurons 2. Prior studies have shown that changes in blood flow are 34 accompanied by comparable changes in glucose metabolism 4,5, making BOLD signal 35 amplitude a reliable proxy for energy demand4. 36 In contrast, techniques like positron emission tomography (PET) and magnetic resonance 37 spectroscopy provide more direct insights into metabolism, but they lack the temporal 38 resolution to capture rapid neural dynamics . Moreover, BOLD signal interpretation is 39 complicated by variability in hemodynamics and neurovascular coupling (NVC), which 40 introduce temporal lags and distortions6, as well as the varying temporal scales of neural 41 processing7. Commonly used measures such as the amplitude of low -frequency 42 fluctuations (ALFF)8, fractional ALFF9, and dynamic ALFF10 quantify BOLD signal power 43 but do not account for temporal deformations, which may mask amplitude variations that 44 are critical to assessing neural energy consumption . Similarly, prevalent dynamic 45 functional connectivity measures—such as sliding window Pearson correlation11, Hidden 46 Markov Models 12, phase synchronization 13, wavelet coherence 14, and lagged cross -47 correlation15— primarily focus on synchronization and are amplitude scale invariant due 48 to normalization procedures. Even many advanced deep learning models 16,17, including 49 spatiotemporal models18 struggle to fully capture these temporal dynamics. 50 There is therefore a pressing need for methodologies that can address both temporal 51 misalignments and amplitude variations in BOLD signals to accurately reflect energy 52 demand mismatches between neural networks. To this end, we propose leveraging 53 dynamic time warping (DTW) to realign BOLD signal pairs, effectively mitigating temporal 54 misalignments caused by hemodynamic variability and variable temporal scales. DTW is 55 a time-series alignment technique that non-linearly warps the temporal dimension to align 56 signals19, allowing us to focus on the amplitude differences between the aligned signals. 57 By analyzing BOLD signals in a way that compensates for temporal variability and 58 hemodynamic delays, these amplitude differences indirectly estimate metabolic energy 59 mismatch between neural networks. Specifically, DTW sums the aligned amplitude 60 differences, and our recently proposed normalized DTW (nDTW)20 computes the average 61 of these differences, aligning with traditional signal energy and power computations 62 respectively, as demonstrated in our study. 63 DTW has been introduced in fMRI studies as an alternative measure of functional 64 connectivity21, demonstrating benefits such as increased sensitivity to motor brain 65 function22, robustness to noise 21, higher test -retest reliability 21,22, and global signal 66 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 3 regression robustness21,22 compared to conventional correlation methods. DTW has also 67 been effective in detecting group differences, including sex differences 21 and clinical 68 disorders like schizophrenia20 and autism spectrum disorders22. In our previous work, we 69 expanded DTW’s application by developing a novel warp elasticity measure, which 70 leverages the alignment function from DTW to capture relative BOLD signal activation 71 speeds (stretching and shrinking temporal dynamics) between networks 23,24. These 72 findings underscore the value of DTW in neuroimaging analyses. 73 In this study, we build on these findings by demonstrating that DTW -aligned amplitude 74 differences serve as an effective metric for capturing activation amplitude discrepancies 75 in the presence of temporal deformations. We incorporate a parameterization of the DTW 76 cost function with a gamma parameter that emphasizes different scales of amplitude 77 disparity—where higher gamma values accentuate larger amplitude differences and 78 lower values highlight smaller ones 25—to create a flexible framework for assessing 79 metabolic energy demand mismatches across varying amplitude scales in the BOLD 80 signal. Our focus is on resting -state fMRI (rs -fMRI) data, recognizing that significant 81 energy is consumed even at rest and that resting -state activity is linked to metabolic 82 processes26. 83 The proposed two-cell model presents a unifying hypothesis that links schizophrenia to 84 dysregulated glucose metabolism in astrocytes and neurons, suggesting that cognitive 85 and psychiatric symptoms arise through a two -hit mechanism 27. This hypothesis is 86 supported by numerous studies connecting schizophrenia to disruptions in metabolic 87 energy28,29. Based on this and recognizing that blood flow—carrying glucose—increases 88 in response to higher energy demand, we investigate imbalances in energy demand 89 between neural networks in individuals with schizophrenia compared to healthy controls 90 using blood flow fMRI data. Consistent with previous findings, we focus on nDTW due to 91 its increased sensitivity to group differences 20. We hypothesize that energy mismatches 92 between networks are more pronounced in broadband signals; therefore, we extend our 93 analysis beyond the classic low -frequency BOLD signals to include higher frequency 94 bands—specifically slow -5 (0.01 –0.027 Hz), slow -4 (0.027 –0.073 Hz), and slow -3 95 (0.073–0.198 Hz)—which have been shown to exhibit functional relevance30. We refer to 96 the broader frequency ranges of 0.01 –0.198 Hz and 0.01 –0.15 Hz as F1 and F2, 97 respectively. In addition, we develop a time -resolved nDTW measure to track how brain 98 networks converge or diverge in their metabolic energy demands over time during rest. 99 By capturing these nuances, our approach offers a novel framework to understand neural 100 network specialization and energy demand dynamics, with the potential to advance the 101 characterization of neuropsychiatric disorders linked to energy metabolism abnormalities. 102 This methodology not only enhances our ability to detect subtle energy imbalances within 103 and between neural networks but also provides new insights into the pathophysiology of 104 conditions such as schizophrenia with the potential to inform more targeted therapeutic 105 interventions. 106 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 4 2. RESULTS 107 2.1. Simulation of DTW sensitivity non-stationarity and signal bandwidth 108 We evaluated DTW as a potential measure of metabolic energy demand mismatch 109 between brain networks through simulations assessing its ability to detect amplitude 110 disparities and its invariance to phase shifts. Sinusoids with varying amplitude and phase 111 modulations were generated, and the performance of DTW, nDTW, and correlation was 112 compared across different 𝛾 values and signal bandwidths (Supplementary Material). 113 DTW aligned original and temporally warped signals, demonstrating its capability to 114 capture temporal deformations (Fig. 1a). Both DTW and nDTW remained invariant to non-115 stationary phase perturbations across increasing modulation indexes . In contrast, 116 correlation decreased with higher phase modulation (Fig. 1b). In contrast, under non -117 stationary amplitude perturbations, correlation remained stable, while DTW and nDTW 118 detected increasing amplitude disparities as modulation indexes rose (Fig. 1b). 119 The sensitivity of DTW and nDTW to amplitude scales was further examined by varying 120 the 𝛾 parameter. Simulations involving amplitude non -stationarity across several 121 modulation indexes revealed that smaller 𝛾 values resulted in higher rates of change in 122 DTW and nDTW similarity scores at lower modulation indexes, indicating enhanced 123 detection of minor amplitude disparities (Fig. 1c). Conversely, larger 𝛾 values exhibited 124 higher rates of change at higher modulation indexes, reflecting improved tracking of 125 substantial amplitude disparities (Fig. 1c). Additionally, while DTW generally maintains 126 phase invariance, nDTW demonstrated slight sensitivity to phase differences when 𝛾 ≥ 1 127 (Fig. 1c). This suggests that nDTW is better capable of responding to phase variations 128 alongside amplitude tracking, offering a more nuanced assessment when absolute phase 129 invariance is not essential. 130 In bandwidth analyses, DTW and nDTW sensitivity to randomness increased with broader 131 full frequency spectrum , approaching the Nyquist frequency, unlike correlation, which 132 remained unaffected (Fig. 1d). This indicates that amplitude differences become more 133 prominent in broader full frequency bands. 134 Figure 1e illustrates nDTW ( 𝛾 = 1) in capturing average aligned amplitude differences 135 between brain networks. Lower nDTW values indicate convergence, potentially reflecting 136 balanced metabolic energy demands, while higher values signify divergence, potentially 137 highlighting greater energy demand discrepancies between networks (Fig. 1e). 138 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 5 139 Fig. 1: Simulation to analyze non-stationarity and signal bandwidth a. Illustration of DTW’s ability to extract temporal deformations. A signal is warped by applying a modulation function (warp function) to its time index, and DTW successfully aligns the original and warped signals, extracting the induced temporal deformatio n. b. Behavior of DTW , nDTW and correlation under amplitude and phase non - stationarity. Time-varying amplitude perturbations (mi=100) and phase perturbations (mi=10) are shown. Similarity scores are derived from DTW and nDTW to highlight high distances (approaching 0) and low distances (approaching 1). Correlation is used as a benchmark. DTW and nDTW tracks amplitude disparities but remains invariant to phase shifts, unlike correlation. Both DTW and nDTW track amplitude disparities as mi increases (Time- varying amplitude sensitivity) but remain invariant to phase shifts, unli ke correlation (Time-varying phase sensitivity). c. Sensitivity of DTW and nDTW to 𝛾 values in amplitude and phase non-stationarity. The derivative of similarity scores under amplitude non -stationarity reveals that for lower mi, smaller 𝛾 values yield higher scores compared to larger 𝛾 values, whereas for higher mi, larger 𝛾 values yield higher scores (Time-varying amplitude tracking). This indicates that smaller 𝛾 values are more adept at tracking minor amplitude disparities, while larger 𝛾 values are more suited to tracking larger amplitude disparities. Although DTW is inherently phase -invariant, nDTW exhibits slightly higher derivatives indication slight sensitivity to phase differences when 𝛾 ≥ 1 (Time-varying phase tracking). d. DTW sensitivity across signal bandwidths. Random signal pairs (i.i.d Gaussian, 𝒩(0,1)) are filtered with increasing bandwidths up to the Nyquist frequency. DTW alignment and distances demonstrate increasing sensitivity to disparities in broadband signals, unlike correlation, which remains unaffected by the increasing frequency spectrum . This highlights DTW’s suitability for capturing disparities in full spectrum signals against slow narrowband ones. e. Illustration of convergence and divergence in brain networks. Using nDTW with 𝛾 = 1, the average amplitude difference between two pairs of aligned brain networks is shown. A small nDTW indicates convergence, while a large nDTW signifies divergence in signal power differences, providing insights into imbalances in metabolic energy demand. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 6 2.2. Signal power imbalances between schizophrenia and control in 140 broader fMRI bandwidth using nDTW 141 The intrinsic brain networks identified using the NeuroMark pipeline were categorized into 142 seven domains: subcortical (SC), auditory (Aud), sensory motor (SM), visual (Vis), 143 cognitive control (CC), default mode (DM), and cerebellum (Cb). 144 To validate the use of nDTW in the context of higher frequency fMRI data, we extended 145 our analysis to include broader frequency bands (slow -5, slow -4, and slow -3; 0.01 –146 0.198 Hz). We also examined the impact of varying 𝛾 values on test–retest reliability and 147 sensitivity to group differences. 148 Focusing on the F1 frequency range, we observed that test –retest reliability improved 149 with higher 𝛾 values, as indicated by lower mean ranks (Fig. 2a). 𝛾 values of ≥1.5 150 demonstrated significantly higher reliability compared to correlation after false discovery 151 rate (FDR) (Wilcoxon signed-rank test, FDR corrected p < 0.05). Given that physiological 152 noise correction in the fMRI data reduces within-subject variability31, thus increasing test-153 retest reliability, the higher test -retest score for higher 𝛾 values may be indicative of the 154 reduced sensitivity of higher 𝛾 values to physiological noise in the data. 155 Figure 2b illustrates that 𝛾 values of 1.5 and 1.25 exhibit the highest sensitivity for 156 detecting group differences across both frequency ranges, outperforming other 𝛾 values 157 and correlation. These 𝛾 values outperform other parameters and the correlation 158 benchmark, as indicated by the greater number of directional arrows pointing towards 159 them and the significant sensitivity enhancements marked by yellow asterisks. However, 160 since 𝛾 = 1.25 has significantly lower test–retest reliability than the correlation benchmark 161 (Fig. 2a), 𝛾 = 1.5 was selected as the optimal parameter for further analysis. Additionally, 162 nDTW consistently detected more group differences in the broader F1 frequency range 163 than in F2 across all 𝛾 values, whereas correlation did not. This indicates that nDTW is 164 more sensitive to capturing clinical functional relevance at higher BOLD signal 165 frequencies while maintaining high test–retest reliability, particularly at 𝛾 = 1.5. 166 Figure 2c shows that controls (lower triangle) consistently exhibit higher convergence in 167 signal power differences, whereas individuals with schizophrenia (upper triangle) display 168 greater divergence across several brain network pairs. This suggests that controls have 169 more balanced metabolic energy demands between brain networks, while schizophrenia 170 patients exhibit imbalances. Controls display significantly more convergent power 171 differences than schizophrenia patients across multiple network pairs (positive red cells 172 in the upper triangle of the FDR-corrected group difference matrix, Fig. 2c). Additionally, 173 Cohen’s d effect sizes (min -max = 0.19 − 0.89; median = 0.49) highlight the practical 174 significance of these group differences. 175 Furthermore, significant associations between nDTW scores and PANSS scores were 176 identified. Positive PANSS associations were primarily observed between Cb networks 177 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 7 and the DM, CC, and superior temporal gyrus ( Aud), Vis and SC networks. In contrast, 178 negative PANSS was mainly associated with connections between the CC, DM, and Vis 179 brain networks. 180 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 8 181 Fig. 2: Reliability and group difference sensitivity across 𝜸 values and frequency bands .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 9 182 2.3. Signal power imbalance effects of including higher frequency fMRI 183 We evaluated the effect of nDTW on expanding the frequency range from F2 to F1 by 184 conducting paired sample t -tests for controls and individuals with schizophrenia 185 separately. Normality of nDTW scores was confirmed using Shapiro -Wilk tests with t-186 values close to 1 across brain network pairs (Fig. 3a). After FDR correction, controls 187 exhibited significantly increased convergence in F1 compared to F2, whereas 188 schizophrenia patients showed greater divergence in F1 versus F2 across specific brain 189 network pairs. These results suggest that broader F1 frequencies enhance metabolic 190 energy balance in controls but exacerbate imbalances in schizophrenia, indicating that 191 higher frequency fMRI enhances the identification of schizophrenia characteristics. 192 Cohen’s d effect sizes were modest (median d = 0.18 for controls and d = 0.15 for 193 schizophrenia). Nonetheless, higher frequencies amplified group distinctions, as reflected 194 by larger effect size differences between controls and schizophrenia across several brain 195 network pairs (Fig. 3a). Upset plot analysis identified 914 significant pairs across both 196 frequencies, with 196 unique to F1 and 10 unique to F2, highlighting the enhanced group 197 sensitivity of nDTW in F1 (Fig. 3b). 198 Figure 3b displays FDR -corrected p -values for significant pairs shared by both 199 frequencies (left panel), with F1 (lower triangle) showing higher practical relevance than 200 a. Test-retest reliability of nDTW was evaluated across multiple 𝛾 values using the HCP dataset, with correlation serving as a benchmark. The mean ranks for all twelve 𝛾 values and correlation are presented for the F1 frequency range. Lower mean rank values indicate higher test-retest reliability. 𝛾 values ≥ 1.5 and 𝛾 = 0.25 show significantly higher test-retest than correlation, as determined by a Wilcoxon signed -rank test (FDR-corrected p < 0.05) on the test–retest ranks. b. Gamma sensitivity in distinguishing schizophrenia from controls. We applied McNemar tests across all 𝛾 values for the F1 and F2 frequency ranges, as well as between F1 and F2. In the resulting matrices, arrows point to the parameter with higher sensitivity, and yellow asterisks highlight significant results with FDR - corrected p-values < 0.05. Results showed that 𝛾 values of 1.5 and 1.25 are the most effective in distinguishing schizophrenia patients from controls compared to other 𝛾 values and correlation for both F1 and F2 frequency ranges. Furthermore, the comparison between F1 and F2 demonstrates that all 𝛾 values increase group difference sensitivity when transitioning from F1 to F2, whereas correlation (rho) does not. These findings suggest that nDTW is more sensitive to group differences when higher frequencies are incorporated. c. Group difference between schizophrenia and controls using a 𝛾 = 1.5. nDTW values were z-scored across both groups and multiplied by -1, converting smaller nDTW values to positive scores (indicating convergence) and larger values to negative scores (indicating divergence). Group averages of nDTW scores reveal that controls exhibit convergence (lower triangle) while individuals with schizophrenia show divergence (upper triangle) across several brain network pairs. FDR - corrected p-values from a generalized linear model, controlling for sex, age, site, and mean frame displacem ent, are displayed in the lower triangle, with significant results (FDR < 0.05) highlighted in the upper triangle. The findings indicate that multiple brain network pairs show significant differences, with controls consistently demonstrating greater convergence than schizophrenia patients. Additionally, Cohen’s d effect sizes exceeding 0.2 underscore the practical relevance of these group differences. d Association of nDTW with PANSS scores. Significant associations between brain network pairs and both positive and negative PANSS scores (FDR-corrected p < 0.05) were identified using a generalized linear model with a Poisson distribution. The analysis controlled for sex, age, site, and mean frame displacement. These significant associations are displayed for both positive and negative PANSS, highlighting the relationship between nDTW scores and symptom severity in schizophrenia. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 10 F2 (upper triangle) based on effect sizes. Pairs unique to F1 had Cohen’s d  > 0.2, 201 underscoring nDTW's increased sensitivity and clinical relevance at higher frequencies. 202 To validate F1-unique pairs, we conducted null hypothesis testing using 1,000 surrogate 203 datasets generated through phase randomization within the 0.15 –0.198 Hz range 204 (Fig. 3c). Several network pairs remained significant, confirming that higher frequency 205 fMRI results are statistically robust and potentially derived from non-stationary amplitude 206 processes. 207 Finally, we examined PANSS associations exclusively linked to F1 and not F2. Positive 208 PANSS correlated with connections between the Cb with SC, DM, CC, Aud and Vis 209 networks. Negative PANSS was primarily associated with CC -DM connections. These 210 associations were absent in F2, suggesting that higher frequency fMRI enhances 211 symptom severity detection. 212 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 11 213 Fig. 3: Clinical significance of high frequency fMRI data inclusion. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 12 214 2.4. Dynamics of time-resolved nDTW 215 To identify recurring states of power imbalance dynamics, we applied k-means clustering 216 with three clusters on the time-resolved nDTW, including both controls and individuals 217 with schizophrenia. An elbow plot (Fig. 4c) confirmed that three clusters were optimal, 218 determined by the derivative of the within-sum of squared distances. For visualization, we 219 subtracted the median of all clusters and multiplied the distribution by −1, ensuring 220 positive values represent convergence and negative values indicate divergence (Fig. 4a). 221 The clusters were classified as follows: state 1—a convergent state, state 3—a divergent 222 state, and state 2—a mixed state showing neither strong convergence nor divergence. 223 Figure 4b illustrates brain network representations for each state. In convergent state 1, 224 the top 5% most convergent pairs predominantly involved sensory networks, including 225 Vis, SM, and Aud networks. The divergent state 3 featured the top 5% most divergent 226 pairs, mainly between the Cb and other networks. The mixed state 2 showed the top 2.5% 227 most convergent and top 2.5% most divergent pairs, combining both dynamics. 228 A Markov chain was constructed to model state transitions, with transition probabilities 229 shown in Fig ure 4d. Mixed state 2 had the highest self -retention probability (0.9021), 230 followed by the convergent state 1 (0.8485) and divergent state 3 (0.8133). No direct 231 transitions occurred between divergent and convergent states, likely due to the 232 monotonicity constraint in the DTW alignment process, which prevents abrupt shifts 233 without intermediate states. 234 The Markov chain’s “stable point” (stationary distribution: Fig. 4e) showed that mixed state 235 2 accounted for 57.36% of system dynamics, the convergent state 1 for 37.21%, and the 236 divergent state 3 for 5.42%. This highlights the predominance of mixed and convergent 237 states in stabilizing power imbalance dynamics. 238 We further analyzed these dynamics by calculating the entropy rate and the convergence 239 rate to the stationary distribution. The entropy rate was 36.6%, where 100% represents 240 complete randomness and 0% indicates full order. This suggests a balance between 241 randomness and order, with a slight preference for order (Fig. 4f) . The rate of 242 a. Statistical analysis of power imbalances between F1 and F2 frequencies. Shapiro-Wilk tests confirmed normality (𝑇-𝑣𝑎𝑙𝑢𝑒𝑠 ≈ 1) for all groups, validating the use of paired t-tests. FDR-corrected p-values of paired t-test indicated that controls exhibit increased signal power convergence, whereas individuals with schizophrenia show power divergence when transitioning from F2 to F1. This suggests that controls demonstrate greater convergence, while schizophrenia patients exhibit increased divergence with the inclusion of higher frequencies. Cohen’s d effect sizes (~0.4) highlight enhanced separability between groups. b. Group differences from a generalized linear model controlling for sex, age, mean frame displacement, and s ite. FDR-corrected p-values across brain network pairs revealed 914 common significant pairs, with F1 showing larger effect sizes and greater practical relevance (left panel). Additionally, 196 significant pairs were uniquely identified in F1 with effect s izes above 0.2 (right panel), indicating modest practical relevance, while only 10 pairs were unique to F2. c. Null hypothesis testing using phase randomization confirmed that the significant brain network pairs uniquely identified in F1 are statistically relevant. d. Associations between normalized nDTW and PANSS scores were exclusively observed in the F1 frequency range, with no significant associations found in F2. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 13 convergence computed by the spectral gap was 0.1851, and the mixing time—defined as 243 the number of timesteps required to reach a tolerance level of 10−3 from the stationary 244 distribution—was 36 timesteps (72 seconds ; TR=2). The mixing time was determined 245 from the farthest distribution point in the Markov chain’s state space, specifically the state 246 [0, 0, 1] (completely divergent), as shown in Fig. 4e. 247 We explored the relationship between entropy and convergence rate and investigated 248 potential interventions to accelerate convergence and enhance cognitive ability by 249 introducing perturbations to each transition probability and measuring the relative 250 changes in convergence and entropy rates (Fig. 4g–h). Increasing self -retention in any 251 brain state led to reductions in both convergence rate and entropy rate. While adding 252 transitions between the mixed state and other states (excluding mixed to divergent) 253 enhanced both rates, transitions between divergent and convergent states were marked 254 with red crosses due to their impracticality as no such transitions were identified in the 255 dynamics (Fig. 4d). 256 Notably, increasing transition probability from mixed to divergent states extend s 257 trajectories away from the stable point (Fig. 4e), reducing convergence speed. While this 258 increased entropy, our findings suggest that higher entropy does not universally 259 accelerate convergence; instead, the effect depends on the specific dynamic transitions 260 within the Markov chain’s state space. 261 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 14 262 Fig. 4: Dynamics of time-resolved nDTW .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 15 263 2.5. Group differences between controls and schizophrenia 264 We assessed group differences in metrics derived from time -resolved nDTW, including 265 mean dwell time, occupancy rate, state transition probabilities, stationary distribution, 266 spectral gap, mixing time, and entropy rate. Schizophrenia patients showed greater 267 association with mixed state 2 and divergent state 3, while controls predominantly 268 engaged with the convergent state 1. Specifically, schizophrenia patients had longer 269 mean dwell times in mixed and divergent states, whereas controls spent more time in the 270 convergent state (Fig. 5a). Similarly, occupancy rates (Fig. 5b) reflected higher mixed and 271 divergent state occupancy in schizophrenia, and greater convergent state occupancy in 272 controls. 273 The stationary distribution (Fig. 5d) revealed that schizophrenia dynamics are more 274 influenced by mixed and divergent states than controls , while controls are more 275 influenced by the convergent state than schizophrenia . State transition probabilities 276 showed controls were more likely to transition from mixed to convergent states (Fig. 5c). 277 In contrast, schizophrenia patients had increased transitions from convergent to mixed, 278 mixed to divergent, and higher self-retention in both mixed and divergent states. 279 280 a. Distribution of recurring brain states identified using K -means clustering on time -resolved nDTW. The median across all 3 cluster centroids were subtracted and multiplied by –1 to ensure positive values represent convergence and negative values represent divergence. Three states were identified: convergent, divergent, and mixed. b. Brain state matrices with corresponding axial brain thumbnails. For convergent and divergent states, the top 5% most convergent and 5% most divergent brain connections are displa yed, respectively. The mixed state shows the top 2.5% most convergent and top 2.5% most divergent connections. c. Elbow plot illustrating the within -sum of squared distances (WSS) for different cluster numbers. The first inflection point in the first derivative of WSS indicates that three clusters are optimal. d. State transition probabilities between the identified brain states, highlighting the likelihood of moving from one state to another. e. Stationary distribution within the Markov chain state space, demonstrating that the stable point of the dynamics is predominantly influenced by the mixed and convergent states. f. Mixing time of the dynamics, indicating the number of timesteps required to converge to the stationary distribution within a tolerance of 10−3. Additionally, the convergence rate is shown through the spectral gap and entropy rate. g. Effects of perturbing state transitions on the spectral gap. Small perturbations were added to each transition probability, and the relative change in spectral gap was calculated to assess changes in convergence rate. h. Effects of perturbing state transitions on the entropy rate. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 16 281 Fig. 5: Group differences and cognitive score associations .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 17 Regarding convergence to the stationary distribution, controls exhibited a higher spectral 282 gap, indicating a faster convergence rate (Fig. 5e). Using a tolerance of 10−3 and starting 283 from the completely divergent state ([0, 0, 1]), controls also had a lower mixing time, 284 reflecting quicker stabilization. No significant difference in entropy rate was observed, 285 suggesting that faster convergence rates in controls may stem from dynamic transition 286 efficiency rather than entropy. 287 2.6. Controls and schizophrenia cognitive score associations 288 We examined associations between time -resolved nDTW metrics and Computerized 289 Multiphasic Interactive Neurocognitive System (CMINDS) scores using generalized linear 290 models (GLM) for controls and schizophrenia patients separately to avoid group bias (Fig. 291 5h). 292 In controls, working memory was positively associated with transitions from convergent 293 to mixed states and negatively associated with convergent self -retention and mixed -to-294 divergent transitions, indicating that transitions into the mixed state support memory 295 performance. Verbal learning was positively associated with transitions from mixed to 296 convergent states and negatively with convergent -to-mixed transitions, suggesting that 297 the convergent state facilitates verbal learning. Additionally, the spectral gap was 298 positively associated with reasoning and problem -solving, highlighting the link between 299 faster convergence rates and higher cognitive performance. 300 In schizophrenia, reasoning and problem -solving were negatively associated with 301 transitions from mixed to convergent states but positively with self-retention in the mixed 302 state. Attention vigilance was positively associated with transitions from divergent to 303 mixed states and negatively with self -retention in the divergent state. These results 304 suggest that specific dynamic transitions, such as moving into the mixed state, enhance 305 a. Mean dwell time in each state for schizophrenia patients and controls, highlighting increased dwell time in mixed and divergent states for schizophrenia and convergent state for controls. b. Occupancy rate of each state, showing longer occupancy in mixed and divergent states for schizophrenia and in convergent state for controls. c. State transition probabilities between states for schizophrenia patients and controls, illustrating higher transition rates from convergent to mixed and mixed to divergent st ates in schizophrenia. d. Group difference in s tationary distribution of the Markov chain. e. Spectral gap, reflecting the convergence rate, is higher in controls compared to schizophrenia patients, suggesting faster convergence to the stable point in controls. f. Mixing time, representing the number of timesteps required to reach the stationary distribution within a tolerance of 10−3, is shorter in controls than in schizophrenia patients, indicating quicker convergence in controls. g. Entropy rate, showing the balance between order and flexibility, does not differ significantly between groups. h. Associations between cognitive scores and nDTW metrics, analyzed separately for each group to account for group bias. In controls, working memory and verbal learning are associated with specific state transitions and spectral gap is associated with reasoning and problem-solving. I n schizophren ia, reasoning, problem -solving, and attention vigilance are linked to different transitions and state retention. All analyses were conducted using generalized linear models controlling for sex, age, site, and mean frame displacement. Results were FDR-corrected with a significance threshold of 0.05. One asterisk (*) indicates p-values between 0.05 and 0.01, two asterisks ( **) indicate p -values between 0.01 and 0.001, and three asterisks ( ***) indicate p-values below 0.001. .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 18 cognitive abilities in schizophrenia, whereas others, like remaining in the divergent state, 306 may impede them. 307 3. DISCUSSION 308 This study developed a comprehensive framework to assess power imbalances between 309 brain networks using fMRI data, offering insights into metabolic energy demands while 310 accounting for varying temporal scales of brain activity. The physiological basis of the 311 BOLD signal is not entirely straightforward. Studies combining fMRI with PET imaging 312 have shown that the magnitude of the BOLD signal can track glucose utilization in the 313 brain4,5. However, other studies indicate that even under conditions of hyperemia (excess 314 amount of blood) and hyperoxia (excess amount of oxygen) , the BOLD signal still 315 accompanies neuronal activation, suggesting that it does not strictly reflect a deficiency 316 of glucose or oxygen 32,33. Research further supports this by proposing that due to the 317 slow response of NVC, neurons initially rely on astrocytic glycogen and lactate stores 318 before blood flow compensates for energy use 34,35. This indicates that blood flow does 319 not simply match real -time neuronal energy consumption 36. Nevertheless, neuronal 320 activation consistently triggers vasoactive signaling, leading to vasodilation through NVC, 321 regardless of brain baseline oxygen or glucose levels 33. Therefore, while the immediate 322 energy source may vary, the magnitude of the BOLD signal consistently reflects the 323 overall energy demand driven by neuronal activation, whether the energy is supplied by 324 blood-derived glucose and oxygen or by intracellular reserves . Tracking amplitude 325 changes in BOLD signals using DTW provides valuable insights into metabolic energy 326 demands, offering robustness against temporal variability introduced by neural 327 processing and hemodynamic modulations. 328 By utilizing our recently proposed nDTW metric 20, we explored functional relevance 329 beyond the typical 0.01–0.15 Hz range. Our time-resolved nDTW identified three intuitive 330 brain states: convergence, divergence, and a mixed state, indicating flexibility in energy 331 distribution. Simulations confirmed the robustness of nDTW in capturing amplitude and 332 power disparities between signals, tracking non -stationary amplitude changes, and 333 remaining invariant to timing and phase deformations. Importantly, our framework 334 demonstrated functional relevance by associating brain dynamics with clinical and 335 cognitive scores. 336 Static results revealed that control subjects exhibited stronger convergence, while 337 schizophrenia showed relatively greater divergence in brain network power disparities. 338 Our static analysis suggests that healthy controls maintain more uniform energy demand 339 across brain networks globally, whereas individuals with schizophrenia exhibit significant 340 imbalances. These imbalances reflect abnormal glucose metabolic demands and align 341 with existing studies demonstrating that schizophrenia is associated with glucose 342 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 19 metabolism abnormalities, including an increased risk of diabetes and insulin 343 resistance37,38. 344 A key finding was the increased convergence in controls and increased divergence in 345 schizophrenia when higher-frequency signals were included. These results indicate that 346 high-frequency BOLD signals—indicating rapid glucose utilization—play a crucial role in 347 supporting healthy cognition and enhancing the distinction between control subjects and 348 those with schizophrenia. In controls, the inclusion of these frequencies appeared to 349 increase the balance of metabolic energy demands across brain networks , while in 350 schizophrenia, it exacerbates imbalances . Furthermore, several brain network pairs 351 showed significant differences between controls and schizophrenia exclusively when 352 higher frequencies were included. The phase-randomized null model confirmed that these 353

were statistically significant and likely attributable to non -stationary amplitude 354 processes detected by nDTW. While synchronization studies have established the 355 functional relevance of higher -frequency BOLD signals 30, our findings provide a novel 356 perspective. We demonstrate that higher-frequency BOLD signals improve the distinction 357 between control individuals and those with schizophrenia by capturing non -stationary 358 energy demand disparity during rapid neural processing. 359 We also found positive associations between nDTW and positive symptoms, particularly 360 in connections involving the cerebellum with the Aud, DM, CC, SC, and Vis domains. The 361 cerebellum, recognized as a key predictive hub in the brain 39, not only regulates rapid 362 cognitive processing—including speed, capacity, and consistency 40— but also plays a 363 central role in orchestrating mental coordination. This observation aligns with the 364 “cognitive dysmetria” theory, which attributes schizophrenia to impaired coordination 365 among prefrontal regions, thalamic nuclei, and the cerebellum41. 366 Our findings further support this theory, revealing that altered cerebellar coordination of 367 neural processing is linked with positive symptoms 41. Specifically, lower divergence—or 368 relatively higher convergence—in schizophrenia was correlated with more severe positive 369 symptoms such as hallucinations, delusions, and disorganized thinking. Considering that 370 our results show significantly stronger convergence in controls compared to individuals 371 with schizophrenia, the observed positive associations with symptom severity may 372 indicate maladaptive or compensatory 42 cerebellar coordination of glucose utilization, 373 inadvertently exacerbating positive symptomatology in the affected networks. 374 Negative associations were observed between brain networks and negative symptoms, 375 indicating that lower convergence (or higher divergence) is linked to more severe negative 376 symptoms. These associations were most prominent between higher -order processing 377 networks such as DM and CC. Given that negative symptoms—characterized by reduced 378 motivation, pleasure, and emotional expression —are linked to high -level cognitive 379 processing43, our findings suggest that greater imbalances in metabolism demands within 380 and between the DM and CC are associated with increased negative symptoms in 381 schizophrenia. Prior studies have connected impaired metabolism to negative symptoms 382 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 20 like depression44,45, and a recent study found that improvements in glucose metabolism—383 achieved through prefrontal cortex stimulation —were associated with reductions in 384 negative symptoms46. 385 Notably, most of the PANSS associations identified were due to the inclusion of higher 386 frequencies, emphasizing the importance of higher frequency BOLD signal . Although 387 higher frequency bands have been associated with various sources of noise 47, work has 388 also identified high frequency fluctuations as a source of non -noise variability 30,48,49. 389 Alterations within such frequency bands may therefore reflect alterations to rapid neural 390 processing and rapid glucose utilization in schizophrenia. 391 Our framework identified three dynamic states: convergent, divergent, and mixed. The 392 convergent state reflects relatively balanced energy demands, suggesting that neural 393 circuits are operating within optimal ranges, with consistent glucose delivery supporting 394 stable metabolism and neural energy homeostasis 50. Conversely, the divergent state 395 represents imbalanced energy demands, where certain brain networks receive either 396 insufficient or excessive blood flow . This imbalance can destabilize neural activity and 397 disrupt the homeostatic regulation of brain energy, ultimately leading to reduced cognitive 398 performance45,50,51. The mixed state lies between these extremes, reflecting neither 399 strong convergence nor strong divergence. Transitions were observed only between the 400 mixed and either the convergent or divergent states, with no direct shifts between 401 convergent and divergent states (Fig. 4d). Thus, the mixed state appears to function as 402 an adaptive intermediary, enabling the brain to switch flexibly between functional states 403 as needed. 404 The brain criticality theory posits that the brain functions near a critical point —a delicate 405 balance between order and disorder52. The mixed state in our framework exemplifies this 406 balance, contributing 57.36% to the stationary distribution. Combined with the Markov 407 chain’s entropy rate of 36.6%, this finding suggests that the brain operates closer to the 408 critical entropy point of 50% rather than being too ordered or too disordered, with a slight 409 bias toward order that remains within the “zone of criticality ”53. Indeed, critical entropy 410 serves as an indicator of brain criticality54—a state that helps the brain maintain flexibility 411 and adaptability in its functions. 412 A key dynamic metric in our study is the convergence rate, which measures how quickly 413 a system reaches its stable state. Our perturbation analysis revealed that increasing self-414 retention probabilities reduced both convergence and entropy rates. This is because 415 higher self-retention enforces predictability, which runs counter to the flexibility required 416 by criticality theory. Generally, higher convergence rates correlated with increased 417 entropy, except for transitions from the mixed to divergent state—highlighting the 418 significant role of transition dynamics in convergence. State space analysis further 419 demonstrated that trajectories moving from the divergent state takes the longest to reach 420 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 21 the stable point (Fig. 4e), thereby validating why increasing mixed to divergent transitions 421 reduces the rate of convergence to the stable point . This pattern underscores the need 422 to balance order and disorder (i.e., to increase entropy) while optimizing the trajectory 423 toward stability. 424 Importantly, control subjects exhibited significantly higher convergence rates—evidenced 425 by metrics such as the spectral gap and mixing time —compared to patients with 426 schizophrenia, despite both groups maintaining similar entropy levels. This suggests that 427 although both groups operate near critical entropy, the dynamics of schizophrenia hinder 428 efficient convergence to stability, likely due to altered transition pathways. These findings 429 further imply that high entropy does not necessarily equate to enhanced information 430 integration or superior cognitive performance, challenging the notion that richer conscious 431 experiences always result from increased entropy55. 432 Cognitive flexibility was further underscored by the associations observed with the mixed 433 state. In control subjects, working memory performance was positively linked to 434 transitions from convergent to mixed states and negatively associated with both self -435 retention in the convergent state and transitions from mixed to divergent states. These 436 relationships suggest that the mixed state —reflecting critical dynamics —facilitates the 437 flexibility necessary for effective higher-order functions such as working memory 56. In 438 individuals with schizophrenia, reasoning and problem -solving abilities were negatively 439 associated with transitions from mixed to convergent states, while showing positive 440 associations with self-retention in the mixed state and transitions from divergent to mixed 441 states. Additionally, attention vigilance was negatively related to self -retention in the 442 divergent state but positively correlated with transitions from divergent to mixed states. 443 Overall, these findings emphasize the pivotal role of the mixed state in supporting 444 cognitive flexibility and indicate that disruptions in transitions to and from this state may 445 contribute to the cognitive impairments observed in schizophrenia. 446 However, the mixed state may not be optimal for every task. For example, in control 447 subjects, transitions from the mixed to convergent state were positively associated with 448 verbal learning, whereas extended self -retention in the mixed state showed negative 449 associations. Given findings that verbal learning benefits from structured cognitive 450 training57—consistent with the ordered nature of the convergent state—it appears that 451 while the mixed state promotes flexibility, tasks that require structured cognition may rely 452 more on the consistency of the convergent state. Overall, these results underscore the 453 importance of balancing flexibility and stability to optimize cognitive performance. 454 Group-level analysis revealed that controls predominantly engaged in the convergent 455 state, while schizophrenia subjects showed greater involvement in the mixed and 456 divergent states. Metrics such as mean dwell time, fraction rate, state transitions and 457 stationary distribution confirmed these patterns. The increased presence of divergent 458 states in schizophrenia aligns with the disorder ’s characteristic disruption in neural 459 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 22 dynamics58. Furthermore, as transitions toward the mixed state were found to decrease 460 the convergen ce rate to stability, the excessive engagement of the mixed state in 461 schizophrenia may reflect an overly sustained flexibility that potentially undermines 462 cognitive processes such as working memory. 463 Although our study provides a unifying platform connecting various brain mechanistic 464 concepts, the use of BOLD fMRI inherently reflects relative rather than absolute metabolic 465 energy demands. Since BOLD is a hemodynamic signal rather than a direct measure of 466 energy consumption, its amplitude does not correspond to absolute metabolic 467 differences. Additionally, preprocessing steps such as z -scoring remove mean and 468 variance information, making the signal inherently relative. Performing ICA further 469 emphasizes this, as ICA focuses on maximizing statistical independence and extracting 470 components based on relative differences rather than absolute amplitudes. As a result, 471 while we cannot directly map differences in network amplitudes onto absolute 472 physiological energy demands, relative changes over time in ICA signals remain 473 meaningful and provide valuable insights into how different brain networks engage with 474 varying levels of neural activity and associated metabolic demand. 475 In summary, we provide evidence to suggest that schizophrenia is characterized by a 476 pronounced disparity in metabolic energy distribution across the brain, both throughout 477 the scanning session and dynamically, becoming more evident during rapid neural 478 processing. Additionally, we observed maladaptive metabolic energy demands in a hub -479 like interaction between the cerebellum and other brain networks linked to positive 480 symptoms, while stronger heterogeneity in energy demands across high -order 481 processing networks, such as the DM and CC, was associated with negative symptoms. 482 Our findings provide a novel and intuitive perspective on brain criticality using blood flow, 483 revealing that cognitive flexibility—a key marker of brain criticality—is linked to improved 484 working memory and fluid intelligence in both controls and individuals with schizophrenia. 485 However, we also establish that cognitive flexibility may hinder performance on structured 486 tasks like verbal learning. Furthermore, transitions into a state of greater energy demand 487 heterogeneity are associated with declines in attention and vigilance in schizophrenia. 488 Finally, we propose a therapeutic intervention in Supplementary material using 489 photobiomodulation and the developed directional DTW. This comprehensive framework 490 linking blood flow to cognition and schizophrenia symptoms offers a unifying perspective 491 for various hypotheses of brain function in neuroscience research. 492 4. METHODS 493 4.1. Dynamic time warping as a measure of energy/power disparity 494 Dynamic time warping is a widely used signal processing algorithm that aligns signals 495 with nonlinear temporal deformations by applying adaptive “elastic” transformations to 496 achieve a meaningful similarity measure 19,59. When temporal deformations complicate 497 phase or timing relationships, it becomes logical to isolate these aspects and concentrate 498 solely on the amplitude of the signals . DTW tackles this by resolving temporal 499 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 23 dissimilarities, optimally aligning signal pairs to minimize a distance metric and bring them 500 to a shared temporal scale, leaving amplitude as the primary source of disparity. 501 The algorithm computes a distance cost matrix that captures the cumulative cost of 502 aligning pairwise elements of the sequences and identifies the optimal warping function 503 to adjust their time axes, minimizing a distance metric. To enhance computational 504 efficiency and ensure meaningful alignments, particularly when guided by domain 505 knowledge, a window size constraint is applied to limit the maximum allowable warping 506 between indices in the time series19. 507 For discrete signals 𝑥 and 𝑦 of lengths 𝑁 and 𝑀 respectively, the distance cost matrix, 508 𝐶𝑀𝐷𝑇𝑊(𝑖,𝑗), is initialized as follows: 509 𝐶𝑀𝐷𝑇𝑊(𝑖,𝑗)𝜖(𝑥,𝑦) = {0, 𝑖𝑓 (𝑖, 𝑗) = (0,0) +∞, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 510 The cells in 𝐶𝑀𝐷𝑇𝑊(𝑖,𝑗) that are within the allowable window, 𝑤, can be defined as: 511 𝐶𝑀𝐷𝑇𝑊(𝑖𝑤,𝑗𝑤)𝜖{|𝑖−𝑗|≤𝑤} = 𝜆𝛾(𝑥𝑖𝑤, 𝑦𝑗𝑤) + 𝑚𝑖𝑛{ 𝐶𝑀𝐷𝑇𝑊(𝑖𝑤−1,𝑗𝑤−1) 𝐶𝑀𝐷𝑇𝑊(𝑖𝑤−1,𝑗𝑤) 𝐶𝑀𝐷𝑇𝑊(𝑖𝑤,𝑗𝑤−1) (1) Following the parameterization of the DTW algorithm recently proposed 25, let 𝛾𝜖ℝ+ 512 denote that 𝛾 is a positive real number. Then, the generalized distance metric is given as: 513 𝜆𝛾(𝑥𝑖𝑤, 𝑦𝑗𝑤) = |𝑥𝑖𝑤 − 𝑦𝑗𝑤| 𝛾 (2) Where 𝛾 > 1 emphasizes large amplitude disparity, 𝛾 < 1 highlights small amplitude 514 disparity and 𝛾 = 1 provides a balance between both25. 515 The optimal warping path or function 𝜑, with length 𝐿, derived from DTW, is the path that 516 belongs to the set of all possible warping functions Π, where the cumulative distance 517 metric between corresponding points in signals 𝑥 and 𝑦 is minimized. It can be expressed 518 as: 519 𝜑 = arg 𝑚𝑖𝑛𝜑𝜖Π ∑ 𝐶𝑀𝐷𝑇𝑊(𝑖𝑤,𝑗𝑤) (𝑖𝑤,𝑗𝑤)𝜖{|𝑖−𝑗|≤𝑤} (3) Where 𝐿 ≥ max (𝑀, 𝑁) 520 The distance cost 𝐷 in DTW, which measures the phase/timing -invariant amplitude 521 disparity between the signal pair, is computed as the sum of distance metric along the 522 warping path. This cost is expressed as: 523 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 24 𝐷 = ∑ 𝐶𝑀𝐷𝑇𝑊(𝜑𝑥(𝜏),𝜑𝑦(𝜏)) 𝐿 𝜏=1 (4) Where 𝜏 is the time index for the warp function, 𝜑𝑥(𝜏) is the warping path along signal 𝑥, 524 and 𝜑𝑦(𝜏) is the warping path along signal 𝑦. 525 Notably, when the distance metric with 𝛾 = 2 is chosen the DTW distance cost becomes: 526 𝐷 = ∑ |𝑥(𝜑𝑥(𝜏)) − 𝑦(𝜑𝑦(𝜏))| 𝛾=2 𝐿 𝜏=1 (5) This represents the energy of the difference signal 𝑥(𝜑𝑥(𝜏)) − 𝑦(𝜑𝑦(𝜏)), with the 527 normalized DTW 20 distance cost 𝐷𝑛, representing the power of the difference signal, as 528 given below: 529 𝐷𝑛 = 1 𝐿 ∑ |𝑥(𝜑𝑥(𝜏)) − 𝑦(𝜑𝑦(𝜏))| 𝛾=2 𝐿 𝜏=1 (6) This positions the algorithm as a measure of t emporal/phase-invariant energy/power 530 disparity between two signals. Leveraging the 𝛾 −parameterization, DTW and nDTW 531 effectively capture energy and power disparities between signals, respectively. This 532 flexible framework adjusts to different amplitude scales based on the value of 𝛾. 533 We define convergence as the scenario where the energy or power differences —534 quantified by DTW and nDTW respectively —are relatively low, indicating that the signal 535 levels are close together. Conversely, when these differences are relatively high, we refer 536 to it as divergence. 537 In the context of fMRI signals, convergence between two brain networks suggests that 538 the metabolic energy demands, or consumption of the corresponding neural networks are 539 homogenous. Conversely, divergence suggests an imbalance in energy demands 540 between these networks. These patterns provide insights into the metabolic energy 541 specificity and functional interactions of brain networks. 542 4.2. Allowable window size for fMRI application 543 The two main challenges of DTW are its computational cost, with a time complexity of 544 𝑂(𝑀 × 𝑁), and the issue of singularities 19, where a single time point in one signal is 545 unnaturally matched to a large portion of the other signal. Both challenges are mitigated 546 by constraining the alignment to a window size that limits the maximum window of 547 allowable alignment60. Given the focus on temporal deformations such as leads/lags and 548 nonlinear dynamics like stretching and shrinking, selecting an appropriate window size to 549 capture relevant temporal deformations in fMRI data is crucial. In our previous study using 550 the proposed warp elasticity approach, we demonstrated that stretching and shrinking 551 dynamics could introduce anti-correlations24. Considering that anti-correlations may also 552 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 25 arise from lead/lag relationships, it is intuitive to choose a window size that spans the 553 length of the signal required to produce maximum anti-correlations (−1). This approach 554 has also been validated as optimal in prior fMRI studies 21. The longest wavelength 555 yielding an anti -correlation of −1 corresponds to the lowest frequency component 556 (0.01𝐻𝑧), which equates to a 100 -second window. Achieving an anti -correlation of −1, 557 requires a time shift of half the wavelength (50 seconds). For the DTW algorithm to 558 capture relevant temporal deformations in fMRI data, two 50 -second windows are 559 sufficient to account for deformations in both directions, resulting in a total window length 560 of 100 seconds 21, corresponding back to the low -cutoff frequency of 0.01𝐻𝑧. However, 561 since the filter design for the data —commonly Butterworth —is typically nonideal, we 562 instead recommend selecting a window size based on the spectral characteristics of the 563 activity signal, determined by the −3𝑑𝑏 point of the low -frequency cutoff of the fMRI 564 signal13,23,24,61. The equation for determining the window size is: 565 𝑓3𝑑𝑏 = 0.88 √𝑁2 − 1 𝐹𝑠 (7) Where 𝑁 represents the number of points in the window, 𝐹𝑠 denotes the sampling 566 frequency, and 𝑓3𝑑𝑏 is the low cutoff frequency. 567 4.3. Time-resolved nDTW 568 To derive a time-resolved evolution of power imbalances between signals, we modify the 569 nDTW measure by eliminating its summation and normalization components. This yields 570 a time-resolved power mismatch metric, 𝐷𝑡𝑟(𝜏), for the aligned signals, ensuring that its 571 average corresponds to the original nDTW value. The formulations are as follows: 572 𝐷𝑡𝑟(𝜏) = |𝑥(𝜑𝑥(𝜏)) − 𝑦(𝜑𝑦(𝜏))| 𝛾 𝐷𝑛 = 1 𝐿 ∑ 𝐷𝑡𝑟(𝜏) 𝐿 𝜏=1 (8) However, the DTW algorithm results in aligned signal lengths 𝐿 that vary across different 573 signal pairs, complicating comparisons between brain network pairs. To ensure that the 574 length of the time-resolved metric equals the signal length and averages to the nDTW, it 575 is necessary to interpolate 𝐷𝑡𝑟(𝜏) with the warped time scale (𝜏) to have the same length 576 as the original time index (𝑡), 𝑁 or 𝑀 (noting that 𝑁 = 𝑀). This interpolation must preserve 577 the integrity of the metric, particularly the first derivatives of the warping functions, which 578 are crucial for capturing the varying temporal scales of fMRI time series as established in 579 our previous study24. 580 To address this, we employed the piecewise cubic Hermite interpolating polynomial 581 (PCHIP) method for interpolating 𝐷𝑡𝑟(𝜏) to the original signal length 𝑁. PCHIP was 582 selected for its ability to preserve the shape and monotonicity 62 of the original warping 583 functions, thereby maintaining the temporal dynamics of the fMRI signals. Unlike standard 584 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 26 cubic splines, PCHIP avoids introducing artificial oscillations and overshoots 62, ensuring 585 the fidelity of the power mismatch measurements. Additionally, PCHIP guarantees 586 continuous first derivatives62, aligning with our requirement to maintain the derivatives of 587 the warping functions. This interpolation ensures that the interpolated metric 𝐷𝑡𝑟 𝑝𝑐ℎ𝑖𝑝(𝜏) 588 accurately averages to the 𝐷𝑛, facilitating consistent and comparable analyses across all 589 brain network pairs. 590 The interpolated time-resolved nDTW metric is thus defined as: 591 𝐷𝑡𝑟 𝑝𝑐ℎ𝑖𝑝(𝜏) = 𝑃𝐶𝐻𝐼𝑃(𝐷𝑡𝑟(𝜏)), 𝑤ℎ𝑒𝑟𝑒 |𝜏| = 𝑁 = 𝑀 𝐷𝑛 ≈ 1 𝑁 ∑ 𝐷𝑡𝑟 𝑝𝑐ℎ𝑖𝑝(𝜏) 𝑁 𝜏=1 (9) To validate our approach, we generated 10,000 independent and identically distributed 592 (i.i.d.) Gaussian signal pairs, each sampled from a standard normal distribution 𝒩(0,1). 593 Calculating the percentage relative difference between the nDTW and the average time-594 resolved nDTW yielded a median percentage difference of 0.20% (Supplementary 595 Material). 596 This methodology ensures that the time -resolved metric retains its temporal dynamics 597 and derivative continuity while standardizing its length. Consequently, preserving 𝐷𝑛 598 within the time -resolved measure maintains its reliability, thereby enhancing the 599 robustness of our analysis. 600 4.4. Resting-state fMRI data 601 The first dataset used is from the Function Biomedical Informatics Research Network 602 (fBIRN) study including resting fMRI data collected from schizophrenia patients and 603 controls63. Scans were collected at a repetition time (TR) of 2 seconds. The study led to 604 a total of 160 controls with an average age of 37.04 ± 10.86 years, ranging from 19 to 59 605 years. Among these, 45 were female and 115 were male. Additionally, there were 151 606 patients diagnosed with schizophrenia, with an average age of 38.77 ± 11.63 years, 607 ranging from 18 to 62 years. In this group, 36 were female and 115 were male. The typical 608 controls and schizophrenic patients were meticulously matched in terms of age, gender 609 distribution, and mean framewise displacement during scans (age: p = 0.18; gender: 610 p = 0.39; mean framewise displacement: p = 0.97). Schizophrenic patients were in a 611 clinically stable condition during the time of their scans. The fBIRN dataset is used in all 612 the clinical analysis performed in this study. 613 The second dataset used in our study is the resting-state fMRI dataset collected from 827 614 subjects via the Human Connectome Project (HCP) database 64,65. We utilized both 615 session scans acquired using a Siemens Skyra 3T scanner with a multiband accelerated, 616 gradient-echo echo-planar imaging (EPI) sequence. The scanning parameters were set 617 to a TR of 0.72s, 72 slices, an echo time (TE) of 58ms, and a flip angle of 90°. A voxel 618 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 27 size of 2 × 2 × 2 mm was used to acquire 1200 time points of fMRI data for each subject. 619 The HCP dataset was used for our test-retest analysis. 620 4.5. Resting-state fMRI pre-processing 621 fMRI data requires extensive preprocessing to correct for various sources of noise and 622 artifacts before analysis. The preprocessing steps commonly applied in fMRI studies 623 include slice timing correction, realignment, spatial normalization, and spatial 624 smoothing66. Following typical preprocessing67, we implemented the NeuroMark pipeline, 625 a fully automated spatially constrained ICA on the preprocessed fMRI data 68. Using the 626 neuromark_fMRI_1.0 template (available in GIFT 627 at http://trendscenter.org/software/gift or http://trendscenter.org/data) we generated 53 628 intrinsic connectivity networks (ICNs) for each subject. These ICNs are grouped into brain 629 functional domains, including subcortical, auditory, sensorimotor, visual, cognitive control, 630 default mode, and cerebellum. 631 4.6. Resting-state fMRI post-processing 632 To enhance the quality of these ICNs, we implemented detrending and despiking 633

to eliminate drifts, sudden fluctuations, and significant artifacts that may have 634 persisted after the initial preprocessing stages. The time series of the ICNs were 635 bandpass-filtered within two frequency ranges: F1 (0.01-0.198Hz) and F2(0.01-0.15Hz). 636 An infinite impulse response filter was designed using the butter function in MATLAB and 637 applied via the filtfilt function to ensure zero phase shifts and preserve phase information. 638 The optimal filter order was estimated using the Buttord function in MATLAB, which 639 returns the lowest order of the Butterworth filter with no more than 3dB ripple in the 640 passband and at least 30dB of attenuation in the stopband. Finally, we performed z -641 scoring on the ICNs. The MATLAB version used was MATLAB R2022a for all steps of the 642 analysis. 643 4.7. Test-retest reliability 644 Typically, test -retest reliability is evaluated using the intra -class correlation coefficient 645 (ICC). However, the ICC relies on a linear parametric ANOVA model that assumes 646 separability and additivity 69. This makes ICC unsuitable for DTW applications, as DTW 647 generates a distance metric based on nonlinear warping, where pairwise distances are 648 interdependent20,21 and thus lack decomposable variance components. We utilize the 649 HCP dataset, which includes two sessions and a large sample size of 827 subjects, 650 facilitating robust test-retest analysis. 651 4.7.1. Within-subject variability 652 We adopt the approach of a previous study by employing a straightforward method to 653 assess within-subject variability: a paired sample t -test between sessions for each brain 654 network pair21. The absolute t-value is reported as a measure of how much the differences 655 between the two sessions deviate from zero, where lower absolute t -values indicate 656 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 28 higher within -subject consistency. The formula for calculating the absolute t -value 657 between two sessions of a single network pair across all subjects is provided below: 658 |𝑡|𝑤𝑖𝑡ℎ𝑖𝑛 = | 𝑑̅ 𝑠𝑑/√𝑛| (10) Where 𝑑̅ is the mean difference between two sessions, 𝑠𝑑 is the standard deviation of the 659 differences, and 𝑛 is the total number of subjects. 660 4.7.2. Between-subject variability 661 To assess between-subject variability, we developed a non-parametric equivalent of the 662 absolute t-value. First, we calculate the differences between the means of each subject’s 663 distribution to capture variability across the group. Since the direction of these differences 664 is not of primary interest, we take their absolute values, resulting in a distribution skewed 665 toward zero. This transformation allows us to perform a non -parametric test to evaluate 666 how a non-parametric central tendency measure, such as the median deviates from zero. 667 Unlike the parametric approach, which scales the central tendency measure by the 668 standard deviation, our non -parametric method scales it by the interquartile range to 669 account for the skewed distribution. The formula for the non -parametric absolute 670 variability measure for a single session is as follows: 671 𝑑𝑖 = |𝑠1,𝑖 − 𝑠1̅ | |𝑡|𝑏𝑒𝑡𝑤𝑒𝑒𝑛 = |𝑚𝑒𝑑𝑖𝑎𝑛(𝑑) 𝐼𝑄𝑅(𝑑)/√𝑛| (11) Where 𝑑𝑖 is the absolute difference between subject 𝑖’s measure in session 1 ( 𝑠1,𝑖) and 672 the group mean for session 1 ( 𝑠1̅ ), 𝐼𝑄𝑅(𝑑) is the interquartile range of the differences 𝑑, 673 and 𝑛 is the total number of subjects. 674 This non-parametric approach provides a robust measure of between -subject variability 675 that is less influenced by outliers and non -normal distributions, thereby enhancing the 676 reliability of our variability assessments. 677 4.7.3. Composite test-retest reliability score 678 To further summarize these two approaches to give a wholistic view of test -retest 679 reliability, we compute a composite score of test-retest as follows: 680 𝑡𝑟𝑡 = 1 − |𝑡|𝑤𝑖𝑡ℎ𝑖𝑛 |𝑡|𝑏𝑒𝑡𝑤𝑒𝑒𝑛 (12) Scores approaching 1 indicate higher test-retest reliability. 681 In our study, we compared 𝑡𝑟𝑡 scores across various gamma values of nDTW against 682 correlation as a benchmark. To achieve this, we ranked the 𝑡𝑟𝑡 scores and performed 683 Wilcoxon signed-rank tests between each gamma value’s 𝑡𝑟𝑡 scores and the correlation 684 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 29 benchmark across all brain network pairs , following a critical difference computation 685 approach. 686 4.8. Group difference and clinical association analyses 687 We employed generalized linear models for both group difference analyses and 688 associations with clinical scores. All GLM analyses controlled for sex, age, site, and mean 689 frame displacement. P-values were adjusted using the false discovery rate method, with 690 significance set at a threshold of 0.05. For group difference analyses and computerized 691 multiphasic interactive neurocognitive system scores, a normal distribution was assumed, 692 whereas a Poisson distribution was applied for PANSS score associations due to the 693 positive integer and skewed nature of the fBIRN dataset’s PANSS scores. Effect sizes 694 were calculated using Cohen’s d. 695 4.9. Group difference sensitivity test 696 To evaluate the sensitivity of various gamma values of nDTW and correlation in 697 distinguishing between controls and individuals with schizophrenia, we conducted 698 McNemar’s sensitivity tests comparing each gamma value and correlation against each 699 other. Following the GLM group difference analyses, we counted the number of significant 700 brain network pairs ( FDR corrected p < 0.05) for each gamma value (nDTW) and for 701 correlation. We then compared the unique and shared significant pairs between each 702 gamma and correlation using McNemar’s test, with p-values adjusted via FDR correction 703 across all comparisons. 704 These sensitivity tests were performed separately for each frequency range, F1 and F2. 705 Additionally, we compared F1 and F2 across all gamma values and correlation metric to 706 assess how transitioning from the lower frequency range F2 to the broader F1 range 707 affects group difference sensitivity. 708 4.10. Validation of high fMRI frequency results. 709 High-frequency fMRI signals are vulnerable to various sources of noise and artifacts, 710 including physiological factors like respiratory and cardiac activity, as well as non -711 neuronal influences such as head motion 70. Comprehensive preprocessing steps 712 followed by our use of ICA help mitigate these concerns. ICA effectively separates 713 independent components, allowing us to identify and remove artifact-related components 714 from the data. Additionally, we included head motion parameters as covariates in our 715 GLM analyses to ensure that detected group differences in high -frequency signals are 716 not attributable to head motion effects. 717 We conducted several additional analyses further to validate our findings within the F1 718 frequency range. 719 4.10.1. Statistical effects on groups with the inclusion of higher frequency 720 We aimed to evaluate the impact of transitioning from the F2 to the broader F1 frequency 721 range on both control and schizophrenia groups. Since F1 and F2 comprise the same 722 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 30 dataset differentiated only by frequency range, a paired sample t-test was appropriate for 723 this analysis. To validate the use of the paired t -test, we first assessed the normality of 724 nDTW scores for each group between F1 and F2 using the Shapiro -Wilk test, ensuring 725 that t-values were close to 1. 726 Subsequently, we conducted paired t -tests separately for the control and schizophrenia 727 groups to examine the effects of transitioning from F2 to F1, with p-values adjusted using 728 the FDR correction. Effect sizes for these comparisons were calculated using Cohen’s d. 729 Additionally, to highlight the differential impact of transitioning from F2 to F1 between 730 groups, we computed the differences in Cohen’s d effect sizes between controls and 731 schizophrenia. This comparison aimed to determine whether the transition to the broader 732 F1 frequency range enhances or diminishes group-specific distinctions. 733 Furthermore, we identified statistical group differences between schizophrenia and 734 controls that were present in F1 but not in F2 by analyzing the shared and unique 735 significant brain network pairs from our GLM analysis (see Section 3.7). We also 736 highlighted their corresponding effect sizes, demonstrating that significant group 737 differences are more pronounced in the F1 frequency range. 738 4.10.2. Null hypothesis test 739 The phase randomization (PR) technique generates surrogate fMRI data by preserving 740 the correlation between brain network pairs 71. This method is suitable for fMRI null 741 hypothesis testing as it assesses non-stationarity by capturing information beyond linear 742 and static correlations 71. In our previous study 20, we employed PR to evaluate the 743 statistical relevance of both DTW and nDTW in fMRI data, demonstrating that both 744 metrics rejected the null hypothesis for several brain region pairs. This validation confirms 745 that the DTW and nDTW methods are sensitive to capturing information beyond 746 traditional correlation measures. Since DTW and nDTW are amplitude-sensitive methods 747 and the intrinsic brain network timecourses are z -scored, they likely capture amplitude 748 non-stationarity, potentially providing insights into temporal variations in metabolic energy 749 demands between networks. 750 In this study, we utilize PR to assess the statistical relevance of nDTW in capturing non -751 stationary amplitude changes specifically due to the inclusion of higher frequency data 752 (0.15–0.198 Hz). The following steps were followed to generate 1,000 surrogate fMRI 753 datasets: 754 I. Phase Randomization: Apply PR to the F1 frequency range data (0.01–0.198 Hz) 755 to create surrogate data, denoted as 𝑍. 756 II. Filtering 𝑍: Use the designed filter applied to achieve original F2 fMRI data (see 757 section 3.4.) on 𝑍 to extract the 0.01 –0.15 Hz frequency range of the surrogate 758 data, resulting in surrogate data 𝑌. 759 III. Isolating Higher Frequencies: Subtract 𝑌 from 𝑍 to obtain the surrogate data 760 residual high frequency 0.15–0.198 Hz, denoted as 𝑋. 761 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 31 IV. Full surrogate data: Add 𝑋 to the original F2 fMRI data (0.01–0.15 Hz) to generate 762 surrogate data with phase randomization affecting only the residual 0.15–0.198 Hz 763 frequency range. 764 This process generates surrogate datasets where randomization impacts solely the 765 residual 0.15–0.198 Hz frequency range. A high-pass filter was not applied to 𝑍 to obtain 766 the 0.15 –0.198 Hz surrogate data because non -ideal filters will artificially create 767 frequency multiples at the edges (~0.15 Hz). 768 We performed a Wilcoxon rank -sum test (FDR -corrected p < 0.05) to compare nDTW 769 scores from the original fMRI data with those from surrogate data. This test assessed 770 whether the significant brain network pairs distinguishing controls from schizophrenia, 771 identified exclusively in the F1 range, also rejected the null hypothesis. 772 4.10.3. PANSS association exclusive to F1 773 We identified significant brain network pairs associated with symptom scores that were 774 exclusive to the F1 frequency range and not present in the F2 range. 775 4.11. Power imbalance dynamics 776 4.11.1. Brain state estimation 777 In time-resolved fMRI studies, k-means clustering is frequently utilized to identify recurring 778 brain states within the temporal dynamics of the data 72. The DTW algorithm inherently 779 anchors the first and last time points during alignment, potentially introducing bias by 780 disproportionately influencing brain states at the temporal extremes. To mitigate this bias, 781 we excluded the first and last 5 time points from the time-resolved nDTW metrics prior to 782 clustering. This exclusion reduces alignment -induced artifacts, enabling the clustering 783 algorithm to more accurately capture intrinsic transitions between brain states. 784 We concatenated the time-resolved nDTW metrics from all subjects in the fBIRN dataset 785 into a single feature matrix, where each time point from every subject serves as an 786 individual sample and each brain network pair constitutes a distinct feature . We 787 implemented the k -means algorithm with cluster numbers ranging from 1 to 10, setting 788 the maximum number of iterations to 10,000 and using 20 random initializations to ensure 789 robust convergence. The city -block (Manhattan) distance metric was selected for its 790 demonstrated robustness in handling high-dimensional data73. To determine the optimal 791 number of clusters, we generated an elbow plot based on the within -cluster sum of 792 squares (WSS). The optimal cluster number was identified as 3, corresponding to the first 793 significant inflection point in the WSS curve’s first derivative. 794 4.11.2. Markov chain & stationary distribution 795 We quantified the transitions between brain states by counting the number of times each 796 subject moved from one state to another and normalizing these counts across each state 797 to construct individual transition probability matrices. This represents the probability of a 798 subject transitioning from one state to another and further provides insights into the 799 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 32 Markov chain of state dynamics. To gain a comprehensive understanding of the general 800 transition patterns independent of specific groups, such as patients and controls, we 801 aggregate the probability transitions across all subjects in the fBIRN dataset. 802 Furthermore, the stationary distribution of the Markov chain represents the long -term 803 behavior of the system, indicating the proportion of time the system spends in each state 804 once it has converged to equilibrium 74. This distribution serves as an indication of a 805 “stable point” when analyzing systems from the perspective of their Markov chains rather 806 than traditional state-space representations. For the stationary distribution to be unique 807 and ensure convergence from any initial state, the Markov chain must be ergodic 75. A 808 Markov chain is considered ergodic if it is irreducible (it is possible to get from any state 809 to any other state, aperiodic (the system does not cycle in a fixed pattern), and positive 810 recurrent (the expected return time to each state is finite). 811 The stationary distribution 𝜋 is computed by solving the following system of linear 812 equations: 813 𝜋𝑃 = 𝜋 ∑ 𝜋𝑘 𝐾 𝑘=1 = 1 (13) Where 𝜋 = [𝜋1, 𝜋2, … , 𝜋𝑘] is the stationary distribution vector, 𝑃 is the transition probability 814 matrix of the Markov chain, 𝐾 is the number of distinct states from K -means (3 in our 815 study). 816 4.11.3. Convergence analysis 817 How fast a chain takes to reach the stationary distribution gives a similar equivalent 818 meaning to how fast a system takes to reach its stable point from any given state of the 819 system. The spectral gap gives the rate of this convergence to equilibrium 74. A larger 820 spectral gap implies faster convergence rate to the stationary distribution. It is computed 821 as follows: 822 𝛾𝑠 = 𝜆1 − 𝜆2 (14) Where 𝜆1, and 𝜆2 are the largest and second largest eigen values of the Markov chain 823 respectively. When the chain is ergodic, 𝜆1 = 1. 824 Beyond the spectral gap, the mixing time provides a more comprehensive measure of 825 how quickly the chain approaches equilibrium 74. The mixing time quantifies the number 826 of steps required for the chain's state distribution to approach the stationary distribution 827 within a specified tolerance level. A prevalent method for determining the mixing time 828 leverages the total variation (TV) distance, which offers a rigorous metric for assessing 829 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 33 the divergence between two probability distributions. The total variation distance at time 830 step 𝑡 is defined as: 831 ‖𝜇𝑃𝑡 − 𝜋‖𝑇𝑉 = 1 2 ∑|𝜇𝑃𝑡(𝑘) − 𝜋(𝑘)| 𝐾 𝑘=1 (14) Where 𝜇 is the initial probability distribution, 𝑃𝑡 represents its 𝑡-th power of transition 832 matrix, indicating the state probabilities after 𝑡 transitions, 𝜋(𝑘) is the stationary probability 833 of state 𝑘 and 𝐾 is the total number of distinct states. 834 The mixing time, denoted as 𝑡𝑀𝑇, is the smallest positive integer 𝑡 for which the TV 835 distance falls below a predefined tolerance level 𝑡𝑜𝑙: 836 𝑡𝑀𝑇 = 𝑚𝑖𝑛{𝑡 ∈ ℝ+ | ‖𝜇𝑃𝑡 − 𝜋‖𝑇𝑉 ≤ 𝑡𝑜𝑙} (15) For our analysis, we set the tolerance level to 𝑡𝑜𝑙 = 10−3. We also find the mixing time 837 using the initial probability distribution 𝜇 farthest from the stationary distribution. 838 These metrics —spectral gap and mixing time —are instrumental in understanding the 839 dynamic behavior of brain network states. The spectral gap provides insight into the 840 inherent convergence rate of the system, while the mixing time offers a tangible measure 841 of how quickly the system reaches equilibrium. Together, they facilitate a comprehensive 842 analysis of power imbalances and their temporal evolution within the fMRI data, 843 enhancing the robustness and reliability of our findings. 844 4.11.4. Entropy rate 845 We investigated how power imbalance dynamics, indicative of the metabolic energy 846 demands between brain networks, inform the flexibility or rigidity of the brain’s metabolic 847 consumption by analyzing the entropy rate. The entropy rate is defined as the limit of the 848 average entropy per step as the number of steps approaches infinity. For an ergodic 849 Markov chain, the entropy rate in bits can be mathematically expressed as: 850 𝐻 = − ∑ 𝜋(𝑖) 𝐾 𝑖=1 ∑ 𝑃(𝑖, 𝑗)𝑙𝑜𝑔2𝑃(𝑖, 𝑗) 𝐾 𝑗=1 (16) Where 𝜋(𝑖) is the stationary probability of being in state 𝑖 and 𝑃(𝑖, 𝑗) is the transition 851 probability from state 𝑖 to state 𝑗. 852 The maximum entropy rate for 𝐾 = 3 states is: 853 𝐻𝑚𝑎𝑥 = log2(3) = 1.585 𝑏𝑖𝑡𝑠 (17) In our study, we present entropy rate as a percentage of 𝐻𝑚𝑎𝑥 to enhance interpretability. 854 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 34 The entropy rate provides insight into the balance between order and randomness within 855 the power imbalance dynamics. A low percentage entropy rate indicates strong order in 856 the dynamics, potentially reducing metabolic flexibility. Conversely, a high percentage 857 entropy rate signifies a high degree of randomness, suggesting unpredictable metabolic 858 energy demands between brain networks. 859 4.11.5. Impact of transition probabilities on Markov chain c onvergence 860 and entropy rate. 861 To evaluate how specific transitions within the Markov chain influence the convergence 862 and entropy rates of the dynamics, we systematically perturbed the aggregated transition 863 matrix. This process involved incrementally increasing the probability of individual 864 transition cells and subsequently renormalizing the corresponding rows to ensure that 865 each row’s probabilities sum to one. Consequently, the targeted transition’s probability 866 increased while the probabilities of other transitions within the same row decreased. 867 For each perturbed transition, we recalculated both the spectral gap and the entropy rate 868 to determine their percentage relative changes compared to the original values. This 869 perturbation procedure was applied to every transition cell across 10,000 random 870 samples, with perturbation magnitudes ranging from 10−1 to 10−3 on a logarithmic scale 871 to ensure uniform sampling. To maintain comparability between transitions, each sample 872 employed a single perturbation magnitude uniformly across all transition cells. 873 4.12. Time-resolved nDTW group difference: schizophrenia vs controls 874 To assess the clinical relevance of the time -resolved nDTW measure, we conducted 875 group difference analyses on several derived metrics: state transitions, stationary 876 distribution, spectral gap, mixing time, and entropy rate (see Section 3.11). Unlike the 877 aggregated approach in Section 3.11, these metrics were calculated individually for each 878 subject. To ensure the uniqueness of these results, we excluded subjects with non -879 ergodic state transitions. 880 Additionally, we included mean dwell time and occupancy rate in our group comparisons. 881 Mean dwell time quantifies the average duration a subject remains in a particular state 882 once entered76, reflecting the persistence of states within the Markov chain. Occupancy 883 rate represents the percentage of time a subject spends in each state 76, indicating the 884 relative prevalence of states during the observation period. 885 Using these subject -specific metrics, we distinguished schizophrenia patients from 886 control subjects and performed group difference analyses employing a GLM, as 887 described in Section 3.8. FDR correction is applied to each metric across all three 888 identified states. 889 4.13. CMINDS score association 890 Similarly, we employed a GLM, as detailed in Section 3.8, to assess associations 891 between each metric and CMINDS scores —including speed of processing, attention 892 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 35 vigilance, working memory, verbal learning, visual learning, and reasoning/problem 893 solving—within the fBIRN dataset. To account for potential group effects, separate GLM 894 analyses were conducted for each group. For transition probabilities, we applied FDR 895 correction to control for multiple comparisons across all transitions. 896 4.14. Directional time-resolved nDTW 897 The derivation of DTW, nDTW, and our developed time -resolved nDTW focuses on the 898 absolute differences between aligned signal pairs, revealing power imbalances between 899 brain networks. However, using only absolute differences obscures directional 900 information—specifically, it does not indicate which brain network has higher metabolic 901 demand. Understanding this directionality is essential for designing targeted interventions 902 that address specific metabolic flexibilities or rigidities within neural circuits. 903 To retain directional information, we modified the time -resolved nDTW by modulating it 904 with directional information, resulting in directional time-resolved nDTW: 905 𝐷𝑡𝑟 𝑑 (𝜏) = 𝒔𝒊𝒈𝒏 (|𝑥(𝜑𝑥(𝜏))| − |𝑦(𝜑𝑦(𝜏))|) × 𝐷𝑡𝑟(𝜏) (18) This modification allows us to detect both the magnitude and direction of power 906 imbalances, identifying which brain network exerts greater metabolic demand relative to 907 its counterpart. We then averaged the directional time -resolved nDTW within each 908 identified cluster, enabling the determination of brain networks that consistently display 909 higher or lower power across states. This refined metric provides a more nuanced 910 understanding of metabolic energy demands, offering precise insights to inform 911 interventions aimed at modulating metabolic flexibility or addressing rigidities within 912 specific neural pathways, thereby enhancing the efficacy of clinical treatments. 913 To evaluate the similarity of directional information across states, we calculated 914 correlations between the directional time -resolved nDTW metrics of each state. 915 Additionally, we performed a Kruskal-Wallis test to determine whether the distributions of 916 directional information differ among the three states. This dual analysis serves as a sanity 917 check, as we anticipate distinct distributions even if their correlations are comparable. 918 5. DECLARATION OF COMPETING INTEREST. 919 None. 920 6. AUTHORS. 921 Sir-Lord Wiafe: Conceptualization, Formal analysis, Methodology, Visualization, Writing 922 –original draft. Spencer Kinsey : Methodology, Validation, Writing – review & editing. 923 Najme Soleimani : Validation, Writing – review & editing. Raymond O. Nsafoa : 924 Validation, Supervision . Nigar Khasayeva : Validation, Writing – review & editing. 925 Amritha Harikumar : Validation, Writing – review & editing. Robyn Miller: Validation, 926 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 36 Methodology, Supervision, Writing – review & editing . Vince Calhoun: Funding 927 acquisition, Validation, Methodology, Resources, Supervision, Writing – review & editing. 928 7. DATA & CODE AVAILABILITY STATEMENT. 929 The codes for all our analyses in MATLAB language can be accessed through GitHub 930 (study code). The data was not collected by us and was provided in a deidentified manner. 931 The IRB will not allow sharing of data or individual derivatives as a data reuse agreement 932 was not signed by the subjects during the original acquisition. 933 8. ACKNOWLEDGMENT. 934 This work was supported by the National Institutes of Health (NIH) grant (R01MH123610) 935 and the National Science Foundation (NSF) grant #2112455. 936 9. REFERENCES. 937 1. Rolfe D, Brown GC. Cellular energy utilization and molecular origin of standard metabolic rate in 938 mammals. Physiological reviews. 1997;77(3):731-758. 939 2. Magistretti Pierre J, Allaman I. A Cellular Perspective on Brain Energy Metabolism and Functional 940 Imaging. Neuron. 2015;86(4):883-901. 941 3. Logothetis NK, Pauls J, Augath M, Trinath T, Oeltermann A. Neurophysiological investigation of 942 the basis of the fMRI signal. nature. 2001;412(6843):150-157. 943 4. Raichle ME. Behind the scenes of functional brain imaging: a historical and physiological 944 perspective. Proceedings of the National Academy of Sciences. 1998;95(3):765-772. 945 5. Tomasi D, Wang G -J, Volkow ND. Energetic cost of brain functional connectivity. Proceedings of 946 the National Academy of Sciences. 2013;110(33):13642-13647. 947 6. Buxton RB, Uludağ K, Dubowitz DJ, Liu TT. Modeling the hemodynamic response to brain 948 activation. Neuroimage. 2004;23 Suppl 1:S220-233. 949 7. Golesorkhi M, Gomez-Pilar J, Zilio F, et al. The brain and its time: intrinsic neural timescales are 950 key for input processing. Communications Biology. 2021;4(1):970. 951 8. Zang YF, He Y, Zhu CZ, et al. Altered baseline brain activity in children with ADHD revealed by 952 resting-state functional MRI. Brain Dev. 2007;29(2):83-91. 953 9. Zou QH, Zhu CZ, Yang Y, et al. An improved approach to detection of amplitude of low-frequency 954 fluctuation (ALFF) for resting -state fMRI: fractional ALFF. J Neurosci Methods. 2008;172(1):137-955 141. 956 10. Fu Z, Tu Y, Di X, et al. Characterizing dynamic amplitude of low -frequency fluctuation and its 957 relationship with dynamic functional connectivity: An application to schizophrenia. NeuroImage. 958 2018;180:619-631. 959 11. Allen EA, Damaraju E, Plis SM, Erhardt EB, Eichele T, Calhoun VD. Tracking whole -brain 960 connectivity dynamics in the resting state. Cereb Cortex. 2014;24(3):663-676. 961 12. Zhang G, Cai B, Zhang A, et al. Estimating Dynamic Functional Brain Connectivity With a Sparse 962 Hidden Markov Model. IEEE Trans Med Imaging. 2020;39(2):488-498. 963 13. Wiafe S-L, Asante NO, Calhoun VD, Faghiri A. Studying time -resolved functional connectivity via 964 communication theory: on the complementary nature of phase synchronization and sliding 965 window Pearson correlation. bioRxiv. 2024:2024.2006. 2012.598720. 966 14. Chang C, Glover GH. Time–frequency dynamics of resting-state brain connectivity measured with 967 fMRI. NeuroImage. 2010;50(1):81-98. 968 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 37 15. Mitra A, Snyder AZ, Hacker CD, Raichle ME. Lag structure in resting -state fMRI. J Neurophysiol. 969 2014;111(11):2374-2391. 970 16. Thapaliya B, Akbas E, Chen J, et al. Brain networks and intelligence: A graph neural network based 971 approach to resting state fmri data. Medical Image Analysis. 2025;101:103433. 972 17. Thapaliya B, Wu Z, Sapkota R, et al. Graph-based deep learning models in the prediction of early-973 stage Alzheimers. Paper presented at: 2024 46th Annual International Conference of the IEEE 974 Engineering in Medicine and Biology Society (EMBC)2024. 975 18. Thapaliya B, Miller R, Chen J, et al. DSAM: A deep learning framework for analyzing temporal and 976 spatial dynamics in brain networks. Medical Image Analysis. 2025:103462. 977 19. Sakoe H, Chiba S. Dynamic programming algorithm optimization for spoken word recognition. 978 IEEE Transactions on Acoustics, Speech, and Signal Processing. 1978;26(1):43-49. 979 20. Wiafe S-L, Kinsey S, Iraji A, Miller R, Calhoun VD. Normalized Dynamic Time Warping Increases 980 Sensitivity In Differentiating Functional Network Connectivity In Schizophrenia. bioRxiv. 981 2024:2024.2010.2031.621415. 982 21. Meszlényi RJ, Hermann P, Buza K, Gál V, Vidnyánszky Z. Resting State fMRI Functional Connectivity 983 Analysis Using Dynamic Time Warping. Frontiers in Neuroscience. 2017;11. 984 22. Linke AC, Mash LE, Fong CH, et al. Dynamic time warping outperforms Pearson correlation in 985 detecting atypical functional connectivity in autism spectrum disorders. NeuroImage. 986 2020;223:117383. 987 23. Wiafe S-L, Faghiri A, Fu Z, Miller R, Calhoun V. Capturing Stretching and Shrinking of Inter-Network 988 Temporal Coupling in FMRI Via WARP Elasticity. 2024. 989 24. Wiafe S-L, Faghiri A, Fu Z, Miller R, Preda A, Calhoun VD. The dynamics of dynamic time warping 990 in fMRI data: a method to capture inter -network stretching and shrinking via warp elasticity. 991 Imaging Neuroscience. 2024. 992 25. Herrmann M, Tan CW, Webb GI. Parameterizing the cost function of dynamic time warping with 993 application to time series classification. Data Mining and Knowledge Discovery. 2023;37(5):2024-994 2045. 995 26. Fukunaga M, Horovitz SG, de Zwart JA, et al. Metabolic origin of BOLD signal fluctuations in the 996 absence of stimuli. J Cereb Blood Flow Metab. 2008;28(7):1377-1387. 997 27. Roosterman D, Cottrell GS. The two -cell model of glucose metabolism: a hypothesis of 998 schizophrenia. Molecular Psychiatry. 2021;26(6):1738-1747. 999 28. Sarnyai Z, Ben -Shachar D. Schizophrenia, a disease of impaired dynamic metabolic flexibility: A 1000 new mechanistic framework. Psychiatry Research. 2024;342:116220. 1001 29. Helaly AMN, Ghorab D. Schizophrenia as metabolic disease. What are the causes? Metab Brain 1002 Dis. 2023;38(3):795-804. 1003 30. Gohel SR, Biswal BB. Functional integration between brain regions at rest occurs in multiple -1004 frequency bands. Brain Connect. 2015;5(1):23-34. 1005 31. Birn RM, Cornejo MD, Molloy EK, et al. The influence of physiological noise correction on test -1006 retest reliability of resting-state functional connectivity. Brain Connect. 2014;4(7):511-522. 1007 32. Lindauer U, Leithner C, Kaasch H, et al. Neurovascular coupling in rat brain operates independent 1008 of hemoglobin deoxygenation. Journal of cerebral blood flow & metabolism. 2010;30(4):757-768. 1009 33. Wolf T, Lindauer U, Villringer A, Dirnagl U. Excessive oxygen or glucose supply does not alter the 1010 blood flow response to somatosensory stimulation or spreading depression in rats. Brain research. 1011 1997;761(2):290-299. 1012 34. Brown AM, Ransom BR. Astrocyte glycogen and brain energy metabolism. Glia. 2007;55(12):1263-1013 1271. 1014 35. Pellerin L, Bouzier‐Sore AK, Aubert A, et al. Activity‐dependent regulation of energy metabolism 1015 by astrocytes: an update. Glia. 2007;55(12):1251-1262. 1016 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 38 36. Hillman EM. Coupling mechanism and significance of the BOLD signal: a status report. Annu Rev 1017 Neurosci. 2014;37:161-181. 1018 37. Dwyer DS, Bradley RJ, Kablinger AS, Freeman AM, 3rd. Glucose metabolism in relation to 1019 schizophrenia and antipsychotic drug treatment. Ann Clin Psychiatry. 2001;13(2):103-113. 1020 38. Newcomer JW, Haupt DW, Fucetola R, et al. Abnormalities in Glucose Regulation During 1021 Antipsychotic Treatment of Schizophrenia. Archives of General Psychiatry. 2002;59(4):337-345. 1022 39. Gatti D, Rinaldi L, Ferreri L, Vecchi T. The Human Cerebellum as a Hub of the Predictive Brain. 1023 Brain Sci. 2021;11(11). 1024 40. Schmahmann JD. An emerging concept: the cerebellar contribution to higher function. Archives 1025 of neurology. 1991;48(11):1178-1187. 1026 41. Andreasen NC, Paradiso S, O'Leary DS. “Cognitive dysmetria” as an integrative theory of 1027 schizophrenia: a dysfunction in cortical -subcortical-cerebellar circuitry? Schizophrenia bulletin. 1028 1998;24(2):203-218. 1029 42. Speers LJ, Bilkey DK. Maladaptive explore/exploit trade -offs in schizophrenia. Trends in 1030 Neurosciences. 2023;46(5):341-354. 1031 43. Millan MJ, Fone K, Steckler T, Horan WP. Negative symptoms of schizophrenia: Clinical 1032 characteristics, pathophysiological substrates, experimental models and prospects for improved 1033 treatment. European Neuropsychopharmacology. 2014;24(5):645-692. 1034 44. Gu X, Ke S, Wang Q, et al. Energy metabolism in major depressive disorder: Recent advances from 1035 omics technologies and imaging. Biomedicine & Pharmacotherapy. 2021;141:111869. 1036 45. Meng F, Wang J, Wang L, Zou W. Glucose metabolism impairment in major depressive disorder. 1037 Brain Research Bulletin. 2025;221:111191. 1038 46. Wu Q, Long Y, Peng X, et al. Prefrontal cortical dopamine deficit may cause impaired glucose 1039 metabolism in schizophrenia. Translational Psychiatry. 2024;14(1):79. 1040 47. Biswal B, Yetkin FZ, Haughton VM, Hyde JS. Functional connectivity in the motor cortex of resting 1041 human brain using echo-planar MRI. Magn Reson Med. 1995;34(4):537-541. 1042 48. DeRamus T, Faghiri A, Iraji A, et al. Modular and state-relevant functional network connectivity in 1043 high-frequency eyes open vs eyes closed resting fMRI data. Journal of Neuroscience Methods. 1044 2021;358:109202. 1045 49. Faghiri A, Iraji A, Damaraju E, Turner J, Calhoun VD. A unified approach for characterizing 1046 static/dynamic connectivity frequency profiles using filter banks. Network Neuroscience. 1047 2021;5(1):56-82. 1048 50. Jamadar SD, Behler A, Deery H, Breakspear M. The metabolic costs of cognition. Trends in 1049 Cognitive Sciences. 1050 51. Gonder-Frederick LA, Zrebiec JF, Bauchowitz AU, et al. Cognitive function is disrupted by both 1051 hypo- and hyperglycemia in school -aged children with type 1 diabetes: a field study. Diabetes 1052 Care. 2009;32(6):1001-1006. 1053 52. Cocchi L, Gollo LL, Zalesky A, Breakspear M. Criticality in the brain: A synthesis of neurobiology, 1054 models and cognition. Progress in Neurobiology. 2017;158:132-152. 1055 53. Moretti P, Muñoz MA. Griffiths phases and the stretching of criticality in brain networks. Nature 1056 communications. 2013;4(1):2521. 1057 54. Shi J, Kirihara K, Tada M, et al. Criticality in the Healthy Brain. Front Netw Physiol. 2021;1:755685. 1058 55. Carhart-Harris RL. The entropic brain - revisited. Neuropharmacology. 2018;142:167-178. 1059 56. Xu L, Feng J, Yu L. Avalanche criticality in individuals, fluid intelligence, and working memory. Hum 1060 Brain Mapp. 2022;43(8):2534-2553. 1061 57. Roheger M, Kessler J, Kalbe E. Structured Cognitive Training Yields Best Results in Healthy Older 1062 Adults, and Their ApoE4 State and Baseline Cognitive Level Predict Training Benefits. Cogn Behav 1063 Neurol. 2019;32(2):76-86. 1064 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint 39 58. Friston KJ. The disconnection hypothesis. Schizophrenia research. 1998;30(2):115-125. 1065 59. Paliwal KK, Agarwal A, Sinha SS. A modification over Sakoe and Chiba's dynamic time warping 1066 algorithm for isolated word recognition. Signal Processing. 1982;4(4):329-333. 1067 60. Berndt DJ, Clifford J. Using Dynamic Time Warping to Find Patterns in Time Series. Paper 1068 presented at: KDD Workshop1994. 1069 61. Faghiri A, Yang K, Ishizuka K, Sawa A, Adali T, Calhoun V. Frequency modulation increases the 1070 specificity of time-resolved connectivity: A resting-state fMRI study. 2023. 1071 62. Fritsch FN, Carlson RE. Monotone piecewise cubic interpolation. SIAM Journal on Numerical 1072 Analysis. 1980;17(2):238-246. 1073 63. Keator DB, van Erp TGM, Turner JA, et al. The Function Biomedical Informatics Research Network 1074 Data Repository. Neuroimage. 2016;124(Pt B):1074-1079. 1075 64. Van Essen DC, Smith SM, Barch DM, Behrens TEJ, Yacoub E, Ugurbil K. The WU -Minn Human 1076 Connectome Project: An overview. NeuroImage. 2013;80:62-79. 1077 65. Van Essen DC, Ugurbil K, Auerbach E, et al. The Human Connectome Project: a data acquisition 1078 perspective. Neuroimage. 2012;62(4):2222-2231. 1079 66. Penny WD, Friston KJ, Ashburner JT, Kiebel SJ, Nichols TE. Statistical parametric mapping: the 1080 analysis of functional brain images. Elsevier; 2011. 1081 67. Fu Z, Batta I, Wu L, et al. Searching Reproducible Brain Features using NeuroMark: Templates for 1082 Different Age Populations and Imaging Modalities. Neuroimage. 2024;292:120617. 1083 68. Du Y, Fu Z, Sui J, et al. NeuroMark: An automated and adaptive ICA based pipeline to identify 1084 reproducible fMRI markers of brain disorders. Neuroimage Clin. 2020;28:102375. 1085 69. Wang Z, Bridgeford E, Wang S, Vogelstein JT, Caffo B. Statistical analysis of data repeatability 1086 measures. International Statistical Review. 2024. 1087 70. Murphy K, Birn RM, Bandettini PA. Resting -state fMRI confounds and cleanup. Neuroimage. 1088 2013;80:349-359. 1089 71. Liégeois R, Yeo BTT, Van De Ville D. Interpreting null models of resting -state functional MRI 1090 dynamics: not throwing the model out with the hypothesis. NeuroImage. 2021;243:118518. 1091 72. Damaraju E, Allen EA, Belger A, et al. Dynamic functional connectivity analysis reveals transient 1092 states of dysconnectivity in schizophrenia. NeuroImage: Clinical. 2014;5:298-308. 1093 73. Aggarwal CC, Hinneburg A, Keim DA. On the surprising behavior of distance metrics in high 1094 dimensional space. Paper presented at: Database Theory —ICDT 2001: 8th International 1095 Conference London, UK, January 4–6, 2001 Proceedings 82001. 1096 74. Levin DA, Peres Y. Markov chains and mixing times. Vol 107: American Mathematical Soc.; 2017. 1097 75. Basharin GP, Langville AN, Naumov VA. The life and work of A.A. Markov. Linear Algebra and its 1098 Applications. 2004;386:3-26. 1099 76. Iraji A, Faghiri A, Lewis N, Fu Z, Rachakonda S, Calhoun VD. Tools of the trade: estimating time -1100 varying connectivity patterns from fMRI data. Social Cognitive and Affective Neuroscience. 1101 2020;16(8):849-874. 1102 1103 1104 .CC-BY-NC-ND 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted March 21, 2025. ; https://doi.org/10.1101/2025.03.20.644399doi: bioRxiv preprint