Intrinsic timescale is a commonly used measure of spontaneous neural dynamics that quantifies the temporal window of processing of neuronal populations. Intrinsic timescale displays a hierarchical cortical organization across multiple species and imaging modalities, with shorter timescales in sensorimotor cortex compared to association cortex. However, less is known about how intrinsic timescale evolves during human brain development and whether its cortical maturation patterns generalize to independent developmental samples. Here we estimate the intrinsic timescale in two independent cross-sectional datasets of youth (HCPD: n = 565; HBN: n = 729; age range 8–22 years) and investigate its neurodevelopmental patterns. We find that developmental patterns in the intrinsic timescale follow a hierarchical pattern that recapitulates an axis spanning sensorimotor to association cortices (S–A axis). Our analysis of an independent healthy young adult dataset (HCPYA: n = 973, age range 22–37 years) underscores the specificity of these developmental findings, suggesting that the intrinsic timescale develops along the S–A axis in youth and stabilizes in adulthood. Together, these results reveal convergence between major axes of cortical organization and development, highlighting intrinsic timescale as a principled marker of hierarchical brain maturation in youth.

Fig 1. Quantifying intrinsic timescale using autocorrelation function.

Fully processed parcellated fMRI time-series were used to estimate intrinsic timescale for each brain region (i.e., parcel) and participant. We estimated an autocorrelation function (ACF) for the normalized (i.e., z-scored) time series. Intrinsic timescale was then quantified as the sum of positive autocorrelation values (ac  >  0; colored parts of the ACFs above) multiplied by TR for each cortical region and individual. This is approximately equivalent to calculating the area under the curve of the positive part of the ACF curve (i.e., shaded areas). As demonstrated in the example plots, a fast-decaying ACF will have shorter memory, and thus a shorter timescale (orange region and ACF plot), compared to a slow-decaying ACF (purple region and ACF plot). Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.g001

Intrinsic timescale is longer in association cortex compared to sensorimotor cortex in youth

We first sought to characterize intrinsic timescale patterns in the two developmental cohorts. Intrinsic timescale varied across the cortex in both developmental datasets: the Human Connectome Project–Development (HCPD; n = 565; [38] and the Healthy Brain Network (HBN; n = 729; [39] (Fig 2). Visual inspection suggested that the cortical distribution of fMRI intrinsic timescale was broadly segregated into lower-order unimodal regions and higher-order transmodal areas (Fig 2a). The mean timescale maps, averaged across individuals, were positively correlated between the two developmental datasets, demonstrating consistent patterns (Fig 2a; r_s = 0.76, p_spin = 0.0002). To assess whether timescale maps in youth reflect the cortical hierarchy, we directly compared the average timescale to the S–A axis [12]; Fig 2b). The cortical distribution of timescale was positively associated with the S–A axis in both datasets. However, this association was statistically significant only in HCPD (Fig 2b scatter plots: HCPD: r_s = 0.37, p_spin = 0.03; HBN: r_s = 0.27, p_spin = 0.12). Further analysis demonstrated that intrinsic timescale was generally significantly longer in association cortex compared to sensorimotor cortex in both datasets (Fig 2b boxplots: HCPD: t = 7.8, p_spin = 0.01; HBN: t = 6.6, p_spin = 0.04; two-tailed). Collectively, these findings indicate that the intrinsic timescale tends to be longer in the association cortex compared to the sensorimotor cortex in youth.

Neurodevelopmental patterns in intrinsic timescale align with the cortical hierarchy

To characterize the developmental patterns of intrinsic timescale in youth, we examined the relationship between participant age and both whole-brain and regional timescale in HCPD, while controlling for participant sex and in-scanner motion. We found that the whole-brain average timescale increases during development in youth (Fig 3a; partial R² = 0.047, p_anova = 1.25 × 10⁻⁷). Furthermore, we found that regional neurodevelopmental patterns were heterogeneous across the cortex, with smaller age effects in the sensorimotor cortex and larger age effects in association cortex (Fig 3b). Notably, the effect size of the developmental associations aligned with the cortical hierarchy defined by the S–A axis (Fig 3c; r_s = 0.41, p_spin = 0.001). Consistent with this hierarchical pattern, supplementary within-subject analyses showed that the categorical separation between sensorimotor and association timescales increased with age (S1 Fig). Further inspection of regional developmental model fits revealed that timescale increased in association cortex but was stable in sensorimotor regions (Fig 3d). Overall, these findings demonstrate that fMRI intrinsic timescale increases with age in youth, mainly driven by increases in timescale in association cortex. Timescale development reflected a hierarchical pattern, recapitulating the S–A axis of cortical organization.

Findings replicate in an independent developmental cohort

We next examined the extent to which these developmental findings were generalizable by replicating our analyses in the HBN dataset. Consistent with HCPD, we found that the average intrinsic timescale increased during development in HBN (Fig 4a; partial R² = 0.037, p_anova = 1.13 × 10⁻⁷). Regional analysis similarly identified heterogeneous developmental effects across the cortex (Fig 4b). As in HCPD, these developmental effects were aligned with the S–A axis (Fig 4c; r_s = 0.19, p_spin = 0.02). Developmental fits showed that overall, timescale remained stable in sensorimotor regions but increased in association cortex (Fig 4d). To directly evaluate the degree to which results converged across the two datasets, we compared developmental patterns of intrinsic timescale between HBN and HCPD (Fig 5). We found that maps of associations with age (i.e., partial R² maps) were similar for the two datasets (Fig 5; r_s = 0.30, p_spin = 0.0002). Overall, 46.5% of the regions in HCPD and 22.5% of the regions in HBN demonstrated significant age effects after FDR correction. The developmental patterns in these regions had a significant spatial overlap between HCPD and HBN (Fig 5; Dice score = 0.42, p_spin = 0.0008). For regions with at least one significant derivative, we summarized the timing of maximal development as the age at which the regional model derivative reached its maximum absolute value (S2a Fig). These estimates should be interpreted cautiously because age of maximal development is only defined for regions with significant derivative-based change; regions with weak, stable, or uncertain age trajectories do not have a well-defined peak-change age. To complement this analysis, we also examined the age at which whole-cortex rates of timescale differences were maximally aligned with the S–A axis, which peaked in early-to-mid adolescence in both datasets (S2b Fig; 14.6 years in HCPD and 12.8 years in HBN).

Timescale remains relatively stable in adulthood

To assess whether the observed age-related differences in timescale were specific to childhood and adolescence or extended into young adulthood, we repeated all analyses in an independent sample of young adults (i.e., 22–37 years old) from the Human Connectome Project–Young Adults (HCPYA; [40]; n = 973) (Fig 6). As expected, intrinsic timescale in young adults reflected the hierarchical organization captured by the S–A axis (Fig 6a; r_s = 0.4, p_spin = 0.01), such that sensorimotor regions had shorter timescales while association regions had longer timescales (t = 7.5, p_spin = 0.02, two-tailed). However, unlike the developmental findings, age analysis demonstrated that the whole-brain average timescale remains relatively stable in adulthood (Fig 6b; partial R² = −0.004, p_anova = 0.04). Regional age effects estimated using GAMs were not statistically significant for the majority of brain regions (11.7% of regions with significant age effects; Fig 6c) and did not exhibit significant variability along the S–A axis (Fig 6d; r_s = −0.1, p_spin = 0.58). Overall, these results emphasize the specificity of the developmental findings in HCPD and HBN to the developmental period of 8–22 years old, suggesting that fMRI intrinsic timescale develops along the S–A axis in youth and stabilizes in adulthood. Given the differences in acquisition parameters and cohort characteristics, the young adult dataset was used to assess the specificity of the developmental findings to the developmental cohorts and the stability of age-related effects within young adulthood, rather than as a direct continuation of the developmental age trajectories.

Sensitivity analysis

To evaluate whether our findings were influenced by confounding factors and specific analytical choices, we performed a series of sensitivity analyses. We first examined whether the findings were independent from confounding factors such as brain volume and parcel size. Given that brain size changes during development, total brain volume might impact the average whole-brain intrinsic timescale. However, we only found a weak association between estimated Total Intracranial Volume (eTIV) and intrinsic timescale (S3a Fig). To ensure that the developmental patterns of intrinsic timescale were independent from this weak correlation, we repeated the age analysis while including eTIV as a model covariate. The results were consistent with the original findings (S3b Fig). Additionally, given that parcel-wise time series are calculated as the average time series across vertices within a given parcel, parcel size might directly impact smoothness of parcel time series and their ACFs. To test this, we compared parcel size and intrinsic timescale for each parcel and each participant in HCPD and found no associations (mean r_s = 0.01). Moreover, given that HBN is a more heterogeneous sample compared to HCPD and includes help-seeking individuals with higher levels of psychopathology, we repeated all analyses using two subsets of HBN individuals (each with N = 600; this approximately matches the HCPD sample size) with low motion or low psychopathology (i.e., p-factor). The findings were consistent with the full HBN sample (S4 Fig). Furthermore, to assess the extent to which the results were influenced by the choice of atlas and number of parcels, we repeated the analysis with a lower resolution Schaefer atlas with 200 parcels in HCPD. The results were consistent with the original findings with the Schaefer 400 atlas (S5 Fig). Moreover, to ensure that findings were independent from signal-to-noise ratio (SNR) of the fMRI signal, we estimated temporal SNR (tSNR) as the ratio of the time-series mean to standard deviation for each region and participant in HCPD. Neither the intrinsic timescale nor the age effects (partial R²) were significantly associated with tSNR (S6 Fig). Additionally, to ensure that findings were independent from data acquisition sites for HCPD and HBN, we repeated the analyses after harmonizing intrinsic timescale across sites. Findings were almost perfectly correlated before and after harmonization for both HCPD and HBN datasets (S7 Fig).

To test whether intrinsic timescale estimates were sensitive to methodological choices and to better contextualize our findings with prior work, we compared our ACF-sum measure with several alternative approaches for quantifying intrinsic timescale in HCPD. First, we quantified intrinsic timescale using two other approaches based on full ACF: (1) fitting an exponential function to ACF and using the exponential decay constant to estimate intrinsic timescale [4]; (2) directly using the first zero-crossing point of ACF as intrinsic timescale. Both approaches generated intrinsic timescale maps and developmental patterns consistent with the original method (S8 Fig). We further repeated the HCPD analyses using lag-1 temporal autocorrelation (lag-1 TA), following prior work [14,31]. We found that lag-1 TA and ACF-sum timescale differed in their mean spatial distribution but showed partially overlapping developmental effects (S9 and S10 Figs). We also decomposed the ACF-sum timescale into short- and long-lag components and found that lag-1 TA was most closely related to short-lag autocorrelation, whereas the full ACF-sum measure was more strongly influenced by longer-lag positive autocorrelation (S11 Fig). Together, these analyses suggest that our core findings are robust across several full-ACF-based timescale definitions, while also showing that different metrics emphasize different features of temporal persistence: lag-1 TA primarily reflects short-lag autocorrelation, whereas the ACF-sum measure is more sensitive to longer-lag positive autocorrelation across the full ACF. This suggests that the developmental effects observed with our ACF-sum measure were driven primarily by long-lag ACF structure.

Participants

We used resting-state fMRI data from 3 independent cohorts: (1) a cross-sectional developmental dataset used as a discovery sample; (2) an independent cross-sectional developmental dataset used to replicate the findings; (3) a young adult dataset to evaluate the specificity of our results.

Neuroimaging data processing

All three datasets were preprocessed and post-processed with similar image processing pipelines. The detailed descriptions are provided based on reports automatically generated by the implemented processing tools, and are released under the CC0 1.0 Universal license.

Intrinsic timescale estimated using the sum of positive autocorrelation

Fully processed parcellated fMRI time series were used to estimate intrinsic timescale for each brain region (i.e., parcel) and participant. Specifically, we first normalized the fMRI time series at each parcel using a z-score normalization; note that this standardization does not affect autocorrelation-based estimates. We then estimated an ACF for the normalized time series using the Time Series Analysis (tsa) tool from statsmodels [65] in Python. Next, we identified the time point where ACF crosses zero autocorrelation value for the first time (i.e., first zero-crossing time point) and calculated the sum of all positive ACF values up to the first zero-crossing point. To obtain an estimate of the intrinsic timescale in seconds, we multiplied the sum of positive ACF values by the corresponding TR for each dataset (HCPD TR = 800 ms; HBN TR = 800 ms; HCPYA TR = 720 ms; [25,66,67]. Note that the fMRI acquisition TR corresponds to the time step between consecutive time points in ACF. This procedure is highly similar to calculating the area under the curve of the positive part of the ACF curve (Fig 1). We repeated this procedure for each parcel and individual, generating a parcellated map of intrinsic timescale for each participant.

We restricted our analyses to low-motion individuals (mean FD < 0.2 mm; see “Participants” for details), given that measures of intrinsic timescale are sensitive to in-scanner motion and physiological noise and are more reliable in higher-quality data with low motion [68]. In addition, we included data quality (i.e., mean FD) as a covariate in subsequent models to account for any residual confounding effects of motion (see “Generalized additive model (GAM)” for details).

Intrinsic timescale estimated using lag-1 temporal autocorrelation

As a complementary estimate of intrinsic timescale, we calculated lag-1 temporal autocorrelation (lag-1 TA) for each brain region and participant, following the approach used by Shinn and colleagues [31] and Bero and colleagues [14]. Lag-1 TA quantifies the similarity between a time series and a copy of itself shifted by one time point. In the context of fMRI, this corresponds to the correlation between the BOLD signal at time point t and the BOLD signal at the immediately following time point, t + 1, within the same brain region. Higher lag-1 TA values therefore indicate slower decorrelation, or greater temporal persistence, across consecutive time points, whereas lower values indicate faster decorrelation from one time point to the next.

Fully processed parcellated fMRI time series were first z-score normalized within each parcel. We then estimated lag-1 TA using the temporal_autocorrelation(x) function from the “spatiotemporal” Python package [31]. This function computes TA-Δ₁, defined as the first-order temporal autocorrelation between consecutive time points. For each participant, the input was a parcel-by-time matrix, and the output was a vector of lag-1 TA values across parcels. Unlike the ACF-sum timescale measure, lag-1 TA captures only the first nonzero lag of the ACF rather than the area under the positive portion of the ACF. Accordingly, lag-1 TA was not integrated across multiple lags and was not multiplied by TR; it is a dimensionless autocorrelation coefficient. We repeated subsequent analyses using these parcellated lag-1 TA maps to compare regional developmental patterns with those estimated using the original ACF-sum timescale metric.

Short- and long-lag contributions to intrinsic timescale

To examine whether the ACF-sum measure of intrinsic timescale was differentially driven by short- versus longer-lag autocorrelation structure, we decomposed the positive ACF area into short- and long-lag components. First, using the normalized parcellated fMRI time series, we estimated the ACF for each parcel and participant as described above. We then computed a group-average ACF by averaging ACF curves across all parcels and participants. To identify a data-driven lag separating the steep early decay of the ACF from the slower positive tail, we fit a two-piece linear model to the positive portion of the group-average ACF. Specifically, for each candidate lag k, we fit one linear function to the early-lag portion of the group-average ACF and a second linear function to the later-lag portion. We selected the lag k that minimized the total sum of squared errors (SSE) across the two fitted line segments. Lag 0 was excluded from this elbow estimation because the ACF at lag 0 is fixed at 1. We searched candidate elbow lags from k = 3–47 within the positive group-average ACF, corresponding to 2.4–37.6 s in HCPD (TR = 0.8 s). The two-piece linear fit identified an elbow at lag k = 5, corresponding to 4.0 s in HCPD.

Using this single group-level elbow lag k, we then returned to the participant- and parcel-level ACFs and calculated two complementary timescale estimates. The short-lag component, Int_short, was calculated as the sum of ACF values from lag 0 through lag k, bounded by the first zero-crossing point when the ACF crossed zero before lag k. The long-lag component, Int_long, was calculated as the sum of positive ACF values after lag k and up to the first zero-crossing point. Both components were multiplied by the dataset-specific TR to express them in seconds. This decomposition ensured that Int_short and Int_long summed to yield the original positive ACF-sum intrinsic timescale estimate for each parcel and participant. We used this approach to test whether short- and long-lag ACF components captured distinct aspects of development, and specifically whether lag-1 autocorrelation preferentially reflected the short-lag component of the ACF.

Generalized additive model (GAM)

We used GAMs to assess linear and nonlinear associations with age in fMRI intrinsic timescale [13,42,69]. The analysis was performed using the “mgcv” package in R 4.2.2. We included age as a smooth term and sex and data quality (i.e., in-scanner motion quantified by mean FD) as linear covariates in each GAM. We used Restricted Maximum Likelihood (REML) to estimate the smoothing parameters and set the maximum basis complexity to 3 for the smooth age term to avoid overfitting. Specifically, GAMs were formulated as ‘mgcv::gam(feature ~ s(age, k = 3, fx = F) + factor(sex) + data_quality)’. We performed two sets of analysis: (1) whole-brain GAMs, where we assessed developmental patterns of the whole-brain average intrinsic timescale; (2) regional GAMs, where we examined developmental patterns of regional intrinsic timescale.

The effect size of associations between intrinsic timescale and age were estimated as partial R². Partial R² was quantified as the difference in SSE between the full model that included the smooth age term and a reduced model with no age term (i.e., only including model covariates) normalized by SSE of the reduced model. Normalization by SSE of the reduced model highlights the relative contribution of the predictor (e.g., age) to the reduced model. We used a signed version of partial R², such that the sign reflected the directionality of observed effects (increase or decrease in intrinsic timescale with increasing age). To obtain the directionality of partial R², we calculated the mean derivative of the model fit for the smooth term (i.e., age). A positive sign was assigned to partial R² if the mean derivative was positive, reflecting an overall increasing trend between the intrinsic timescale and age, whereas a negative sign was assigned to partial R² for negative mean derivative. The statistical significance of the associations between intrinsic timescale and age was assessed using analysis of variance (ANOVA) between the full model and reduced model that excluded the smooth age term. Results were corrected for multiple comparisons by controlling for the false discovery rate (FDR correction; Q < 0.05).

Age of maximum alignment between age-related timescale differences and the S–A axis

To assess when age-related differences in intrinsic timescale were most strongly aligned with the cortical hierarchy, we estimated age-specific rates of timescale differences from the first derivative of each region’s GAM smooth function for age [42]. These derivatives represent the estimated rate of timescale differences per one-year age interval. For each age from 8 to 22 years, we quantified the spatial association between regional derivative values and regional S–A axis ranks using Spearman’s rank correlation. To account for uncertainty in derivative estimates, we sampled 10,000 draws from the posterior distribution of the derivative for each regional GAM smooth and computed the age-specific correlation between derivative values and S–A axis ranks for each draw. The median correlation across draws was used to summarize the age-specific alignment between timescale differences and the S–A axis, and the 95% credible interval was calculated from the distribution of sampled correlations. We then identified the age of maximum alignment as the age at which the derivative-to-S–A axis correlation was maximal for each posterior draw. The median and 95% credible interval of these ages were used to summarize the timing of maximum alignment between age-related timescale differences and the S–A axis.

Data acquisition site harmonization

To evaluate whether the findings were influenced by acquisition site, we harmonized intrinsic timescale estimates across data acquisition sites using CovBat-GAM [69,70,71,72,73]. Harmonization was performed in R 4.2.2 using the covfam function from the ComBatFamily package. This approach extends ComBat/CovBat harmonization to account for both mean and covariance differences across batches. In this analysis, the study site was treated as the batch variable. Age was included as a smooth term to preserve linear and nonlinear age-related effects, and sex was included as a linear covariate. To assess whether harmonization influenced our findings, we repeated the main analyses before and after site harmonization and compared the resulting whole-brain age associations, regional timescale maps, and regional age-effect maps.

Sensorimotor-association timescale differentiation

To test whether timescale differentiation between sensorimotor and association cortex increased with age, we performed two supplementary analyses in HCPD. First, we conducted a within-subject distance-based analysis. Cortical regions were divided into lower-ranking sensorimotor and higher-ranking association regions based on their position along the S–A axis. For each participant, we computed pairwise absolute timescale differences among regions within the same end of the hierarchy and across sensorimotor-to-association regions. We then compared across-hierarchy differences to within-hierarchy differences using a two-sample t test, yielding a subject-level t-value indexing sensorimotor-association timescale differentiation. We then tested whether this differentiation measure varied with age using a GAM with the same model structure as the main analysis.

Second, we performed a complementary analysis in which regional timescale values were directly compared between sensorimotor and association cortex for each participant using a two-sample t test. The resulting subject-level t-values were then modeled as a function of age using the same GAM framework. For both analyses, effect size was quantified using partial R² and significance was assessed by comparing the full model including age to a reduced model containing only covariates.

Spatial null model

We used spatial autocorrelation-preserving null models (i.e., “spin” tests) to compare cortical patterns between different features or datasets [74,75]. Specifically, we generated null brain maps with a randomized spatial distribution while preserving the spatial autocorrelation inherent to the data using the neuromaps toolbox [76]. The randomization was implemented by applying 10,000 randomly sampled rotations to the spherical projections of the data for the Schaefer atlas [74,76]. The rotations were applied to one hemisphere and then mirrored to the other hemisphere.

The resulting randomized (“spun”) brain maps were used for two types of significance testing. First, to evaluate the significance of associations between two cortical maps (e.g., the correlation between age effects and the S–A axis), we generated a null distribution of association statistics (e.g., Spearman’s rank correlation coefficient) by correlating one empirical map with 10,000 spun versions of the other map. The statistical significance of the observed association was determined by comparing the empirical statistic to this null distribution using a two-tailed test. Second, for analyses testing differences in mean values between two distributions (e.g., t-tests), spun brain maps were used to generate a null distribution of mean differences. The empirical mean difference was then compared against this null distribution to assess statistical significance (two-tailed test).

S1 Fig. Sensorimotor-association timescale differentiation in HCPD.

We tested whether timescale differentiation between sensorimotor and association regions increased with age. (a) In a within-subject distance-based analysis, cortical regions were divided into lower-ranking sensorimotor and higher-ranking association regions based on their position along the S–A axis. For each participant, pairwise timescale differences were computed within sensorimotor/association regions and across sensorimotor-to-association regions. Across-hierarchy differences were compared to within-hierarchy differences using a two-sample t test, generating a subject-level estimate of timescale differentiation. This measure showed a significant positive association with age using a GAM (partial R² = 0.025, p_anova = 0.0007). (b) In a complementary analysis, sensorimotor and association timescale values were directly compared within each participant using a two-sample t test. The resulting subject-level t-values also showed a significant positive association with age using a GAM (partial R² = 0.026, p_anova = 9.36 × 10⁻⁵). Each data point represents an individual participant. GAM fits are shown with black lines. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s001

(TIFF)

S2 Fig. Timing of age-related differences in intrinsic timescale across two developmental cohorts.

(a) Derivatives of the intrinsic timescale regional developmental patterns (Figs 3d and 4d) were calculated as finite differences between consecutive time points to estimate the age at which timescale development demonstrated maximal change. Specifically, the age of peak change was estimated as the age where the derivative was greatest in absolute magnitude. (b) To estimate the age at which age-related timescale differences were maximally aligned with the S–A axis, we calculated the first derivative of each region’s GAM smooth function for age. These derivatives represent the estimated rate of timescale differences per one-year age interval. At each age, we quantified the spatial correlation between regional derivative values and regional S–A axis ranks. To account for uncertainty in derivative estimates, we sampled 10,000 draws from the posterior derivative distribution and calculated the derivative-to-SA axis correlation for each draw. The black line shows the median correlation across draws, and the gray band shows the 95% credible interval. The pink band shows the 95% credible interval for the age of maximum alignment. Maximal alignment occurred at 14.6 years in HCPD and 12.8 years in HBN. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s002

(TIFF)

S3 Fig. Effects of total brain volume on intrinsic timescale and its developmental patterns.

(a) To ensure that findings were independent from changes in brain size during development, we directly examined the relationship between estimated Total Intracranial Volume (eTIV) and intrinsic timescale in the HCPD dataset and found a weak association between the two. Each data point in the scatter plot represents an individual participant. (b) To assess whether developmental patterns of intrinsic timescale were independent from changes in brain size, we repeated the age analysis (i.e., GAMs) while controlling for eTIV as a model covariate. The age effects (i.e., partial R²) were consistent with the original analysis. Each data point in the scatter plot represents a brain region from the Schaefer 400 atlas. r_s denotes Spearman’s rank correlation coefficient. Linear regression lines are added for visualization purposes only. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s003

(TIFF)

S4 Fig. Consistent findings in two subsets of the HBN dataset.

Given that HBN is a more heterogeneous sample compared to HCPD and includes help-seeking individuals with higher levels of psychopathology, we repeated all the analyses using two subsets of HBN individuals (each with N = 600; this approximately matches the HCPD sample size) with (a) low motion and (b) low psychopathology (i.e., p-factor). In both cases, the developmental patterns in intrinsic timescale (i.e., partial R²) were consistent with the full HBN sample. Each data point in the scatter plot represents a brain region from the Schaefer 400 atlas. r_s denotes Spearman’s rank correlation coefficient. Linear regression lines are added for visualization purposes only. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s004

(TIFF)

S5 Fig. Consistent findings using the Schaefer 200 atlas.

Developmental patterns of intrinsic timescale were replicated in the Schaefer 200 atlas in HCPD. (a) Consistent with the finding with the Schaefer 400 atlas (Fig 3), GAM results demonstrated that average timescale increases during development in youth (partial R² = 0.041, p_anova = 9.86 × 10⁻⁷). (b) Region-wise GAMs identified heterogeneous age effects (i.e., partial R²) across the cortex. (c) Age effects on intrinsic timescale were hierarchically organized along the S–A axis. Significance of the association between age effects and S–A axis rank was assessed using 10,000 spin tests. r_s denotes Spearman’s rank correlation coefficient. Linear regression line is added for visualization purposes only. (d) Regional trajectories of developmental patterns in intrinsic timescale were obtained from region-wise GAM results (as shown in panel b). Each line corresponds to the model fit of each cortical region and is colored based on the region’s rank along the S–A axis. Similar to the findings with the Schaefer 400 atlas, the results demonstrate that intrinsic timescale in association regions increases during development in youth while it remains relatively stable in sensorimotor regions. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s005

(TIFF)

S6 Fig. Findings were independent from temporal signal-to-noise ratio (tSNR).

To test whether the findings were sensitive to signal-to-noise ratio (SNR) of the fMRI signal, we estimated temporal SNR (tSNR) as the ratio of the time-series mean to standard deviation for each region and participant in HCPD. Neither (a) the intrinsic timescale nor (b) the age effects (partial R²) were significantly associated with tSNR. Each data point in the scatter plot represents a brain region from the Schaefer 400 atlas. The data points are colored based on their rankings on the S–A axis. r_s denotes Spearman’s rank correlation coefficient. Linear regression lines are added for visualization purposes only. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s006

(TIFF)

S7 Fig. Findings were robust to site harmonization.

To evaluate whether the findings were influenced by acquisition sites, we harmonized intrinsic timescale estimates across study sites using CovBat-GAM. Age effects (i.e., partial R²) were highly consistent before and after site harmonization. Each data point in the scatter plots represents a brain region from the Schaefer 400 atlas. The data points are colored based on their rankings on the S–A axis. r_s denotes Spearman’s rank correlation coefficient. Significance was assessed using 10,000 spin tests. Linear regression lines are added for visualization purposes only. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s007

(TIFF)

S8 Fig. Convergent results with methodological sensitivity analysis in quantifying intrinsic timescale.

To test whether intrinsic timescale estimates were sensitive to methodological choices, we quantified the intrinsic timescale in HCPD using two other approaches. (a) The intrinsic timescale was quantified as the first zero-crossing time point of the fMRI time series autocorrelation functions (ACF). (b) The intrinsic timescale was quantified as the exponential decay constant from an exponential fit to ACF. Both approaches generated intrinsic timescale maps consistent with the original method. Additionally, the developmental patterns were consistent with the original analysis. Each data point in the scatter plot represents a brain region from the Schaefer 400 atlas. The data points are colored based on their rankings on the S–A axis. r_s denotes Spearman’s rank correlation coefficient. Linear regression lines are added for visualization purposes only. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s008

(TIFF)

S9 Fig. Developmental patterns of lag-1 temporal autocorrelation in HCPD.

Lag-1 temporal autocorrelation (lag-1 TA) was estimated in HCPD as a complementary measure of intrinsic timescale. (a) Cortical map of lag-1 TA and relationship between whole-brain average lag-1 TA and age. Each data point in the scatter plot represents an individual participant (partial R² = 0.008, p_anova = 0.03). (b) Region-wise GAMs identified heterogeneous age effects (i.e., partial R²) across cortex. FDR-corrected maps show regions with significant age effects (p_FDR < 0.05; 28.2% of regions with significant age effects). (c) Age effects on lag-1 TA were not significantly associated with the S–A axis. Significance of the association between age effects and S–A axis rank was assessed using 10,000 spin tests. r_s denotes Spearman’s rank correlation coefficient. Linear regression line is added for visualization purposes only. (d) Regional trajectories of developmental patterns of lag-1 TA were obtained from region-wise GAM results. Each line corresponds to the model fit of each cortical region and is colored based on the region’s rank along the S–A axis. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s009

(TIFF)

S10 Fig. Relationship between ACF-sum intrinsic timescale and lag-1 temporal autocorrelation.

We compared regional maps of intrinsic timescale estimated using the original ACF-sum method and lag-1 temporal autocorrelation (lag-1 TA) in HCPD. (a) Regional intrinsic timescale maps estimated using the two approaches were not significantly correlated. (b) In contrast, regional age effects (i.e., partial R²) estimated using the two approaches were significantly correlated. Each data point in the scatter plot represents a brain region from the Schaefer 400 atlas. The data points are colored based on their rankings on the S–A axis. r_s denotes Spearman’s rank correlation coefficient. Significance was assessed using 10,000 spin tests. Linear regression lines are added for visualization purposes only. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s010

(TIFF)

S11 Fig. Short- and long-lag contributions to ACF-sum intrinsic timescale.

To examine whether ACF-sum intrinsic timescale was differentially driven by short- versus longer-lag autocorrelation structure, we decomposed the positive ACF area into short-lag and long-lag components. (a) The short-lag component was modestly associated with the intrinsic timescale estimated using the original ACF-sum method. The regional age effects from the short-lag component were moderately associated with age effects from the original method. (b) In contrast, the short-lag component was highly associated with lag-1 temporal autocorrelation (lag-1 TA), both for regional timescale maps and regional age effects. (c) The long-lag component and its age effects were strongly associated with the original findings estimated using the ACF-sum method. (d) In contrast, the long-lag component was not significantly associated with lag-1 TA in regional timescale maps and was only weakly associated with lag-1 TA regional age effects. Each data point in the scatter plot represents a brain region from the Schaefer 400 atlas. The data points are colored based on their rankings on the S–A axis. r_s denotes Spearman’s rank correlation coefficient. Significance was assessed using 10,000 spin tests. Linear regression lines are added for visualization purposes only. Data and code needed to regenerate this figure can be found at https://doi.org/10.5281/zenodo.22799416.

https://doi.org/10.1371/journal.pbio.3004037.s011

(TIFF)