Domain: neuroscience
Editorial outcome: major_revision

Abstract
The human cortical functional hierarchy, spanning from primary sensorimotor to transmodal association regions, represents a fundamental principle of brain organization. Here, we show lifespan changes in the sensorimotor-association (S-A) gradient in the cortical functional hierarchy using multimodal neuroimaging data from 33,247 participants aged 32 postmenstrual weeks to 80 years. We identify three critical neurodevelopmental milestones: initiation (third trimester to perinatal period), establishment (infancy to early childhood), and expansion-stabilization (late childhood to adulthood). Pronounced gradient changes are predominantly observed during the first decade, with continued refinement extending into mid-adulthood. Spatiotemporally heterogeneous growth patterns in functional gradients align with evolutionary hierarchies, segregation–integration dynamics, structural maturation, and cognitive spectrum development, proceeding along a dominant S-A growth axis. These findings establish a unified neurodevelopmental framework that links connectome gradient dynamics to multifaceted functional and structural properties, advancing our understanding of cortical hierarchy maturation across the lifespan.

Introduction
The functional hierarchy in the human cerebral cortex serves as a fundamental principle of brain organization. At the macroscale, functional connectome gradients provide a critical framework for delineating the intrinsic spatial topography of the brain 1,2. Among these gradients, the sensorimotor-association (S-A) gradient represents the principal axis of the cortical functional hierarchy, spanning from primary and unimodal regions to transmodal association regions 3. This gradient axis is anchored at its extremes by the most distal cortical regions, with multimodal regions occupying intermediate positions along the continuum 1,3,4. This cortex-wide organization ensures the spatial continuity of information transfer across adjacent areas, facilitating a continuous cognitive spectrum ranging from basic sensory processing, perception, and motor action to abstract cognition along the cortical mantle 3,5. Notably, the S-A gradient axis aligns with spatial variations in key neurobiological properties 2, including cortical thickness 6,7, intracortical myelination 8–10, and gene expression signatures 11,12. Perturbations in this hierarchical organization along the S-A axis have been linked to the pathophysiology of neuropsychiatric disorders 13–17. However, despite its central role in cortical organization, the spatiotemporal dynamics of S-A axis across the lifespan, as well as its coevolution with multifaceted structural properties, remain poorly understood. Addressing these knowledge gaps is essential for advancing our understanding of the principles governing cortical functional organization and their implications for both normative cognition and the pathogenesis of neuropsychiatric disorders.

The development patterns of the S-A axis have garnered considerable interests 16,18–23. Emerging evidence suggests that the anterior–posterior gradient configuration begins to form during prenatal development 18,19, manifesting as an immature S-A hierarchy. Although the S-A axis becomes discernible during preschool years 16, its global and regional topography undergoes gradual refinement throughout childhood and adolescence 20,21, continuing into adulthood 22,23. This protracted process facilitates progressive segregation and integration, supporting the emergence of advanced cognitive abilities 19,20. Despite these advances, critical gaps remain in the lifespan maturation of the functional hierarchy. First, most existing studies of the functional hierarchy have focused on discrete developmental stages (e.g., the third trimester 19, neonatal period 18, adolescence 20,21 and adulthood 3,22), resulting in a fragmented perspective that obscures nonlinear gradient dynamics across the full age continuum. Moreover, prior work has often been limited by small (e.g., N < 100) 18,19 or modest (N = hundreds to thousands) 16,20–22 sample sizes, along with a predominant reliance on linear modelling for gradient development 19,20,22. Although a recent study 23 established a continuous, nonlinear mapping of functional gradients from birth to 100 years of age, the developmental patterns remain constrained by a modest sample size (N = 3,556) and particularly sparse data coverage during childhood. Second, previous studies have failed to delineate how functional hierarchies coevolve with specific structural and geometric features, such as cortical thickness 2,6,7, intracortical myelination 2,8,9, and intrinsic geometry 3,24, across the lifespan. While the S-A axis aligns spatially with these structural attributes in adult brains, the developmental functional-structural synchrony remains unclear. A previous study 23 linked lifespan changes in functional gradients with a composite structural gradient, but its analysis strategy that integrates multiple structural features into a morphometric similarity network, did not directly assess how specific structural properties shape functional gradients. Therefore, addressing these fundamental questions will require large-scale, lifespan-spanning multimodal neuroimaging datasets, harmonized processing pipelines and non-linear analytical frameworks capable of modelling lifespan trajectories of functional gradients, structural properties, and their dynamic interactions.

To address these challenges, we utilized multimodal neuroimaging data from 33,247 participants aged between 32 postmenstrual weeks (PMW) and 80 years, collected at 132 global sites. By employing a diffusion map embedding approach 3,25, we identified the S-A gradient of the functional connectome. Specifically, we sought to i) delineate neurodevelopmental milestones of the spatial organization and extent of the S-A gradient axis, ii) map nonlinear trajectories of the S-A axis at the global, system, and regional levels via generalized additive models for location, scale, and shape (GAMLSS) 26, iii) investigate how lifespan changes in the functional gradient axis correlate with multifaceted structural and geometric hierarchies, and iv) conduct a NeuroSynth-based meta-analysis 27 to examine the alignment of developing S-A gradient with adult cognitive spectrum. These integrative multimodal and multifaceted analyses provide insights into the principles governing the lifespan maturation of the cortical hierarchy.

Results
Following a rigorous quality control protocol 28, we analyzed structural MRI and task-free functional MRI data from 33,247 participants aged 32 PMW and 80 years. Detailed demographic information, acquisition parameters, quality control procedures, and preprocessing protocols are described in the Methods section and in our prior work 28. Preprocessed functional data were mapped onto cortical surfaces, and functional connectomes were constructed by computing pairwise correlation coefficients between the time series of all cortical surface vertices (4,609 × 4,609) 28. Using well-established diffusion map embedding techniques 3,25, we decomposed both age-specific, group-level connectomes and person-specific connectomes into multiple gradients representing functional connectivity variations across the brain. Our primary focus was on the lifespan growth of the principal S-A gradient, given its pivotal role in cognition 3, development 19–21,29, disease 13–17, and evolution 30,31.

Topographically distinct growth stages of the S-A gradient axis
We investigated whether the S-A gradient undergoes distinct stages across the lifespan using both group- and individual-level analyses. We first examined how the hierarchical S-A axis emerges during early development and matures into a canonical gradient of cortical functional organization. To this end, we divided the lifespan into 26 distinct age groups 28 and generated 26 age-specific, group-based functional connectomes. Age ranges were defined with finer intervals for early developmental stages, reflecting rapid changes, and broader intervals for adulthood, reflecting relative stability. Initially, the connectome gradient manifests as an anterior‒posterior axis, which gradually evolves into the canonical S-A axis that spans from primary and unimodal regions at one end to transmodal frontal, temporal, and parietal association regions at the other end (Fig. 1a). This maturation process highlights substantial reorganization of the functional topography of the S-A axis over the lifespan. To delineate the temporal sequence of S-A axis establishment and maturation, we performed K-means clustering analysis on the 26 group-based gradient maps (Fig. 1b & Supplementary Fig. 1). This analysis revealed that functional gradients across age groups are optimally clustered into three distinct stages (Fig. 1b, c).

i) Phase I: Initiation. This phase covers the third trimester and perinatal period (approximately 32 PMW to 1 month). During this stage, cortical topography begins to differentiate functionally between anterior and posterior regions, accompanied by a gradual broadening of the gradient range (Fig. 1a, d).

ii) Phase II: Establishment. Spanning infancy to early childhood (approximately 3 months to 4 years), this phase is characterized by a remarkable transition in functional topography, as the anterior‒posterior gradient is replaced by the S-A gradient (Fig. 1a, d). Early in infancy (approximately 3‒9 months), the anterior‒posterior gradient range narrows, and multiple peaks emerge along this axis (Fig. 1a, d), reflecting increased functional differentiation in cortical organization. By late infancy and early childhood (approximately 12 months to 4 years), the canonical S-A axis is established, as evidenced by two distinct peaks in the ridge distribution corresponding to the sensorimotor and association cortices, respectively (Fig. 1a, d).

iii) Phase III: Expansion and stability. This phase extends from late childhood through adolescence and adulthood (5 years to 80 years). During late childhood and adolescence, the cortical topography of the S-A gradient becomes clearly discernible and undergoes progressive expansion in the gradient range, followed by a prolonged period of relative stability in cortical topography over several decades (Fig. 1a, c, and d). Notably, it is during adolescence that the S-A gradient solidifies as the dominant principal gradient (Supplementary Fig. 2).

We also mapped the visual-sensorimotor (V-S) and the modulation-representation (M-R) gradients maps across the lifespan. The V-S gradient emerges from an early anterior-posterior pattern before birth and refines progressively from infancy through adulthood. In contrast, the M-R gradient is preliminarily established by early childhood and gradually stabilizes into its mature configuration during adolescence and adulthood (Supplementary Fig. 3-4). We further computed the spatial similarity between these group-specific gradients and the canonical S-A axis2 and observed that it increased steadily with age in Phase I, slightly decreased in Phase II, and eventually stabilized at a high level in Phase III (Supplementary Fig. 5). Notably, the slight decrease observed during Phase II is accompanied by a period of transient instability in cortical gradient organization, during which the dorsal attention (DA) and ventral attention (VA) systems show temporary changes in their alignment with large-scale gradients, consistent with broader shifts in gradient organization (Supplementary Fig. 6). This trend supports the view that the S-A axis is not fully established at birth but rather emerges through progressive reorganization. The three developmental phases were also identified through individual-level gradient clustering without predefined age bins, within-age-group-harmonized analysis, and within-dataset replication analyses (Supplementary Fig. 7-11). Nonetheless, the definitions and age boundaries of developmental phases could be influenced by methodological choices such as site effects and clustering analysis.

Normative growth of the S-A gradient axis across lifespan
We systematically delineated the lifespan growth trajectories of the S-A gradient at the global, system and regional levels. Lifespan growth curves were modelled using GAMLSS 26, with age included as a smoothing term, sex and in-scanner head motion as fixed effects, and site as a random effect. The growth rates were derived as the first derivatives of these growth curves.

i) Global level. Three global gradient measures 14 were evaluated for the S-A axis: the explanation ratio, which measures the percentage of variance in the connectivity profile accounted for by the gradient; the gradient range, which measures the difference between the two extremes of the gradient; and the standard deviation, which measures the heterogeneity of gradient values across the cortex. The explanation ratio of the S-A gradient exhibited a prolonged period of nonlinear increase, peaking in the late fifth decade (46.8 years, 95% bootstrap confidence interval (CI): 43.4-49.9), followed by a gradual decline (Fig. 2a). Similarly, both the gradient range and standard deviation of the S-A axis demonstrated nonlinear increases, peaking in the fourth decade (gradient range: 32.0 years, 95% bootstrap CI: 29.2-36.1; standard deviation: 36.2 years, 95% bootstrap CI: 33.9-38.8; Fig. 2a).

ii) System level. We examined the lifespan growth trajectories of the gradient range, standard deviation, and gradient values for each system within the S-A axis. Here, individual-level functional systems were derived by integrating age-specific Yeo 7-network atlases with a person-specific iterative algorithm 28,32,33. Generally, the growth curves for the gradient range and standard deviation across all brain systems (Fig. 2b) mirrored the patterns observed at the global level (Fig. 2a). Notably, the default-mode network (DMN) and frontoparietal network (FPN) exhibited higher gradient ranges and standard deviations, whereas the sensorimotor systems presented lower values throughout the lifespan, reflecting greater functional differentiation in association cortices than in primary cortices. Gradient value analysis revealed distinct nonlinear growth patterns across systems (Supplementary Fig. 12). Specifically, the visual, ventral attention, and dorsal attention systems reach their extreme values during adolescence, whereas the DMN, FPN, and somatomotor systems exhibited more protracted trajectories, reaching their extreme values in the fourth decade of life (Fig. 2c). To further characterize the distribution of gradient values across systems along the S‑ A axis, we performed a descriptive exclusion analysis by iteratively removing each system's vertices and quantifying the shift in mean gradient value normalized by vertex count. This analysis revealed system-specific differences in gradient distribution across developmental stages (repeated-measures analysis of variance (ANOVA), F (12,199460) = 341.8, p < 0.001; Fig. 2d and Supplementary Tables 1-4). During Phase I and Phase II, excluding somatomotor-associated vertices produced the largest shifts in the mean gradient value (all post-hoc ps < 0.001), whereas during Phase III, the largest shifts resulted from excluding DMN-associated vertices (all post-hoc ps < 0.001).

iii) Regional level. We characterized the regional maturation processes of gradient values along the S-A axis. The most pronounced changes occurring in the first decade of life and continued into at least early adulthood (Fig. 2e). Notably, substantial changes in regional gradient values were observed not only at the two extremes of the S-A axis (e.g., posterior DMN and lateral visual regions), but also at intermediate positions (e.g., ventral attention regions such as the insula and dorsal anterior cingulate, and dorsal attention regions such as the superior frontal and parietal cortices). To elucidate the spatiotemporal growth principles across the cortex, we applied principal component analysis (PCA) to the fitted growth curves of all vertices. The first principal component accounted for 65.4% of the variance in the regional gradient profiles (Fig. 2f), representing the dominant lifespan growth axis. Along this axis, distinct growth patterns of gradient values were evident in the sensorimotor and association cortices (Fig. 2f). This dominant lifespan growth axis closely aligned with the S-A gradient axis (Spearman’s ρ = 0.91, p < 0.001; Supplementary Fig. 13) and the macaque-to-human evolutionary hierarchy (Spearman’s ρ = 0.32, p = 0.02; Fig. 2g).

Lifespan association between the S-A axis and segregation/integration
The hierarchical organization of cortical architecture facilitates efficient information transfer between sensorimotor and association regions 4, which is critically dependent on the dynamic equilibrium between functional segregation and integration 34. Here, we investigated the lifespan coupling between S-A gradient and segregation-integration dynamics via graph theoretical analysis. Our findings revealed that the global clustering coefficient (C , segregation index) and inverse characteristic shortest path length (1/L , integration index) exhibited distinct nonlinear maturation trajectories: the C demonstrated an initial increase, peaking in early adulthood (32.7 years, 95% bootstrap CI: 30.7–35.0), followed by a gradual decline (Fig. 3a), whereas 1/L p showed an initial decease, reaching its minimum in mid-adulthood (41.1 years, 95% bootstrap CI: 39.6–42.7), before subsequently increasing (Fig. 3b).

Regional clustering coefficient (C ) decreased from primary to transmodal association regions, independent of age (Fig. 3c). While sensorimotor regions presented the highest levels of local clustering throughout the cortex, these primary regions still presented a nonlinear increase in functional specialization from infancy through early adulthood, followed by a decline in later aging stages (Fig. 3c). Regional efficiency (1/L ) demonstrated substantial developmental changes primarily during the first decade of life: sensorimotor regions exhibited peak efficiency during infancy, which gradually decreased with development, whereas association regions showed marked increases in efficiency, becoming dominant from childhood through adolescence and adulthood (Fig. 3d). This pattern reflects a developmental trajectory of reduced integration in primary regions alongside enhanced integration in association regions as functional specialization progresses. Spatial correlation analyses revealed significant associations between clustering coefficient maps and functional gradient maps across the lifespan (all p < 0.05, FDR-corrected; Supplementary Fig. 14a), whereas regional efficiency showed no significant correlations (all p > 0.05, FDR-corrected; Supplementary Fig. 14b). This finding suggests that functional gradient maturation is more strongly coupled with functional segregation than with integration.

Finally, we applied PCA to the fitted growth curves of regional clustering (C ) and efficiency (1/L ) values across all vertices. The first PCA components accounted for 74.9% and 57.6% of the variance in the regional clustering and efficiency growth profiles, respectively (Fig. 3e, f). Notably, both regional clustering and efficiency demonstrated spatially heterogeneous maturation patterns across the lifespan. The lifespan S-A growth axes of regional clustering and efficiency were correlated with those of the canonical connectome gradient (regional clustering: Spearman’s ρ= -0.78, p < 0.001; Fig. 3g; regional efficiency: Spearman’s ρ= 0.19, p = 0.02; Fig. 3h).

To complement the graph-theoretical results, we additionally examined averaged within- and between-system FC across vertices 35,36. Both measures showed systematic variation along the S-A axis and were significantly coupled with S-A gradient maturation across the lifespan (Supplementary Fig. 15-16). PCA further revealed that developmental axes of both measures were correlated with the growth axis of the S-A gradient (within-system FC: Spearman’s ρ= -0.51, p < 0.001; between-system FC: Spearman’s ρ= -0.56, p < 0.001; Supplementary Fig. 15). These findings highlight that the lifespan maturation of the S-A gradient and functional segregation/integration proceed in a convergent manner along an overarching, hierarchical growth axis.

Lifespan association between the S-A axis and structural hierarchies
We investigated how the S-A gradient axis coevolves with multifaceted geometric and structural features across the lifespan. Our investigation focused on three features: geometric distance 3,24, macrostructural cortical thickness 7, and intracortical T1w/T2w contrast ratio 9.

For geometric distance analysis, we first identified vertices within the top and bottom 5% of S-A gradient values (Fig. 4a), representing regions of maximal functional differentiation. For each individual, we calculated the mean normalized geometric distances between paired vertices with extreme gradient values. This geometric distance exhibited a nonlinear trajectory across the lifespan, peaking during young adulthood (20.7 years, 95% bootstrap CI: 17.2–23.2), followed by a gradual decline (Fig. 4b). Cosine similarity analysis between geometric distances and gradient value differences revealed robust geometry–hierarchy coupling (r > 0.88 across all ages), with similarity decreasing until late adolescence (peak at 17.0 years, 95% bootstrap CI: 14.3-19.5) and stabilizing thereafter (Fig. 4c). These findings suggest that while geometric constraints strongly influence the functional hierarchy, their influence gradually diminishes during development.

For structural analyses, our investigation revealed the following key findings: First, transmodal association cortices exhibited greater cortical thickness and reduced T1w/T2w contrast ratio than did sensorimotor cortices throughout the lifespan (Fig. 4d, e), which is consistent with established S-A axis characteristics in adults7,10. Structural changes were most pronounced during the first decade of life and continued into the third decade, with distinct regional patterns: cortical thinning predominated in transmodal association cortices, whereas the primary visual cortex showed thickening (Fig. 4d), and T1w/T2w contrast ratio patterns revealed decreased levels in paralimbic association cortices (e.g., insula and dorsal anterior cingulate) alongside increased levels in sensorimotor regions (Fig. 4e). Second, spatial correlation analyses revealed significant associations between functional gradient maps and both cortical thickness and intracortical T1w/T2w contrast ratio maps across the lifespan (all p < 0.05, FDR-corrected; Supplementary Fig. 17), with the strongest correlations occurring during the first two decades. Third, PCA of the structural growth curves revealed that the first principal components explained 97.3% and 89.0% of the variance in the cortical thickness and intracortical T1w/T2w contrast ratio growth profiles, respectively (Fig. 4f, g). The cortical thickness growth axis exhibited a regionally heterogeneous pattern spanning from sensorimotor to association cortices (Fig. 4f), with a positive correlation with the functionally defined S-A growth axis (Spearman’s ρ= 0.36, p < 0.001; Fig. 4h). In contrast, the T1w/T2w contrast ratio-based growth axis extended from the sensorimotor cortices to paralimbic association cortices (e.g., the insula and cingulate cortex) (Fig. 4g) and was not correlated with either the lifespan growth axis of the functional gradient (Spearman’s ρ = 0.04, p = 0.42; Fig. 4i) or the cortical thickness (Spearman’s ρ = 0.01, p = 0.25). These findings collectively demonstrate that lifespan maturation of functional gradients is primarily coupled with cytoarchitectural development rather than myeloarchitectural development.

Lifespan S-A gradient alignment with adult cognitive spectrum
In adults, the S-A gradient axis supports a cognitive spectrum that spans from perception and action to abstract cognition 3. However, how lifespan changes in the S–A axis align with this adult cognitive spectrum remains unknown. To address this issue, we conducted a cognitive decoding meta-analysis on both group- and person-specific S-A gradients in the NeuroSynth database 27. The dominant S-A gradient was first identified from young adult Human connectome project data (n=897) 37. We then constructed a reference spectrum map comprising 24 cognitive terms distributed along the S-A axis (for details, see Methods). This reference map showed a clear continuum: perception and action terms anchored the primary/unimodal end, while abstract cognitive terms dominated the transmodal end (Fig. 5a), replicating earlier findings 3.

For group-level analysis, we generated cognitive spectrum maps along the S-A axis for each of the 26 age-specific groups via identical procedures to those used for reference map construction (Fig. 5a). We quantified two key metrics: i) the mean gradient bin width associated with cognitive terms, reflecting the degree of clustering of the corresponding brain areas along the S-A axis, and ii) the correlation between the median gradient bin of each term and its corresponding gradient percentile, reflecting the covariation between the cognitive and S-A axes (Fig. 5b, left panel). Our analyses revealed that compared with those of earlier developmental stages, Phase III (expansion and stability stage) presented narrower mean gradient widths (Kruskal‒Wallis test, η2 = 0.62, p < 0.001; post-hoc tests, both p ≤ 0.01; Fig. 5b, middle panel) and stronger cognitive‒gradient correlations (Kruskal‒Wallis test, η2 = 0.73, p < 0.001; post-hoc both p ≤ 0.002; Fig. 5b, right panel). Individual-level analysis of the cognitive spectrum maps revealed that the mean term width decreased rapidly during early development and stabilized after adolescence (peak at 15.3 years, 95% bootstrap CI: 14.5‒16.6; Fig. 5c), whereas the cognitive term‒gradient correlations increased progressively, reaching maximal values in the third decade of life (29.3 years, 95% bootstrap CI: 26.1‒31.8; Fig. 5c). Finally, quantitative comparison of group-based cognitive spectrum maps with the reference map via Dice coefficients revealed a progressive increase in similarity from the perinatal period through early adulthood, followed by a modest decline (Fig. 5d). Individual-level analysis revealed a similar developmental trajectory, with Dice coefficient peaking in the third decade (26.1 years, 95% bootstrap CI: 25.4‒28.7; Fig. 5e). These results reveal a progressive strengthening of alignment between age-specific functional gradients and the adult cognitive spectrum, reflecting the developmental convergence of cortical organization toward adult cognitive-functional topography.

Sex differences in the lifespan S-A gradient
To investigate sex differences in the lifespan trajectory of the S-A gradient axis, we fitted GAMLSS models with a sex-by-age interaction term at the global, system, and regional levels. Globally, females showed higher explanation ratio and greater gradient range and variation across the lifespan than males (all p < 0.05, FDR-corrected, Supplementary Fig. 18a). At the system level, females also showed greater gradient ranges and variations than males (all p < 0.001, FDR-corrected, Supplementary Fig. 18b). Notably, in several systems, including visual, limbic, FPN, and DMN, females exhibited later peak ages than males (all p < 0.001, FDR-corrected, Supplementary Fig. 18c and Supplementary Table 5). Regionally, females showed higher gradient values in medial prefrontal and parietal regions and later temporal regions, whereas males exhibited higher values in the occipital regions (p < 0.001, FDR-corrected, Supplementary Fig. 18d). These findings indicate that sex differences in cortical hierarchy manifest not only in magnitude but also in developmental timing.

Sensitivity analyses
To validate the lifespan growth patterns of the S-A gradient, we conducted a comprehensive sensitivity analysis via the following five strategies 28 (for details, see Methods and Supplementary Fig 19-24): i) To minimize the influence of head motion, we replicated the main analysis via data from 24,494 participants with a stricter motion threshold (mFD < 0.2 mm). ii) To assess the potential effects of sampling, we performed 1,000 bootstrap resampling iterations. iii) To evaluate the reproducibility of our results, we implemented a split-half replication strategy. iv) To examine the potential influence of specific sites, we performed leave-one-site-out analyses. v) To address biases from unbalanced sample sizes and site distributions across age groups, we employed a balanced resampling approach (1,000 times). The results of these sensitivity analyses were statistically compared with the main results. At the global and system levels, all the growth curves strongly correlated with the main results (global measures: Pearson’s r =0.99±0.01, MSE=9.5×10-5±1.7×10-4; Fig. 6a; system-level measures: Pearson’s r =0.98±0.03, MSE=3.5×10-5±5.5×10-4; Fig. 6b). The peak ages of gradient values for each system identified in the validation results closely aligned with those observed in the main results (mean standard deviation age = 0.91 years; Fig. 6c). All sensitivity analyses replicated system-based exclusion results (Fig. 6d and Supplementary Table 6). The regional lifespan growth rate of the S-A gradient was highly consistent across analytical strategies: 95.4% of vertices showed intraclass correlation coefficient values greater than 0.5 and the overall minimum mean square error (MMSE) remained low (< 5 × 10-6) (Fig. 6e). Furthermore, the lifespan growth axis showed strong consistency across strategies (r = 0.97 ± 0.01), with each strategy clearly differentiating the S-A axis (Fig. 6f).

To evaluate the influence of sample size, we conducted a resampling-based sensitivity analysis across sample sizes ranging from 200 to 33,000 participants. Both mean absolute error (MAE) and root mean squared error (RMSE) decreased as sample size increased (Supplementary Fig. 25). Notably, errors at the global and system levels declined sharply until approximately 5,000 participants, after which further reductions plateaued. Regionally, a sample size of ~6,000 participants was required to achieve stable error estimates for more than half of cortical vertices, while a quarter of vertices demanded even larger samples. This analysis indicates that larger samples are essential for reliable gradient growth curves and offers practical guidance for future lifespan neuroimaging studies.

Methods
Data acquisition, quality control, and preprocessing
In this study, we collected 3T structural and task-free functional MRI data, comprising 44,576 imaging scans from 42,428 participants ranging in age from 32 PMW to 80 years across 172 distinct sites. For participants with multiple imaging sessions, only the first scan was selected for inclusion, ensuring that each individual contributed a single observation to the analyses. Each dataset received approval from its respective local ethics committee, and written informed consent was obtained from all participants or their legal guardians. The detailed participant demographics and imaging scan parameters can be found in our prior work 28. To ensure the accuracy and reliability of the neuroimaging data, we implemented a rigorous multistep quality control (QC) framework that combines automated assessments and expert manual inspections (https://github.com/sunlianglong/BrainChart-FC-Lifespan/blob/main/QC/README.md). The final dataset comprised 33,247 participants aged between 32 PMW and 80 years with cross-sectional high-quality structural and functional MRI from 132 sites. The structural MRIs were preprocessed via cortical surface-based processing pipelines 58–60, including skull stripping, tissue segmentation, cortical reconstruction, and cortical surface measurement calculation. The functional MRIs were initially preprocessed in volumetric space with slice timing, motion correction, echo planar imaging (EPI) distortion correction, registration to anatomical images, intensity normalization, and subsequently in surface space with linear trend removal, nuisance regression, temporal bandpass filtering, scrubbing, surface-based smoothing, and resampling to fsaverage4 space (2562 vertices for each hemisphere). More detailed QC and preprocessing procedures can be found in our prior study 28.

Connectome gradient analysis
Following functional MRI preprocessing, we computed pairwise Pearson’s correlations between the timeseries of each vertex pair in fsaverage4 (excluding the medial wall) to generate a 4609 × 4609 functional connectome matrix for each participant 28. The functional connectome matrices underwent Fisher r-to-z transformation to enhance the normal distribution of correlation coefficients.

We first investigated how the S-A axis is established during early development and when it becomes a canonical gradient of cortical functional organization at the group level. Given that functional gradients can be influenced by scanning parameters 61, we selected several major datasets, including the Developing Human Connectome Project (dHCP), Baby Connectome Project (BCP), Children Brain Development Project (CBDP), Human Connectome Project-Development (HCP-D), Human Connectome Project-Young (HCP-Y), and Human Connectome Project-aging (HCP-A), that cover the entire age range and share similar scanning parameters, such as similar spatial and temporal resolutions, and longer scan durations, which are preferable for obtaining reliable functional gradients (the age range covered by each dataset and scanning parameters are shown in Supplementary Table 7). Among the datasets meeting these criteria, the selected datasets also have the largest sample sizes within their corresponding age groups. The participants from these datasets were divided into 26 distinct age groups with finer intervals during the early years and broader intervals during adulthood. To calculate the group-level mean functional connectome, individual connectivity matrices were weighted using Gaussian distribution probability values corresponding to the ages within each group. The group mean connectome for each age group was then computed by averaging the weighted connectivity matrices, with normalization by the sum of the Gaussian weights. The age-specific group mean connectomes were then converted back to r-values via a hyperbolic tangent function, ensuring scaling between -1 and 1.

We identified cortical functional connectome gradients via the BrainSpace toolbox (https://github.com/MICA-MNI/BrainSpace ) 62. To enhance the reliability of gradient estimation while preserving meaningful individual variability, we applied a row-wise top 10% threshold to each group mean connectome, retaining the strongest connections for subsequent gradient analysis 61. This method was validated by repeating the analyses with a stricter top 5% threshold, which yielded highly consistent results (Supplementary Fig. 30). We subsequently computed the row-wise cosine similarity to generate the affinity matrix, which represents the similarity of the vertex-wise functional connectivity profiles (Supplementary Fig. 31). The well-validated diffusion map embedding 3,25 was applied to identify principal eigenvectors, known as functional gradients, which explain spatial variations in functional connectivity. Applying this algorithm to the affinity matrix allows us to identify multiple low-dimensional gradients that explain the variance in functional connectivity in descending order. For each eigenvector, a value was assigned to each cortical vertex, yielding a cortical map reflecting the gradient topography and visualizing macroscale continuous transitions in overall connectivity patterns. The diffusion embedding is controlled by a single parameter α, which influences the density of sampling points on the underlying manifold. On the basis of prior work, we set α = 0.5, which is well suited for brain connectivity data analysis 3.

To ensure accurate comparisons between groups, we applied Procrustes rotation 62 to align age-specific group-level functional gradients to iterative group-level reference gradients. Initially, all group-level functional gradients were aligned to reference gradients obtained from young adults (aged 23-35 years, n = 897) in the Human Connectome Project dataset 37. We then averaged the aligned gradients of all the groups to obtain new reference gradients. Subsequent alignment iterations refined these reference gradients by aligning all group-level gradients to the new reference gradients, repeating this process 10 times. The final aligned gradients resulted from a combination of the original gradients. We identified the correspondence between the original gradients and the aligned ones on the basis of the largest values of the final transformation matrices 14. Compared with constructing the template directly from individual gradients, this multi-step alignment strategy reduces biases from sample size imbalances and development variations in gradient topography.

Individual gradients were computed from each participant’s functional connectome using the same procedure and aligned with their corresponding age-specific group-level gradients, which were iteratively aligned to the reference gradient as described above. Following Procrustes alignment, we observed a high spatial correspondence between group-level gradients aligned using different template construction strategies, indicating that the core large-scale topography of the S-A gradient was preserved across alignment procedures (Supplementary Fig. 32a). At the individual level, cosine similarity analyses revealed highly comparable lifespan trajectories across all template choices (Supplementary Fig. 32b), demonstrating that the broad developmental pattern of increasing alignment toward an adult-like S-A organization is robust to reasonable variations in template selection. Within this interpretive framework, template-based similarity metrics primarily reflect alignment to a reference framework rather than independent evidence of biological invariance. In line with recent developmental and comparative gradient studies, we applied no scaling during alignment to retain biologically meaningful differences in gradient range across individuals and age groups. Notably, applying scaling substantially attenuated the developmental trajectory of gradient range (Supplementary Fig. 33), despite maintaining high spatial correspondence (average Spearman’s ρ = 0.99).

Our subsequent analyses focused on the S-A gradient, given its central role in cortical hierarchical organization encompassing cognition 3, development 19–21,29, disease 13–17, and evolution 30,31.

Identifying distinct growth stages of the S-A gradient axis
To elucidate the temporal sequence of the S-A gradient development, we categorized the age-specific group-level functional gradients into different stages. We first computed the similarity between any two age-specific group-level gradients via cosine correlation coefficients. The similarity matrix was then transformed into a dissimilarity matrix by subtracting the correlation values from 1. Classical multidimensional scaling was applied to reduce the dimensionality of the dissimilarity matrix to two dimensions, explaining 87.9% of the variance in total. This step simplifies the data while preserving its dissimilarity structure, making it more suitable for clustering. We then performed K-means clustering on the two-dimensional dissimilarity matrix to define clusters (phases) of functional gradients from all age groups. The number of clusters was evaluated from 1 to 10, and the optimal number of clusters was determined via a winner-take-all approach on the basis of thirty indices via the NbClust package 63.

Characterizing the growth pattern of the S-A gradient axis
To quantify age-related changes in functional connectome gradients, we comprehensively calculated measures at the global, system and regional levels. At the global level, we computed the explanation ratio, range, and standard deviation for each individual functional gradient. The explanation ratio reflects the percentage of variance in the functional connectivity profile explained by the gradient. The gradient range measures the difference between the two extremes of the gradient with a larger range indicating a greater difference in functional profiles between vertices at the extremes of the gradient. The standard deviation of the functional gradients represents the heterogeneity of functional connectivity profiles across the cortex. The lifespan growth curves of these measures were fitted via the GAMLSS (for a detailed description, see Modelling growth curves across the lifespan section).

At the system level, we used the well-validated age-specific Yeo 7-network atlas to assign cortical vertices to functional systems 28,32,33. We calculated the gradient range, standard deviation, and gradient values for each functional system and each individual. We then constructed the growth curves of these system-level measures across the lifespan via the GAMLSS. To describe the distribution of gradient values across functional systems along the S-A axis, we performed a descriptive exclusion analysis in which vertices belong to individual system were removed from the precomputed gradients and the resulting changes in mean gradient values were quantified (Δ mean gradient value) for each individual. The Δ mean gradient value for each system was normalized by vertex number of that system. We then performed repeated-measures -ANOVA on the Δ mean gradient value to assess variations in system contribution across different developmental phases. In this analysis, the functional system was regarded as a within-subject factor, and the developmental phase was regarded as a between-subject factor. Interaction effects between the functional system and developmental phase indicate whether the development of the Δ mean gradient value differs among functional systems. The main effects of the functional system revealed variation in the Δ mean gradient value across different systems. The main effects of developmental phase indicate overall changes in the Δ mean gradient value across different phases. Post-hoc pairwise t-test comparisons with Bonferroni correction were performed to examine specific differences between functional systems and phases. All statistical tests were two-sided.

We then investigated the spatiotemporal growth patterns of the functional gradients at the regional level. We first fitted the growth curves of the functional gradient for each vertex via GAMLSS, and computed the first derivative to estimate the growth rate. Next, to identify the spatial axis that explained the largest variance in the age-related changes in gradient values across the cortex, we applied data-driven PCA to the fitted growth curves of functional gradient across all vertices. The first PC is referred to as the principal lifespan growth axis, capturing variability in lifespan growth profiles within a low-dimensional representation. To characterize the variation in growth curves and growth rate curves along this principal lifespan growth axis, we divided the axis into 20 decile bins and calculated the average gradient values and growth rates for all vertices within each bin. To explore the relationships between this lifespan growth axis and other neurobiological properties, we computed Spearman’s correlation coefficients between the principal lifespan growth axis and the S-A axis derived from the functional connectome of HCP young adults. To evaluate potential nonlinearity in spatial correlation analyses, we compared first-order (linear) and second-order (quadratic) models using Akaike Information Criterion (AIC)-based model selection. A quadratic term was retained only when it yielded a ≥5% reduction in AIC relative to the linear model, otherwise linear models were retained. Model order was limited to second order to balance flexibility and interpretability. Additionally, we assessed the correlation between the lifespan growth axis and the cortical evolutionary hierarchy, which was quantified by mapping macaque-to-human cortical expansion 30. This cortical evolutionary hierarchy was obtained from the neuromaps dataset 64 and downsampled to the fsaverage4 space via the neuromaps toolbox (https://github.com/netneurolab/neuromaps). The significance of these correlation analyses were determined using non-parametric 1,000 spin permutation tests that control for spatial autocorrelations 65.

Modelling growth curves of the S-A gradient axis
We used the Generalized Additive Models for Location, Scale, and Shape (GAMLSS) model to estimate age-related changes in functional gradients across the lifespan from 33,247 functional MRI scans. GAMLSS is a flexible statistical framework that models not only the mean of a brain phenotype but also its variability and distribution shape 26. The model was optimized to estimate nonlinear lifespan growth curves, their confidence intervals, and the first derivatives of these curves. The reliability of the models was assessed via several sensitivity analysis procedures, including bootstrap resampling, leave-one-site-out, balanced resampling, and split-half replication procedures (see the Sensitivity analyses section).

i) Distribution selection
We employed the gamlss package (version 5.0-6) in R to establish growth curves for each global, system-level and regional measure. The process begins by evaluating a range of three- or four-parameter GAMLSS distribution families to identify the best-fitting distribution, with GAMLSS fitted to global brain phenotypes, including the explanation ratio, gradient range, and standard deviation for each distribution type. The selection process started with fitting a baseline normal distribution model, which served as a reference. We then systematically evaluated a broad range of potential distributions via the chooseDist function, which ranks distributions on the basis of the Bayesian Information Criterion (BIC). The optimal distribution was selected as the one with the lowest BIC.

We found that while the optimal model varied across different brain phenotypes, the skew exponential power type 4 (SEP4) distribution generally yielded the lowest BIC compared with that of other distribution families (Supplementary Fig. 34). Thus, we used the SEP4 distribution in GAMLSS to estimate the age-related changes in functional gradients, with age, sex, mFD and site as predictors. To align the dHCP dataset with other cohorts, we converted gestational age to chronological age. The conversion was based on the assumption that 40 gestational weeks corresponds to birth (0 years of chronological age). The specific formula applied was: Chronological age (years) = (Gestational age at scan – 40) × 7 / 365. This procedure ensured all ages were referenced to a consistent postnatal timeline.

ii) Model fitting
SEP4 is a four-parameter distribution: median (𝜇), coefficient of variation (𝜎), skewness (𝜈) and kurtosis (𝜏). The parameters were specified as follows: The location parameter (𝜇) included age as a smoothing term, sex and mFD as fixed effects, and site as a random effect. The scale parameter (𝜎) included age as a smoothing term, and sex as a fixed effect. The skewness (𝜐) and kurtosis (𝜏) parameters were fixed at 1. The model equation is as follows:
𝑦 = 𝑆𝐸𝑃4(𝜇,𝜎,𝜈,𝜏) (1)
𝜇 = 𝑓 (𝑎𝑔𝑒)+𝛽1(𝑠𝑒𝑥)+𝛽2(𝑚𝐹𝐷)+𝑧 (𝑠𝑖𝑡𝑒) (2)
𝜇 𝜇 𝜇 𝜇
𝜎 = 𝑓 (𝑎𝑔𝑒)+𝛽 (𝑠𝑒𝑥) (3)
𝜎 𝜎
𝜈 = 𝛽 (4)
𝜈
𝜏 = 𝛽 (5)
𝜏
Considering the complexity of lifespan growth trends, we explored a range of model parameters to capture age-related changes. Specifically, we fitted three models with varying degrees of freedom for the age spline in the location parameter: 3, 4, or 5 degrees of freedom. The degree of freedom for the age spline in the scale parameter was set to 3. The convergence criterion was set to a default log-likelihood value of 0.001 between iterations, with a maximum of 200 iteration cycles. The best-fitting model was selected on the basis of having the lowest BIC value.

Once the optimal model was selected, we used it to generate predictions and compute the median growth curve, along with the model parameters (μ, σ, ν, and τ). These parameters were then used to calculate quantiles. To explore sex-specific growth curves, quantile curves for ages ranging from 32 PMW to 80 years were generated separately for males and females across different sites. For each site and sex combination, we computed quantiles at the 0.05, 0.25, 0.5, 0.75, and 0.95 percentiles, resulting in age-specific growth curves.

iii) Model evaluation
To comprehensively evaluate the fit of our GAMLSS, we employed a combination of residual analysis and predictive performance metrics. We first analyzed normalized quantile residuals, which are recommended for GAMLSS diagnostics because of their suitability for complex distributional models. Specifically, we examined residuals plots against the fitted values of the location parameter (μ), residuals against observation indices, kernel density estimation of residuals and the normal quantile‒quantile (Q‒Q). Our results revealed a well-behaved residual pattern (Supplementary Fig. 35), with residuals distributed uniformly around the horizontal zero line in the residuals versus fitted values or observation indices plots. Kernel density plots suggested indicating an approximately normal distribution, and Q‒Q plot displayed a straight-line relationship with a near-zero intercept and a slope close to one. The detrended transformed Owen’s plots confirmed that residuals remained within the nonparametric confidence intervals (Supplementary Fig. 36). These diagnostic results indicate that the GAMLSS provided an adequate fit, with residual patterns supporting the assumption of normality and no evident systematic deviations.

To complement the residual diagnostics and evaluate the predictive accuracy, we conducted a cross-validation procedure by dividing the dataset into training and testing sets. We split all the participants into two groups, balanced for age and site. The two groups consisted of 16,661 and 16,586 participants respectively. We computed the R-squared (R2) values to quantify the explained variance in the central tendency in the testing set. Furthermore, centile calibration was evaluated via randomized z-scores, with normality assessed via the Shapiro‒Wilk test. A W statistic close to 1 indicated a good fit. We also examined skewness and kurtosis to measure distribution symmetry and tail behavior. A skewness near zero reflected a symmetric distribution, whereas a kurtosis close to zero indicated well-behaved residual distributions without heavy tails (Supplementary Fig. 37 and Supplementary Table 8).

Relating network segregation/integration to the S-A gradient axis
The segregation and integration of the functional connectome reflect local information specialization and global information integrity within the brain network, respectively 34. In our prior work, we demonstrated that functional segregation and integration mediate the development of global S-A gradient metrics in children 20. Here, we aimed to investigate whether the maturation of the S-A gradient is coupled with the development of functional segregation and integration across the lifespan. We first applied a 10% density threshold to each individual’s functional connectivity matrix and binarized the network. We then computed graph theoretical metrics via our Parallel Graph-theoretical Analysis (PAGANI) toolbox 66 (https://www.nitrc.org/projects/pagani_toolkit/) to evaluate functional segregation and integration. For each node, we calculated the clustering coefficient (C ) and path length (1/L ). The clustering coefficient is the fraction of triangles around a node and quantifies the local clustering and closeness of a network. The path length is the average shortest path length from a given node to all other nodes, providing a measure of network integration and global routing efficiency. The global C and 1/L were derived by averaging these metrics across all nodes. To mitigate potential issues caused by extreme path lengths, we used the inverse of the path length to ensure more accurate model estimation by reducing the influence of outliers. To verify the robustness of the findings, we additionally performed weighted network analyses (Supplementary Fig. 38).

The GAMLSS was applied to fit the lifespan growth trajectories of the global and regional measures. After obtaining the growth curves for the C and 1/L at each vertex, we computed Spearman’s correlation coefficients between the fitted cortex-wide S-A gradient maps and the fitted C and 1/L maps at the corresponding ages. This approach allowed us to examine the growth pattern of the spatial correspondence between the S-A gradient and functional segregation/integration across the lifespan. In addition, we performed PCA on the fitted growth profiles of C and 1/L to identify the principal lifespan growth axis for both segregation and integration. We then assessed the correlation between the lifespan growth axis of the S-A gradient and the lifespan growth axes of segregation and integration, respectively. Statistical significance was determined via 1,000 non-parametric spin tests that control for spatial autocorrelations 65. All p-values of the spatial correlation between the S-A gradient and functional segregation/integration were corrected with multiple comparisons controlling for a false discovery rate (FDR) 67 of less than 0.05.

Relating structural and geometric hierarchies to the S-A gradient axis
To explore the structural basis of functional organization across the human lifespan, we assessed the spatial correlation between the S-A gradient and various multifaceted structural attributes, including: i) the geometric distance, defined as the shortest path between two vertices along the cortical surface, reflecting intrinsic cortical geometry 3,24; ii) macrostructural cortical thickness, representing the distance between white matter and pial surfaces, which correlates with cortical cytoarchitecture 7; and iii) intracortical T1w/T2w contrast ratio, derived from the signal intensity ratio of T1-T2-weigthted images, providing an indirect measure of regional myeloarchitecture 9.

i) Geometric distance
Considering the geometric constraints on the functional dynamics 24 and spatial arrangement 3, we explored whether the lifespan growth of the S-A gradient correlates with changes in the macroscale geometry. We analyzed the lifespan changes of the geometric distance and its correlation with the S-A gradient. We considered only vertices with top and bottom 5% S-A gradient values because their differences are indicative of cortical functional differentiation and computational constraints. The geometric distance between any two vertices is computed as the shortest path between them on the triangle surface mesh via the wb_command -surface-geodesic-distance implemented in the workbench command (https://www.humanconnectome.org/software/workbench-command/-surface-geodesic-distance). We computed the geometric distance on the mid-thickness surface in the fsaverage_LR32k native space. To account for inter-individual differences in brain size, geodesic distances were normalized by the square root of the cortical surface area. We then obtained the geometric distance matrix between the paired top and bottom gradient vertices. The growth curve and its growth rate curve for the mean geodesic distance were charted via GAMLSS to evaluate the lifespan changes in the geometry between most functionally differentiated regions. Then, the cosine similarity between paired geometric distances and gradient value differences was computed to measure the correlation between the geometric distance and S-A gradient. The growth curve and its growth rate curve for the cosine similarity were charted via GAMLSS to evaluate the lifespan changes in the geometric constraints on functional differentiation.

ii) Cortical thickness
Cortical thickness is one of the most commonly used cortical morphological measures and is related to cortical cytoarchitecture 7. We aimed to investigate how lifespan changes in cortical thickness are related to the increase in functional differentiation, as indicated by the S-A gradient across the lifespan. The cortical thickness computation was implemented in the HCP structural preprocessing pipeline 60. The cortical thickness was computed between the white (gray‒white interface) and pial surfaces in native space. For each vertex, the shortest distance of that vertex to any other vertex on the other surface was found. The shortest distances from the pial to the white surface and from the white to the pial surface were then averaged to compute the cortical thickness. The cortical thickness for each participant was resampled from native space to fsaverage_LR32k space and then to fsavareage4 space. The final cortical thickness data were derived from 24,992 extensive scans collected from 77 sites. The age and site distributions can be found in Supplementary Fig. 39a.

iii) Intracortical T1w/T2w contrast ratio
Intracortical myelination, a key marker of brain plasticity, develops in stages and different regions across the lifespan with primary motor and sensory regions reaching all milestones earlier than association regions do 49. We explored the impact of intracortical myelination growth on functional S-A gradient growth across the human lifespan. Intracortical myelination was assessed via a previously validated T1-weighted (T1w) /T2-weighted (T2w) contrast ratio approach 9,60 and subsequent improvements 68,69. To mitigate inter-scanner variability and enhance cross-site comparability, all T1w and T2w images underwent bias correction and intensity normalization prior to the computation of the T1w/T2w ratio. We estimated intracortical T1w/T2w contrast ratio within the accurate, high-resolution cortical ribbon volume produced during the preprocessing pipeline. The T1w image was divided by the aligned T2w image within voxels between the white and pial surfaces. The division helps cancel out signal intensity bias related to the sensitivity profile of radio frequency receiver coils, enhancing myelin contrast. This ratio was then mapped onto the mid-thickness surface in fsaverage_LR32k space, which was created by averaging the white and pial surfaces, minimizing partial volume effects. The T1w/T2w ratio maps were resampled to fsaverage4 space. The final intracortical T1w/T2w contrast ratio data were extracted from 18,502 scans from 11 sites. The age and site distribution can be found in Supplementary Fig. 39b.

To quantify the structural substrates of S-A gradient growth, we examined its relationship with morphometric (cortical thickness) and microstructural metrics (T1w/T2w contrast ratio). First, we fitted the growth curves of the structural attributes for each vertex via GAMLSS with age and sex as fixed effects, and site as a random effect. We obtained the fitted structural attribute maps and their corresponding growth rate maps across the lifespan. We computed the Spearman correlation coefficients between the fitted cortex-wide S-A gradient maps and the fitted structural attribute maps at corresponding ages. We obtained the growth curve of the spatial correspondence between the gradient and structural attributes across the lifespan. The same PCA analyses were performed on the fitted growth profiles of cortical thickness and intracortical T1w/T2w contrast ratio. We then assessed the correlation between the lifespan growth axis of the S-A gradient and the lifespan growth axes of cortical thickness and intracortical T1w/T2w contrast ratio, respectively. Significance was determined using non-parametric spin tests that control for spatial autocorrelations 65. All p-values of the spatial correlation between the S-A gradient and structural features were corrected, with multiple comparisons controlling the FDR 67 to be less than 0.05.

Mapping the S-A gradient maturation to the adult cognitive spectrum
To investigate how the organization of the cognitive spectrum along the cortical S-A gradient varies across the lifespan, we conducted a functional decoding meta-analysis using both group-based and person-specific S-A gradients derived from the NeuroSynth database 27. First, we created a reference map of the cognitive spectrum along the S-A axis using data from young adults (n=897) in the HCP database 37. The cortical regions were divided into 20 bins along the S-A axis and projected to the Montreal Neurological Institute (MNI) volume space via the BrainStat toolbox (https://github.com/MICA-MNI/Brainstat). Following the approach outlined by Margulies et al., 2016 3, we obtained 24 cognitive terms describing the cognitive spectrum, spanning from perception and action to abstract cognition. We computed the associations between each of the 20 regions and each of the term maps respectively. These cognitive terms were ranked by their weighted average position along the gradient, as illustrated in Fig. 5a. Notably, the HCP database was excluded from the following cognitive spectrum analyses.

For the group-based S-A gradients, we computed the cognitive spectrum maps for each of the 26 age-specific groups, following the same procedure used to generate the reference map. For each map, we calculated the width of the distribution of each cognitive term within the gradient bins to measure the degree of clustering of the corresponding brain regions along the S-A axis. Additionally, we computed Spearman’s correlation coefficient between the median position of each term in the gradient bin and the gradient percentile, reflecting the covariation between the cognitive and S-A axes. The mean width across terms and Spearman’s correlation differences among stages were tested via the non-parametric Kruskal‒Wallis test followed by post-hoc pairwise t-test comparisons with Bonferroni multiple comparison correction. All statistical tests were two-sided. To investigate the temporal growth of the organization for each cognitive term, we modelled their growth curves via the GAMLSS. We also fitted the growth curves of the mean width across terms and Spearman’s correlation coefficient using the GAMLSS.

To quantify the similarity between each group-based cognitive spectrum map and the reference map, we computed Dice coefficients, with both the S-A and cognitive axes fixed. We also generated person-specific cognitive spectrum maps along the cortical gradients and assessed their similarity to the reference map via the same procedure. Finally, to characterize changes in the organization of the cognitive spectrum, we fitted a growth curve and its growth rate curve of the Dice coefficient via the GAMLSS.

Sensitivity analysis
Although the GAMLSS provides automated parameter optimization to best fit the included data, covariates, and random effects, we conducted several sensitivity analyses 28 to test the robustness and reliability of the optimized models.

i) Stricter head motion threshold analysis
To reduce the potential influence of head motion on functional MRI data, we applied a more stringent head motion threshold by excluding participants with a mean framewise displacement (FD) greater than 0.2 mm. After applying this stricter threshold, we reanalyzed the lifespan growth curves for our main findings. Specifically, we reran the GAMLSS using data from the 24,494 participants who met the stricter motion criteria for S-A gradient measures at the global, system and regional levels.

ii) Bootstrap resampling analysis
To examine the lifespan growth curves for a broader population and compute the confidence intervals (CIs), we performed 1,000 bootstrap iterations via stratified sampling with replacement. In each iteration, the relative proportions of age and sex within each interval were maintained. The lifespan age range (from 32 PMW to 80 years) was divided into 10 equal intervals. We then fitted 1,000 growth curves and calculated 95% CIs for the median (50th centile) curve of S-A measures at the global, system and regional levels. These CIs were derived by computing the mean and 1.96 times the standard deviation of the growth curves and growth rates across all 1,000 bootstrap iterations.

iii) Split-half replication analysis
To evaluate the replicability of the lifespan growth curves in independent datasets, we implemented a split-half replication strategy. All participants were randomly assigned to two equally sized groups, with stratification by site. The two groups consisted of 16,661 and 16,586 participants, respectively. We estimated lifespan growth curves via the GAMLSS for each group independently.

iv) Leave-one-site-out analysis
To evaluate the robustness and consistency of the lifespan growth curves, we performed leave-one-site-out analyses. For each of the 132 sites, we excluded data from that site and refitted the GAMLSS to the remaining data from the 131 other sites. We compared the median centile curves derived from the remaining sites with those from the full dataset. The 95% CIs were again computed from the mean and 1.96 times the standard deviation of the growth curves and growth rates across all leave-one-site-out iterations.

v) Balanced resampling analysis
To address potential biases due to unbalanced sample sizes and sites across different age groups, we employed a balanced resampling strategy. This approach ensured an equal number of participants and sites across each age group. The lifespan age range was divided into 16 age bins, each spanning five years. For each bin, the number of participants and sites was determined on the basis of the bin with the fewest participants and sites, excluding the age groups under 5 years and over 70 years. The smallest participant group had 457 participants (in the age bin 35-40 years), and the smallest site group had 23 sites (in the age bin 40-45 years). For each resampling iteration, 457 participants from 23 sites were randomly selected for all age bins, resulting in 6,770 participants per resample. This process was repeated 1,000 times, and the GAMLSS was refitted for all the global and system-level measures. The consistency of the resampling results was evaluated by computing the correlation between the median growth curves from each resampling and the complete dataset. We repeated the above process once and refitted the GAMLSS for all gradient values of cortical surface vertices.

We compared the results of the sensitivity analyses with the main findings. At the global and system levels, we computed Pearson’s correlation coefficients and the MSE between the growth curves from the main results and those from the sensitivity analyses. To validate the lifespan growth of the gradient values for the functional systems, we calculated the standard deviation of the peak age across the validation strategies within each system. Additionally, we conducted repeated-measures-ANOVA on the Δ mean gradient value for each validation strategy. At the regional level, we assessed the consistency of the growth rate curves across the validation strategies by calculating the ICC and MMSE. The ICC was determined by comparing the variance between different time points and the variance within each time point across validations. The MMSE score was calculated as the minimal value of the average squared differences between the growth rates from the validation and those from the main results. Finally, for the lifespan growth axis of the S-A gradient, we computed Spearman’s correlation coefficient between the main results and sensitivity analyses.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this paper, the author maps the development of the sensory-association (SA) gradients across the human lifespan and explores how these gradients relate to structural brain measures. While the characterization of the SA gradient is a valuable contribution to understanding brain organization and maturation, I have major concerns regarding several aspects of the study. These include questions about the novelty of the work in light of existing literature, the appropriateness and rigor of the methodological approaches, the clarity and robustness of the reported results, and the strength and validity of the conclusions drawn from the findings.\nComments:\n1. Introduction - “First, existing studies focus mainly on discrete developmental stages (e.g., the third trimester 19 , neonatal period, adolescence, and adulthood), without provided a continuous spatiotemporal perspective from prenatal development through late adulthood.” This is not true and misleading in light of Ref 23.\n2. Introduction - “Furthermore, methodological limitations, such as small (less than one hundred) or modest (hundreds to thousands) sample sizes and an overreliance on simplistic linear modelling approaches, fundamentally constrain our ability to capture nonlinear neurodevelopmental dynamics.” This is also not true. The methods used for data processing and analysis (FreeSurfer, GAMLSS, etc.) in the current work are mostly publicly available and are quite standard. There are no methodological hurdles preventing other groups from carrying out similar work.\n3. Introduction - “Second, current conceptual frameworks predominantly emphasize functional gradient remodelling during development while largely overlooking its structural foundations. In mature brains, the S-A axis aligns spatially with structural features indexed by cortical thickness, intracortical myelination, and intrinsic geometry.” This is not true. Association between primary gradients and structural gradients has already been investigated for example in reference 23.\n4. Why 26 age groups? This creates artificial temporal demarcations between different life periods that do not necessarily reflect true biological development. This is particularly problematic when different brain regions have distinct development timelines that do not coincide with these artificial temporal partitions. It is uncertain if age grouping will result in artificial disconunity in the analysis, not just at the global whole-brain level but also local vertex level.\n5. How does the choice of 26 age groups influence the clustering results? The patterns shown in Figure 1b may vary not only with the number of age groups but also with how these groups are distributed over time. Moreover, the age grouping adopted from the authors’ previous publication may not be appropriate here, as this study focuses on gradients rather than simple features of functional connectivity. The authors should demonstrate that this age grouping does not introduce bias by comparing it against clustering results obtained from the entire dataset without imposing predefined age groups. 6. How is the difference in sample size across age groups dealt with? The authors did a “balanced resampling analysis” by equal resampling across age bins. However, it is unclear why the 16 “age bins” are different from the 26 “age groups”. Also, does the fact that resampling uniformly with a smaller number of samples per age bin still replicates the main results suggest that the large sample size used in the study may not be necessary?\n7. How are differences in data variances between datasets and across time dealt with?\nThere can be nonlinear variance changes across age groups and within age groups and also across the different datasets and within the datasets.\n8. A large sample size like the one used in this paper is not always accessible to all researchers. Please show empirical evidence that such a large sample size gives conclusions that cannot be drawn using a smaller sample size.\n9. Anterior-posterior gradient during early infancy – This is contradictory to previous observation that the SA gradient already exists at birth. The SA gradient does not start to occur at 3 months but before that.\n10. Grouping the broad age range from 5 years to 80 years as a single phase (Phase III) is too coarse and does not reflect the fact that growth inflection points happen quite often during adolescence.\n11. Datasets like the developing HCP (dHCP) use gestational age. Postnatal datasets such as the BCP use chronological age. How is the discrepancy between gestational age and chronological age resolved?\n12. “Collectively, these results revealed previously unreported spatiotemporal patterns and critical neurodevelopmental milestones of the functional hierarchy, as characterized by the initiation, establishment, and expansion/stability of the S-A axis of cortical functional organization.” Again, this is untrue in light of ref 23.\n13. Range and variance are measuring almost the same thing, i.e., the spread of the gradient values, and are therefore redundant. Please pick one to report.\n14. In Figure 2c and throughout the paper, it is unclear what “gradient score” is referring to.\n15. The results shown in Figure 2e seem to indicate that the SA gradient changes very little after 20 years. This contradicts Figs 2a and 2b and other results shown in the paper.\n16. T1-weighted and T2-weighted images are not quantitative and can have arbitrary signal gains associated with different scanners, sequences, and choices of acquisition parameters. It is unclear how this is dealt with in this study, which involves datasets from quite a number of different sources.\n17. Were the inter-vertex geometric distances normalized with respect to the “total” inter-vertex distance. Since the brain is changing in size across the lifespan, it is unclear how the brain size will affect the analysis on geometric distance.\n18. Why multi-step Procrustes alignment using age groups, instead of constructing the gradient template directly from the gradients of the individuals? 19. Cortical surfaces were generated using different pipelines and hence cortical thickness measurements are not consistent. For example, iBEAT surfaces and FreeSurfer might not be entirely consistent. The conclusions might be affected by this inconsistency.\n20. The fMRI data is not necessarily scanned under the same condition, even in the lifespan HCP datasets. For example, the BCP data for the first three years are scanned with the children asleep. How is this discrepancy handled to ensure the validity of findings?\n21. Phase I ends at one month postnatally, and Phase III begins at age five. These timepoints coincide with the end of the dHCP and the end of the BCP, respectively. The dHCP, BCP, and the remaining lifespan HCP studies use different imaging protocols. Therefore, it is likely that the demarcation of Phases I, II, and III might be an artifact of differences in image acquisition.\n22. How is longitudinal information, repeated observations of each subject, dealt with? Not doing so might bias the results.\n23. The effects of sex may vary nonlinearly across age, especially given the long developmental time span covered in this study. However, this nonlinearity does not appear to have been modeled.\n24. Section “Connectome gradient analysis”: It appears that group-level and lifespan reference gradients are constructed without first performing harmonization. Is this approach not affected by site effects?\n25. The analysis fundamentally relies on the averaging of connectivity matrices. The computation of gradients from a raw FC matrix is highly nonlinear. The output of this process on the mean of a collection of FC matrices is not the same as the mean of the output of the process on the collection of FC matrices. Is there worry that inter-individual variation is lost via this template creation strategy? Don’t these templates emphasize the central tendency of the group FC rather than the architectural features that diffusion embedding extracts on an individual level?\n26. Along those lines, averaging and then thresholding is not the same as thresholding then averaging. Did you examine the difference? Why did you choose the former?\n27. Did you examine the FC degree distribution before and after thresholding?\n28. Why are gradient extremes so large during infancy and then smaller during childhood?\n29. The authors mention that they use Procrustes alignment to achieve correspondence across gradients, but they do not mention whether they use scaling during alignment or not. This choice is important and can strongly affect downstream gradient metrics. Did you compare and if so, what did you find?\n30. Authors report 3 distinct phases identified via K means. K has to be chosen and the number of age bins across which you average the gradients would have a profound impact on the clusters identified. Is there any analysis of how cohorts intersect with these clusters?\n31. Cognitive “reference spectrum map” incorrectly applies adult-derived cognitive associations across the entire lifespan, ignoring the fact that cognition and its substrates are dynamic wrt age. 32. Cognitive analysis is circular.\n33. No significant correlation between myelin growth axis and SA growth axis undermines claims of unified multiscale framework.\n34. Knockout analysis doesn’t really tell us much. Removing somatomotor during early life and default mode during adulthood both having the largest effect is entirely expected given the topography of these gradients at these timepoints. Examining degree distributions alongside this knockout analysis may be informative.\n35. For the knockout analysis, the mean change in gradient value after removing entire resting state networks from the FC matrices is computed. This seems problematic, as the gradients must look entirely different after the removal of whole networks. How is alignment then carried out if the gradients look different? It would be more transparent to also show these network-ablated gradients. It’s not entirely clear that this analysis is meaningful in any way, as the networks that have the largest impact on the gradient values are the networks that anchor the gradient at opposite poles.\n36. Knockout analysis could be influenced by vertex counts / FC strength across networks.\nNeed to control for this somehow.\n37. The claim that the findings \"underscore the complex interplay between functional and structural changes\" is somewhat vague and generic. The actual evidence for this is a set of modest cross-sectional correlations.\n38. In line 808-809 you state “Our subsequent analyses focused mainly on the S-A gradient”, as far as I can see, only SA is analyzed?\n39. Did you try alternative template construction schemes? It would also be informative to see a plot and GAMLSS plot of cosine similarity to the template versus age.\n40. Analysis of gradients appears to be superficial, mostly based on summary statistics such as range and variance, although gradients provide rich spatial information of functional hierarchy. Interpretation and analysis of gradient topography and metrics are in general lacking.\n41. Any analysis on how fMRI time series length impacts gradients?\n42. Using only about 2000 vertices per hemisphere is too coarse. Please show example surfaces with this amount of vertices. Do the conclusions change with a greater number of vertices per hemisphere? It would be highly informative to see how these results change with higher resolution surfaces and with lower FWHM smoothing kernels.\n43. Please report the individual subject eigenspectra (specifically the gaps between lambda 1 and lambda 2, etc) across age. Large shifts in eigengap would suggest that one axis is not adequate to capture multi-scale changes.\n44. It would be informative to see what happens if you perform harmonization on the vertex-wise gradients themselves.\n45. Compute clustering and characteristic path length on weighted, signed networks (no binarization) to reduce circular dependency with the S-A axis. 46. The big conclusion is that the brain's functional organization develops from sensorimotor to association areas over the lifespan. Didn't we already know this? What is the profound new insight here, beyond putting a specific (and likely model-dependent) number on it?\n47. The authors use totally different, complex preprocessing pipelines for the infant and adult data. How do we know that the differences we see in gradients during infancy are not simply artifacts of this? The infant and adult cohorts were preprocessed with different slice-timing, despiking, and ICA-denoise settings. Where is the direct evidence that these pipeline disparities don’t dominate the apparent “infant inflation” of gradient extremes?\n48. NeuroSynth is based on meta-analysis and does not necessarily reflect direct cognitive measurements. How did you choose the specific 24 neurosynth terms you landed on out of hundreds of available terms? Did you do the analysis with multiple different sets of terms?\n49. Motion is markedly higher in children, yet gradients are built before regressing mean FD.\nCan you show scatterplots of gradient range versus FD by age band, or rerun the pipeline after scrubbing to equalize motion?\n50. The spline in GAMLSS is capped at 5 df regardless of vertex or metric. Sharper infant transitions could be smoothed away. What happens if you raise df or use penalized thin-plate splines with GCV?\n51. The K-means “lifespan phases” are accepted at face value. Where are the silhouette widths, gap statistic, or stability analysis across bin widths (e.g., 1-year vs 3-year bins)?\n52. The entire developmental narrative is built on gradients extracted from age-bin mean FC matrices. Yet diffusion embedding is a highly nonlinear transform, meaning that the gradients of the mean FC do not necessarily reflect the “mean” gradients computed from individual FC. To carry out any meaningful study of a particular gradient axis, a template must be constructed, but little justification is given in the paper on the choice of this particular scheme for template derivation. Some comparison with other template schemes, perhaps using a naive full-lifespan template may help justify the strategy used in the manuscript.\n53. The graph theoretic metrics are computed on the binarized row-wise 10% thresholded FC matrices while the gradients are computed with the 10% thresholded weighted FC matrices. Why were weighted implementations of characteristic path length and clustering coefficient not used and how does this impact the result?\n54. In Figure 2c, it is initially unclear that what is being reported is actually the mean magnitude of the gradient within each network. This ambiguity could be resolved by plotting the trajectories of the mean gradient values (not magnitudes) and marking the maxima (and minima for networks on the negative pole).\n55. Certain countries are restricted from accessing the Lifespan HCP datasets (see https://grants.nih.gov/grants/guide/notice-files/NOT-OD-25-083.html). Please confirm that the current study is compliant with the data use agreement of the NIMH Data Archive (NDA).", "reviewer_2": "This manuscript integrates a vast amount of data to investigate the development of the principal cortical hierarchy over the lifespan, from 32 PMW through 80 years of age. The primary finding is that functional organization proceeds in 3 phases over the lifespan, ultimately establishing then refining the S-A axis. The authors investigate how different parameters of the gradient develop globally as well as at the network and region level as well as how the gradient development relates to functional properties (graph theoretical metrics of FC) and structural properties (distance, thickness, and “myelin”). I commend the authors for their very impressive effort in pulling together this large lifespan dataset and for the many informative analyses (both primary and sensitivity) that were included in the paper. This is truly a feat and emphasizes what is possible with modern collaborative imaging data sharing and analysis workflows. On the whole, I think this is an informative and impactful work that will be of broad interest. That said, while the dataset and breadth of analyses are impressive, there are several methodological and conceptual choices that may limit interpretability and should be clarified prior to publication. I summarize these in the comments below.\n1. The authors continually refer the primary functional gradient as the S-A axis. However, the gradient itself looks quite different in Phase I vs the other phases. How do the authors determine that this principal gradient is indeed the S-A axis in all cases? Further, at least one prominent definition of the S-A axis comes from Sydnor et al. (Neuron, 2021) which defines an S-A axis based a variety of spatial maps. This definition of the axis has since been commonly used in published work. I would agree that the principal gradient derived here in the adolescents and adults looks very similar to that gradient, but the authors should comment on their criteria for determining if the result of their gradient analysis did indeed produce the “S-A axis” or something else. It would also be prudent to compare their derived axis to this other publicly available S-A axis as well as discuss the different manners in which this axis can be defined.\n2. Related, there were methodological choices made in the connectome gradient analysis that require more rationale. For example, the authors set all negative correlations in the FC matrices to zero. Why was this done and why are these correlations not deemed meaningful? How does this impact the outcome? Additionally, seemingly arbitrary decisions are made about which edges to include, for example thresholding at the top 10% of positive correlations without explanation. The authors should must describe how and why these choices were made and how these choices impact the outcomes. This is critical as these choices may bias results in ways that differ substantially by age—especially when the age-range spans from a newborn brain and an 80 year old brain. For example, it is possible that the distribution of positive and negative correlations varies substantially across the lifespan which would feed forward into downstream analyses.\n3. The analysis of how functional connectivity itself influences gradient development is an important one, particularly since the gradients are themselves based on the same FC data matrices. The focus on segregation and integration makes sense conceptually, but the reliance on graph theoretical properties to define segregation and integration abstracts the data quite a bit, and those results were the least clear. The graph theoretical approach taken here to look at segregation and integration imposes analysis decisions that may obscure informative patterns in the FC data (e.g. threshold FC matrices, binarizing values to remove variability in FC strength). This approach may miss interesting patterns that are dependent on the nuances of segregation and integration among different pairs of connections and communities (for example different combinations of changes in within and between network connectivity could potentially give rise to a similar overall integration and integration score). Recent work (e.g. Pines et al., 2022; Luo et al., 2024, both in Nature Communications) highlights how a more fine grained approach that looks at differences in FC strength more directly can yield important insights into the establishment and embedding of the cortical hierarchy. Given the prominence and highly related nature of these recent studies, the authors should consider adding a parallel approach to their investigation of FC differences underlying gradient development in the present work. At minimum, they should contextualize their findings in relation to this recent work, which is notably absent from the current manuscript references.\n4. The examination of the structural properties that may contribute to gradient development is a strength of the paper. I have a minor comment that the authors are overly narrow in their labelling of the T1/T2 ratio as “myelin” or “myeloarchitecture” per se. Of course, this ratio measure is sensitive to myelin differences, but it is also sensitive to anything else that may influence T1 and T2 relaxation times, which could be a myriad of things over this wide developmental span. As such, the authors should avoid labelling differences in T1/T2 effects as differences in “myelin” specifically. This work would be equally as impactful if they simply refer to these effects as differences in T1/T2 contrast ratio and then discuss the factors that could drive this signal.\n5. Related to the findings on geometric distance, the authors claim that “These findings suggest that while geometric constraints strongly influence the functional hierarchy, their influence gradually diminishes during development, potentially facilitating increased functional flexibility.” The last part of this conclusion “facilitating increased functional flexibility” bears more scrutiny and explanation. How was this conclusion drawn and how might diminishing influence of geometric distance influence functional flexibility?\n6. Minor: Lines 309-310. The text states “Regional clustering analysis (Cp-node) revealed a progressive decline in spatial topography from primary to transmodal association regions. What does a decline in spatial topography mean in this context? This phrasing could be clarified.", "reviewer_3": "This is a well-written and timely study that leverages a large neuroimaging dataset across a wide age range (32 postmenstrual weeks to 80 years), addressing an important question regarding the development of the cortical functional hierarchy (sensorimotor-association axis). The sample size and cross-site design are impressive. The identification of developmental shifts in the S-A gradient across the lifespan nicely summarizes recent findings and further extends them by investigating associations with cortical thickness, area expansion, myelination, functional integration and segregation, as well as cognition. I have several comments for the authors to consider:\n• The S-A gradient has been shown to account for an increasing proportion of variance in the functional connectome from childhood to adulthood, and it is not consistently the dominant gradient across development (Lariviere et al., 2019; Xia et al., 2021; Dong et al., 2021). Since the current study spans the full lifespan, it would be valuable to examine how the explained variance of the S-A gradient (and the ordering of dominant vs. secondary gradients) shifts across developmental stages.\n• This study identified age-related shifts in the S-A gradient consistent with previous findings, especially during early development. The authors propose an interesting idea that this shift is likely driven by increasing functional specialization and the strengthening of long-range connections. Could the authors provide direct evidence supporting this hypothesis?\n• Given that network configuration (e.g., size and distribution) changes over development (Sun et al., 2025, Nature Neuroscience, from the same research group), it is unclear why the authors opted to use the adult Yeo 7-network configuration to calculate segregation and integration. Would age-specific network definitions yield more accurate estimates? I noticed that the functional segregation growth rates reported here (Fig. 3c) and those in Sun et al. (2025) appear to differ. Could the use of age-specific network configurations help explain this discrepancy?\n• A related point: the spatial distribution of the S-A gradient changes across the lifespan, transitioning from an anterior–posterior configuration in prenatal/neonatal stages to a mature axis spanning somatomotor/visual to default mode regions. Could the authors elaborate further on how the topographic organization of the S-A gradient relates to age-related changes in network topology?\n• Fig. 3d appears to show that integration growth also aligns with the S-A axis (e.g., increased integration in the default mode network and decreased integration in somatomotor regions). However, only a weak association was reported between the growth axis of integration and the S-A gradient (r = 0.19). Could the authors discuss this finding further?\n• The authors performed Procrustes rotation for cross-subject/cross-age comparisons. Since the raw gradient may reflect individual or age-specific patterns, it would be helpful if the authors could provide the extent and effect of the rotation procedure? For example, how much rotation was applied to align the S-A gradient across individuals and age groups?\nMinor\n• Fig. 1c: The meaning of “dimension 1” and “dimension 2” is unclear and may be confused with the gradient dimensions. Please define these explicitly.\n• Fig. 5a: The meaning of “z-stat % (max = 15)” and the labels “25, 50, 100” are unclear.\n• Some sites appear to include longitudinal data. How were longitudinal data handled in the analysis?\n• Please clarify the abbreviation “rm-ANOVA” as “repeated-measures ANOVA.”"}, "author_response": "We sincerely thank the editor and all the reviewers for their comments on our manuscript, titled “Spatiotemporal dynamics of the human cortical functional hierarchy across the lifespan” (ID: NCOMMS-25-33712-T). We found the comments to be constructive and highly valuable for improving our manuscript, which we revised based on the comments received. Our point-by-point responses to the reviewers’ comments are presented below.\nReviewer #1 (Remarks to the Author):\nIn this paper, the author maps the development of the sensory-association (SA) gradients across the human lifespan and explores how these gradients relate to structural brain measures. While the characterization of the SA gradient is a valuable contribution to understanding brain organization and maturation, I have major concerns regarding several aspects of the study. These include questions about the novelty of the work in light of existing literature, the appropriateness and rigor of the methodological approaches, the clarity and robustness of the reported results, and the strength and validity of the conclusions drawn from the findings.\nR: We thank this reviewer for all positive comments and insightful comments/suggestions.\nComments:\nQ1.1 Introduction - “First, existing studies focus mainly on discrete developmental stages (e.g., the third trimester 19 , neonatal period, adolescence, and adulthood), without provided a continuous spatiotemporal perspective from prenatal development through late adulthood.” This is not true and misleading in light of Ref 23.\nR: We have rephrased the relevant descriptions and hopefully it is clearer now (page 4):\n“First, most existing studies of the functional hierarchy have focused on discrete developmental stages (e.g., the third trimester 19, neonatal period 18, adolescence 20,21 and adulthood 3,22), resulting in a fragmented perspective that obscures fine-grained transitions and nonlinear dynamics of cortical gradients across the full age continuum. Moreover, prior work has often been limited by small (e.g., N < 100) 18,19 or modest (N = hundreds to thousands) 16,20–22 sample sizes, along with a predominant reliance on linear modelling to characterize gradient development 19,20,22. Although a recent study 23 established a continuous, nonlinear mapping of functional gradients from birth to 100 years of age, the inferred developmental patterns remain constrained by a modest sample size (N = 3,556) and particularly sparse data coverage during childhood.”\nQ1.2 Introduction - “Furthermore, methodological limitations, such as small (less than one hundred) or modest (hundreds to thousands) sample sizes and an overreliance on simplistic linear modelling approaches, fundamentally constrain our ability to capture nonlinear neurodevelopmental dynamics.” This is also not true. The methods used for data processing and analysis (FreeSurfer, GAMLSS, etc.) in the current work are mostly publicly available and are quite standard. There are no methodological hurdles preventing other groups from carrying out similar work.\nR: We would like to clarify this point raised by this reviewer. We agree that both standard preprocessing pipelines (e.g., FreeSurfer) and advanced nonlinear modeling frameworks (e.g., GAMLSS) are well-established within the community. Our aim was not to imply that these tools are non-standard, but to highlight that many prior studies have yet to fully integrate large-scale datasets with nonlinear modeling to characterize developmental trajectories in the cortical hierarchy. We have revised the sentence accordingly (page 4):\n“First, most existing studies of the functional hierarchy have focused on discrete developmental stages (e.g., the third trimester 19, neonatal period 18, adolescence 20,21 and adulthood 3,22), resulting in a 1 fragmented perspective that obscures fine-grained transitions and nonlinear dynamics of cortical gradients across the full age continuum. Moreover, prior work has often been limited by small (e.g., N < 100) 18,19 or modest (N = hundreds to thousands) 16,20–22 sample sizes, along with a predominant reliance on linear modelling to characterize gradient development 19,20,22.”\nQ1.3 Introduction - “Second, current conceptual frameworks predominantly emphasize functional gradient remodelling during development while largely overlooking its structural foundations. In mature brains, the S-A axis aligns spatially with structural features indexed by cortical thickness, intracortical myelination, and intrinsic geometry.” This is not true. Association between primary gradients and structural gradients has already been investigated for example in reference 23.\nR: We thank the reviewer for raising this point. We have revised the text to clarify this (page 4-5):\n“Second, current research frameworks lack the capacity to delineate how functional hierarchies coevolve with specific structural features, such as cortical thickness 2,6,24, intracortical myelination 2,8,9, and intrinsic geometry 3,25, across the lifespan. While the S-A axis aligns spatially with these structural attributes in mature brains, the developmental synchrony between functional and structural maturation remains largely unclear, particularly given their distinct maturational timelines in early life 2. A previous study 23 reported coupling between lifespan changes in functional gradients and a composite structural gradient derived from multiple features; however, its analysis strategy that integrates multiple structural properties prior to dimensionality reduction does not directly assess the contribution of specific structural properties to the emergence of functional gradients.”\nQ1.4 Why 26 age groups? This creates artificial temporal demarcations between different life periods that do not necessarily reflect true biological development. This is particularly problematic when different brain regions have distinct development timelines that do not coincide with these artificial temporal partitions. It is uncertain if age grouping will result in artificial disconunity in the analysis, not just at the global whole-brain level but also local vertex level.\nR: We thank the reviewer for raising this point. Currently, there is no universally accepted standard for defining age groups. Our selection of 26 age groups was informed by developmental principles indicating that brain maturation occurs rapidly during infancy and childhood but slows progressively through adolescence and adulthood. Therefore, we applied finer age partitions at younger stages and broader intervals at later ages to better capture these nonlinear developmental dynamics.\nTo further address the concern regarding potential artificial discontinuities, we also conducted complementary analyses at the individual level without imposing predefined age bins. Specifically, functional gradients were computed for each participant and then clustered across all individuals, using chronological age as a quasi-continuous variable. This individual-level approach confirmed the main findings of three phases of functional gradient development: (Phase I) early anterior-posterior differentiation during infancy, (Phase II) emergence of the S-A gradient in early childhood, and (Phase III) its gradual expansion, refinement, and stabilization throughout adolescence and adulthood. The convergence between the group-based and individual-level analyses indicate that our conclusions are not affected by the chosen age-grouping scheme.\nThe Supplementary analyses were added in the revised manuscript (page 7):\nResults: “The three developmental phases were also identified through both individual-level gradient clustering without predefined age bins and within-age-group-harmonized analysis (Supplementary Fig. 3, 4). This indicates that the observed patterns were robust and not driven by specific age-grouping schemes or site related confounds.”\nSupplementary Results: “To validate the robustness of our findings, we first verified that developmental phases identified in the group-based analyses were consistently reproduced in individual-level clustering without predefined age bins, confirming that the observed patterns were not driven by the specific age-grouping scheme (Supplementary Fig. 3).”\nSupplementary Fig. 3 | Distinct growth stages of the sensorimotor-association gradient across the lifespan identified from individual gradients. a, K-means clustering analysis based on the pairwise correlations between the individual gradient maps, resulting in three clusters. b, Three distinct phases of the S-A gradient are visualized via multidimensional scaling (MDS). “Dimension 1” and “dimension 2” are defined as the first and second axes obtained through MDS of the functional distance matrix.\nQ1.5 How does the choice of 26 age groups influence the clustering results? The patterns shown in Figure 1b may vary not only with the number of age groups but also with how these groups are distributed over time. Moreover, the age grouping adopted from the authors’ previous publication may not be appropriate here, as this study focuses on gradients rather than simple features of functional connectivity. The authors should demonstrate that this age grouping does not introduce bias by comparing it against clustering results obtained from the entire dataset without imposing predefined age groups.\nR: To evaluate whether the clustering results were influenced by the choice of age-group definitions, we conducted an additional analysis at the individual level across all participants without imposing any predefined age bins. Functional gradients were derived separately for each participant, and clustering was then performed on the complete dataset using chronological age as a continuous variable.\nThis individual-level analysis reproduced the same three major developmental phases identified in the age-binned analysis: (Phase I) early anterior-posterior differentiation during infancy, (Phase II) emergence of the sensorimotor-association (S-A) gradient in early childhood, and (Phase III) progressive expansion and stabilization of this gradient through adolescence and adulthood. The high consistency between the group-based and individual-level results demonstrates that the choice of 26 age groups does not bias the clustering outcomes. These findings confirm the robustness of our developmental trajectory characterization and validate that the observed patterns are not artifacts of grouping strategy.\nQ1.6 How is the difference in sample size across age groups dealt with? The authors did a “balanced resampling analysis” by equal resampling across age bins. However, it is unclear why the 16 “age bins” are different from the 26 “age groups”. Also, does the fact that resampling uniformly with a smaller number of samples per age bin still replicates the main results suggest that the large sample size used in the study may not be necessary?\nR: We complemented the age-binned analysis with individual-level clustering that does not depend on predefined bins. This analysis again confirmed the same three developmental phases, demonstrating that the main findings are robust to differences in group size and binning scheme.\nThe use of 16 “age bins” in the balanced resampling analysis reflects a methodological adjustment aimed at ensuring an adequate and comparable number of participants in each bin for statistical resampling, rather than a conceptual change from the 26 developmental age groups used in the main analysis.\nRegarding the necessity of a large sample size, as detailed in our response to Q1.8, larger samples substantially improve the stability and reliability of lifespan growth curve estimates, particularly at the regional and vertex levels, by reducing variance and enhancing the precision of developmental trajectory characterization. Thus, while smaller, balanced subsamples supports the validity of our conclusion, the full sample remains essential for achieving optimal reliability and generalizability of brain gradient trajectories.\nQ1.7 How are differences in data variances between datasets and across time dealt with?\nThere can be nonlinear variance changes across age groups and within age groups and also across the different datasets and within the datasets.\nR: We took several steps to minimize potential biases introduced by variance differences across datasets and age groups in the original analyses.\nFirst, at the preprocessing stage, all datasets were processed using identical pipelines: projection to fsaverage_LR32k, nuisance regression with 24 motion parameters plus WM, CSF, and global signal, linear detrending, temporal bandpass filtering (0.01-0.08 Hz), scrubbing (FD > 0.5 mm with one preceding and two subsequent volumes interpolated), surface smoothing (FWHM = 6 mm), and resampling to fsaverage4 space. These harmonized pipelines ensured comparable signal properties across datasets.\nSecond, at the modelling stage, we explicitly accounted for age- and site-related heteroscedasticity using the GAMLSS framework with a flexible skew exponential power type 4 (SEP4) distribution. Specifically, we modeled the variance (scale parameter, σ) as a smooth function of age and sex (sigma.fo=~bs(Age) + Sex), which accommodates nonlinear variance changes across the lifespan. Site was further included as a random effect to capture both mean and variance differences across datasets.\nThird, model diagnostics confirmed the adequacy of our approach. We examined normalized quantile residuals, residuals-fitted plots, kernel density estimates, and normal Q-Q plots (Supplementary Fig. 28), all of which indicated well-behaved residuals without systematic trends. Detrended Owen’s plots further showed that residuals remained within nonparametric confidence intervals (Supplementary Fig. 29). Together, these diagnostics suggest that the GAMLSS effectively captured heteroscedasticity.\nFinally, sensitivity analyses using a leave-one-site-out strategy demonstrated that our main findings were highly consistent across datasets, supporting the robustness of our results against potential variance differences.\nQ1.8 A large sample size like the one used in this paper is not always accessible to all researchers. Please show empirical evidence that such a large sample size gives conclusions that cannot be drawn using a smaller sample size.\nR: We thank the reviewer for raising this important point. To empirically evaluate the benefits of large sample sizes, we performed a resampling-based analysis while controlling for both age and sex distributions. Specifically, the lifespan range (32 PMW to 80 years) was divided into 10 equal intervals, and subsamples were drawn while preserving the original proportions of age and sex within each interval. For global- and system-level analyses, sample sizes were varied from 200 to 33,000 participants, with finer increments at smaller sizes. For regional-level analyses, 13 predefined sampling conditions were used. For each sample size, we performed 50 independent resampling iterations, refitted the lifespan growth curves, and computed the mean absolute error (MAE) and root mean squared error (RMSE) relative to the reference curve estimated from the full dataset. The minimum sample size was defined as the point at which the observed MAE/RMSE decreased to within 1.5 times the error estimated at the maximum sample size for regional-level gradient values.\nThe results showed that both MAE and RMSE decreased as sample size increased. At the global and system levels, error dropped sharply when sample size grew from a few hundred to ~5,000 participants, after which gains plateaued. At the regional level, required sample sizes varied across cortical vertices. A sample size of at least ~6,000 participants was needed for more than 50% of cortical vertices to achieve stable error levels, while ~25% of vertices required even larger samples.\nIn sum, these analyses demonstrate that larger sample sizes substantially enhance the stability and reliability of lifespan growth curve estimates, particularly at the regional scale. These results provide empirical guidance on the minimum sample sizes required for robust conclusions in future lifespan neuroimaging studies. We have added relevant results in the revised manuscript (page 16):\n“To evaluate the influence of sample size on our findings, we conducted a resampling-based sensitivity analysis, varying sample sizes from 200 to 33,000 participants. We observed that both mean absolute error (MAE) and root mean squared error (RMSE) decreased with increasing sample size (Supplementary Fig. 17). Notably, errors at the global and system levels dropped sharply until about 5,000 participants, after which gains plateaued. Furthermore, at the regional level, a sample size of ~6,000 participants is necessary to achieve stable error estimates for more than half of the cortical vertices, while a quarter of vertices required even larger samples. This work confirms that larger samples are crucial for reliable gradient growth curves, especially at the regional scale, and provide valuable guidance for future lifespan neuroimaging studies.”\nSupplementary Fig. 17 | Relationship between sample size and fitting error for gradient measures. a, Relationship between sample size and fitting error for global-level gradient measures (MAE shown in blue, RMSE in red). b, Relationship between sample size and fitting error for system-level gradient measures (MAE in the left panel, RMSE in the right panel). c, Relationship between sample size and fitting error for regional-level gradient values. The minimum sample size at which the MAE (left) /RMSE (right) reached a stable plateau for each vertex is mapped onto the cortical surface. 5 The pie charts summarize the distribution of these minimum sample size across all vertices. Sample sizes ranged from 200 to 33,000 participants, with 50 resampling iterations per size, and age/sex distributions were preserved across 10 lifespan intervals (32 PMW to 80 years).\nQ1.9 Anterior-posterior gradient during early infancy – This is contradictory to previous observation that the SA gradient already exists at birth. The SA gradient does not start to occur at 3 months but before that.\nR: We thank the reviewer for this valuable comment and the opportunity to clarify. In our data, we observed a prominent anterior-posterior (A-P) gradient during the prenatal and neonatal periods, which progressively transitions into a more differentiated S-A gradient after approximately 3 months of age. This observation of an A-P gradient in early infancy is consistent with prior studies. For instance, Xia et al. (2024) reported the emergence of an anterior-posterior gradient as early as 34.6 postmenstrual weeks 1. Lariviere et al. (2020) described anterior-posterior variations in term-born neonates, interpreting these as an immature precursor of the differentiation between unimodal and transmodal areas 2.\nThe apparent inconsistency with Ref23 likely reflects methodological differences. Specifically, Ref23 employed a generalized additive mixed model (GAMM) to fit gradient development across the entire age range. Given the brevity of the early developmental window and the relatively small sample size available, such modeling may smooth or shape the less differentiated early-stage patterns into a more continuous S-A trajectory. In contrast, our approach used group-level analyses within discrete age bins, enabling us to directly characterize gradient architecture at each developmental stage without imposing model-based continuity across age. We have revised the text by adding a small discussion (page 18):\n“Notably, we found that the functional topography at birth is not characterized by the S-A gradient but is instead dominated by an anterior-posterior differentiation. This finding is aligned with previous studies of preterm, term-born, and neonatal cohorts 18,19, suggesting that this pattern represents a developmental precursor to the mature S-A gradient. In contrast, one recent study 23 suggests that the S-A gradient is already present at birth. The discrepancy may partly arise from methodological differences: the previous study 23 applied lifespan generalized additive mixed modeling, which could smooth over less differentiated early-stage patterns. Our stage-specific analyses, by comparison, capture gradient organization directly at each period, offering a more fine-grained characterization of early functional topography.”\nQ1.10 Grouping the broad age range from 5 years to 80 years as a single phase (Phase III) is too coarse and does not reflect the fact that growth inflection points happen quite often during adolescence.\nR: We would like to further clarify this point. Our analysis for defining Phase III was not to suggest that age span from 5 to 80 years constitutes a homogenous developmental stage. Rather, Phase III was conceptualized as a broader period of gradient stabilization and refinement-a stage in which the S-A gradient, once established, undergoes progressive fine-tuned throughout childhood, adolescence, and adulthood, and then continues to evolve into late life. We have clarified this point in the revised manuscript (page 19):\n“our definition of Phase III (5-80 years) reflects the prolonged refinement and stabilization of the S-A gradient from childhood, adolescence to adulthood, rather than developmental uniformity through this period.”\nQ1.11 Datasets like the developing HCP (dHCP) use gestational age. Postnatal datasets such as the BCP use chronological age. How is the discrepancy between gestational age and chronological age resolved?\nR: To ensure consistency across datasets, we converted gestational age in the dHCP dataset into chronological age by treating 40 gestational weeks as the time of birth (i.e., 0 years). Chronological age was then calculated as:\n(𝐺𝑒𝑠𝑡𝑎𝑡𝑖𝑜𝑛𝑎𝑙 𝑎𝑔𝑒 𝑎𝑡 𝑠𝑐𝑎𝑛−40)×7\n𝐶ℎ𝑟𝑜𝑛𝑜𝑙𝑜𝑔𝑖𝑐𝑎𝑙 𝑎𝑔𝑒 (𝑦𝑒𝑎𝑟𝑠)=\n“To align the dHCP dataset with other cohorts, we converted gestational age to chronological age. The conversion was based on the assumption that 40 gestational weeks corresponds to birth (0 years of chronological age). The specific formula applied was: Chronological age (years) = (Gestational age at scan – 40) × 7 / 365. This procedure ensured all ages were referenced to a consistent postnatal timeline.”\nQ1.12 “Collectively, these results revealed previously unreported spatiotemporal patterns and critical neurodevelopmental milestones of the functional hierarchy, as characterized by the initiation, establishment, and expansion/stability of the S-A axis of cortical functional organization.” Again, this is untrue in light of ref 23.\nR: We appreciate the reviewer’s comment. The related statements have been removed from the revised version.\nQ1.13 Range and variance are measuring almost the same thing, i.e., the spread of the gradient values, and are therefore redundant. Please pick one to report.\nR: While both range and standard deviation quantify aspects of distribution spread, they emphasize distinct properties. The range captures the full extent of gradient differentiation (maximum-minimum), highlighting the span of the axis and sensitivity to extreme values. In contrast, the standard deviation reflects the average dispersion of all gradient values (not just the maximum and minimum) around the mean, providing a more robust, outlier-resistant measure of variability. For example, consider two groups of numbers with the same range: one group is tightly clustered around the mean, and the other is spread out more evenly between the minimum and maximum values. Both groups have the same range, but the first group will have a smaller standard deviation, indicating less variability around the mean, whereas the second group will have a larger standard deviation, reflecting greater dispersion. Because these two indices offer complementary perspectives‒one emphasizing the maximal differentiation and the other overall dispersion‒we report both to provide a fuller characterization of developmental and topographical features of the gradient.\nQ1.14 In Figure 2c and throughout the paper, it is unclear what “gradient score” is referring to.\nR: To improve clarity, we have revised the terminology throughout the manuscript and figure legends, replacing “gradient score” with “gradient value”. This change more accurately reflects the numerical values derived from the gradient decomposition.\nQ1.15 The results shown in Figure 2e seem to indicate that the SA gradient changes very little after 20 years. This contradicts Figs 2a and 2b and other results shown in the paper.\nR: We would like to make a clarification about this point. First, the apparent stability of the S-A gradient after age 20 in Fig. 2e was partly due to the original colormap design rather than the underlying results. Specifically, the broad white center of the colormap compressed the visual contrast around zero values, making age-related changes less discernible. In the revised figure (see updated Fig. 2e), we applied a colormap with more balanced contrast to better capture subtle variations across adulthood.\nSecond, our analyses do indicate that the S-A gradient continues to change after early adulthood, although at a slower rate compared with the rapid alterations observed in early childhood. These adult changes are consistent with the global- and system-level results shown in Fig 2a and 2b, which demonstrate relative stabilization of gradient measures from the third decade onward to early aging, while still exhibiting gradual refinements at the regional levels.\nFig. 2 | Lifespan growth of the sensorimotor-association gradient. a, Lifespan growth curves (top panel) and growth rate curves (bottom panel) for global gradient measures. The solid line represents the 50th centile (median), while the dotted lines represent the 5th, 25th, 75th, and 95th centiles. The growth rate was estimated by the first derivative of the median growth curve, with a 95% confidence interval (shaded in gray) obtained from 1,000 bootstrap resamplings. The horizontal dotted line (y=0) marks the transition from growth to decline, with the peak age indicated by a purple dot. b, Lifespan growth curves of gradient range and standard deviation for each functional system defined according to the age-specific Yeo 7-network atlas 30,34,35. c, Lifespan growth of the gradient value for each functional system. The peak age is marked by a black dot. d, Contributions to lifespan growth of the S-A gradient from each functional system, estimated via knockout-based simulation analysis (left panel). A repeated-measures ANOVA tested variations in system contributions across developmental phases (right panel). ***, p<0.001. e, Lifespan growth rate maps of the S-A gradient at representative ages. f, Lifespan growth axis of the S-A gradient. The left panel shows regional gradient values across age. PCA identified the first component as the dominant lifespan growth axis, explaining 65.4% of the variance (middle panel). The growth axis was divided into 20 decile bins, and mean gradient values were computed for all vertices within each bin (right panel). g, Spatial correlation between the lifespan growth axis and the evolutionary hierarchy. Evolutionary hierarchy was quantified using macaque-to-human cortical expansion 32. Statistical significance was assessed via 1,000 non-parametric spin permutation tests 36. CI, confidence interval; VIS, visual; SM, somatomotor; DA, dorsal attention; VA, ventral attention; LIM, limbic; FP, frontoparietal; DM, default-mode; w, week; m, month; yr, year.\nQ1.16 T1-weighted and T2-weighted images are not quantitative and can have arbitrary signal gains associated with different scanners, sequences, and choices of acquisition parameters. It is unclear how this is dealt with in this study, which involves datasets from quite a number of different sources.\nR: We appreciate the reviewer’s insightful comment. It is indeed true that T1-weighted and T2-weighted images are not inherently quantitative, and their signal intensities can vary across scanners, sequences, and acquisition protocols.\nNevertheless, the T1w/T2w ratio provides a reliable proxy for relative cortical myelin content because myelin consistently shortens T1 relaxation times (increasing T1w signal) and T2 relaxation times (decreasing T2w signal), producing a systematic contrast that is amplified when taking their ratio 6,7. This approach has been widely adopted in developmental, clinical and comparative studies 6,8–11. Importantly, several studies have demonstrated that the T1w/T2w ratio produces robust and anatomically meaningful spatial patterns that remain consistent across different scanners and acquisition settings. For example, Glasser et al. reported that group-average datasets acquired on different scanners with different sequences at independent sites yielded qualitatively very similar results 6,7.\nIn the present study, we took several steps to further mitigate site-related variability. First, during preprocessing we applied biases and intensity normalization, which reduce transmit field biases and minimize intensity inhomogeneities. These steps are critical for improving cross-site comparability of T1w/T2w-derived myelin maps. Second, in our statistical modeling, we explicitly included site as a random effect in the GAMLSS framework, thereby adjusting for potential residual site-specific differences when estimating growth trajectories.\nTogether, these strategies substantially reduce concerns regarding inter-scanner variability and support the robustness of our T1w/T2w-based myelin estimates across multiple datasets. In the revised manuscript, we have added relevant description and hopefully it is clearer now (page 29).\n“To mitigate inter-scanner variability and enhance cross-site comparability, all T1w and T2w images underwent bias correction and intensity normalization prior to the computation of the T1w/T2w ratio.……. we fitted the growth curves of the structural attributes for each vertex via GAMLSS with age and sex as fixed effects, and site as a random effect.”\nQ1.17 Were the inter-vertex geometric distances normalized with respect to the “total” inter-vertex distance. Since the brain is changing in size across the lifespan, it is unclear how the brain size will affect the analysis on geometric distance.\nR: We thank the reviewer for raising this point. To account for inter-individual differences in brain size and their potential impact on geodesic distance measures, we normalized geodesic distances by the square root of the cortical surface area. This was done when analyzing distances between vertices within the top and bottom 5% of S-A gradient values. We then re-fitted the lifespan growth curves using these normalized distances, and the results remained highly consistent with our original findings. We have rephrased the relevant sections in Methods and Results as follows:\nMethods (page 28): “To account for inter-individual differences in brain size, geodesic distances were normalized by the square root of the cortical surface area.”\nResults (page 12): “For each individual, we calculated the mean normalized geometric distances between paired vertices with extreme gradient values. This geometric distance exhibited a nonlinear trajectory across the lifespan, peaking during young adulthood (20.7 years, 95% bootstrap CI: 17.2–23.2), followed by a gradual decline (Fig. 4b).”\nQ1.18 Why multi-step Procrustes alignment using age groups, instead of constructing the gradient template directly from the gradients of the individuals?\nR: We adopted a multi-step Procrustes alignment based on age groups to reduce potential biases arising from sample size imbalances and age-related variations in gradient topography. If a group-level gradient template were constructed directly from all individual embeddings, the resulting template would likely be disproportionately influenced by age groups with larger sample size and might obscure meaningful developmental stage-specific differences. By first aligning gradients within narrower age groups and then iteratively converging toward a common template, our approach ensures that the reference template captures both shared organizational features and age-specific variations in gradient topography. This strategy provides a more balanced and representative framework for comparing gradients across the lifespan. We have rephrased the relevant sections (page 23):\n“To ensure accurate comparisons between groups, we applied Procrustes rotation 79 to align age-specific group-level functional gradients to iterative group-level reference gradients……Compared with constructing the template directly from individual gradients, this multi-step alignment strategy reduces biases from sample size imbalances and development variations in gradient topography.”\nQ1.19 Cortical surfaces were generated using different pipelines and hence cortical thickness measurements are not consistent. For example, iBEAT surfaces and FreeSurfer might not be entirely consistent. The conclusions might be affected by this inconsistency.\nR: We used different cortical surface reconstruction pipelines for different age groups (iBEAT for infants and FreeSurfer for older participants). This decision was intentional 12–14, as each tool has been specifically optimized for the developmental stage it targets. Applying a single pipeline across the entire age range would likely introduce greater bias, since adult-oriented tools such as FreeSurfer are suboptimal for infant brains, while infant-tailored pipelines are not designed for adult anatomy.\nPrevious validation studies support the appropriateness of this strategy 12–14. For example, Wang et al. (2023) demonstrated that iBEAT V2.0, a deep learning-based pipeline optimized for infant cortical reconstruction, produces accurate and reproducible surfaces across scanners and imaging protocols 12. Similarly, Makropoulos et al. (2018) established a neonatal minimal processing pipeline designed to capture the unique anatomical features of early development 13. For older participants, the widely adopted Human Connectome Project pipeline and FreeSurfer have been rigorously validated in multisite settings and are considered standards for adult cortical surface analysis 14.\nBy using pipelines tailored to the anatomical characteristics of each age group, we aimed to minimize reconstruction errors and enhance the reliability of cortical thickness estimation across the lifespan. All subsequent analyses were conducted on surfaces registered to a common template space, further mitigating concerns about cross-pipeline inconsistency. We have rephrased the relevant sections in the Supplementary Methods section (page 2):\n“This age-appropriate pipeline selection minimizes cortical reconstruction bias across the lifespan, and all surfaces were subsequently registered to a common template to ensure cross-pipeline comparability.”\nQ1.20 The fMRI data is not necessarily scanned under the same condition, even in the lifespan HCP datasets. For example, the BCP data for the first three years are scanned with the children asleep. How is this discrepancy handled to ensure the validity of findings?\nR: We thank the reviewer for raising this point. Indeed, scanning conditions necessarily vary across lifespan datasets due to practical and ethical considerations. For example, in the Baby Connectome Project (BCP), infants aged 0-35 months were typically scanned during natural sleep, whereas older children could be scanned either awake or asleep depending on feasibility 15. Similarly, in the Developing Human Connectome Project (dHCP), virtually all infants were scanned during natural sleep. In contrast, large-scale adolescent and adult cohorts such as ABCD, HCP, ADNI, CamCAN, and UK Biobank explicitly required participants to remain awake during resting-state scans (eyes open or eyes closed, depending on protocol). Unfortunately, precise per-scan records of sleep/awake state are not publicly available for dHCP and BCP, preventing direct stratification by brain state.\nWe fully acknowledge that brain state can influence functional activity patterns 16 and may affect gradient estimation. However, several factors mitigate this concern. First, the gradients observed show robust and biologically plausible age-related trajectories that are consistent across datasets and scales, suggesting that they are not solely driven by state differences. Second, all datasets underwent harmonized preprocessing pipelines, which help minimize systematic biases. Third, the vary large sample size provides statistical stability, reducing the influence of individual-level variation in scan conditions. We now explicitly note this as a limitation in the revised manuscript (page 21):\n“Finally, the large fMRI datasets exhibited heterogeneity in scanning states: neonate and infant data were primarily acquired during natural sleep, whereas children and adults were scanned during wakefulness. Given differences in functional brain activity between sleep and wake states 73, we cannot rule out the potential influence of this acquisition discrepancy on the observed patterns. Future studies that control for this variable by using identical scanning states across age groups will be important to clarify this issue.”\nQ1.21 Phase I ends at one month postnatally, and Phase III begins at age five. These timepoints coincide with the end of the dHCP and the end of the BCP, respectively. The dHCP, BCP, and the remaining lifespan HCP studies use different imaging protocols. Therefore, it is likely that the demarcation of Phases I, II, and III might be an artifact of differences in image acquisition.\nR: We thank the reviewer for pointing out this point. To assess whether the phase boundaries coinciding with the transitions between the dHCP, BCP, and HCP datasets could be driven by differences in acquisition protocols, we explicitly evaluated the potential influence of such site and protocol effects. A single harmonization across the entire lifespan would not be appropriate in this context, as the topographic architecture of functional connectivity systematically changes with age. Applying one global harmonization procedure could inadvertently distort age-specific topographic features or impose artificial similarity across distinct developmental stages.\nInstead, we minimized acquisition-related variability by removing site-specific effects within each age group before gradient estimation. Specifically, for each age group, site was modeled as a random effect in the vertex-wise FC data and regressed out prior to constructing group-level connectivity matrices and computing gradients. When gradients and clustering were derived from these within-age-group-harmonized data, the same three developmental phases (early anterior-posterior differentiation; emergence of the S-A axis in early childhood; and subsequent expansion, refinement, and stabilization) were reproduced and closely matched the results reported in the main text. Moreover, when determining optimal number of clusters, both 3- and 4-cluster solutions were identified with equal frequency across clustering indices. We reported the 3-cluster solution in the main text, while the 4-cluster solution only further subdivided the earliest phase (32w to 3 months). This concordance indicates that the observed phase boundaries reflect genuine developmental transitions rather than artifacts arising from dataset boundaries or acquisition differences. We have added a corresponding description in the revised Results (page 7) and Supplementary Results (page 3):\nResults: “The three developmental phases were also identified through both individual-level gradient clustering without predefined age bins and within-age-group-harmonized analysis (Supplementary Fig. 3, 4). This indicates that the observed patterns were robust and not driven by specific age-grouping schemes or site related confounds.”\nSupplementary Results: “To ensure that the phase boundaries were not artifacts of differences in acquisition protocols across datasets, site was modeled as random effect and regressed out within each age group prior to computing group-level connectivity matrices and gradients. Analyses using these harmonized data reproduced the same three phases, indicating that the observed phase boundaries reflect genuine developmental transitions rather than dataset or acquisition differences (Supplementary Fig. 4).”\nSupplementary Fig. 4 | Distinct growth stages of the sensorimotor-association gradient across the lifespan from the harmonized functional connectome. a, S-A gradient maps generated for 26 distinct age groups derived from the harmonized functional connectome. b, K-means clustering analysis based on the pairwise correlations between the 26 group-based gradient maps, resulting in three clusters (main text) and an alternative four-cluster solution that further subdivides the earliest phase. c, Three distinct phases of the S-A gradient are visualized via multidimensional scaling (MDS). “Dimension 1” and “dimension 2” are defined as the first and second axes obtained through MDS of the functional distance matrix. d, Distribution of gradient values in the S-A axis across lifespan. w, week; m, month; yr, year.\nQ1.22 How is longitudinal information, repeated observations of each subject, dealt with? Not doing so might bias the results.\nR: To minimize potential bias introduced by repeated measurements from the same participants, we restricted our analyses to a single imaging session per individual across all longitudinal datasets. This ensure that each participant contributed only once to the cross-sectional lifespan analyses, thereby preventing inflated sample sizes or undue weighting of individuals with repeated scans. This procedure is described in the Methods section of the revised manuscript (page 22):\n“For participants with multiple imaging sessions, only the first scan was selected for inclusion, ensuring that each individual contributed a single observation to the analyses.”\nQ1.23 The effects of sex may vary nonlinearly across age, especially given the long developmental time span covered in this study. However, this nonlinearity does not appear to have been modeled.\nR: To test whether the effects of sex on gradient measures vary nonlinearly across the lifespan, in the revised manuscript we performed additional analyses at the global, system, and regional levels.\nSpecifically, we used the GAMLSS framework, which flexibly captures nonlinear trajectories across age. Two nested models were compared: a baseline model including main effects of age and sex, and a full model including an age-by-sex interaction term. Both models controlled for mean FD as a covariate, and site was modeled as a random effect. Significance of the interaction was assessed using the generalized likelihood ratio tests (GLRT), with FDR correction applied (q < 0.05).\nAt the global level, females showed significantly higher explanation ratio, gradient range, and variation compared with males across the lifespan (all p<0.05, FDR-corrected; Supplementary Fig. 35a). At the system level, the overall hierarchical organization was preserved, with the default-mode and frontoparietal systems consistently showed higher gradient ranges and standard deviations, while the sensorimotor system displayed lower values throughout life. Importantly, females exhibited significantly greater gradient ranges and variations than males (all p<0.001, FDR-corrected; Supplementary Fig. 35b). Moreover, peak age analyses revealed that in several functional systems, including visual, limbic, frontoparietal, and default mode networks, females demonstrated significantly later maturational peaks than males (all p<0.001, FDR-corrected; Supplementary Fig. 35c & Supplementary Table 8), suggesting delayed maturational timing in females. At the regional level, females exhibited higher gradient values in temporal, precuneus, medial prefrontal and somatomotor regions, whereas males showed higher values in occipital and lateral prefrontal regions (p<0.001, FDR-corrected; Supplementary Fig. 35d).\nThese findings have been incorporated into the Supplementary Methods (page 3) and Results sections (page 7):\nSupplementary Methods:\n“To test whether the effects of sex on gradient measures vary nonlinearly across the lifespan, we performed additional analyses at the global, system, and regional levels. Specifically, we used the GAMLSS framework, which flexibly captures nonlinear trajectories across age. Two nested models were compared: a baseline model including main effects of age and sex, and a full model including an age-by-sex interaction term. Both models, controlled for mean FD as a covariate, and site was modeled as a random effect. Significance of the interaction was assessed by comparing the model fits of the baseline and full models using the generalized likelihood ratio test (GLRT). All p-values were corrected for multiple comparisons using the FDR (q < 0.05). The direction of sex differences was determined based on the mean differences between the fitted sex-specific growth curves.”\nSupplementary Results:\n“To investigate whether sex differences in cortical gradient measures vary across age"}, {"round": 2, "reviewer_comments": {"reviewer_1": "General comments:\n- Upon inspecting the code for the group-level gradient analysis, “Gradient_Analysis.m” in the publicly available repo, there are a few things that seem worth flagging for clarification. First, I see a potential flaw with the way that age-specific group-FC is computed: they apply Gaussian age weight within each bin (good), but then appear to average by the number of subjects rather than by the sum of the weights. That effectively makes the scale of the group FC depend on the age dispersion and sample size in each bin, and it means subjects contribute very unevenly within a bin in a way that isn’t clearly justified. Second, the use of a single reference age to anchor the alignment, followed by reordering the gradients back into chronological age, assumes that component identities and eigenvalue ordering are stable across the lifespan, something that has been measured elsewhere to not hold.\n- In figure S3, which shows the clustering into lifespan stages based on individual gradients, the clusters in the MDS plot do not really look like clearly defined clusters but instead it looks like there are arbitrarily oriented cluster divisions. This is more suggestive of a “gradual transitions” narrative rather than a “distinct growth phases” narrative.\nFurther concerns are outlined below, using numbering consistent with the previous iteration.\nResponse 1.1: The current study however does not provide a “continuous” mapping of functional gradients because the temporal information is divided in the 26 clusters.\nResponse 1.3: It is unclear what is meant by “ its analysis strategy that integrates multiple structural properties prior to dimensionality reduction”. Gradients is a dimensionality reduction technique.\nResponse 1.4: If the results do not depend on binning, why is binning needed in the first place. The developmental phases are a somewhat artificial demarcation of the temporal developmental patterns since regional developmental patterns might not conform to these phases.\nResponse 1.5: Again, phases are very coarse groupings. Regional results might change with the age binning pattern.\nResponse 1.6: Although it is in general true that a larger sample size improves stability and reliability (provided that harmonization can be confirmed to be done correctly), it is unclear here what new neuroscientific findings exactly can be found here compared with studies of smaller sample sizes.\nResponse 1.7: The data, other than the fMRI data, was not processed using identical pipelines. There are additional acquisition and study design related variability that cannot be taken care of by using identical pipelines.\nResponse 1.9: The authors attribute the discrepancy with Ref 23 to methodological differences, suggesting that generalized additive mixed modeling (GAMM) 'smoothed over' early-stage patterns. However, a more parsimonious explanation lies in the dimensionality and rank-ordering of the eigenvectors.\nIt is well-documented that in the neonatal brain, the eigengap between the first and second gradients is often narrow, and the variance explained by the VS axis and the SA axis can be comparable. The 'emergence' of the S-A gradient at 3 months may simply represent a rank-swap between AP and SA gradients. Without showing both the first and second gradients and their corresponding eigenvalues before 3 months, it is very difficult to tell whether the observed AP gradient indeed evolves into the SA axis or if there are higher dimensional dynamics that underlie this transition. Furthermore, without explicitly considering temporal continuity in the modeling process, it is also difficult to understand how the SA axis gradually comes into existence from the AP axis.\nBy restricting the analysis to the primary gradient and aligning it to an adult template, the authors may be failing to detect an SA axis that is present but secondary in the perinatal period. To support the claim that the SA gradient is absent rather than merely sub-dominant, the authors must show Gradient 1, Gradient 2, and the corresponding scree plots (eigenvalues) for the <3 month age bins.\nResponse 1.13: Range can be estimated more robustly using robust statistics for the maximum and minimum, for example by using respectively the 99-th and 1-th percentiles. The examples on two groups with different clustering around the mean sound hypothetical and might not reflect the context of the current work. Range and variance are highly correlated, more so than being complementary.\nResponse 1.16: The response regarding the T1w/T2w ratio misses the biological constraint. While intensity normalization handles scanner variance, it does not account for the biophysics of neonatal tissue. In the first year of life, extremely high water content and low iron/myelin concentration render the T1w/T2w ratio a poor proxy for myeloarchitecture compared to adults. The lack of correlation between the structural and functional gradient growth axes (Q1.33) is likely a Type II error driven by this metric failure, rather than a true biological decoupling. The authors should discuss this biophysical limitation rather than attributing the null result solely to developmental asynchrony.\nResponse 1.19: Pipelines customized for specific time periods are not necessarily consistent across time between different cohorts. Cortical meshes can be defined differently and can cause inter-pipeline discrepancies particularly when the meshes are mapped to the sphere for downstream analysis.\nFor lifespan pipelines, see for example https://doi.org/10.1038/s41551-024-01283-7\nResponse 1.20: The authors acknowledge the confound that neonates/infants were scanned during sleep while older subjects were awake, but suggest that harmonized preprocessing and large sample sizes mitigate this. This reasoning is flawed. Preprocessing removes noise, but it cannot correct for the fundamental electrophysiological state difference between NREM sleep and wakefulness. The attenuation of the DMN and the alteration of long-range integration are biological signatures of sleep, not artifacts.\nIf the 'Establishment' phase (Phase II) correlates with the age at which subjects begin scanning while awake, the study may simply be tracking the 'Establishment of Wakefulness.' Large sample sizes do not average out a systematic bias present in ~100% of the cohort at a specific age. The authors must cite literature or provide analysis demonstrating that the S-A gradient topology is invariant to sleep state in infants, or explicitly re-frame Phase II as potentially state-dependent.\nResponse 1.21: The response does not fully address the concern. Regressing out site effects within age bins removes only mean shifts, not the protocol-specific covariance structure that differs substantially between dHCP, BCP, and HCP. Because the clustering is performed on gradients derived from datasets with different TRs, preprocessing pipelines, SNR profiles, and temporal degrees of freedom, the phase boundaries may still reflect acquisition differences rather than neurodevelopmental transitions.\nMoreover, the provided analyses remain qualitative. No quantitative cluster stability metrics (ARI/Rand/NMI, silhouette scores, perturbation tests) are reported, and cohort identity is not tested as a predictor of cluster membership. The equal frequency of 3- vs 4-cluster solutions further indicates instability. As written, the evidence is insufficient to conclude that the three phases are independent of dataset boundaries.\nAdditional quantitative validation is needed to rule out acquisition-driven structure.\nResponse 1.25: A Spearman correlation of 0.85 implies that roughly 28% of the variance is unexplained/ inconsistent between the two methods. In the context of mapping subtle developmental shifts in high dimensional topology is not trivial. The authors should not rely solely on the correlation coefficient between these two cases. They must demonstrate that the dimensionality and shape of the embeddings (preferably multiple dimensions) are equivalent between these methods. If the gradient-then-average result yields a more variable result across individuals, that is a relevant biological finding that should be reported and not smoothed over.\nResponse 1.30: The additional analyses described in the rebuttal help, but they do not yet demonstrate that the three identified “phases” are independent of age-binning choices or cohort boundaries.\nFirst, the individual-level clustering analysis is presented qualitatively, but no quantitative stability metrics (e.g., ARI/Rand/NMI between binned and unbinned cluster assignments) are provided. Without these, it is not possible to evaluate whether the same three clusters truly emerge consistently or whether visual similarity reflects smoothing and alignment decisions.\nSecond, showing that phases “partially overlap” with cohort inclusion does not address the core concern. The relevant question is whether cohort identity predicts cluster membership once age is controlled. This requires explicit testing (e.g., multinomial regression, χ²/MI tests, or adjusted associations). Simply harmonizing by modeling site as a random effect does not rule out cohort-specific variance structure influencing gradient topology.\nThird, the clustering procedure itself remains highly sensitive to K-choice, scaling, and the construction of the age-averaged templates. Without examining cluster stability across K, across random initialization, and across alternative harmonization strategies, the conclusion of three robust phases remains premature.\nIn short, while the added analyses are useful, the evidence provided does not yet demonstrate that the three clusters reflect genuine neurodevelopmental transitions rather than methodological or cohort-driven structure. Additional quantitative validation is required.\nResponse 1.31: Thank you for the clarification and for adding a limitation about the adult-derived cognitive reference map. However, I do not feel the concern is fully addressed, and the current text of the Results/Discussion still seems inconsistent with how you now describe this analysis.\nYou state that the goal was not to infer cognitive functions in younger cohorts, but the manuscript repeatedly interprets the findings as tracking the development and maturation of cognition itself (e.g., “formation and maturation of the cognitive spectrum from basic perception to abstract cognition,” “lifespan trajectory and maturation patterns of the cognitive spectrum along the S-A gradient”). Methodologically, however, what you do is: (i) construct an adult “cognitive spectrum” by ranking NeuroSynth term maps along an adult HCP S-A gradient, and then (ii) project those same adult term maps onto age-specific gradients that have been aligned to an adult S-A template. This can only quantify how much developing gradients converge on an adult functional/cognitive ontology; it cannot support inferences about the intrinsic “cognitive spectrum” in infancy and childhood, where both cognitive repertoire and task structure differ.\nRelatedly, I still see a degree of circularity. The adult S-A gradient is used both to define the cognitive reference spectrum (by ordering terms along that axis) and to align all age-specific gradients. Measures such as Dice similarity to the reference map are therefore, to a considerable extent, re-expressions of how closely each age group’s gradients have been forced into an adult S-A framework, rather than an independent test of cognitive organization across development.\nIn my view, the analysis does not justify strong claims about “lifespan maturation of the cognitive spectrum” or the “formation” of a perception-to-abstract-cognition spectrum across the entire age range. At most, it shows that brain regions which, in adults, are preferentially associated with certain cognitive terms become progressively more organized along an adult S-A axis.\nResponse 1.35: The authors’ clarification—that they excluded gradient values after computation rather than re-embedding network-ablated matrices—confirms that this analysis is circular. Since the Somatomotor and Default Mode networks are mathematically defined as the anchors (the extremes) of the S-A gradient, \"knocking them out\" and observing a shift in the mean is a tautological arithmetic certainty, not a biological insight into system \"contributions\" or \"interplay.\" The results describe the definition of the gradient, not its developmental mechanism. We strongly suggest the authors remove the causal language (e.g., \"somatomotor system contributed most significantly\") and frame this strictly as a descriptive property of the distribution.\nResponse 1.39: Thank you for adding the cosine similarity analysis. However, I do not feel this fully addresses my question about alternative template construction, nor the interpretational issues around template choice. First, you explicitly acknowledge that you did not implement alternative template schemes, yet you present the individual→template cosine similarity as evidence that your current scheme “effectively captures age-specific functional topographies.” Because individual gradients are Procrustes-aligned to the age-specific templates before similarity is computed, this metric is not an independent validation of the template but is partly tautological: it will naturally increase with age as SNR and motion improve and as the alignment itself enforces correspondence. Second, I see a tension between this new analysis and your own canonical S-A axis results. In the main text, similarity between group-level gradients and the canonical S-A axis from Sydnor et al. increases in Phase I, decreases in Phase II, and then stabilizes at a high level in Phase III (Supplementary Fig. 2). Yet the individual→age-template similarity in Supplementary Fig. 25 increases through childhood and adolescence and peaks in early adulthood. This divergence suggests that the age-binned template plus multi-step Procrustes procedure may be absorbing the substantial reconfiguration of gradients in infancy–early childhood into the templates themselves, smoothing over the mismatch to an independent canonical SA axis. In my view, this new analysis does not substitute for testing alternative template construction strategies (e.g., different reference ages, pooled templates, or reduced reliance on the adult HCP anchor), nor does it resolve the concern that the developmental story might depend sensitively on the current template choice.\nResponse 1.40: I appreciate the clarification. However, my concern was less about the number of reported metrics and more about how fully the analysis leverages the rich spatial information in the gradients. Most of the results you describe are based on 1D summary measures (e.g., global range, network-averaged values, bin widths, centile trajectories) rather than analyses of spatial patterning per se (e.g., within-network reconfigurations, local clusters, or multivariate spatial relationships among gradients). In my view, the current work still primarily treats gradients as scalar scores to be summarized across regions and fitted with age-curves, which is different from a deeper topographic analysis of gradient organization. At minimum, I would encourage tempering statements that the S–A gradient provides a full “developmental blueprint,” to reflect this limitation.\nResponse 1.42: Is there evidence that this is true in infant brains too? If not, some comparison between high resolution and lower resolution gradients should be done during early life to verify.\nResponse 1.43: In your reply, you state that analyzing eigenspectra is beyond the scope of the current work. I find this problematic. The manuscript repeatedly claims to characterize hierarchical organization across the lifespan, yet there is no quantitative analysis of how the hierarchical series of embedding axes changes with age. Reporting that the first gradient explains ~20% of variance at each age does not address whether the relative importance of the leading axes (e.g., λ₁, λ₂, λ₃, …) shifts over development, nor whether a single axis is adequate to summarize multi-scale reorganization.\nWhile I appreciate that this study focuses on the S–A axis, the dynamics of that axis—particularly in infancy—are tightly coupled to the behavior of the remaining spectrum. Large changes in eigengaps (e.g., λ₁–λ₂, λ₂–λ₃) would directly bear on whether the S–A axis remains a uniquely dominant dimension at all stages, or whether other axes transiently become equally or more informative. At minimum, plotting age-related trajectories of eigengaps at the subject or group level would provide an important sanity check on the assumption that one axis is sufficient.\nGiven that the eigenvalues underlying the gradients must already be computed in your existing pipeline, this does not seem like an outlandish or technically burdensome request. If such analyses are not performed, I would strongly encourage you to temper claims about “hierarchical organization” being captured by a single S–A gradient across the lifespan, and to state clearly that the multiscale structure of the eigenspectrum remains uncharacterized.\nResponse 1.47: The appeal to convergence with prior infant studies is reassuring but does not provide direct evidence that pipeline disparities are not contributing to the inflated gradient range in your data, since those studies also rely on infant-specific acquisition and preprocessing choices. Ideally, one would see at least a limited sensitivity analysis (e.g., a subset of infant and/or adult data reprocessed with a common pipeline, or a demonstration that the infant gradient range is robust across reasonable denoising choices or after controlling for motion and temporal SNR). If such analyses are not feasible within the scope of the current revision, I would encourage you to soften the claim that preprocessing differences cannot systematically bias the gradients and to more explicitly acknowledge this as a limitation when interpreting the inflated infant extremes.\nResponse 1.51: I appreciate the clarification about the use of NbClust and the suite of clustering indices. However, my original concern was not only how K was chosen, but also about (i) seeing at least some of the underlying metrics (e.g., silhouette widths or gap statistic across candidate K) and (ii) the stability of the “lifespan phases” across different age binning schemes. The current response does not address the bin-width stability, and the NbClust results are only summarized qualitatively.\nResponse 1.55: Please include a statement in the manuscript to reflect that this study is performed in compliance with the NDA data use agreement.\n(Remarks on code availability)\nThe code suggests a potential flaw in the implementation; see the “General Comments” section of my report for details.\nReviewer #2\n(Remarks to the Author)\nI commend the authors on their efforts to address all reviewer comments through substantial and insightful additional analyses. The inclusion of these new analyses has provided a much more thorough description of the data and has added confidence that the overall patterns of results are not substantively impacted by different analytic choices. This leaves the authors with the difficult task of synthesizing and interpreting the expanded corpus of results into a cohesive story. I feel this part is somewhat incomplete, leaving important food for thought on the table. Sharpening the story in light of all the new analyses could help to improve the impact of this huge body of work. I mention a few of these areas below.\nSex effects observed in the new analyses are surprising and bear further scrutiny. There should be some discussion of these findings. There is currently no reference to this analysis or finding in the main text of the manuscript. Rather than burying these results in the supplement, the paper would be strengthened by referencing them directly in the text and offering an interpretation or discussion of their potential meaning in the discussion section\nThe flip in gradient dominance during phase 2 (e.g. Supplementary Figure 1) and the corresponding dip in correlation of the primary gradient and the canonical SA axis (Supplementary Figure 2) is quite interesting. It seems that something unique and important going on during phase II that may be developmentally reorganizing the A-P gradient into what will eventually become the S-A gradient by means of a dominance of gradient 2. The authors discuss possible reasons for the S-A gradient emergence in general in the discussion, but some additional interpretation of what is going on during this unique period, and what it means for “principal” vs secondary gradients during development—and these sorts of PCA approaches in general—would be useful to explore more.\nSome of the spatial correlation plots appear to actually have non-linear relationships, yet they are always analyzed and overlaid with a linear fit (e.g Fig 3g, Supplementary Fig 9a). I am left wondering if appreciating this non-linearity would be meaningful for understanding how all the pieces fit together.\n(Remarks on code availability)\nThe code appears to be documented and complete; however, I did not try to execute it so cannot comment on whether it can be run successfully.\nReviewer #3\n(Remarks to the Author)\nMy comments have been addressed. I have no further comments.\n(Remarks on code availability)"}, "author_response": ""}, {"round": 3, "reviewer_comments": {"reviewer_1": "I would like to commend the authors on a large and– for the most part–extremely thorough revision. I appreciate their patience and diligence in responding to my many highly technical and specific comments. Their doing so has dramatically improved the manuscript and its scientific defensibility.\nFigure S12 is quite interesting, however it is difficult to get a clear sense for the relative rank-ordering of each resting-state network when the gradient value axis is scaled and positioned differently for each network. I appreciate that for this figure, the scatterplots of individual subject data make the axis ranges vary significantly, and seeing the variance/ distribution is interesting and informative. Maybe, in addition to showing these plots, the network centroid trajectories could also be shown on the same axis, as the relationships between trajectory values and timings are potentially more interesting than each trajectory in isolation. Also on this note, in section (ii) of the results section “Normative growth of the S-A gradient axis across lifespan”, you note “Specifically, the visual, ventral attention, and dorsal attention systems peaking in adolescence”. I think your meaning here is closer to “reach their extreme value” rather than “peak” as these networks reach peaks on the sensory side of the axis.\n1.6: The concern raised is not fully addressed for several reasons. First, a smaller sample size does not automatically imply that the data are fragmentary or restricted to a limited age range; meaningful and generalizable insights can still be derived from well-curated datasets. Second, prior studies have consistently demonstrated that gradient patterns exhibit systematic changes across multiple scales, including global, system-level, and regional levels, suggesting that these patterns are robust and detectable even in smaller cohorts. Third, relationships between structural and functional gradients have been successfully investigated in studies with relatively modest sample sizes that span the entire human lifespan, indicating that comprehensive lifespan trends can be captured without extremely large cohorts. Fourth, attempts have already been made to characterize lifespan trajectories of gradients in relation to cognitive functions using smaller datasets, further supporting the feasibility of lifespan analyses even when the sample size is limited. Collectively, these points suggest that smaller sample sizes, when appropriately designed, can still provide meaningful insights into gradient-based brain organization across development and aging.\n1.9: Thank you for displaying these results. They are reassuring and quite interesting.\n1.25: Thank you for carrying out these analyses. It remains slightly unclear how the 0.77 value across the first two gradient dimensions is obtained. Is this the mean similarity across the two dimensions? If so, what does the breakdown by axis? Although the 0.77 value indicates broad spatial similarity, it may be slightly misleading to claim that the approaches are fully comparable. This is not a fatal flaw and there is no strong argument that one way is more correct than the other, however it bears noting that the approaches do yield substantively different gradients when it comes to topography.\n1.30: Seeing the clustering metrics is indeed informative.\n1.31: Thank you for adding these caveats, I believe they make the claims more defensible.\n1.39: This is very reassuring, thank you for carrying out this analysis. This appears to be a significant result in and of itself– that template construction method has very little effect on topography.\n1.43: This is highly informative to see. To clarify, the values depicted in Fig S2 are the explanation ratios associated with the group-level template-aligned gradients computed from the bin-average connectivity matrices, is that correct?", "reviewer_2": "I appreciate the additional analyses added by the authors to further explore the developmental trade-offs between the first and second principal gradients. It appears to highlight an interesting tension between and gradient based perspective on brain organizational structure and a network based perspective. I also appreciate the additional analysis in response to the very thoughtful comments of reviewer 1. I have no further concerns to be addressed prior to publication.", "reviewer_3": "My comments have been addressed. I have no further comments."}, "author_response": "Comments by the reviewers and our point-by-point responses:\nWe sincerely thank the editor and all the reviewers for their comments on our manuscript, titled\n“Spatiotemporal dynamics of the human cortical functional hierarchy across the lifespan” (ID:\nNCOMMS-25-33712-T). We found the comments to be constructive and highly valuable for improving\nour manuscript, which we revised based on the comments received. Our point-by-point responses to the\nreviewers’ comments are presented below.\nReviewer #1 (Remarks to the Author):\nIn this paper, the author maps the development of the sensory-association (SA) gradients across the human\nlifespan and explores how these gradients relate to structural brain measures. While the characterization of\nthe SA gradient is a valuable contribution to understanding brain organization and maturation, I have\nmajor concerns regarding several aspects of the study. These include questions about the novelty of the\nwork in light of existing literature, the appropriateness and rigor of the methodological approaches, the\nclarity and robustness of the reported results, and the strength and validity of the conclusions drawn from\nthe findings.\nR: We thank this reviewer for all positive comments and insightful comments/suggestions.\nComments:\nQ1.1 Introduction - “First, existing studies focus mainly on discrete developmental stages (e.g., the\nthird trimester 19 , neonatal period, adolescence, and adulthood), without provided a continuous\nspatiotemporal perspective from prenatal development through late adulthood.” This is not true and\nmisleading in light of Ref 23.\nR: We have rephrased the relevant descriptions and hopefully it is clearer now (page 4):\n“First, most existing studies of the functional hierarchy have focused on discrete developmental stages\n(e.g., the third trimester 19, neonatal period 18, adolescence 20,21 and adulthood 3,22), resulting in a\nfragmented perspective that obscures fine-grained transitions and nonlinear dynamics of cortical\ngradients across the full age continuum. Moreover, prior work has often been limited by small (e.g., N <\n100) 18,19 or modest (N = hundreds to thousands) 16,20–22 sample sizes, along with a predominant reliance\non linear modelling to characterize gradient development 19,20,22. Although a recent study 23 established a\ncontinuous, nonlinear mapping of functional gradients from birth to 100 years of age, the inferred\ndevelopmental patterns remain constrained by a modest sample size (N = 3,556) and particularly sparse\ndata coverage during childhood.”\nQ1.2 Introduction - “Furthermore, methodological limitations, such as small (less than one hundred)\nor modest (hundreds to thousands) sample sizes and an overreliance on simplistic linear modelling\napproaches, fundamentally constrain our ability to capture nonlinear neurodevelopmental dynamics.” This\nis also not true. The methods used for data processing and analysis (FreeSurfer, GAMLSS, etc.) in the\ncurrent work are mostly publicly available and are quite standard. There are no methodological hurdles\npreventing other groups from carrying out similar work.\nR: We would like to clarify this point raised by this reviewer. We agree that both standard preprocessing\npipelines (e.g., FreeSurfer) and advanced nonlinear modeling frameworks (e.g., GAMLSS) are well-\nestablished within the community. Our aim was not to imply that these tools are non-standard, but to\nhighlight that many prior studies have yet to fully integrate large-scale datasets with nonlinear modeling to\ncharacterize developmental trajectories in the cortical hierarchy. We have revised the sentence accordingly\n(page 4):\n“First, most existing studies of the functional hierarchy have focused on discrete developmental stages\n(e.g., the third trimester 19, neonatal period 18, adolescence 20,21 and adulthood 3,22), resulting in a\n1\nfragmented perspective that obscures fine-grained transitions and nonlinear dynamics of cortical\ngradients across the full age continuum. Moreover, prior work has often been limited by small (e.g., N <\n100) 18,19 or modest (N = hundreds to thousands) 16,20–22 sample sizes, along with a predominant reliance\non linear modelling to characterize gradient development 19,20,22.”\nQ1.3 Introduction - “Second, current conceptual frameworks predominantly emphasize functional\ngradient remodelling during development while largely overlooking its structural foundations. In mature\nbrains, the S-A axis aligns spatially with structural features indexed by cortical thickness, intracortical\nmyelination, and intrinsic geometry.” This is not true. Association between primary gradients and\nstructural gradients has already been investigated for example in reference 23.\nR: We thank the reviewer for raising this point. We have revised the text to clarify this (page 4-5):\n“Second, current research frameworks lack the capacity to delineate how functional hierarchies coevolve\nwith specific structural features, such as cortical thickness 2,6,24, intracortical myelination 2,8,9, and\nintrinsic geometry 3,25, across the lifespan. While the S-A axis aligns spatially with these structural\nattributes in mature brains, the developmental synchrony between functional and structural maturation\nremains largely unclear, particularly given their distinct maturational timelines in early life 2. A previous\nstudy 23 reported coupling between lifespan changes in functional gradients and a composite structural\ngradient derived from multiple features; however, its analysis strategy that integrates multiple structural\nproperties prior to dimensionality reduction does not directly assess the contribution of specific structural\nproperties to the emergence of functional gradients.”\nQ1.4 Why 26 age groups? This creates artificial temporal demarcations between different life periods\nthat do not necessarily reflect true biological development. This is particularly problematic when different\nbrain regions have distinct development timelines that do not coincide with these artificial temporal\npartitions. It is uncertain if age grouping will result in artificial disconunity in the analysis, not just at the\nglobal whole-brain level but also local vertex level.\nR: We thank the reviewer for raising this point. Currently, there is no universally accepted standard for\ndefining age groups. Our selection of 26 age groups was informed by developmental principles indicating\nthat brain maturation occurs rapidly during infancy and childhood but slows progressively through\nadolescence and adulthood. Therefore, we applied finer age partitions at younger stages and broader\nintervals at later ages to better capture these nonlinear developmental dynamics.\nTo further address the concern regarding potential artificial discontinuities, we also conducted\ncomplementary analyses at the individual level without imposing predefined age bins. Specifically,\nfunctional gradients were computed for each participant and then clustered across all individuals, using\nchronological age as a quasi-continuous variable. This individual-level approach confirmed the main\nfindings of three phases of functional gradient development: (Phase I) early anterior-posterior\ndifferentiation during infancy, (Phase II) emergence of the S-A gradient in early childhood, and (Phase III)\nits gradual expansion, refinement, and stabilization throughout adolescence and adulthood. The\nconvergence between the group-based and individual-level analyses indicate that our conclusions are not\naffected by the chosen age-grouping scheme.\nThe Supplementary analyses were added in the revised manuscript (page 7):\nResults: “The three developmental phases were also identified through both individual-level gradient\nclustering without predefined age bins and within-age-group-harmonized analysis (Supplementary Fig. 3,\n4). This indicates that the observed patterns were robust and not driven by specific age-grouping schemes\nor site related confounds.”\nSupplementary Results: “To validate the robustness of our findings, we first verified that developmental\nphases identified in the group-based analyses were consistently reproduced in individual-level clustering\n2\nwithout predefined age bins, confirming that the observed patterns were not driven by the specific age-\ngrouping scheme (Supplementary Fig. 3).”\nSupplementary Fig. 3 | Distinct growth stages of the sensorimotor-association gradient across the lifespan identified\nfrom individual gradients. a, K-means clustering analysis based on the pairwise correlations between the individual\ngradient maps, resulting in three clusters. b, Three distinct phases of the S-A gradient are visualized via multidimensional\nscaling (MDS). “Dimension 1” and “dimension 2” are defined as the first and second axes obtained through MDS of the\nfunctional distance matrix.\nQ1.5 How does the choice of 26 age groups influence the clustering results? The patterns shown in\nFigure 1b may vary not only with the number of age groups but also with how these groups are distributed\nover time. Moreover, the age grouping adopted from the authors’ previous publication may not be\nappropriate here, as this study focuses on gradients rather than simple features of functional connectivity.\nThe authors should demonstrate that this age grouping does not introduce bias by comparing it against\nclustering results obtained from the entire dataset without imposing predefined age groups.\nR: To evaluate whether the clustering results were influenced by the choice of age-group definitions, we\nconducted an additional analysis at the individual level across all participants without imposing any\npredefined age bins. Functional gradients were derived separately for each participant, and clustering was\nthen performed on the complete dataset using chronological age as a continuous variable.\nThis individual-level analysis reproduced the same three major developmental phases identified in the\nage-binned analysis: (Phase I) early anterior-posterior differentiation during infancy, (Phase II) emergence\nof the sensorimotor-association (S-A) gradient in early childhood, and (Phase III) progressive expansion\nand stabilization of this gradient through adolescence and adulthood. The high consistency between the\ngroup-based and individual-level results demonstrates that the choice of 26 age groups does not bias the\nclustering outcomes. These findings confirm the robustness of our developmental trajectory\ncharacterization and validate that the observed patterns are not artifacts of grouping strategy.\nQ1.6 How is the difference in sample size across age groups dealt with? The authors did a “balanced\nresampling analysis” by equal resampling across age bins. However, it is unclear why the 16 “age bins”\nare different from the 26 “age groups”. Also, does the fact that resampling uniformly with a smaller\nnumber of samples per age bin still replicates the main results suggest that the large sample size used in\nthe study may not be necessary?\nR: We complemented the age-binned analysis with individual-level clustering that does not depend on\npredefined bins. This analysis again confirmed the same three developmental phases, demonstrating that\nthe main findings are robust to differences in group size and binning scheme.\nThe use of 16 “age bins” in the balanced resampling analysis reflects a methodological adjustment aimed\nat ensuring an adequate and comparable number of participants in each bin for statistical resampling,\nrather than a conceptual change from the 26 developmental age groups used in the main analysis.\nRegarding the necessity of a large sample size, as detailed in our response to Q1.8, larger samples\nsubstantially improve the stability and reliability of lifespan growth curve estimates, particularly at the\nregional and vertex levels, by reducing variance and enhancing the precision of developmental trajectory\ncharacterization. Thus, while smaller, balanced subsamples supports the validity of our conclusion, the\nfull sample remains essential for achieving optimal reliability and generalizability of brain gradient\ntrajectories.\nQ1.7 How are differences in data variances between datasets and across time dealt with?\nThere can be nonlinear variance changes across age groups and within age groups and also across the\ndifferent datasets and within the datasets.\nR: We took several steps to minimize potential biases introduced by variance differences across datasets\nand age groups in the original analyses.\nFirst, at the preprocessing stage, all datasets were processed using identical pipelines: projection to\nfsaverage_LR32k, nuisance regression with 24 motion parameters plus WM, CSF, and global signal,\nlinear detrending, temporal bandpass filtering (0.01-0.08 Hz), scrubbing (FD > 0.5 mm with one\npreceding and two subsequent volumes interpolated), surface smoothing (FWHM = 6 mm), and\nresampling to fsaverage4 space. These harmonized pipelines ensured comparable signal properties across\ndatasets.\nSecond, at the modelling stage, we explicitly accounted for age- and site-related heteroscedasticity using\nthe GAMLSS framework with a flexible skew exponential power type 4 (SEP4) distribution. Specifically,\nwe modeled the variance (scale parameter, σ) as a smooth function of age and sex (sigma.fo=~bs(Age) +\nSex), which accommodates nonlinear variance changes across the lifespan. Site was further included as a\nrandom effect to capture both mean and variance differences across datasets.\nThird, model diagnostics confirmed the adequacy of our approach. We examined normalized quantile\nresiduals, residuals-fitted plots, kernel density estimates, and normal Q-Q plots (Supplementary Fig. 28),\nall of which indicated well-behaved residuals without systematic trends. Detrended Owen’s plots further\nshowed that residuals remained within nonparametric confidence intervals (Supplementary Fig. 29).\nTogether, these diagnostics suggest that the GAMLSS effectively captured heteroscedasticity.\nFinally, sensitivity analyses using a leave-one-site-out strategy demonstrated that our main findings were\nhighly consistent across datasets, supporting the robustness of our results against potential variance\ndifferences.\nQ1.8 A large sample size like the one used in this paper is not always accessible to all researchers.\nPlease show empirical evidence that such a large sample size gives conclusions that cannot be drawn\nusing a smaller sample size.\nR: We thank the reviewer for raising this important point. To empirically evaluate the benefits of large\nsample sizes, we performed a resampling-based analysis while controlling for both age and sex\ndistributions. Specifically, the lifespan range (32 PMW to 80 years) was divided into 10 equal intervals,\nand subsamples were drawn while preserving the original proportions of age and sex within each interval.\nFor global- and system-level analyses, sample sizes were varied from 200 to 33,000 participants, with\nfiner increments at smaller sizes. For regional-level analyses, 13 predefined sampling conditions were\nused. For each sample size, we performed 50 independent resampling iterations, refitted the lifespan\ngrowth curves, and computed the mean absolute error (MAE) and root mean squared error (RMSE)\nrelative to the reference curve estimated from the full dataset. The minimum sample size was defined as\nthe point at which the observed MAE/RMSE decreased to within 1.5 times the error estimated at the\nmaximum sample size for regional-level gradient values.\nThe results showed that both MAE and RMSE decreased as sample size increased. At the global and\nsystem levels, error dropped sharply when sample size grew from a few hundred to ~5,000 participants,\nafter which gains plateaued. At the regional level, required sample sizes varied across cortical vertices.\nA sample size of at least ~6,000 participants was needed for more than 50% of cortical vertices to achieve\nstable error levels, while ~25% of vertices required even larger samples.\nIn sum, these analyses demonstrate that larger sample sizes substantially enhance the stability and\nreliability of lifespan growth curve estimates, particularly at the regional scale. These results provide\nempirical guidance on the minimum sample sizes required for robust conclusions in future lifespan\nneuroimaging studies. We have added relevant results in the revised manuscript (page 16):\n“To evaluate the influence of sample size on our findings, we conducted a resampling-based sensitivity\nanalysis, varying sample sizes from 200 to 33,000 participants. We observed that both mean absolute\nerror (MAE) and root mean squared error (RMSE) decreased with increasing sample size (Supplementary\nFig. 17). Notably, errors at the global and system levels dropped sharply until about 5,000 participants,\nafter which gains plateaued. Furthermore, at the regional level, a sample size of ~6,000 participants is\nnecessary to achieve stable error estimates for more than half of the cortical vertices, while a quarter of\nvertices required even larger samples. This work confirms that larger samples are crucial for reliable\ngradient growth curves, especially at the regional scale, and provide valuable guidance for future lifespan\nneuroimaging studies.”\nSupplementary Fig. 17 | Relationship between sample size and fitting error for gradient measures. a, Relationship\nbetween sample size and fitting error for global-level gradient measures (MAE shown in blue, RMSE in red). b,\nRelationship between sample size and fitting error for system-level gradient measures (MAE in the left panel, RMSE in the\nright panel). c, Relationship between sample size and fitting error for regional-level gradient values. The minimum sample\nsize at which the MAE (left) /RMSE (right) reached a stable plateau for each vertex is mapped onto the cortical surface.\nto 33,000 participants, with 50 resampling iterations per size, and age/sex distributions were preserved across 10 lifespan\nintervals (32 PMW to 80 years).\nQ1.9 Anterior-posterior gradient during early infancy – This is contradictory to previous observation\nthat the SA gradient already exists at birth. The SA gradient does not start to occur at 3 months but before\nthat.\nR: We thank the reviewer for this valuable comment and the opportunity to clarify. In our data, we\nobserved a prominent anterior-posterior (A-P) gradient during the prenatal and neonatal periods, which\nprogressively transitions into a more differentiated S-A gradient after approximately 3 months of age. This\nobservation of an A-P gradient in early infancy is consistent with prior studies. For instance, Xia et al.\n(2024) reported the emergence of an anterior-posterior gradient as early as 34.6 postmenstrual weeks 1.\nLariviere et al. (2020) described anterior-posterior variations in term-born neonates, interpreting these as\nan immature precursor of the differentiation between unimodal and transmodal areas 2.\nThe apparent inconsistency with Ref23 likely reflects methodological differences. Specifically, Ref23\nemployed a generalized additive mixed model (GAMM) to fit gradient development across the entire age\nrange. Given the brevity of the early developmental window and the relatively small sample size available,\nsuch modeling may smooth or shape the less differentiated early-stage patterns into a more continuous S-\nA trajectory. In contrast, our approach used group-level analyses within discrete age bins, enabling us to\ndirectly characterize gradient architecture at each developmental stage without imposing model-based\ncontinuity across age. We have revised the text by adding a small discussion (page 18):\n“Notably, we found that the functional topography at birth is not characterized by the S-A gradient but is\ninstead dominated by an anterior-posterior differentiation. This finding is aligned with previous studies of\npreterm, term-born, and neonatal cohorts 18,19, suggesting that this pattern represents a developmental\nprecursor to the mature S-A gradient. In contrast, one recent study 23 suggests that the S-A gradient is\nalready present at birth. The discrepancy may partly arise from methodological differences: the previous\nstudy 23 applied lifespan generalized additive mixed modeling, which could smooth over less differentiated\nearly-stage patterns. Our stage-specific analyses, by comparison, capture gradient organization directly\nat each period, offering a more fine-grained characterization of early functional topography.”\nQ1.10 Grouping the broad age range from 5 years to 80 years as a single phase (Phase III) is too coarse\nand does not reflect the fact that growth inflection points happen quite often during adolescence.\nR: We would like to further clarify this point. Our analysis for defining Phase III was not to suggest that\nage span from 5 to 80 years constitutes a homogenous developmental stage. Rather, Phase III was\nconceptualized as a broader period of gradient stabilization and refinement-a stage in which the S-A\ngradient, once established, undergoes progressive fine-tuned throughout childhood, adolescence, and\nadulthood, and then continues to evolve into late life. We have clarified this point in the revised\nmanuscript (page 19):\n“our definition of Phase III (5-80 years) reflects the prolonged refinement and stabilization of the S-A\ngradient from childhood, adolescence to adulthood, rather than developmental uniformity through this\nperiod.”\nQ1.11 Datasets like the developing HCP (dHCP) use gestational age. Postnatal datasets such as the\nBCP use chronological age. How is the discrepancy between gestational age and chronological age\nresolved?\nR: To ensure consistency across datasets, we converted gestational age in the dHCP dataset into\nchronological age by treating 40 gestational weeks as the time of birth (i.e., 0 years). Chronological age\nwas then calculated as:\n(𝐺𝑒𝑠𝑡𝑎𝑡𝑖𝑜𝑛𝑎𝑙 𝑎𝑔𝑒 𝑎𝑡 𝑠𝑐𝑎𝑛−40)×7\n𝐶ℎ𝑟𝑜𝑛𝑜𝑙𝑜𝑔𝑖𝑐𝑎𝑙 𝑎𝑔𝑒 (𝑦𝑒𝑎𝑟𝑠)=\npostnatal datasets within a common developmental timeline. Importantly, this step was applied solely for\nmodeling smooth lifespan trajectories in the GAMLSS framework and does not influence other analyses\nor alter the overall conclusions of the study. In this revised manuscript, we made a further clarification.\n“To align the dHCP dataset with other cohorts, we converted gestational age to chronological age. The\nconversion was based on the assumption that 40 gestational weeks corresponds to birth (0 years of\nchronological age). The specific formula applied was: Chronological age (years) = (Gestational age at\nscan – 40) × 7 / 365. This procedure ensured all ages were referenced to a consistent postnatal\ntimeline.”\nQ1.12 “Collectively, these results revealed previously unreported spatiotemporal patterns and critical\nneurodevelopmental milestones of the functional hierarchy, as characterized by the initiation,\nestablishment, and expansion/stability of the S-A axis of cortical functional organization.” Again, this is\nuntrue in light of ref 23.\nR: We appreciate the reviewer’s comment. The related statements have been removed from the revised\nversion.\nQ1.13 Range and variance are measuring almost the same thing, i.e., the spread of the gradient values,\nand are therefore redundant. Please pick one to report.\nR: While both range and standard deviation quantify aspects of distribution spread, they emphasize\ndistinct properties. The range captures the full extent of gradient differentiation (maximum-minimum),\nhighlighting the span of the axis and sensitivity to extreme values. In contrast, the standard deviation\nreflects the average dispersion of all gradient values (not just the maximum and minimum) around the\nmean, providing a more robust, outlier-resistant measure of variability. For example, consider two groups\nof numbers with the same range: one group is tightly clustered around the mean, and the other is spread\nout more evenly between the minimum and maximum values. Both groups have the same range, but the\nfirst group will have a smaller standard deviation, indicating less variability around the mean, whereas the\nsecond group will have a larger standard deviation, reflecting greater dispersion. Because these two\nindices offer complementary perspectives‒one emphasizing the maximal differentiation and the other\noverall dispersion‒we report both to provide a fuller characterization of developmental and topographical\nfeatures of the gradient.\nQ1.14 In Figure 2c and throughout the paper, it is unclear what “gradient score” is referring to.\nR: To improve clarity, we have revised the terminology throughout the manuscript and figure legends,\nreplacing “gradient score” with “gradient value”. This change more accurately reflects the numerical\nvalues derived from the gradient decomposition.\nQ1.15 The results shown in Figure 2e seem to indicate that the SA gradient changes very little after 20\nyears. This contradicts Figs 2a and 2b and other results shown in the paper.\nR: We would like to make a clarification about this point. First, the apparent stability of the S-A gradient\nafter age 20 in Fig. 2e was partly due to the original colormap design rather than the underlying results.\nSpecifically, the broad white center of the colormap compressed the visual contrast around zero values,\nmaking age-related changes less discernible. In the revised figure (see updated Fig. 2e), we applied a\ncolormap with more balanced contrast to better capture subtle variations across adulthood.\nSecond, our analyses do indicate that the S-A gradient continues to change after early adulthood, although\nat a slower rate compared with the rapid alterations observed in early childhood. These adult changes are\nconsistent with the global- and system-level results shown in Fig 2a and 2b, which demonstrate relative\nstabilization of gradient measures from the third decade onward to early aging, while still exhibiting\ngradual refinements at the regional levels.\nFig. 2 | Lifespan growth of the sensorimotor-association gradient. a, Lifespan growth curves (top panel) and growth\nrate curves (bottom panel) for global gradient measures. The solid line represents the 50th centile (median), while the\ndotted lines represent the 5th, 25th, 75th, and 95th centiles. The growth rate was estimated by the first derivative of the\nmedian growth curve, with a 95% confidence interval (shaded in gray) obtained from 1,000 bootstrap resamplings. The\nhorizontal dotted line (y=0) marks the transition from growth to decline, with the peak age indicated by a purple dot. b,\nLifespan growth curves of gradient range and standard deviation for each functional system defined according to the age-\nspecific Yeo 7-network atlas 30,34,35. c, Lifespan growth of the gradient value for each functional system. The peak age is\nmarked by a black dot. d, Contributions to lifespan growth of the S-A gradient from each functional system, estimated via\nknockout-based simulation analysis (left panel). A repeated-measures ANOVA tested variations in system contributions\nacross developmental phases (right panel). ***, p<0.001. e, Lifespan growth rate maps of the S-A gradient at representative\nages. f, Lifespan growth axis of the S-A gradient. The left panel shows regional gradient values across age. PCA identified\nthe first component as the dominant lifespan growth axis, explaining 65.4% of the variance (middle panel). The growth\naxis was divided into 20 decile bins, and mean gradient values were computed for all vertices within each bin (right panel).\ng, Spatial correlation between the lifespan growth axis and the evolutionary hierarchy. Evolutionary hierarchy was\nquantified using macaque-to-human cortical expansion 32. Statistical significance was assessed via 1,000 non-parametric\nspin permutation tests 36. CI, confidence interval; VIS, visual; SM, somatomotor; DA, dorsal attention; VA, ventral\n8\nattention; LIM, limbic; FP, frontoparietal; DM, default-mode; w, week; m, month; yr, year.\nQ1.16 T1-weighted and T2-weighted images are not quantitative and can have arbitrary signal gains\nassociated with different scanners, sequences, and choices of acquisition parameters. It is unclear how this\nis dealt with in this study, which involves datasets from quite a number of different sources.\nR: We appreciate the reviewer’s insightful comment. It is indeed true that T1-weighted and T2-weighted\nimages are not inherently quantitative, and their signal intensities can vary across scanners, sequences, and\nacquisition protocols.\nNevertheless, the T1w/T2w ratio provides a reliable proxy for relative cortical myelin content because\nmyelin consistently shortens T1 relaxation times (increasing T1w signal) and T2 relaxation times\n(decreasing T2w signal), producing a systematic contrast that is amplified when taking their ratio 6,7. This\napproach has been widely adopted in developmental, clinical and comparative studies 6,8–11. Importantly,\nseveral studies have demonstrated that the T1w/T2w ratio produces robust and anatomically meaningful\nspatial patterns that remain consistent across different scanners and acquisition settings. For example,\nGlasser et al. reported that group-average datasets acquired on different scanners with different sequences\nat independent sites yielded qualitatively very similar results 6,7.\nIn the present study, we took several steps to further mitigate site-related variability. First, during\npreprocessing we applied biases and intensity normalization, which reduce transmit field biases and\nminimize intensity inhomogeneities. These steps are critical for improving cross-site comparability of\nT1w/T2w-derived myelin maps. Second, in our statistical modeling, we explicitly included site as a\nrandom effect in the GAMLSS framework, thereby adjusting for potential residual site-specific\ndifferences when estimating growth trajectories.\nTogether, these strategies substantially reduce concerns regarding inter-scanner variability and support the\nrobustness of our T1w/T2w-based myelin estimates across multiple datasets. In the revised manuscript, we\nhave added relevant description and hopefully it is clearer now (page 29).\n“To mitigate inter-scanner variability and enhance cross-site comparability, all T1w and T2w images\nunderwent bias correction and intensity normalization prior to the computation of the T1w/T2w\nratio.……. we fitted the growth curves of the structural attributes for each vertex via GAMLSS with age\nand sex as fixed effects, and site as a random effect.”\nQ1.17 Were the inter-vertex geometric distances normalized with respect to the “total” inter-vertex\ndistance. Since the brain is changing in size across the lifespan, it is unclear how the brain size will affect\nthe analysis on geometric distance.\nR: We thank the reviewer for raising this point. To account for inter-individual differences in brain size\nand their potential impact on geodesic distance measures, we normalized geodesic distances by the square\nroot of the cortical surface area. This was done when analyzing distances between vertices within the top\nand bottom 5% of S-A gradient values. We then re-fitted the lifespan growth curves using these\nnormalized distances, and the results remained highly consistent with our original findings. We have\nrephrased the relevant sections in Methods and Results as follows:\nMethods (page 28): “To account for inter-individual differences in brain size, geodesic distances were\nnormalized by the square root of the cortical surface area.”\nResults (page 12): “For each individual, we calculated the mean normalized geometric distances between\npaired vertices with extreme gradient values. This geometric distance exhibited a nonlinear trajectory\nacross the lifespan, peaking during young adulthood (20.7 years, 95% bootstrap CI: 17.2–23.2), followed\nby a gradual decline (Fig. 4b).”\nQ1.18 Why multi-step Procrustes alignment using age groups, instead of constructing the gradient\ntemplate directly from the gradients of the individuals?\nR: We adopted a multi-step Procrustes alignment based on age groups to reduce potential biases arising\nfrom sample size imbalances and age-related variations in gradient topography. If a group-level gradient\ntemplate were constructed directly from all individual embeddings, the resulting template would likely be\ndisproportionately influenced by age groups with larger sample size and might obscure meaningful\ndevelopmental stage-specific differences. By first aligning gradients within narrower age groups and then\niteratively converging toward a common template, our approach ensures that the reference template\ncaptures both shared organizational features and age-specific variations in gradient topography. This\nstrategy provides a more balanced and representative framework for comparing gradients across the\nlifespan. We have rephrased the relevant sections (page 23):\n“To ensure accurate comparisons between groups, we applied Procrustes rotation 79 to align age-specific\ngroup-level functional gradients to iterative group-level reference gradients……Compared with\nconstructing the template directly from individual gradients, this multi-step alignment strategy reduces\nbiases from sample size imbalances and development variations in gradient topography.”\nQ1.19 Cortical surfaces were generated using different pipelines and hence cortical thickness\nmeasurements are not consistent. For example, iBEAT surfaces and FreeSurfer might not be entirely\nconsistent. The conclusions might be affected by this inconsistency.\nR: We used different cortical surface reconstruction pipelines for different age groups (iBEAT for infants\nand FreeSurfer for older participants). This decision was intentional 12–14, as each tool has been\nspecifically optimized for the developmental stage it targets. Applying a single pipeline across the entire\nage range would likely introduce greater bias, since adult-oriented tools such as FreeSurfer are suboptimal\nfor infant brains, while infant-tailored pipelines are not designed for adult anatomy.\nPrevious validation studies support the appropriateness of this strategy 12–14. For example, Wang et al.\n(2023) demonstrated that iBEAT V2.0, a deep learning-based pipeline optimized for infant cortical\nreconstruction, produces accurate and reproducible surfaces across scanners and imaging protocols 12.\nSimilarly, Makropoulos et al. (2018) established a neonatal minimal processing pipeline designed to\ncapture the unique anatomical features of early development 13. For older participants, the widely adopted\nHuman Connectome Project pipeline and FreeSurfer have been rigorously validated in multisite settings\nand are considered standards for adult cortical surface analysis 14.\nBy using pipelines tailored to the anatomical characteristics of each age group, we aimed to minimize\nreconstruction errors and enhance the reliability of cortical thickness estimation across the lifespan. All\nsubsequent analyses were conducted on surfaces registered to a common template space, further\nmitigating concerns about cross-pipeline inconsistency. We have rephrased the relevant sections in the\nSupplementary Methods section (page 2):\n“This age-appropriate pipeline selection minimizes cortical reconstruction bias across the lifespan, and\nall surfaces were subsequently registered to a common template to ensure cross-pipeline comparability.”\nQ1.20 The fMRI data is not necessarily scanned under the same condition, even in the lifespan HCP\ndatasets. For example, the BCP data for the first three years are scanned with the children asleep. How is\nthis discrepancy handled to ensure the validity of findings?\nR: We thank the reviewer for raising this point. Indeed, scanning conditions necessarily vary across\nlifespan datasets due to practical and ethical considerations. For example, in the Baby Connectome Project\n(BCP), infants aged 0-35 months were typically scanned during natural sleep, whereas older children\ncould be scanned either awake or asleep depending on feasibility 15. Similarly, in the Developing Human\nConnectome Project (dHCP), virtually all infants were scanned during natural sleep. In contrast, large-\nscale adolescent and adult cohorts such as ABCD, HCP, ADNI, CamCAN, and UK Biobank explicitly\n10\nrequired participants to remain awake during resting-state scans (eyes open or eyes closed, depending on\nprotocol). Unfortunately, precise per-scan records of sleep/awake state are not publicly available for dHCP\nand BCP, preventing direct stratification by brain state.\nWe fully acknowledge that brain state can influence functional activity patterns 16 and may affect gradient\nestimation. However, several factors mitigate this concern. First, the gradients observed show robust and\nbiologically plausible age-related trajectories that are consistent across datasets and scales, suggesting that\nthey are not solely driven by state differences. Second, all datasets underwent harmonized preprocessing\npipelines, which help minimize systematic biases. Third, the vary large sample size provides statistical\nstability, reducing the influence of individual-level variation in scan conditions. We now explicitly note\nthis as a limitation in the revised manuscript (page 21):\n“Finally, the large fMRI datasets exhibited heterogeneity in scanning states: neonate and infant data were\nprimarily acquired during natural sleep, whereas children and adults were scanned during wakefulness.\nGiven differences in functional brain activity between sleep and wake states 73, we cannot rule out the\npotential influence of this acquisition discrepancy on the observed patterns. Future studies that control for\nthis variable by using identical scanning states across age groups will be important to clarify this issue.”\nQ1.21 Phase I ends at one month postnatally, and Phase III begins at age five. These timepoints\ncoincide with the end of the dHCP and the end of the BCP, respectively. The dHCP, BCP, and the\nremaining lifespan HCP studies use different imaging protocols. Therefore, it is likely that the\ndemarcation of Phases I, II, and III might be an artifact of differences in image acquisition.\nR: We thank the reviewer for pointing out this point. To assess whether the phase boundaries coinciding\nwith the transitions between the dHCP, BCP, and HCP datasets could be driven by differences in\nacquisition protocols, we explicitly evaluated the potential influence of such site and protocol effects. A\nsingle harmonization across the entire lifespan would not be appropriate in this context, as the topographic\narchitecture of functional connectivity systematically changes with age. Applying one global\nharmonization procedure could inadvertently distort age-specific topographic features or impose artificial\nsimilarity across distinct developmental stages.\nInstead, we minimized acquisition-related variability by removing site-specific effects within each age\ngroup before gradient estimation. Specifically, for each age group, site was modeled as a random effect in\nthe vertex-wise FC data and regressed out prior to constructing group-level connectivity matrices and\ncomputing gradients. When gradients and clustering were derived from these within-age-group-\nharmonized data, the same three developmental phases (early anterior-posterior differentiation; emergence\nof the S-A axis in early childhood; and subsequent expansion, refinement, and stabilization) were\nreproduced and closely matched the results reported in the main text. Moreover, when determining\noptimal number of clusters, both 3- and 4-cluster solutions were identified with equal frequency across\nclustering indices. We reported the 3-cluster solution in the main text, while the 4-cluster solution only\nfurther subdivided the earliest phase (32w to 3 months). This concordance indicates that the observed\nphase boundaries reflect genuine developmental transitions rather than artifacts arising from dataset\nboundaries or acquisition differences. We have added a corresponding description in the revised Results\n(page 7) and Supplementary Results (page 3):\nResults: “The three developmental phases were also identified through both individual-level gradient\nclustering without predefined age bins and within-age-group-harmonized analysis (Supplementary Fig. 3,\n4). This indicates that the observed patterns were robust and not driven by specific age-grouping schemes\nor site related confounds.”\nSupplementary Results: “To ensure that the phase boundaries were not artifacts of differences in\nacquisition protocols across datasets, site was modeled as random effect and regressed out within each\nage group prior to computing group-level connectivity matrices and gradients. Analyses using these\n11\nharmonized data reproduced the same three phases, indicating that the observed phase boundaries reflect\ngenuine developmental transitions rather than dataset or acquisition differences (Supplementary Fig. 4).”\nSupplementary Fig. 4 | Distinct growth stages of the sensorimotor-association gradient across the lifespan from the\nharmonized functional connectome. a, S-A gradient maps generated for 26 distinct age groups derived from the\nharmonized functional connectome. b, K-means clustering analysis based on the pairwise correlations between the 26\ngroup-based gradient maps, resulting in three clusters (main text) and an alternative four-cluster solution that further\nsubdivides the earliest phase. c, Three distinct phases of the S-A gradient are visualized via multidimensional scaling\n(MDS). “Dimension 1” and “dimension 2” are defined as the first and second axes obtained through MDS of the\nfunctional distance matrix. d, Distribution of gradient values in the S-A axis across lifespan. w, week; m, month; yr, year.\nQ1.22 How is longitudinal information, repeated observations of each subject, dealt with? Not doing so\nmight bias the results.\nR: To minimize potential bias introduced by repeated measurements from the same participants, we\nrestricted our analyses to a single imaging session per individual across all longitudinal datasets. This\nensure that each participant contributed only once to the cross-sectional lifespan analyses, thereby\npreventing inflated sample sizes or undue weighting of individuals with repeated scans. This procedure is\ndescribed in the Methods section of the revised manuscript (page 22):\n“For participants with multiple imaging sessions, only the first scan was selected for inclusion, ensuring\n12\nthat each individual contributed a single observation to the analyses.”\nQ1.23 The effects of sex may vary nonlinearly across age, especially given the long developmental\ntime span covered in this study. However, this nonlinearity does not appear to have been modeled.\nR: To test whether the effects of sex on gradient measures vary nonlinearly across the lifespan, in the\nrevised manuscript we performed additional analyses at the global, system, and regional levels.\nSpecifically, we used the GAMLSS framework, which flexibly captures nonlinear trajectories across age.\nTwo nested models were compared: a baseline model including main effects of age and sex, and a full\nmodel including an age-by-sex interaction term. Both models controlled for mean FD as a covariate, and\nsite was modeled as a random effect. Significance of the interaction was assessed using the generalized\nlikelihood ratio tests (GLRT), with FDR correction applied (q < 0.05).\nAt the global level, females showed significantly higher explanation ratio, gradient range, and variation\ncompared with males across the lifespan (all p<0.05, FDR-corrected; Supplementary Fig. 35a). At the\nsystem level, the overall hierarchical organization was preserved, with the default-mode and frontoparietal\nsystems consistently showed higher gradient ranges and standard deviations, while the sensorimotor\nsystem displayed lower values throughout life. Importantly, females exhibited significantly greater\ngradient ranges and variations than males (all p<0.001, FDR-corrected; Supplementary Fig. 35b).\nMoreover, peak age analyses revealed that in several functional systems, including visual, limbic,\nfrontoparietal, and default mode networks, females demonstrated significantly later maturational peaks\nthan males (all p<0.001, FDR-corrected; Supplementary Fig. 35c & Supplementary Table 8), suggesting\ndelayed maturational timing in females. At the regional level, females exhibited higher gradient values in\ntemporal, precuneus, medial prefrontal and somatomotor regions, whereas males showed higher values in\noccipital and lateral prefrontal regions (p<0.001, FDR-corrected; Supplementary Fig. 35d).\nThese findings have been incorporated into the Supplementary Methods (page 3) and Results sections\n(page 7):\nSupplementary Methods:\n“To test whether the effects of sex on gradient measures vary nonlinearly across the lifespan, we\nperformed additional analyses at the global, system, and regional levels. Specifically, we used the\nGAMLSS framework, which flexibly captures nonlinear trajectories across age. Two nested models were\ncompared: a baseline model including main effects of age and sex, and a full model including an age-by-\nsex interaction term. Both models, controlled for mean FD as a covariate, and site was modeled as a\nrandom effect. Significance of the interaction was assessed by comparing the model fits of the baseline\nand full models using the generalized likelihood ratio test (GLRT). All p-values were corrected for\nmultiple comparisons using the FDR (q < 0.05). The direction of sex differences was determined based on\nthe mean differences between the fitted sex-specific growth curves.”\nSupplementary Results:\n“To investigate whether sex differences in cortical gradient measures vary across age"}]

Decision letter