Criticality and increased intrinsic neural timescales in stroke

TL;DR

This study shows that stroke causes lasting changes in brain activity timing, which could help predict recovery and guide treatment.

Contribution

The study introduces intrinsic neural timescales as potential biomarkers for stroke recovery and links them to criticality in brain dynamics.

Findings

Stroke patients showed prolonged intrinsic neural timescales in multiple cortical regions.

Dynamic changes in brain timing were more severe in patients with poor recovery outcomes.

Computational models suggest stroke-induced changes reflect a shift toward criticality in brain dynamics.

Abstract

Stroke disrupts brain function beyond focal lesions, altering multiscale temporal dynamics essential for information processing. We investigated intrinsic neural timescales (INT) and other properties of long-range temporal correlations, using longitudinal fMRI data from 15 ischemic stroke patients across 6 months, and compared them to age-matched controls. Results show that stroke patients exhibited significantly prolonged INT in multiple cortical regions, reflecting slowed temporal dynamics and disrupted hierarchy. These dynamic changes persisted through recovery and were more pronounced in patients with poor outcomes, especially within cognitive control networks. Computational modeling suggested that stroke-induced INT prolongation driven by heightened neuronal excitability reflects a dynamic shift towards criticality. Our findings position long-range temporal correlations and INT as potential biomarkers for monitoring and predicting functional recovery. This framework provides a novel perspective on stroke-induced brain changes and suggests avenues for targeted neurorehabilitation using interventions aiming at restoring intrinsic temporal dynamics.

Figures

One-month post-stroke analysis of intrinsic neural timescales.A Resting-state fMRI data were collected from a cohort of 15 stroke patients and age-matched healthy controls. Recruited stroke patients (*n* = 15) underwent five follow-up fMRI scans approximately every 30–40 days over 6 months post-stroke. B Brain parcellation comprising 32 regions of interest (ROIs) organized into 8 functional networks—Default Mode Network (DMN), Sensorimotor Network (SMN), Visual Network (VIS), Salience Network (SAN), Language Network (LN), Frontoparietal Network (FPN), Dorsal Attention Network (DAN), and Cerebe

A. The corresponding autocorrelation curves, which are used to compute the INT, are shown in Supplementary Material Fig. 2B. As shown in Fig. [4A](#Fig4), INT and overall network activity (see Supplementary Material Fig. 2C) are directly shaped by $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document}$ . Notably, intrinsic timescales increase sharply as the network approaches the critical point ( $

A, B) and the Hurst exponent (Supplementary Fig. 3C, D) from both the raw neuronal dynamics and the HRF-transformed BOLD signals. Although the HRF convolution introduced temporal smoothing and added noise, a peak in both INT and Hurst exponent remained clearly visible near *σ* = 1 (the critical point). This indicates that the criticality-related modulation of timescales survives the nonlinear transformation into the BOLD domain. While this transformation is not perfect, the inclusion of these HRF-based results (Supplementary Fig. 3) provides further support that our model captures relevant tem

Full text

Introduction

The human brain operates through complex dynamics that enable flexible cognition, efficient information processing, and adaptive responses to the environment^1–5^. These dynamics emerge from interactions between neural populations and are crucial for maintaining a balance between stability and flexibility in brain function^6^. As a fundamental property of the dynamics of brain regions, intrinsic neural timescales (INT) reflect the temporal window over which neural activity is integrated^2,3,7,8^. Alterations in INT have been observed across various conditions, including temporal lobe epilepsy^9,10^, Alzheimer’s disease^11^, and Parkinson’s disease^12^, highlighting its potential as a biomarker for disrupted neural dynamics and functional impairments. INT provides valuable insights into how neural architecture shapes functional behavior and information processing^8,13^, and represents a critical dynamical feature linking brain structure and function. Altered INT are often related^14^ or caused by alterations in gray matter volume^15^, reflecting structural and functional changes in the brain, thereby serving as a bridge between anatomical alterations and dynamic neural processes.

Despite these recent advances in our understanding of INT in different conditions, its role in stroke remains unexplored. As stroke typically causes substantial structural damage to the brain, INT can be a relevant tool to reveal critical changes in brain dynamics. Stroke can result in localized lesions and disruptions to white matter integrity, leading to widespread functional impairments, including deficits in motor control, cognitive processes, and sensory integration^16–18^. While significant progress has been made in understanding post-stroke atypical temporal dynamics in neural activity through studies of network reorganization^19–21^, functional connectivity^22–24^, and network efficiency^25^, the specific impact of stroke lesions on INT and their role in brain function and recovery remains unclear.

Furthermore, INT is derived from the autocorrelation function (ACF) of a neural signal and is measured by its decay properties^3,15,26,27^. A slow decay of ACF typically occurs at a critical point and is known in the literature as critical slowing down^6,28–31^. In practical terms, this means that INT is maximized at a critical point^15^ and various computational properties are also optimized around the critical state^6,28^. This concept of criticality provides a powerful framework for linking INT to distinct dynamical states^15^: subcritical, critical, and supercritical. In the subcritical and supercritical states, INT is short due to rapid (exponential) signal decay; at criticality, they persist longer, following a slow decay. Throughout this framework, criticality has been applied to various brain states^32^, including anesthesia^33,34^, epilepsy^35,36^, neurodegeneration^37^, cognition^38,39^, psychiatry^40^, and sleep medicine^41^. However, its relevance to stroke remains unexplored, presenting an opportunity to investigate how stroke lesions disrupt neural criticality and influence functional recovery. Additionally, since long-range temporal correlations (TCs) have been shown to reflect critical brain dynamics^42–44^, we further examined post-stroke brain dynamics by analyzing key TCs properties, including the time delay to reach a specific ACF strength^43,44^ and the Hurst exponent^45^, as complementary measures to INT. We hypothesize that analyzing TCs and INT through the lens of criticality will provide deeper insights into the temporal reorganization of post-stroke brain dynamics and their influence on recovery trajectories.

This longitudinal study examines how stroke disrupts the brain’s multiscale temporal organization by investigating alterations in INT as well as other influential properties of long-range TCs^42–45^ and their hierarchical arrangement across functional networks. Grounded in the criticality framework, we combine longitudinal fMRI with biologically inspired, parsimonious neuronal network modeling to characterize post-stroke temporal reorganization. Empirically, we track INT trajectories during recovery; computationally, we simulate how changes in neuronal excitability reshape brain network dynamics. We test a central hypothesis that stroke generates unique INT signatures that distinguish patients from healthy controls across temporal hierarchies and inform functional recovery outcomes. By bridging empirical findings with computational models of critical dynamics, this work establishes a novel understanding of post-stroke recovery that spans spatial scales, from neuronal excitability to network reconfiguration, while identifying INT as a mechanistic marker of temporal disruption and a potential clinical biomarker for prognosis.

Results

Demographic and clinical characteristics

There were no significant differences in age, sex category, and mean framewise displacement between stroke patients and healthy controls (Table 1). Patients were first scanned on an average of 23.06 days (standard deviation: 4.32 days) after stroke onset. Recruited patients underwent five follow-up scans at evenly spaced intervals of 30–40 days over 6 months post-stroke. Recovery trajectories of patients were quantified using the Brunnstrom stage score^46^. We use 2-stage improvements at the final assessment as a cutoff to distinguish the good (n = 6) and poor (n = 9) recovery subgroups^47,48^. There are no significant differences in age or sex categories between the two groups, while the lesion volume is higher in the poor outcome group (19.43 ± 19.54 cm^3^) than in the good outcome group (3.93 ± 2.24 cm^3^). The details of recovery outcome assessment at each follow-up for both groups can be seen in Supplementary Material Table 2.Table 1. Demographics and clinical informationStroke (n = 15)Healthy controls (n = 15)P-valuesAge, years63.8 (47–81)68.6 (41–81)0.18Sex, Female470.17Mean FD0.033 (0.01–0.05)0.088(0.02-0.11)0.99Days since stroke23.06 (14–42)——Lesion volume (ml)13.01 (0.91, 63.55)——NIHSS7.26 (1–20)——Values are presented as mean (range). Age and mean framewise displacement values were compared by means of two-sided t-tests, and sex (the number of females) by means of a Pearson’s chi-squared test.FD framewise displacement, NIHSS National Institutes of Health Stroke Scale.

Abnormal INT following stroke

One month after the stroke, stroke patients showed significantly longer INTs compared to healthy controls, both globally (Fig. 1D) and within each specific functional network (Fig. 1E). The global average INT was notably higher in stroke patients (t = 8.23, p < 0.0001, Welch’s t-test, FDR corrected), highlighting widespread alterations in the temporal dynamics of brain activity. When examining specific functional networks, stroke patients showed significant increases in INT across all networks, as seen in the statistics and p-values in Table 2.Fig. 1. One-month post-stroke analysis of intrinsic neural timescales.A Resting-state fMRI data were collected from a cohort of 15 stroke patients and age-matched healthy controls. Recruited stroke patients (n = 15) underwent five follow-up fMRI scans approximately every 30–40 days over 6 months post-stroke. B Brain parcellation comprising 32 regions of interest (ROIs) organized into 8 functional networks—Default Mode Network (DMN), Sensorimotor Network (SMN), Visual Network (VIS), Salience Network (SAN), Language Network (LN), Frontoparietal Network (FPN), Dorsal Attention Network (DAN), and Cerebellar Network (CE)—was utilized to extract regional BOLD time series. C Intrinsic neural timescales were estimated by calculating the area under the curve in which the autocorrelation function (ACF) is positive. Compared to healthy controls, 1-month post-stroke patients (PAT1) exhibited significantly longer INT both globally, for the whole brain (D, t = 8.23, p < 0.0001, Welch’s t-test, FDR corrected) and within specific functional networks (E, see Table 1 for detailed statistics). Each dot represents one subject. **** indicates p < 0.0001 (FDR-corrected). F Hierarchy of intrinsic neural timescales are also disrupted after a stroke, compared to healthy controls (Spearman’s rho = −0.90, p < 0.0001). The upper and lower error bars display the largest and smallest values within 1.5 times IR above the 75th percentile and below the 25th percentile, respectively. Abbreviations: INT intrinsic neural timescales, PAT1 15 patients recruited and scanned for the first time (first follow-up, approximately 20–30 days post-stroke). CTR: the 15 healthy controls. Data presented in (D–F) are derived specifically from the first follow-up (PAT1).Table 2. Comparison of INT across all networks between healthy controls and stroke patientsNetworkT-statistics (CTR < PAT)p-valuesDMN5.92<0.0001SMN5.27<0.0001VIS3.93<0.0001DAN6.82<0.0001FPN6.60<0.0001LN6.72<0.0001CE4.28<0.0001P-values are FDR corrected.DMN default mode network, SMN sensorimotor network, VIS visual network, SAN salience network, LN language network, FPN frontoparietal network, DAN dorsal attention network, CE cerebellar network.

Reconfigured hierarchical INT in post-stroke patients

In addition to the global and network-specific abnormalities in INT, stroke patients exhibited a significant disruption in the hierarchical organization of INT. Compared with healthy controls (Fig. 1F), we found the gradients of mean INT of functional networks are significantly disrupted following stroke. In particular, the CE functional network exhibited the shortest INT (Mean  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm \,$$\end{document} Std: 1.422 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm \,$$\end{document} 0.257) in the control group and the longest INT at PAT1, which corresponds to 1 month after stroke (Mean  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm \,$$\end{document} Std: 8.178 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm$$\end{document} 5.901). Furthermore, the default mode network (DMN), which had long INT for controls, exhibited the shortest INT in stroke patients (Mean  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm \,$$\end{document} Std: 6.068, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm \,$$\end{document} 2.875). Despite this relative decrease within the stroke group, the INT of the DMN remained longer than that of the controls. Considering all functional networks, the gradient of INT that represents a hierarchical organization of brain dynamics was severely disrupted in stroke patients at 1-month post-stroke (Spearman’s rho = −0.90, p < 0.0001), as shown in Fig. 1F.

Development of INT within 6 months after stroke

The global average INT was compared between healthy controls and stroke patients across different time points post-stroke to examine the development within a 6-month period post-stroke (Fig. 2A). Repeated-measures one-way ANOVA analysis revealed significant differences in INT (F5,70 = 13.05, p < 0.0001, FDR corrected). The post-hoc t-test indicates that the prolonged INT observed in stroke patients persists for at least 5 months after the stroke attack (tCTR-PAT1 = 8.22, p < 0.0001; tCTR-PAT2 = 10.80, p < 0.0001; tCTR-PAT3 = 5.94, p < 0.0001; tCTR-PAT4 = 7.61, p < 0.0001; tCTR-PAT5 = 5.80, p < 0.0001). This enduring elevation in INT suggests a long-lasting disruption in the brain’s temporal dynamics following a stroke, as depicted in the global INT comparisons. A whole-brain map of the prolonged INT in stroke patients at different time points further highlights the widespread nature of these changes (the example of PAT1 can be seen in Fig. 2B. Other measures of long-range TCs, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.1}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.5}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{0}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{0.1}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{0.5}$$\end{document} , further support analogous effects of stroke lesions on the brain’s temporal dynamics (FHurst (5,70) = 26.11, p < 0.0001, FINT01(5,70) = 9.73, p < 0.0001, FINT05(5,70) = 2.69, p < 0.0001, FACF0(5,70) = 26.98, p < 0.0001, FACF01(5,70) = 15.82, p < 0.0001, FACF05(5,70) = 3.38, p < 0.01, FRD corrected. See Fig. 3).Fig. 2. Impact of stroke on intrinsic neural timescales.A Prolonged INT observed in stroke patients persists up to 5 months post-stroke (tCTR-PAT1 = 8.22, p < 0.0001; tCTR-PAT2 = 10.80, p < 0.0001; tCTR-PAT3 = 5.94, p < 0.0001; tCTR-PAT4 = 7.61, p < 0.0001; tCTR-PAT5 = 5.80, p < 0.0001, FDR corrected). B Whole_-brain map showing prolonged INT in stroke patients at 1-month post-stroke (PAT1). C At 2 months post-stroke (PAT2), significant differences in INT emerged between patients with good and poor recovery outcomes (tpoor-good_ = 2.30, p = 0.03, Welch’s t-test). D Network-level differences in INT between good and poor outcome groups at PAT2. The upper and lower error bars display the largest and smallest values within 1.5 times IR above the 75th percentile and below the 25th percentile, respectively. **** indicates p < 0.0001 (FDR-corrected); ** indicates p < 0.05 (FDR-corrected); * indicates p < 0.05 (FDR-corrected); ns indicates no significance.Fig. 3. Longitudinal effects of long-range temporal correlation measurements.A Hurst exponent, B \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.1}$$\end{document} , C \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.5}$$\end{document} , D \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACW}}_{0}$$\end{document} , E \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACW}}_{0.1}$$\end{document} and F \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACW}}_{0.5}$$\end{document} sustain the stroke lesions’ disruption on the brain’s temporal dynamics. FHurst (5,70) = 26.11, p < 0.0001, FINT01(5,70) = 9.73, p < 0.0001, FINT05(5,70) = 2.69, p < 0.0001, FACF0(5,70) = 26.98, p < 0.0001, FACF01(5,70) = 15.82, p < 0.0001, FACF05(5,70) = 3.38, p < 0.01, FRD corrected. The upper and lower error bars display the largest and smallest values within 1.5 times IR above the 75th percentile and below the 25th percentile, respectively. **** indicates p < 0.0001 (FDR-corrected); *** indicates p < 0.001 (FDR-corrected);** indicates p < 0.05 (FDR-corrected); * indicates p < 0.05 (FDR-corrected).

INT as a prognostic biomarker for recovery

To further explore recovery dynamics, patients were categorized into two groups according to motor function outcomes at approximately six months post-stroke: a good and a poor recovery group. At 2 months post-stroke (PAT2), a significant difference in INT emerged between the two groups (Fig. 2C), with the poor recovery group exhibiting significantly longer INT than the good recovery group (PAT2: tpoor-good = 2.30, p = 0.03, FDR corrected), while at other timepoints, no significant differences were observed. This finding suggests that differences in INT at 2 months post-stroke may serve as an early predictor of recovery outcomes 6 months later, highlighting the potential utility of INT as a prognostic biomarker.

A two-way ANOVA was then used at PAT2 to test whether the recovery outcomes (good vs. poor) and the functional networks (DMN, SMN, DAN, etc.) have an effect on INT. We found a statistically significant difference in INT by both outcomes (Frecovery = 11.99, p < 0.001, FDR corrected) and by various function networks (Fnetworks = 3.83, p < 0.001, FDR corrected), though the interaction between these terms was not significant. Post-hoc t-test analysis revealed network-specific significant differences in INT between the pairwise recovery groups (Fig. 2D). Patients with good recovery outcomes exhibited notably shorter INT compared to those with poor recovery outcomes in the salience network (SAN; tgood-poor = −3.10, p = 0.016), dorsal attention network (DAN; tgood-poor = −5.96, p = 0.006) and language network (LN; tgood-poor = −4.31, p = 0.031). These findings indicate that lower INT within specific functional networks is associated with better recovery trajectories, further emphasizing the role of INT as a network-level biomarker for predicting post-stroke recovery outcomes.

Distance to criticality explains increased INT in stroke

Finally, we interpret the abnormal INT in stroke within the criticality framework. By fixing the network mean degree ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K$$\end{document} ) and varying the propagation probability (λ), we control the branching ratio ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} ), a key metric that governs neural activity spread and serves as an indicator of criticality. The activity is defined as the instantaneous density of active neurons as a function of time. The examples of neuronal network activities simulated with various \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document}  = [0.08, 0.1, 0.12] with fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=10$$\end{document} , which represent the subcritical, critical, and supercritical, can be seen in Supplementary Material Fig. 2A. The corresponding autocorrelation curves, which are used to compute the INT, are shown in Supplementary Material Fig. 2B. As shown in Fig. 4A, INT and overall network activity (see Supplementary Material Fig. 2C) are directly shaped by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} . Notably, intrinsic timescales increase sharply as the network approaches the critical point ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma =1$$\end{document} ), a hallmark of critical slowing down (see Supplementary Material Fig. 2D for the single trial relationship between the INT and criticality).Fig. 4. Computation modeling explains how stroke lesions prolong intrinsic neural timescales and alter network dynamics.A Intrinsic neural timescale is shaped by the branching ratio ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} ), indicating that stroke can shift the brain dynamics from a slightly subcritical state (blue) toward criticality (red), with the potential to enter a supercritical state. Near a phase transition, cortical network dynamics can be modeled as a branching process, where intrinsic neural timescales peak at the critical point^15^. B–G Similarly, complementary measures (see “Methods” for details)—including the Hurst exponent, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.1}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.5}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{0}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{0.1}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{0.5}$$\end{document} —consistently peak at the critical point, reflecting signs of critical slowing down, further demonstrating the robustness and convergence of our findings across multiple temporal metrics. This result is an average of 50 trials. Each dot represents the intrinsic neural timescale (INT) or temporal correlations (TCs) value obtained from the computational model at a given branching ratio ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} ). These values were generated by varying the synaptic propagation probability ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{\lambda }}$$\end{document} ) while keeping network connectivity constant. All networks \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=\mathrm{100,000}$$\end{document} neurons and a mean degree \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=10$$\end{document} with varied values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document} to satisfy the relationship \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{\sigma }}=K\cdot {\rm{\lambda }}$$\end{document} . The external driving is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r={10}^{-5}$$\end{document} . See “Methods” for details.

To further establish the correspondence between our computational model and fMRI signals, we extended the neural simulations and transformed the resulting activity into BOLD-like time series using the canonical hemodynamic response function^49^ (HRF; TR = 2 s). Specifically, we simulated 300 s (5 min; time step = 0.001 s) of neuronal activity at different σ values and obtained 150 BOLD-equivalent samples via HRF convolution and downsampling (1 TR). We then computed (INT; Supplementary Fig. 3A, B) and the Hurst exponent (Supplementary Fig. 3C, D) from both the raw neuronal dynamics and the HRF-transformed BOLD signals. Although the HRF convolution introduced temporal smoothing and added noise, a peak in both INT and Hurst exponent remained clearly visible near σ = 1 (the critical point). This indicates that the criticality-related modulation of timescales survives the nonlinear transformation into the BOLD domain. While this transformation is not perfect, the inclusion of these HRF-based results (Supplementary Fig. 3) provides further support that our model captures relevant temporal features and helps bridge the gap between fast (milliseconds) neuronal-level dynamics and slow (seconds) empirical fMRI observations. Hence, when linked to the prolonged INT observed in stroke patients, these results suggest that stroke increases neural excitability (reflected by λ), shifting brain dynamics closer to criticality and reducing the distance to criticality (DTC). Complementary measures related to long-range TCs—including the Hurst exponent, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.1}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.5},$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{0}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{0.1}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{0.5}$$\end{document} —further support this interpretation, consistently indicating critical slowing down across multiple estimates (Fig. 4B–G). This shift renders post-stroke neural networks more susceptible to prolonged activity fluctuations, potentially pushing the system into a supercritical state. Such a transition may drive maladaptive neuroplasticity, contributing to dysfunctional network reorganization and impairing recovery trajectories in stroke patients.

Discussion

This study investigates the impact of stroke on INT and long-range TCs as well as their dynamics during recovery, by analyzing functional neuroimaging data from ischemic stroke patients (ISP) over a 6-month post-stroke period. Our findings reveal that stroke induces persistent disruptions in INT and other long-range TCs, both globally and within specific functional networks, reflecting altered temporal dynamics and hierarchical organization. Notably, differences in INT between patients with good and poor recovery outcomes emerge at 2 months post-stroke, highlighting the potential of INT as a predictor of long-term recovery. The observation of significantly shorter INT in patients with good outcomes further suggests that post-stroke recovery may involve the normalization of INT. These insights advance our understanding of post-stroke brain dynamics and underscore the potential of INT as a biomarker for recovery trajectories.

Changes in INT have been observed across the lifespan, with longer INT in young adults compared to elderly individuals^15^. Atypical INT patterns have also been widely documented in various brain disorders, including shorter INT in sensory regions in autism^14^, mixed patterns in schizophrenia (shorter in parietal and occipital regions^50^ but longer in self-referential processes^51^, and shorter INT in Alzheimer’s disease^11^ and temporal lobe epilepsy^9,10^. Additionally, significant differences in INT have been reported in Parkinson’s disease, with shorter INT in late-stage and longer INT in early-stage patients^12^. Building on this foundation, this study investigates how stroke affects INT and its dynamics during recovery. Our results reveal that stroke patients exhibit significantly elevated INT and long-range TCs compared to healthy controls, with this increase persisting throughout the 5-month follow-up period. The observed whole-brain prolongation of INT in stroke patients relative to controls was highly robust and consistent across multiple complementary measures of TCs and the Hurst exponent. All of these measures converged on the same pattern of substantial INT increases in stroke patients. These converging results, together with supporting neurophysiological^52,53^ and modeling evidence^54–56^ of post-stroke hyperexcitability^57–60^, indicate that the observed INT prolongation is unlikely to be an artifact of the analysis approach and instead reflects genuine alterations in brain dynamics after stroke. The persistent increase in INT suggests enduring disruptions in the temporal integration of neural activity^7,8^, likely reflecting widespread dysfunction in brain network dynamics^19–21^. The consistent elevation of INT across key functional networks—including the DMN, salience network (SAN), and frontoparietal network (FPN)—indicates a global impact of stroke. This disruption may impair functional specialization^61,62^ and network efficiency^25^, contributing to the cognitive and motor deficits commonly observed in stroke patients^61,63^.

In addition to elevated INT, our study revealed a significant disruption in the hierarchical organization of INT in stroke patients. The hierarchy INT is consistently observed across both small^13,64,65^ and large-scale fMRI datasets^4,5^. Unimodal regions, such as sensory and motor networks, typically exhibit shorter INT, whereas transmodal regions, including higher-order networks like the central-executive network (CEN), DAN, and DMN, tend to display longer INT^5,7,64,66,67^. Our finding of INT in healthy individuals supports and extends the hierarchical organization of INT, demonstrating that temporal dynamics across networks are balanced and optimized for efficient information processing and integration^7,13^. It is worth noting that our network definition was derived from a group ICA map^68,69^, resulting in a relatively coarse parcellation in which large-scale networks such as SMN and VIS may encompass both lower- and higher-order subregions. This spatial averaging across heterogeneous regions can influence the apparent INT hierarchy, potentially explaining why our control group exhibited longer INTs in SMN and VIS compared to the canonical hierarchy reported in previous literature^3,7,70^. Additional methodological factors, such as reliance on resting-state fMRI alone, specific preprocessing choices, and age-matching of the control cohort to the stroke group, may also contribute to these differences. Importantly, in healthy controls, INT hierarchy showed the expected systematic ordering across most networks^5,7,15^, and this organization was strongly altered in stroke patients as indicated by a marked negative correlation in the hierarchy relative to controls (Spearman’s rho = −0.90, p < 0.0001). The cerebellum network, whose function is associated with motor coordination and learning, including balance and posture, as well as other cognitive functions^71^, was affected the most. This likely reflects the importance of the cerebellum network in the reconfiguration of brain dynamics involving motor and other cognitive functions. In a similar vein, during this reconfiguration, the DMN, which typically exhibits long INT^13^ and stable dynamics ^44^^,^^72^, also demonstrated significant changes in its dynamics, showing the shortest INT and the least stability among functional networks in stroke patients^73^.

A key finding of this study is the divergence in INT trajectories based on recovery outcomes^16,48,63,73^. At 2 months post-stroke, patients with poor recovery exhibited significantly longer INT than those with good recovery. This suggests that early-stage INT measurements could serve as a predictive marker for long-term outcomes. In addition, this finding also led to the hypothesis that good recovery could be associated with a normalization of the altered INT, i.e., a network reduction in the magnitude of INT (for example, salient network and DAN), whereas poor recovery could be associated with an inability to promptly recover the abnormal INT to normal levels. These findings align with the potential of the temporal profile of the responses to transcranial magnetic stimulation (TMS) observed in EEG signals following stroke^74^. Nevertheless, noted that even among these patients, INT values did not return to the levels observed in healthy controls. This suggests that full recovery of intrinsic temporal dynamics may unfold over a much longer timescale than 6 months. Stroke-induced disruptions to excitability and network topology may take extended periods to stabilize and reorganize, especially within higher-order transmodal networks. Therefore, longer-term longitudinal studies are essential to fully characterize the trajectory of INT recovery and validate its prognostic utility in clinical practice.

Our findings provide new insights into the mechanisms underlying altered INT after stroke, integrating empirical evidence with computational modeling. Post-stroke increases in INT observed in our data may reflect both pathological hyperexcitability^57–60^ and adaptive network reorganization^19,21^. From a neurophysiological perspective, stroke can induce pronounced changes in synaptic plasticity—such as long-term potentiation (LTP) and homeostatic adjustments—which modulate network excitability and alter the propagation of neural activity^57^. These plastic mechanisms influence the emergence of network hubs and stabilize global activity patterns, both of which are linked to proximity to criticality. In our modeling framework, the synaptic propagation probability (λ) serves as a proxy for these excitability changes. This idea is associated with previous modeling studies, which have demonstrated that resting-state networks emerge naturally near criticality when large-scale dynamics are constrained by empirical structural connectivity^75^. Also, empirical research has demonstrated that the INTs, derived from the BOLD signal, reflect neural activity integrated over long timescales corresponding to the slow fluctuations (<0.1 Hz) captured by resting-state fMRI^15,75^. Our results extend this principle to the stroke context, suggesting that post-stroke increases in excitability can push the system toward a more critical state, thereby lengthening INT and potentially enhancing integration capacity in reorganized networks. This interpretation fits well with infra-slow scale-free dynamics recently discussed by Ao et al.^76^, in which criticality and scale-free dynamics form a background dynamic state that modulates foreground temporal processing, such as INT. From this perspective, stroke-induced changes in excitability may disrupt the background scale-free regime, thereby altering INT.

Increasing λ in our computational model effectively mimics enhanced recurrent excitation, consistent with recent studies showing that longer INT in the cortex is associated with stronger local recurrent connectivity and elevated baseline activity^77,78^. An increase in λ represents an increase in excitability, which functionally parallels the effect of enhanced local recurrent excitation, thereby sustaining activity and prolonging integration windows. In this sense, stroke-related increase in INT may reflect a shift toward stronger effective recurrent excitation in reorganized networks, consistent with the hyperexcitability documented in post-stroke physiology^57–60^. By varying λ while holding network structure constant, we could systematically shift the model between subcritical, critical, and supercritical regimes, observing that INT peaked at criticality and declined with distance to criticality (DTC). This behavior aligns with theoretical predictions and prior modeling work on critical brain dynamics^15,43^. Notably, such excitability changes may arise not only from acute neurochemical alterations but also from progressive intra-regional structural remodeling following stroke^79,80^, including changes in synaptic reorganization^81^, local recurrent connectivity^82^, and circuit rewiring^83^. These processes, which unfold across several weeks to months, may sustain or even exacerbate effective excitability within local networks, contributing to prolonged INT. This interpretation is consistent with recent findings linking INT to recurrent synaptic structure^77,78^, and highlights the value of INT as a dynamic marker of both physiological state and ongoing plasticity in the context of stroke recovery^56,81,84^. Nevertheless, structural remodeling can also alter the effective structural degree K of the network model, which—together with λ—determines the branching ratio (σ = K·λ) and thus the DTC. Such K-driven shifts can move dynamics toward or away from σ ≈ 1 and thereby modulate INT. In the current computational framework, we varied λ while holding K fixed; future work relaxing this constraint and allowing K to vary should capture the network reorganization more explicitly.

The relationship between stroke and INTs can be understood through the framework of criticality^6,28,32,85,86^, which posits that the brain operates near a transition point between ordered and disordered states to benefit from its optimal computational properties. Shifts in the DTC have been applied to various clinical contexts, including depression^30,87^, schizophrenia^88,89^, epilepsy^90,91^, insomnia^91^, Alzheimer’s disease^92^, Parkinson’s disease^93^, and also states like aging^15,94–97^, anesthesia^98^, meditation^99^, cognitive tasks^100^, sleep homeostasis^101^, and sustained wakefulness^102^. INTs and other properties of long-range TCs serve as a proxy for the brain’s proximity to criticality. Stroke disrupts this balance through mechanisms such as excitotoxicity, inflammation, and maladaptive plasticity, often resulting in hyperexcitability and an increased risk of seizures^103–105^. Post-stroke brains exhibit prolonged long-range TCs compared to healthy controls, indicating a shift toward criticality (see Fig. 4) and potentially into a supercritical state, where pathological dynamics with excessive synchronization can emerge. This shift aligns with the criticality model, in which long-range TCs, commonly used to measure the DTC^106^—peak at criticality and decay in both subcritical and supercritical regimes^15^, linking post-stroke changes in temporal dynamics to network instability. To further contextualize these findings within the framework of brain criticality, the complementary metrics of long-range TCs, including the Hurst exponent and autocorrelation-based indices, were examined^42–44^. These analyses consistently revealed elevated TCs in stroke patients, reinforcing the interpretation that stroke shifts brain dynamics toward a more temporally persistent and potentially critical regime.

Notably, INT is unevenly affected across the brain, causing alterations in the hierarchy of timescales, as some brain regions are more affected than others. This disruption in the hierarchical organization of timescales further highlights how stroke-induced changes can lead to an imbalance in the brain’s functional networks. Thus, the criticality framework provides a unifying perspective on stroke-induced disruptions, translating the complexity of molecular and vascular processes into their net effect on neural excitability. However, variability in stroke location and severity complicates the clinical quantification of deviations from criticality in individual patients. Despite these challenges, criticality remains a valuable conceptual tool for understanding how stroke alters brain network dynamics.

Our findings suggest that the prolonged long-range TCs observed in stroke patients may reflect a maladaptive shift in neural dynamics, driven by excessive network persistence and impaired information processing. The observed association between shorter INT and improved recovery highlights the potential of modulating INT as a therapeutic target for functional restoration^8^. Given the strong link between INT, neuronal excitability, and network interactions, our results indicate non-invasive brain stimulation (NIBS) techniques—such as TMS and transcranial electrical stimulation (TES)—as promising strategies to normalize aberrant INT^1,8,107–112^. Specifically, excitatory stimulation (e.g., anodal transcranial direct current stimulation, high-frequency TMS, intermittent theta burst stimulation, and entrainment by transcranial alternating current stimulation^113^) could be used to reduce abnormally prolonged INT, aligning with the dynamics observed in patients with better recovery outcomes^112^. By targeting influential hub regions or larger brain systems with altered INT, NIBS offers a novel therapeutic avenue to restore balanced network dynamics, enhance functional recovery, and improve post-stroke outcomes. Future research integrating real-time neurophysiological monitoring and computational modeling could further refine these interventions^109^, enabling personalized stimulation protocols tailored to individual neurophysiological profiles. This approach holds significant promise for optimizing recovery trajectories and advancing precision medicine in stroke rehabilitation.

This study provides valuable insights into the neural dynamics of stroke recovery and underscores important avenues for future exploration. The longitudinal design, featuring five follow-ups over 6 months, provides a robust framework for investigating INT in stroke patients and could be extended to refine recovery trajectory models. The sample size (N = 15) provides meaningful insights and establishes a strong foundation for future large-scale studies to further validate and expand upon these findings. Besides, the findings of this study are based on stroke patients with heterogeneous lesions; hence, stronger findings are expected in a more homogeneous group. In terms of modeling, while our approach leveraged the Kinouchi–Copelli framework^114^ to simulate shifts in network excitability via the propagation probability λ, it did not explicitly include inhibitory neuron populations. This simplification allows a clear definition of criticality through the branching ratio, but omits the biological specificity of excitation–inhibition (E/I) dynamics. Given that stroke is known to disrupt inhibitory circuits and GABAergic signaling^115^, future models incorporating distinct excitatory and inhibitory units could offer a more mechanistic understanding of how E/I imbalance contributes to altered INT and recovery potential. Additionally, the focus of the study on fMRI data opens the possibility for integrating multimodal approaches, such as EEG or diffusion-weighted imaging, that could deepen our understanding of the temporal and structural changes associated with INT dynamics. The identification of meaningful patterns of INT across functional networks marks a notable advancement. Future research utilizing advanced techniques, such as machine learning or personalized network analyses, holds significant potential to reveal individual-level variability and improve predictive accuracy. While this study focused on motor recovery, future research could explore the impact of INT and TCs on cognitive performance, particularly in functions associated with the salience, language, and DANs, which showed significant predictive value for recovery outcomes.

Methods

Subjects

Fifteen ISP admitted to the First Affiliated Hospital of Shantou University Medical College were recruited for this longitudinal study. The recruited patients were followed up in 6 months across five timepoints. Inclusion criteria were as follows: (i) The first stroke or previous stroke without sequelae; (ii) The diagnosis was in line with the main points of cerebrovascular disease diagnosis approved by the fourth National Cerebrovascular Disease Academic Conference in 1995; (iii) The onset time was 2 weeks to 6 months, and Brunnstrom stage 0–3 of the affected hand; (iv) No metal implants or pacemakers in the body; (v) No previous history of epilepsy; (vi) Informed consent signed by the patient or his immediate family member. Participants were excluded if they met any of the following conditions: (i) Patients with Parkinson’s disease or other neurological diseases; (ii) Patients with serious diseases of the heart, lung, and other organs; (iii) A history of mental illness, drug abuse, and alcohol abuse; (iv) Patients with an unstable condition, such as blood pressure index: low-pressure <90 mmHg or high-pressure >160 mmHg.

The fifteen patients presented with mild to severe motor function deficits (National Institutes of Health Stroke Scale, NIHSS: mean: 7.26, 1–20). Among these patients, 5 had right-sided, and 10 had left-sided strokes, with an average age of 63.81 years (standard deviation: 11.68 years). The corresponding lesion map of all stroke patients has been attached to Supplementary Material Fig. 1. The cohort consisted of 4 males and 11 females, with the first MRI scan performed on average 23.06 days post-stroke (standard deviation: 4.32 days). Additionally, a control group of 15 age-matched healthy individuals with no history of stroke and a normal neurological examination (7 males, 6 females, mean age: 68.61 years, standard deviation: 6.42 years) was included for comparison. No significant differences were observed between the stroke patients and healthy controls in terms of age (p = 0.18, two-sample t-tests) or sex distribution (p = 0.17, Pearson’s chi-squared test). Demographic features of recruited stroke patients and healthy controls are listed in Table 1.

Each patient underwent five resting-state fMRI scans spaced approximately 30–40 days apart over a 6-month period post-stroke. For clarity, we refer to the first, second, third, fourth, and fifth follow-up scans across all patients as PAT1, PAT2, PAT3, PAT4, and PAT5, respectively. In practice, all fifteen recruited patients finished the follow-up plan. Hence, there are 75 scans for 15 patients and 15 scans of 15 healthy controls for this unique longitudinal dataset (90 MRI scans in total). Given the observational nature of the study, treatment regimens were not standardized across participants, allowing for natural variability in clinical care. Recovery trajectories of patients were quantified using the Brunnstrom stage score^46^. This validated measure of upper limb motor recovery function is assessed across four dimensions. To test the association between INT and post-stroke recovery, patients demonstrating at least a 2-stage improvement by the final assessment were classified as having good recovery, while those showing improvement by one or fewer stages were classified as poor recovery^47,48^. The detailed demographic characteristics and the follow-up clinical features of stroke patients can be seen in Supplementary Material Table 2. There was no statistically significant difference in time since stroke (p = 0.59, permutation test) and size of lesion volumes (p = 0.16, permutation test) between subgroups.

fMRI acquisition, pre-processing, and denoising

Resting-state fMRI data were acquired for all patients at five follow-up visits using a 3.0T Discovery MRI scanner with an 8-channel head coil at the SUMC MRI Center. High-resolution T1-weighted anatomical images were acquired using a multi-planar rapid gradient echo sequence with the following parameters: 129 slices, repetition time (TR) = 2250 ms, echo time (TE) = 4.52 ms. Following the anatomical scan, resting-state functional MRI data were collected using a single-shot gradient-echo EPI sequence: TR = 2000 ms, TE = 30 ms, flip angle = 90°, and voxel size = 3.43 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 3.43 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\times$$\end{document} 5.0 mm^3^ (no gap). A total of 210 volumes were obtained over a 7-min duration for each MRI scan (see Fig. 1A).

The fMRI data underwent preprocessing using a tailored pipeline within the CONN functional connectivity toolbox^116^, integrated with the Statistical Parametric Mapping software (SPM12)^117^. For fMRI from stroke patients, the bias of the hemisphere was eliminated by flipping the right hemispheric lesions to the left along the midsagittal plane (see Supplementary Material Fig. 1C for the left lesions of all patients). For each participant, the first 10 dummy volumes were discarded. The subsequent functional images were corrected for slice timing and head motion (no significant difference in framewise displacement detected, p = 0.99, two-sample t-tests). Outlier detection was performed using Artifact Detection Tools^116^ to identify anomalous time points for each participant. These outliers were removed as covariates, and the remaining function images were then normalized to the Montreal Neurological Institute (MNI) space after masking lesioned tissue.

Non-smoothed functional images were further processed through the default denoising pipeline^116^, which included removing confounding effects and applying temporal band-pass filtering. Potential confounding effects in blood oxygenation level-dependent (BOLD) signal include noise components from (i) cerebral white matter and cerebrospinal areas, (ii) estimated motion parameters identified during the realignment step of pre-processing, (iii) outlier scans from outlier identification (i.e., scrubbing), and (iv) constant and (v) first-order linear session effects. While recent ECoG‑fMRI studies have demonstrated that the global signal contains meaningful neuronal information related to arousal and neuro‑vascular coupling^118,119^, global-signal regression (GSR) can introduce artefactual anticorrelations and potentially remove genuine neural fluctuations^120^. Hence, no global signal regression was applied. We regressed out the five confounds with CompCor algorithm^121^ (a component-based noise correction method for BOLD timeseries) implemented in CONN^116^, and then the residual time series was band-pass filtered in the 0.008–0.09 Hz range^116^.

Head motion control

The head motion effect was controlled by calculating the individual mean and maximum framewise displacements (FD)^122^:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\rm{FD}}}_{{\rm{i}}}=\left|\Delta {{\rm{d}}}_{{\rm{ix}}}\right|+\left|\Delta {{\rm{d}}}_{{\rm{iy}}}\right|+\left|\Delta {{\rm{d}}}_{{\rm{iz}}}\right|+\left|\Delta {{\rm{\alpha }}}_{{\rm{i}}}\right|+\left|\Delta {{\rm{\beta }}}_{{\rm{i}}}\right|+\left|\Delta {{\rm{\gamma }}}_{{\rm{i}}}\right|$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {d}_{{ix}}={d}_{\left(i-1\right){x}}-{d}_{{ix}}$$\end{document} and similarly for other rigid body parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[{d}_{{ix}}\,,{d}_{{iy}}\,,{d}_{{iz}},\,{\alpha }_{i},\,{\beta }_{i},\,{\gamma }_{i}]$$\end{document} . Participants with a maximum displacement exceeding 1.5 mm and a maximum rotation above 1.5° were excluded. While in practice, no FD was above the threshold, and no subjects were excluded. Besides, 24 motion parameters calculated from the six original motion parameters using Volterra expansion^123^ were regressed out as nuisance covariates.

Brain parcellation

To analyse changes in network dynamics resulting from stroke lesions, we utilized a functional brain parcellation derived from CONN’s group independent component analysis (ICA) of the HCP dataset (497 subjects). This parcellation includes 32 regions of interest (ROIs) spanning the entire brain, grouped into eight large-scale functional networks or systems: DMN, sensorimotor network (SMN), visual network (VIS), DAN, SAN, FPN, LN, cerebellar network (CE) (see Supplementary Material Table 2 for the details of eight large-scale functional network and 32 ROI as well as their peak coordinates). ICA spatial maps were applied to each participant’s fMRI BOLD data to extract a representative timeseries for each ROI (see Fig. 1B for parcellation map).

Estimation of intrinsic neural timescales

We estimated intrinsic timescale values for all regional BOLD time series. This involves first computing the ACF for each timeseries:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACF}}_{k}\,=\,\frac{{\sum }_{t\,=\,k+1}^{T}{(y}_{t}-\bar{y})({y}_{t-k}-\bar{y})}{{{\sum }_{t=1}^{T}{(y}_{t}-\bar{y})}^{2}}\,$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y$$\end{document} denotes the preprocessed regional BOLD timeseries, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{y}$$\end{document} is the mean value across time points. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${t}$$\end{document} is the length of time bins, which is the time of repetition (TR = 2000 ms) of the MRI scan, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T$$\end{document} is the number of time points (200 in the experiment).

Due to the lower temporal resolution of fMRI data, the intrinsic timescale was calculated as the area under the curve of the ACF from one to the time lag at which the correlation reaches zero (see Fig. 1C for illustration):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${INT}={TR}\cdot {\sum }_{k=1}^{{n}_{{lag}}-1}{AC}{F}_{k},{AC}{F}_{{n}_{{lag}}} < 0$$\end{document}

It has been proven that the method proposed to compute the INTs for fMRI data was highly correlated with those calculated from simultaneously recorded EEG data^14^. This process was applied to all regional BOLD timeseries of patients and healthy controls, resulting in an intrinsic timescale map for the entire brain network for each participant.

To further validate potential biases in the estimation of INT across the literature, we incorporated additional metrics commonly used in studies of long-range TCs. Specifically, we computed alternative INT measures based on the area under the ACF curve up to defined thresholds: the point at which the ACF first decays to 0.5 (denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.5}$$\end{document} ) and 0.1 ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.1}$$\end{document} )^42,43^. Given the well-established relationship between temporal autocorrelation and the Hurst exponent^124^, we also included the Hurst exponent as a complementary measure of TCs. Furthermore, we quantified ACF strength^43^, defined as the time lag at which the ACF curve first crosses zero after linear interpolation (denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACW}}_{0}$$\end{document} ). To enhance robustness, we additionally computed threshold-based ACF indices—0.5 ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACW}}_{0.5}$$\end{document} ) and 0.1 ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACW}}_{0.1}$$\end{document} )—corresponding to the time lag above these respective threshold values. These supplementary TCs-based measures serve to cross-validate and reinforce our INT-based findings, offering a more comprehensive characterization of post-stroke temporal brain dynamics. See Table 3 for the details about INTs and additional measures of TCs.Table 3. Intrinsic neural timescales and additional measures of temporal correlationsTCsDefinitionIllustration \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${INT}$$\end{document} Area under the ACF curve up to 0^14,15^. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.1}$$\end{document} Area under the ACF curve up to the time lag where ACF first decays to 0.1. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{INT}}_{0.5}$$\end{document} Area under the ACF curve up to the time lag where ACF first decays to 0.5^129^. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACW}}_{0}$$\end{document} Time lag at which the ACF first crosses zero (after linear interpolation)^44,129^ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACW}}_{0.1}$$\end{document} Time lag at which the ACF first decays to the threshold value of 0.1. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{ACW}}_{0.5}$$\end{document} Time lag at which the ACF first decays to the threshold value of 0.5^43^.HurstSlope of the linear fit between log of fluctuation magnitude and log of time window size, indicating the strength of long-range temporal correlations^124^.

Computational model of stroke-induced network excitability

To investigate stroke-related alterations in INT, we adopted a computational modeling approach grounded in the neurophysiological principles of synaptic plasticity. As reviewed by Stampanoni Bassi et al.^57^, LTP can promote the emergence of highly connected hub regions that support efficient integration, while homeostatic plasticity helps prevent instability by limiting excessive connectivity in peripheral nodes. Stroke disrupts this balance, often inducing local hyperexcitability^103,105^and large-scale network reorganization^19,21^. We hypothesized that such synaptic changes shift the brain’s operating point relative to its critical state, thereby altering its intrinsic timescale structure.

Hence, we implemented a large-scale excitable neuronal network model based on the Kinouchi–Copelli framework^114^ to simulate and explain stroke-related alterations in INT. In particular, a random network of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=\mathrm{100,000}$$\end{document} excitable spiking neurons and a mean average degree \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K=10$$\end{document} was utilized to initialize brain neuronal networks. The neuronal network dynamics were then modeled using the Kinouchi-Copelli model^114^. Each neuron is a cyclic cellular automaton that follows a discrete-time process defined on the state space (quiescent, spike, and refractory). When at a quiescent state, neurons can spike by (i) an external drive modeled as a Poisson process with rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r=\,{10}^{-5}$$\end{document} , which activates neurons with probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=1-{e}^{-r\cdot {{\rm{\delta }}}_{t}}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\delta }_{t}\,=\,1$$\end{document} ms is the time step, or by (ii) receiving input propagation from a spiking connected neuron with probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document} . The spiking neuron will enter a refractory state at the next time step and then return to the quiescent state after a period of 8 ms. We set the time length as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=5050$$\end{document} (unless otherwise stated), including a 50-ms transient time window, and then the neuronal network activity is defined as the instantaneous density of active neurons as a function of time.

The branching ratio \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{\sigma }}(=K\cdot {\rm{\lambda }}$$\end{document} ) of the neuronal networks represents the average number of spikes generated by each excited neuron in the next time step. The brain is in a critical state when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{\sigma }}=1$$\end{document} , where the network activity propagates in a balanced manner, maximizing information transfer and computational efficiency^6,114,125^. In the subcritical state ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rm{\sigma }} < 1$$\end{document} ) activity decays rapidly, leading to diminished neural interactions and reduced information integration. Conversely, in the supercritical state ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma > 1$$\end{document} ) activity grows excessively. For branching processes, the branching ratio can indicate deviations from criticality, which is known as the DTC, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${DTC}=1-{\rm{\sigma }}$$\end{document} in subcritical states and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${DTC}=\sigma -$$\end{document} 1 in supercritical states. The INT is maximized at the critical state, and it decreases with DTC towards both sub- and supercritical states^15^. Due to the nature of the branching ratio, changes in network degree \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K$$\end{document} or the propagation probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document} will affect \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma$$\end{document} and, thus, cause alterations in INT. The effects of varying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${K}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N$$\end{document} on the relationship between INT and DTC were previously explored in the context of aging^15^ and cognitive performance^43^. Given that stroke often results in hyperexcitability^57–60^ and an increased risk of seizures^126,127^, here we vary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda$$\end{document} while fixing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K$$\end{document} for simplicity to examine how stroke-induced changes in neural excitability alter the relationship between DTC and INT. In addition, to relate simulated neuronal activity to BOLD signals, we convolved a long stimulated neuronal network activity (300 s) with the canonical HRF as implemented in SPM^117^. The simulated neuronal activity was generated at a 1 ms time step and then convolved with the HRF, downsampled to a repetition time (TR) of 2 s to yield BOLD-like time series. We computed INT and the TCs on both the raw neuronal activity and the HRF-transformed BOLD signals at each σ to examine whether criticality-related modulations of timescales are preserved in the fMRI domain.

Our decision to vary synaptic propagation probability (λ) while keeping network connectivity (K) constant was specifically intended to capture changes in neural excitability that are characteristic of stroke-induced hyperexcitability. This is consistent with evidence that stroke alters the balance of excitatory and inhibitory transmission and affects both short-term and long-term synaptic plasticity^57^, which in turn modulates the effective propagation of activity across neural networks. These plastic mechanisms influence the emergence of network hubs and the stabilization of network-wide activity, both of which are tightly linked to criticality and timescale dynamics. Hence, such stroke-induced network reorganization can be effectively abstracted by varying λ, which controls the probability of post-synaptic activation in our model. By linking λ to excitability and plasticity, our model bridges microscale synaptic changes with macroscale functional dynamics, reinforcing the relevance of DTC and INT as biomarkers sensitive to stroke-induced network alterations.

Ethics declarations

This research complies with all relevant ethical regulations across all named affiliations. Written informed consent was obtained from all participants in accordance with the Declaration of Helsinki and procedures established by the Medical Research Ethics Committee of the First Affiliated Hospital of Shantou University Medical College (approval number: SUMC-2021-78-K). All aspects of this study were also approved by the Monash University Human Research Ethics Committee.

Statistical analysis

Welch’s t-test was performed to investigate if there were statistically significant differences in the whole brain, global average, and network-level INT between healthy controls and patients. Whenever it involves multiple comparisons, the p-values were corrected with the false discovery rate (FDR) method at a corrected \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha =0.05$$\end{document} . Repeated-measures one-way ANOVA (level of significance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${p} < \,0.05$$\end{document} , FRD corrected^128^) was used to test if INT significantly differed across different time points. Two-way ANOVA was used to test whether the recovery outcome and functional network have an effect on INT. In the case of significant ANOVA results, post hoc t-tests were performed.

Supplementary information

Supplementary Information