Spatial heterogeneity and subtypes of functional connectivity development in youth

Domain: neuroscience
Published: 2026-01-23
Editorial outcome: major_revision

Abstract
Brain age prediction has been widely utilized to assess functional connectivity (FC) development, but conventional global brain age indices are limited in capturing spatial heterogeneity across the cortex. This study introduces a regional brain development index to characterize fine-grained FC maturation across cortical regions. We examined its spatial variability and stratified individuals into subtypes with distinct region-wise FC developmental patterns. Using data from the Philadelphia Neurodevelopmental Cohort (ages 8-23 years), we identified three distinct subtypes and found that individuals with advanced FC developmental pattern aligning with the sensorimotor-association axis exhibited superior cognitive performance. Robustness was confirmed through replication in the Human Connectome Project Development cohort. Further analyses revealed associations between FC development and gene expression linked to neural differentiation, synaptogenesis, and myelination. These findings suggest that spatial heterogeneity in FC development reflects cortical microstructure and hierarchical organization, underscoring its critical role in neurocognitive maturation during youth.

Introduction
The development of functional networks is a dynamic and complex process closely linked to ongoing maturation of the brain’s functional organization, which underpins the progressive enhancement of cognitive, emotional, and social capacities throughout childhood and adolescence1,2. Research has established the critical role of cortical hierarchy in brain function, and accumulating evidence has demonstrated that the cortical development follows spatial and temporal patterns aligned with this hierarchy from perspectives of macroscale neurobiological properties1,3 and underlying microscale mechanisms in developmental plasticity, such as cortical excitation, inhibition, and glial cell types1,4,5. Recent studies further show that the spatial-temporal patterns of functional network development in youth systematically align with cortical hierarchy organization, particularly along the sensorimotor-association (S-A) axis, with developmental increases in functional connectivity (FC) in sensorimotor cortices and decreases in association cortices6. However, how these FC development patterns vary across individuals remains unclear. Neurocognition undergoes a prolonged development trajectory parallel to the maturation of the brain’s functional organization, with individual differences in youth linked to important outcomes such as academic performance and quality of life2, while deficits in neurocognition are associated with vulnerability to mental illness. Therefore, advancing our understanding of FC development at the individual level may provide valuable insights into typical development processes and deviations from these typical trajectories, with implications for cognitive and psychological functioning as well as mental illness.

Brain age has increasingly been utilized to evaluate brain status and track brain changes at the individual level, leveraging neuroimaging data across various modalities7–22. The difference between brain age and chronological age, known as the brain age gap (BAG), provides a quantitative and interpretable metric for characterizing brain changes associated with development, aging, and brain disorders. BAG measures derived from functional connectome have revealed that individuals with accelerated brain development exhibit enhanced FC in specific functional brain networks, including the default mode network (DMN), during childhood and adolescence12. Existing studies have established a linkage between brain age derived from structural imaging data and cognition, demonstrating that individuals with advanced brain development exhibit superior cognitive processing speed compared to those with delayed development8. While studies have identified associations between FC and neurocognition in youth using predictive modeling23–26, it remains unclear whether and how brain age estimates derived from the functional connectome are associated with neurocognition in youth.

Most existing studies estimate brain age from a global perspective by mapping whole-brain, high-dimensional FC profiles to a single unitary index, which inherently overlook the spatial variability of FC changes across brain regions. To address this limitation, several studies have proposed fine-grained modeling at the regional level10,13,15,17,18,20,27. Distinct and heritable patterns of brain aging, derived from structural imaging features, have been linked to brain disorders, such as schizophrenia, multiple sclerosis, and dementia10. Additionally, heterogenous aging patterns across brain regions in schizophrenia have been revealed using multimodal neuroimaging data20. Advanced brain aging in specific brain regions has also been shown to mediate cognitive outcomes in cerebral small vessel disease, highlighting the potential to capture heterogeneous brain changes beyond a unitary measure of global brain age15. Given that variations in FC profiles, functional topography, and structural-functional coupling are associated with brain development and neurocognition2,6,9,11,12,23,24,26–29, exploring regional patterns of brain development measures derived from FC and their associations with neurocognition could provide valuable insights into the development of brain functional organization and neurocognitive processes.

In this study, we utilized a large sample (total n = 1149) of youths aged 5–23 years from the Philadelphia Neurodevelopmental Cohort (PNC)30 and the Human Connectome Project Development cohort (HCP-D)31 to investigate region-wise FC development using functional MR imaging (fMRI) data. For each cortical region, we developed a region-specific age prediction model based on its regional FC profile. Each model provides an individualized brain age measure, and the difference between this measure and chronological age, referred to as the regional brain development (RBD) index, reflects region-wise FC development. We hypothesized that the RBD index would exhibit spatial variability and inter-subject heterogeneity, linked to developmental differences in neurocognition. To explore the spatial variability of regional FC development and its associations with neurocognition, we stratified individuals into subtypes characterized by distinct FC development patterns using a cluster analysis based on individual subjects’ whole-brain RBD indices. Subsequently, we investigated differences in cognitive measures across these subtypes. Our analysis revealed that individuals with advanced FC development parallel to the hierarchical cortical organization defined by the S-A axis exhibited superior cognitive profiles. To investigate the biological mechanism underlying the spatial variability of FC development, we examined the relationship between regional FC development patterns and gene expression. We predicted that regional FC development patterns aligning with the cortical hierarchy organization would show enrichment in the expression of genes associated with neuron development.

Results
We used fMRI data from the PNC (n = 693; ages 8–23) to examine the spatial variability and inter-individual heterogeneity in regional FC development. The RBD index, defined as brain age minus chronological age and corrected for bias and covariates, was used as a surrogate measure to assess FC development across cortical regions. Specifically, for each cortical region, a dedicated brain age prediction model was built on its FC profile, which was defined as its FC measures to all the other cortical regions, with Pearson’s correlation coefficient serving as the FC measure. Five-fold cross-validation was implemented to predict brain age across all individuals and cortical regions. An overview of the analysis process is illustrated in Fig. 1. Using the RBD index, we began by examining spatial variability in FC development and its inter-individual heterogeneity. Next, we applied clustering to identify covarying patterns of FC development across cortical regions, stratifying individuals into subtypes characterized by distinct regional FC development profiles. We then compared the spatial patterns across subtypes and assessed their alignment with established principles of cortical hierarchy in brain development. To further characterize these subtypes, we analyzed their neurocognitive profiles using cognitive measures and meta-analytic functional decoding. We subsequently explored the biological underpinnings of these FC development patterns by examining their associations with gene expression. Finally, we evaluated the generalizability and robustness of the identified subtypes of FC development using data from an independent cohort.

Spatial variation and inter-individual heterogeneity in functional connectivity development
To characterize the spatial variation of FC development across cortical regions, we analyzed differences in their capacity for age prediction. Two metrics, accuracy (r) and mean absolute error (MAE), were used to quantify each region’s predictive capacity, with accuracy (r) computed as the Pearson’s correlation coefficient between the predicted age (prior to bias correction) and chronological age. Cortical regions exhibited varying capacity for age prediction (Fig. 2A), with certain functional networks showing higher predictive capacity, as indicated by higher r and lower MAE values. These networks included the somatomotor network, ventral attention network (VAN), and frontoparietal network (Fig. 2B). However, heterogeneity was evident within individual functional networks: for example, the medial prefrontal aspect of the DMN displayed higher predictive capacity compared with lateral elements of the DMN.

Building on the observed variation in FC development across cortical regions, we next investigated inter-individual differences in the pace of regional FC development. To assess this, we quantified the brain development for each cortical region and individual using the RBD index, defined as brain age minus chronological age and corrected for bias and covariates32, with positive/negative values indicating advanced/delayed development. Individuals of the same chronological age (9, 13, or 17 years) showed variability in region-wise brain development (Fig. 2C), highlighting inter-individual variability in the pace of regional FC development. Additionally, distinct spatial covariation patterns in RBD were evident at the individual level. For instance, the 13-year-old individual (Fig. 2C, bottom-middle) exhibited consistently advanced development in frontoparietal and DMN regions alongside delayed development in sensorimotor regions, whereas the 9-year-old individual (Fig. 2C, top left) showed a largely opposite pattern.

Subtypes of functional connectivity development and their links to neurocognition
To investigate the observed inter-individual heterogeneity and potential spatial covariation in regional FC development, we identified covarying patterns of FC development and stratified individuals into subgroups based on their whole-brain RBD index maps. Using the K-means clustering algorithm, individuals were grouped into three subtypes, with each subtype characterized by a distinct FC development pattern (P1, P2, and P3) represented by the centroids of their RBD index maps (Fig. 3A). Specifically, P1 showed a global trend of delayed FC development across cortical regions; P2 exhibited advanced FC development localized to the association cortex, prominently within DMN regions; P3 displayed widespread advanced FC development, particularly pronounced in sensorimotor regions.

To further characterize the spatial variation of these distinctive development patterns, we examined their alignment with the S-A axis, which represents a unified cortical organization across diverse neurobiological properties and has been demonstrated to be significantly associated with FC development6. While the development pattern of P1 did not significantly align with the S-A axis (Fig. 3B, Spin-test, r = 0.186, p = 0.188), those of P2 and P3 showed significant correlations, with P2 showing a strong positive correlation (Fig. 3C, Spin-test, r = 0.580, p < 0.0001) and P3 showing a strong negative correlation (Fig. 3D, Spin-test, r = −0.511, p = 0.0001). Notably, the cortical regions with predominantly advanced development in the P2 were located at the high-ranking association end of the S-A axis, while those in P3 were located at the low-ranking sensorimotor end. Given the overall advanced FC development of P2 and P3, we then examined their differences in FC development across regions. Specifically, P2 demonstrated significantly more advanced FC development in DMN regions compared with P3, whereas P3 demonstrated significantly greater advanced development in sensorimotor regions relative to P2 (Fig. 3E).

We next analyzed the cognitive performance of individuals across subtypes in three cognitive domains: executive function, social cognition, and memory, to investigate their behavioral implications. Individuals in P2 demonstrated significantly higher cognitive performance than those in P1 and P3 (Fig. 3F), with respect to executive function (Wilcoxon rank sum test, z = 3.73 and p = 0.0002 for P2 vs. P1, z = 2.59 and p = 0.010 for P2 vs. P3), social cognition (Wilcoxon rank sum test, z = 3.80 and p = 0.0001 for P2 vs. P1, z = 2.36 and p = 0.019 for P2 vs. P3), and memory (Wilcoxon rank sum test, z = 2.17 and p = 0.030 for P2 vs. P1, z = 2.24 and p = 0.025 for P2 vs. P3). Interestingly, no significant differences were observed between P1 and P3, despite globally advanced FC development of individuals in P3 (Fig. S1B). There were no significant differences across subgroups regarding chronological age and sex. These findings highlight the critical role of spatial heterogeneity in FC development during youth, particularly within regions that support neurocognitive functions.

To better understand the relationship between distinctive FC development patterns and properties of the brain connectome, we also analyzed their alignment with the principal functional gradient and cortical myelin content. The P2 showed a strong positive correlation with the principal functional gradient (Fig. 3G, Spin-test, r = 0.604, p < 0.0001), consistent with prior findings that hierarchical functional organization supports neurocognition development. Moreover, the P2 exhibited a significantly negative correlation with cortical myelin content (Fig. 3H, Spin-test, r = −0.512, p < 0.0001), suggesting a relationship with the ongoing process of myelination during brain maturation. On the other hand, the P3 exhibited a reversed alignment pattern, showing a negative correlation with the principal functional gradient and a positive correlation with the cortical myelin content (Fig. S2). These findings underscore the critical interplay between functional organization and structural maturation of brain connectome in shaping neurocognitive outcomes.

In addition, we utilized a meta-analytic functional decoding technique33 based on the NeuroQuery database34 to interpret the functional significance of the identified FC development patterns. The decoding results also revealed significant associations between the P2 and neurocognition-related terms, such as cognitive, decision, and context (Fig. 3I).

We next conducted a series of sensitivity analyses to investigate the robustness of the identified FC development patterns. The subtypes of FC development remained robust in terms of spatial patterns, cognitive profiles, and individual-level subtype assignments, when controlling for additional covariates including socioeconomic status (SES) and parental education (Fig. S5), when using Gaussian Mixture Models (GMM) based clustering for the subtyping (Fig. S6), when using the Schaefer-300 atlas and Schaefer-200 atlas for cortical parcellation (Fig. S7), and when using the Fisher’s r-to-z transformed FC measures to derive the RBD index (Fig. S8).

Associations between regional FC development and gene expression
To investigate the biological mechanism underlying the regional FC development, we analyzed the relationship between the FC development patterns and gene expression data from the AHBA35. We first examined the alignment of FC development patterns with the first principal component (PC1) of gene expression36. The P2 exhibited significant negative correlation (Fig. 4A, Spin-test, r = −0.440, p = 0.0017), while no significant correlations were observed for P1 and P3 (Fig. S3). We then used spatial correlation analysis to examine the association between the P2 pattern and the expression of 4112 genes in the cortex. This analysis identified 735 positively correlated genes and 783 negatively correlated genes (q < 0.05, FDR corrected, Spin-test). Among them, the neuronal differentiation 6 (NEUROD6) gene was the most positively correlated, while the Mitochondrial Ribosomal Protein L48 (MRPL48) gene was the most negatively correlated (Fig. 4B). Using Specific Expression Analysis (SEA)37, we found that the set of positively correlated genes was enriched for neurodevelopmental genes with cortical expression spanning stages from early childhood to young adulthood (Fig. 4C), suggesting that these genes may play a role in supporting the regional FC development in youth.

We further performed Gene Ontology (GO) enrichment analysis on the set of correlated genes using the online tool of GOrilla38, which revealed that the significantly correlated genes were primarily enriched in GO terms associated with “membrane”, “neuron projection”, and “synapse” (q < 0.05, FDR corrected, Fig. 4D). Given that regional differences in cortical gene expression may be influenced by differences in cellular composition of the cortex, we conducted a cell-type-specific enrichment analysis39 to examine the enrichment of correlated genes in adult cortical cell types. This analysis revealed that the correlated genes were enriched across multiple adult cell types, including astrocytes, oligodendrocyte precursor cells, inhibitory and excitatory neurons (Fig. 4E).

The FC development patterns are generalizable to the HCP-Development cohort
To assess the robustness of our findings, we investigated regional FC development using data from the HCP-D cohort (n = 456; ages 5–22) and evaluated whether the identified FC development patterns could stratify the HCP-D individuals into distinct subgroups. Five-fold cross-validation was utilized to predict brain age for all cortical regions and individuals in the HCP-D cohort, as in the previous analyses of the PNC cohort. The spatial variation of FC development across cortical regions exhibited a pattern consistent with that observed in the PNC cohort (Spin-test, r = 0.467, p < 0.0001), with regions showing higher predictive accuracy primarily located in the somatomotor network, VAN, and the anterior part of the DMN (Fig. 5A, B). Next, HCP-D individuals were stratified into three subgroups based on the FC development patterns identified in the PNC cohort. Specifically, each individual’s RBD index map was compared with the PNC cohort’s centroids of RBD index maps (P1, P2, and P3), and the individual was assigned to the closest matching subgroup (Fig. 5C). We then analyzed cognitive performance of individuals across subgroups using a battery of cognitive assessments relevant to executive function, memory, and general cognition. Overall, individuals in subgroup P2 demonstrated the highest cognitive performance, particularly on the Flanker Test (Wilcoxon rank sum test, z = 2.90, p = 0.004) and the List Sort Task (Wilcoxon rank sum test, z = 2.15, p = 0.032) compared to subgroup P1 (Fig. 5D). No significant differences were observed between subgroups P1 and P3 with respect to these two assessments. Individuals in subgroup P2 also demonstrated the highest cognitive performance in terms of other cognitive domains, including fluid cognition, crystallized cognition, and total cognition (Fig. S9). These findings demonstrate that the three FC development patterns identified in the PNC cohort are robust and generalizable to the HCP-D cohort. Moreover, the results reinforce that the spatial heterogeneity in FC development across individuals is linked to neurocognitive performance. Replication analysis applying the same clustering procedure within the HCP-D cohort revealed consistent FC development patterns and subtyping results (Fig. S10), further supporting the reproducibility of the findings.

Methods
Participants and data processing
Data were drawn from two datasets: the PNC (n = 693) and the Human Connectome Project Development (HCP-D; n = 456).

The PNC30 is a community sample of children and adolescents from the greater Philadelphia area, collected for studying brain development, including 1601 participants with neuroimaging data. For this study, demographic, cognitive assessment, and neuroimaging data from 693 participants (301 males and 392 females) aged 8–23 years were included after applying the inclusion criteria as previously described11. Informed consent was obtained from all participants or their parents/guardians, and minors provided assent. All study procedures were approved by the institutional review boards of both the University of Pennsylvania and the Children’s Hospital of Philadelphia. One resting-state and two task-based (n-back and emotion identification) fMRI scans for each participant were used in this study. The neuroimaging preprocessing followed previously described procedures and methods11. Cognitive assessment was conducted using the Penn computerized neurocognitive battery (Penn CNB), and three factor-analysis based composite scores57, including executive and complex cognition, social cognition, and episodic memory, were included for analysis in this study.

The HCP-D31 includes 652 participants recruited across four sites in the U.S. to study healthy brain development in children and adolescents. A sample of 456 participants (214 males and 242 females) aged 5–22 years with complete demographic, cognitive assessment, and neuroimaging data were included in this study after applying the inclusion criteria as previously described6. All study procedures were approved by a central Institutional Review Board at Washington University in St. Louis. Resting-state and task-based (carit, emotion, and guessing tasks) fMRI scans were used for all participants. The neuroimaging acquisition and preprocessing procedures have been previously described6. Participants completed a battery of cognitive assessments, including segments of the NIH toolbox and others31. Cognitive scores from the Dimensional Change Card Sort Task, the Flanker Test, the ListSort Working Memory Task, the Oral Reading Test, the Pattern Comparison Processing Speed Task, the Picture Sequence Memory Task, and the Picture Vocabulary Test were included in this study. Composite scores, including fluid cognition, crystalized cognition, and total cognition from the HCP-D data release31 were also included. Demographic information and variables of interest for both datasets may be found in Tables 1 and 2.

For both samples, preprocessed fMRI time series were parcellated into 400 cortical regions using the Schaefer-400 atlas58, and time series from resting-state and task fMRI scans were concatenated to compute FC measures for the analysis.

Brain development index
For region-specific brain age prediction, a dedicated linear ridge regression model to predict chronological age was built for each cortical region on its FC profile that was defined as its FC values to all the other cortical regions, with each value measured as a Pearson’s correlation coefficient between the two regions. Five-fold cross-validation was utilized to obtain predicted brain age for all individuals and cortical regions.

For global brain age prediction, the whole-brain functional connectome (upper triangle entries of the FC matrix, containing FC values between all pairs of cortical regions) was used as features, and a linear ridge regression model was trained for the prediction of chronological age.

Both the region-specific and global prediction models were built and evaluated with the same five-fold cross-validation, with identical training and evaluation data splits. For each prediction model, the regularization parameter in ridge regression was optimized in the range [2^10, 2^9, …, 2^4, 2^5] using nested cross-validation within the training sample. The prediction models were implemented using the scikit-learn library (v1.4.1) in Python (v3.9.18).

Given the predicted brain age for each cortical region, the difference between this measure and chronological age was referred as the RBD index, reflecting region-specific FC development. To alleviate the known systematic bias in the brain age prediction due to the “regression to the mean” effect, the RBD index was adjusted by a linear bias correction procedure32. In addition, covariates including sex, head motion, and site information (for the HCP-D sample) were regressed out from the RBD index for all subsequent analysis. The same procedure was applied to obtain global BAG measures based on the predicted global brain age.

Subtyping FC development
Given the RBD index maps for all participants of the PNC cohort, the K-means clustering algorithm was utilized to identify subtypes of FC development and stratify the participants into different subgroups. Specifically, the clustering was performed with whole-brain RBD indices as features and cosine distance as the distance metric, as implemented in MATLAB (R2021a).

The number of clusters was set to 3, determined based on both reproducibility and within-cluster point-to-centroid distances of the clustering results. Specifically, all participants were randomly split into two halves, and clustering was performed on each half separately with the number of clusters k set from 2 to 7. This procedure was repeated for 100 times. The clustering results for each half was denoted as C_k,1 and C_k,2. For each k, the clustering centroids of the first half were then applied to samples in the second half to obtain the clustering results C_k,1->2, and vice versa to obtain C_k,2->1. The average of RandIndex(C_k,1, C_k,2->1) and RandIndex(C_k,2, C_k,1->2) was adopted as the reproducibility measure for clustering results with k clusters, where RandIndex(C, C') refers to the Rand index for two sets of clustering results C and C'. Additionally, to evaluate the data-fitting performance of the clustering results, the distance of each data point to its corresponding cluster centroid was calculated and averaged across all data points. The optimal number of clusters was determined by reaching a balance between reproducibility and data-fitting performance, where the two measures interacted as shown in Fig. S4A.

In addition to reproducibility based and data-fitting based measures (Fig. S4A), we further justified the choice of number of clusters using the NbClust package59 and according to the data distribution characteristics. Particularly, NbClust uses a wide variety of indices to determine the optimal number of clusters using Hierarchical clustering or K-means methods. The most frequently selected cluster numbers across clustering methods were 2 and 3 (Fig. S4B), suggesting them as optimal choices for our dataset. To further verify this finding, we examined the distribution of individual RBD patterns in a 2D embedding space using t-SNE projection (Fig. S4C). It was observed that clusters P1 and P3 in the 3-clusters solution aligned closely with clusters P1 and P2 in the 2-cluster solution, while cluster P2 in the 3-cluster solution captured fine-grained inter-individual variability between P1 and P3. Based on all these quantitative and qualitative assessments, the number of clusters was set to 3 in our analysis.

The HCP-D individuals were stratified into distinctives subgroups based on the FC development patterns identified in the PNC cohort for the generalization analyses. Specifically, each individual’s RBD index map was compared to the PNC cohort’s centroid maps, and the individual was assigned to the subgroup with the closest match. For instance, individuals i and i+1 were assigned to subgroups P2 and P1, respectively, as their RBD index maps showed the highest similarity to the corresponding centroids (Fig. 5C). To further assess the robustness of the identified FC development patterns, the clustering procedure was also replicated within the HCP-D cohort.

Meta-analytic functional decoding of FC development pattern
The meta-analytic decoding tool33 was used to link the FC development pattern to functional/behavioral profiles. Specifically, the functional terms and term-based meta-analytic maps from the NeuroQuery database34 was utilized with the correlation decoder33 for the functional decoding. The functional decoding results were illustrated by a word cloud, with larger font size indicating higher functional relevance.

Gene enrichment analysis
Gene expression data was obtained from the Allen Human Brain Atlas (AHBA)35 and processed with the abagen toolbox60 and the Schaefer-400 atlas in the MNI space to obtain regional gene expression. Since not all genes in the AHBA dataset are assumed to be expressed in the brain, only brain-expressed genes were retained for analysis, following an established procedure61. Additionally, the analysis focused on the left hemisphere with 200 brain regions to minimize sample variability across brain regions, as the gene expression data were sampled from both hemispheres of only two of the six donors62. A total of 4112 genes were included in our analysis.

The spatial alignment of regional FC development patterns and the expression patterns of each gene was quantified using Peason’s correlation coefficient, and its significance was tested using the spin test with 1,0000 permutations63,64. The set of significantly correlated genes was included in the subsequent enrichment analysis. A SEA tool37 was first used to examine the enrichment for cortical-expressed genes at different developmental stages through childhood to young adulthood. A rank-based GO enrichment analysis using Gorilla38 was then conducted to test the enrichment for the ontology of cellular component. As the regional differences in gene expression may reflect underlying differences in cellular composition across cortical regions, a median rank-based approach61 was adopted to assess the enrichment for the adult cell types39.

Sensitivity analysis
To investigate whether our findings were robust and consistent across analytic choices, such as inclusion of additional covariates, clustering algorithms, atlas used for cortical parcellation, and computation of FC measures, a series of sensitivity analyses were conducted. Specifically, FC development subtyping was evaluated by further including SES and parental education as covariates to derive the RBD index. Neighborhood SES measure quantified by geocoded variables as in previous studies65,66 was adopted for the PNC cohort, and SES measure quantified as the natural log of the income-to-needs ratio was adopted for the HCP-D cohort as in ref. 67. Furthermore, subtyping results using GMM based clustering were evaluated, as suggested that K-means and GMM perform well in terms of stability and generalizability by a prior study on MRI based subtyping68. Finally, subtyping analyses were evaluated using additional cortical parcellations, including the Schaefer-200 atlas and Schaefer-300 atlas. The subtyping results using Fisher’s r-to-z transformed correlation values as FC measures were also evaluated.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Spatial Heterogeneity and Subtypes of Functional Connectivity Development in Youth——Comments\nLi and colleagues leverage large-sample fMRI data and machine learning-based brain age prediction to propose a regional\nbrain development index (RBD) and identify three individual-level subtypes of functional connectivity (FC) over\ndevelopment. They further demonstrate how these subtypes relate to cognition, cortical hierarchy, myelination, and gene\nexpression—linking multiple biological levels. Notably, the P2 subtype (characterized by advanced association cortex\ndevelopment) shows better cognitive performance across two independent cohorts, highlighting the robustness and\nscientific value of their findings. This work pushes beyond traditional global metrics (like brain age gap, BAG), offering the\nfirst demonstration of FC developmental spatial heterogeneity and its multi-level implications at the individual level. Despite\nthese strengths, I do have several comments and suggestions to improve the manuscript:\n1. It seems that the authors do not adequately account for critical confounders such as socioeconomic status (SES), parental\neducation, or health conditions. The large-scale cohort studies such as ABCD and IMAGEN routinely include these as\ncovariates, since they are tightly linked to brain structure, function, and cognitive development in youth (Jirsaraie et al., 2023\nNat Commun; Noble et al., 2015 Nature). I thus suggest that the authors include these variables as covariates in their main\nor supplementary analyses, or at least stratifying their results accordingly. Doing so would substantially strengthen the\ncredibility and scientific validity of your conclusions.\n2. Although the three FC developmental subtypes differ in executive function, social cognition, and memory, the cognitive\nbattery is relatively narrow and does not capture a broader spectrum of behavioral and psychological abilities, such as\nemotion processing, motivational regulation, or creative thinking. Given the multidimensional nature of adolescent behavior\nand cognition—and their relevance to neurodevelopment (see Casey et al., 2018, Biol Psychiatry)—I encourage the authors\nto more fully acknowledge this limitation and, where possible, connect with additional behavioral or affective dimensions\navailable in open datasets.\n3. The proposed three-subtype spatial clustering is mainly compared to the conventional whole-brain BAG, but there is no\ndirect benchmarking against other established spatial clustering methods—such as structural brain age subtyping (Wang et\nal., 2019, Nat Commun), spectral clustering (Satterthwaite et al., 2015, Neuron), or GMM. Please consider providing a\nsystematic comparison of clustering performance (e.g., ARI, NMI, stability), cognitive associations, and spatial distribution\noverlap to clarify the innovation and boundaries of your approach. It would be also helpful if the authors can provide further\njustification for the chosen number of clusters (k), including cross-cohort reproducibility or biological interpretability.\n4. While the manuscript connects gene expression, myelination, FC, and cognition, most analyses are correlational. It is now\nwell accepted that relationships across biological scales are often complex, with mediation, moderation, and feedback loops\n(see Fornito et al., 2015, Nat Rev Neurosci). I would suggest explicitly discussing potential mechanistic links between the\nmicro, meso, and macro levels (e.g., via multimodal causal modeling, SEM, network mediation analysis) and acknowledging\nthe current limitations of causal inference.\n5. The analyses rely on two large cross-sectional cohort studies (PNC and HCP-D), which are valuable for describing group-\nlevel differences across age but cannot capture within-individual developmental trajectories. The field broadly recognizes\nthat longitudinal follow-up is necessary to truly understand developmental rate, sequence, and causality (see Kraemer et al.,\n2000, Am J Psychiatry; Tamnes et al., 2017, Nat Neurosci). Please discuss this limitation, and emphasize the need for future\nlongitudinal studies (e.g., ABCD, IMAGEN) that can follow the same individuals over time to more accurately characterize\nbrain developmental heterogeneity and its cognitive and molecular correlates.\n6. The introduction does not sufficiently relate classical neurodevelopmental theories (e.g., Felleman & Van Essen, 1991,\ncortical hierarchy; Huttenlocher, 1979, synaptic pruning; Hebb, 1949, plasticity) to the novel aspects of spatial subtyping. I\nwould recommend more explicitly connecting the study’s approach to these foundational models in both the introduction and\ndiscussion, which would help a broader audience appreciate the conceptual advance.\n7. The results are organized as spatial variability/heterogeneity, subtyping, cognition/gene associations, and replication.\nHowever, the scientific rationale and linkage between each step could be more clearly articulated—for instance, why certain\nmolecular analyses were performed post-subtyping, and how each step is motivated by prior findings or hypotheses.\nExplicitly stating the underlying logic for each major analysis will make the paper more cohesive and accessible.\n8. Figure 1 shows the RBD analysis pipeline, but does not clearly label model inputs and outputs, the formula for RBD, or\nhow RBD feeds into downstream clustering and association analyses. Please consider annotating workflow nodes, core\nequations, and analytic steps to improve interpretability, especially for interdisciplinary readers.\n9. Currently, subtyping is based only on brain imaging. There is increasing interest in behavioral (e.g., LCA) or joint brain-\nbehavior subtyping approaches (see Drysdale et al., 2017, Am J Psychiatry; Smith et al., 2022, Trends Cogn Sci). I\nencourage the authors to discuss prospects for integrating behavioral subtyping with brain-based subtyping, and perhaps\neven attempt a preliminary behavioral subtyping analysis with available data to explore overlap and complementarity, which\ncould inform future stratified intervention strategies.\n10. All major analyses are correlational, and there is no direct causal chain established between FC subtype, cognition, or\ngene expression. Please make a clear distinction between correlation and causality in the discussion and avoid over-\ninterpretation. Where possible, supplement with sensitivity analyses—such as subgroup analysis by age, sex, or SES;\nmediation analysis or SEM to explore whether FC subtype mediates links between molecular/structural factors and cognitive\noutcomes. Also, highlight that future work leveraging longitudinal data, intervention, or multi-omic integration will be needed\nto clarify causal mechanisms. Such enhancements would greatly increase the scientific impact and robustness of findings.", "reviewer_2": "In this manuscript, Li and colleagues examined spatial patterns of functional brain development across two normative\nsamples and characterized data-driven subgroups by their cognitive profiles. The study is well motivated: deriving regional\nbrain age estimates addresses an important limitation of traditional brain age studies that rely on a single global metric.\nMy primary concern relates to the authors' approach to replicating the data-driven analysis. The authors applied k-means\nclustering to identify subgroups based on patterns of regional brain development (RBD) indices in their primary cohort (PNC)\nand compared cognitive scores across subgroups. Because clustering will always produce a solution, independent\nreplication is critical to demonstrate that subgroups are meaningful rather than artifacts. To address this, the authors\nclassified individuals from their replication cohort (HCP-D) using PNC-derived centroids and concluded that subgroup\nreplicability was supported by similar cognitive differences across datasets. However, this approach does not validate\nreplicability.\nTrue replication would require rerunning clustering independently in both cohorts and demonstrating that (1) subgroup\npatterns of RBD indices and (2) associated cognitive profiles are similar. Without establishing (1), it remains unclear whether\ncomparable subgroup structures exist. Classifying HCP-D participants based on PNC centroids assumes, rather than tests,\nsubgroup generalizability and risks overestimating replicability. For example, it is theoretically possible that multiple,\ndrastically different RBD patterns could yield similar distributions of cognitive scores without reflecting the same underlying\nneurodevelopmental patterns. To convincingly demonstrate replicability, the authors should rerun clustering independently\nin the HCP-D cohort and compare both subgroup structures and cognitive profiles to those observed in PNC.\nSeparately, the cognitive measures differ substantially between the two datasets (executive functioning, social cognition,\nand memory in PNC, versus fluid cognition, crystallized cognition, and total cognition in HCP-D). As a result, even if\nsubgroup differences were consistent in direction or magnitude, it would be difficult to interpret this as evidence of\nreplication, given that the constructs are not directly comparable. To strengthen their conclusions, the authors should\nharmonize cognitive domains across cohorts or more explicitly acknowledge this limitation.\nAdditional comments for improving the manuscript:\n• The use of cross-sectional data to study development is a major limitation that is not discussed in the manuscript. This\nshould be acknowledged.\n• Normative modeling, which similarly derives individualized regional deviations from expected trajectories (and can\nincorporate non-linearity and covariates, stated limitations of the chosen approach), is not mentioned. While I do not expect a\nchange in approach, a brief rationale for choosing RBD over normative modeling would be helpful.\n• The choice of the Schaefer 400 parcellation should be justified, particularly given evidence suggesting that the 200-parcel\nresolution produces more representative connectivity networks (Luppi et al., 2021). Additionally, the omission of the typically\npaired Melbourne subcortical atlas is a missed opportunity to capture subcortical contributions.\n• Raw correlation values were used as inputs to the ridge regression model. Typically, Fisher r-to-z transformation is applied\nto stabilize variance and normalize correlations; the rationale for deviating from this standard practice should be provided.\n• I appreciate the authors' approach to determining the optimal number of clusters. However, it would strengthen the\nmanuscript to include more standard clustering validity metrics (e.g., Silhouette coefficient) in the supplemental material for\ncomparison.\n• More detail is needed on the clustering methodology: Why was k-means chosen over alternative algorithms (e.g., spectral\nclustering)? Why was the cosine distance metric used?\n• Complete statistical reporting is needed throughout the manuscript, including test statistics, exact p-values, and effect\nsizes. Currently, many results are only partially reported (e.g., indicating p < 0.05 without further detail).\n• The manuscript should specify that Python and scikit-learn were used for the ridge regression modeling, including software\nversions.\n• A table summarizing cohort characteristics (age, age range, sex, race/ethnicity, scanner information, motion parameters,\ncognitive performance, etc.) should be included."}, "author_response": "We sincerely thank the reviewers for their insightful comments and constructive suggestions.\nFollowing their comments and suggestions, we have revised our manuscript and carried out more\nanalyses. All the revisions and new results have been highlighted in the revised manuscript.\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nSpatial Heterogeneity and Subtypes of Functional Connectivity Development in Youth——\nComments\nLi and colleagues leverage large-sample fMRI data and machine learning-based brain age\nprediction to propose a regional brain development index (RBD) and identify three individual-\nlevel subtypes of functional connectivity (FC) over development. They further demonstrate how\nthese subtypes relate to cognition, cortical hierarchy, myelination, and gene\nexpression—linking multiple biological levels. Notably, the P2 subtype (characterized by advanced\nassociation cortex development) shows better cognitive performance across two independent cohorts,\nhighlighting the robustness and scientific value of their findings. This work pushes beyond\ntraditional global metrics (like brain age gap, BAG), offering the first demonstration of FC\ndevelopmental spatial heterogeneity and its multi-level implications at the individual level.\nDespite these strengths, I do have several comments and suggestions to improve the\nmanuscript:\n1. It seems that the authors do not adequately account for critical confounders such as\nsocioeconomic status (SES), parental education, or health conditions. The large-scale cohort\nstudies such as ABCD and IMAGEN routinely include these as covariates, since they are tightly\nlinked to brain structure, function, and cognitive development in youth (Jirsaraie et al., 2023 Nat\nCommun; Noble et al., 2015 Nature). I thus suggest that the authors include these variables as\ncovariates in their main or supplementary analyses, or at least stratifying their results\naccordingly. Doing so would substantially strengthen the credibility and scientific validity of your\nconclusions.\nResponse: We appreciate the reviewer’s insightful comments and constructive suggestions.\nFollowing the suggestions, we have performed additional analyses by including the SES and\nparental education as covariates. Specifically, neighborhood SES measure quantified by\ngeocoded variables as in previous studies [ref r1 and r2] was adopted for the PNC cohort,\nand SES measure quantified as the natural log of the income-to-needs ratio was adopted for\nthe HCP-D cohort as in [ref r3]. As illustrated in Fig. S5 below, the subtyping results after\ncontrolling for these variables remain consistent with our previous findings in terms of spatial\npatterns (Fig. S5A), cognitive profiles (Fig. S5 B and D), and individual-level subtype\nassignments (Fig. S5 C and E). These results have been included in the supplementary\nmaterial of the revised manuscript.\n[r1] Moore, T. M., et al. \"Characterizing social environment's association with neurocognition using\ncensus and crime data linked to the Philadelphia Neurodevelopmental Cohort.\" Psychological\nmedicine 46.3 (2016): 599-610.\n[r2] Murtha, Kristin, et al. \"Associations between neighborhood socioeconomic status, parental\neducation, and executive system activation in youth.\" Cerebral Cortex 33.4 (2023): 1058-1073.\n[r3] Weissman, David G., et al. \"Family income is not significantly associated with T1w/T2w ratio in the\nHuman Connectome Project in Development.\" Imaging Neuroscience 1 (2023): 1-10.\nFig. S5. Subtypes of functional connectivity development identified based on RBD index maps,\ncontrolling for SES and parental education. (A) Centroid RBD index maps of the three subgroups\nidentified in the PNC cohort. (B) Statistical comparisons of cognition scores across three subgroups\nin the PNC cohort. (C) Crosstab of subtype assignments before and after controlling for SES and\nparental education (subtyping solution 1 vs. subtyping solution 2) in the PNC cohort, demonstrating\nconsistency in subtype assignments. (D) Statistical comparisons of cognition scores across three\nsubgroups in the HCP-D cohort. (E) Crosstab of subtype assignments before and after controlling for\nSES and parental education (subtyping solution 1 vs. subtyping solution 2) in the HCP-D cohort.\n2. Although the three FC developmental subtypes differ in executive function, social cognition,\nand memory, the cognitive battery is relatively narrow and does not capture a broader spectrum\nof behavioral and psychological abilities, such as emotion processing, motivational regulation, or\ncreative thinking. Given the multidimensional nature of adolescent behavior and cognition—and\ntheir relevance to neurodevelopment (see Casey et al., 2018, Biol Psychiatry)—I encourage the\nauthors to more fully acknowledge this limitation and, where possible, connect with additional\nbehavioral or affective dimensions available in open datasets.\nResponse: Thanks for this thoughtful comment. We agree that incorporating additional\nbehavioral and psychological measures would provide a more comprehensive\nunderstanding of the relevance of FC development to youth behavior and cognition. We\nhave acknowledged this as a limitation and a direction for future work in the revised\nDiscussion section:\n“Fifth, the associations between FC development and cognition were investigated using a\nrelatively narrow cognition battery. Future studies incorporating a broader spectrum of\nbehavioral and psychological measures are needed to comprehensively capture the\nmultidimensional relevance of FC development to youth behavior and cognition.”\n3. The proposed three-subtype spatial clustering is mainly compared to the conventional whole-\nbrain BAG, but there is no direct benchmarking against other established spatial clustering\nmethods—such as structural brain age subtyping (Wang et al., 2019, Nat Commun), spectral\nclustering (Satterthwaite et al., 2015, Neuron), or GMM. Please consider providing a systematic\ncomparison of clustering performance (e.g., ARI, NMI, stability), cognitive associations, and\nspatial distribution overlap to clarify the innovation and boundaries of your approach. It would be\nhelpful if the authors can provide further justification for the chosen number of clusters (k),\nincluding cross-cohort reproducibility or biological interpretability.\nResponse: Thank you for the constructive suggestion. Imaging data used for subtyping are\ngenerally unlabeled, making it nontrivial to identify the optimal clustering algorithm solely\nbased on clustering performance, as standard clustering metrics are not necessarily aligned\nwith cognitive or behavior differences across subtypes. A prior study [ref r4] on MRI based\nsubtyping suggested that K-means and Gaussian Mixture Models (GMM) perform well in\nterms of stability and generalizability. Building on this study, we also adopted GMM to\nsubtype FC development patterns in our study.\nAs shown in Fig. S6, the subtyping results obtained using both K-means and GMM\nmethods are largely consistent in terms of spatial patterns (Fig. S6A), cognitive profiles\n(Fig. S6B), and individual-level subtype assignments (Fig. S6C), demonstrating the\nrobustness of the identified FC development patterns. These results have been included in\nthe supplementary material of the revised manuscript.\nIn addition to stability based (ARI) and data-fitting based (within-cluster distance) measures\n(Fig. S4A), we further justified the choice of number of clusters using the NbClust package\n[ref r5] and according to the data distribution characteristics. Particularly, NbClust uses a\nwide variety of indices to determine the optimal number of clusters using Hierarchical\nclustering or K-means methods. As shown in Fig. S4B, the most frequently selected cluster\nnumbers across clustering methods were 2 and 3, suggesting them as optimal choices for\nour dataset. To further verify this finding, we examined the distribution of individual RBD\npatterns in a 2D embedding space using t-SNE projection (Fig. S4C). It was observed that\nclusters P1 and P3 in the 3-cluster solution aligned closely with clusters P1 and P2 in the 2-\ncluster solution, while cluster P2 in the 3-cluster solution captured fine-grained inter-\nindividual variability between P1 and P3. Based on all these quantitative and qualitative\nassessments, the number of clusters was set to 3 in our analysis.\nMoreover, we replicated the clustering procedure in the HCP-D cohort and identified\nconsistent RBD subtypes (Fig. S10), further supporting the robustness and reproducibility of\nthe 3-cluster solution. All these results have been added to the revised manuscript and\nsupplemental materials.\n[r4] Zhang, Shirui, et al. \"A Comparative Analysis of Subtyping Methodologies on Cross-sectional\nsMRI Data.\" Neuroscience Bulletin (2025): 1-7.\n[r5] Charrad, Malika, et al. \"NbClust: an R package for determining the relevant number of clusters in\na data set.\" Journal of statistical software 61 (2014): 1-36.\nFig. S6. Subtypes of functional connectivity development identified based on RBD index maps using\nGMM. (A) Centroid RBD index maps of the three subgroups identified in the PNC cohort. (B)\nStatistical comparisons of cognition scores across the three subgroups. (C) Crosstab of subtype\nassignments by K-means (subtyping solution 1) and GMM (subtyping solution 2), demonstrating\nsubtype consistency between clustering methods.\nFig. S4. Determination of the number of subtypes based on clustering performance and data\ndistribution. (A) Reproducibility and data-fitting measures for FC development subtyping across\ndifferent numbers of subtypes (ranging from 2 to 7), evaluated over 100 split-half replications.\nReproducibility was assessed using the Adjusted Rand index, and data-fitting was quantified as\n1−𝑑𝑎𝑣𝑔, where 𝑑𝑎𝑣𝑔 represents the average of all data points' distances to their corresponding\ncluster centroids. (B) Optimal number of clusters suggested by various metrics in the NbClust\npackage, using both hierarchical clustering and K-means methods. Blue-colored blocks indicate the\nnumber of clusters recommended by each metric. (C) Distribution of individual FC development\npatterns projected into a 2D space using t-Distributed Stochastic Neighbor Embedding (t-SNE), with\noverlaid subtype assignments from both 2-cluster and 3-cluster solutions. Each dot represents one\nindividual. A crosstab of the 2-cluster and 3-cluster assignments is shown on the right.\nFig. S10. Subtypes of functional connectivity development identified based on RBD index maps in\nthe HCP-D cohort. (A) Centroid RBD index maps of the three subgroups identified. (B)\nStatistical comparisons of cognition scores across the three subgroups. (C) Crosstab of subtype\nassignments based on matching to the PNC subtype centroids (subtyping solution 1) and the\nHCP-D subtype centroids (subtyping solution 2), demonstrating subtype consistency between\ncohorts.\n4. While the manuscript connects gene expression, myelination, FC, and cognition, most\nanalyses are correlational. It is now well accepted that relationships across biological scales are\noften complex, with mediation, moderation, and feedback loops (see Fornito et al., 2015, Nat\nRev Neurosci). I would suggest explicitly discussing potential mechanistic links between the\nmicro, meso, and macro levels (e.g., via multimodal causal modeling, SEM, network mediation\nanalysis) and acknowledging the current limitations of causal inference.\nResponse: Thank you for this thoughtful comment. We agree that the current analyses\nlinking FC to gene expression and myelination are correlational in nature, and that causal\ninference is a promising direction for uncovering potential mechanistic links, particularly\nwhen multimodal data (e.g., transcriptome data and imaging data) are available from the\nsame cohort. We have acknowledged this limitation in the revised Discussion section:\n“Third, the associations between FC development patterns and gene expression data are\ncorrelational and rely on the Allen Human Brain Atlas, which represents transcriptomic data\nfrom adult brains. Replication studies using transcriptome data from youth populations are\nneeded to confirm and extend these findings. Moreover, future analysis using causal\ninference approaches may help elucidate mechanistic links between FC development,\nmolecular factor, and cognition.”\n5. The analyses rely on two large cross-sectional cohort studies (PNC and HCP-D), which are\nvaluable for describing group-level differences across age but cannot capture within-individual\ndevelopmental trajectories. The field broadly recognizes that longitudinal follow-up is necessary\nto truly understand developmental rate, sequence, and causality (see Kraemer et al., 2000, Am\nJ Psychiatry; Tamnes et al., 2017, Nat Neurosci). Please discuss this limitation, and emphasize\nthe need for future longitudinal studies (e.g., ABCD, IMAGEN) that can follow the same\nindividuals over time to more accurately characterize brain developmental heterogeneity and its\ncognitive and molecular correlates.\nResponse: Thank you for this thoughtful comment. We agree that longitudinal data is\nessential for capturing the associations between intra-individual changes in FC and\ncognition throughout development. We have acknowledged this as a limitation in the\nrevised Discussion section:\n“In addition, the current analysis was based on cross-sectional data to characterize the\ndevelopment differences across individuals. Future research incorporating longitudinal data\nis needed to investigate intra-individual changes in FC and cognition throughout\ndevelopment.”\n6. The introduction does not sufficiently relate classical neurodevelopmental theories (e.g.,\nFelleman & Van Essen, 1991, cortical hierarchy; Huttenlocher, 1979, synaptic pruning; Hebb,\n1949, plasticity) to the novel aspects of spatial subtyping. I would recommend more explicitly\nconnecting the study’s approach to these foundational models in both the introduction and\ndiscussion, which would help a broader audience appreciate the conceptual advance.\nResponse: Thank you for the constructive suggestion. We have elaborated this in the\nrevised manuscript. In the Introduction:\n“Research has established the critical role of cortical hierarchy in brain function, and\naccumulating evidence has demonstrated that the cortical development follows spatial and\ntemporal patterns aligned with this hierarchy from perspectives of macroscale\nneurobiological properties [refs r6 and r7] and underlying microscale mechanisms in\ndevelopmental plasticity, such as cortical excitation, inhibition, and glial cell types [refs r6,\nr8, r9]. Recent studies further show that the spatial-temporal patterns of functional network\ndevelopment in youth systematically align with cortical hierarchy organization, particularly\nalong the sensorimotor-association (S-A) axis, with developmental increases in functional\nconnectivity (FC) in sensorimotor cortices and decreases in association cortices [ref r10].\nHowever, how these FC development patterns vary across individuals remains unclear.\nNeurocognition undergoes a prolonged development trajectory parallel to the maturation of\nthe brain’s functional organization, with individual differences in youth linked to important\noutcomes such as academic performance and quality of life [ref r11], while deficits in\nneurocognition are associated with vulnerability to mental illness. Therefore, advancing our\nunderstanding of FC development at the individual level may provide valuable insights into\ntypical development processes and deviations from these typical trajectories, with\nimplications for cognitive and psychological functioning as well as mental illness.”\n[r6] Sydnor, Valerie J., et al. \"Neurodevelopment of the association cortices: Patterns, mechanisms,\nand implications for psychopathology.\" Neuron 109.18 (2021): 2820-2846.\n[r7] Felleman, Daniel J., and David C. Van Essen. \"Distributed hierarchical processing in the primate\ncerebral cortex.\" Cerebral cortex (New York, NY: 1991) 1.1 (1991): 1-47.\n[r8] Toyoizumi, Taro, et al. \"A theory of the transition to critical period plasticity: inhibition selectively\nsuppresses spontaneous activity.\" Neuron 80.1 (2013): 51-63.\n[r9] Larsen, Bart, and Beatriz Luna. \"Adolescence as a neurobiological critical period for the\ndevelopment of higher-order cognition.\" Neuroscience & Biobehavioral Reviews 94 (2018): 179-195.\n[r10] Luo, Audrey C., et al. \"Functional connectivity development along the sensorimotor-association\naxis enhances the cortical hierarchy.\" Nature communications 15.1 (2024): 3511.\n[r11] Keller, Arielle S., et al. \"Hierarchical functional system development supports executive\nfunction.\" Trends in cognitive sciences 27.2 (2023): 160-174.\n7. The results are organized as spatial variability/heterogeneity, subtyping, cognition/gene\nassociations, and replication. However, the scientific rationale and linkage between each step\ncould be more clearly articulated—for instance, why certain molecular analyses were performed\npost-subtyping, and how each step is motivated by prior findings or hypotheses. Explicitly\nstating the underlying logic for each major analysis will make the paper more cohesive and\naccessible.\nResponse: Thank you for the constructive suggestion. We have further elaborated the\nanalytical workflow in the revised manuscript.\n“We began by examining spatial variability in FC development and its inter-individual\nheterogeneity. Next, we applied clustering to identify covarying patterns of FC development\nacross cortical regions, stratifying individuals into subtypes characterized by distinct\nregional FC development profiles. We then compared the spatial patterns across subtypes\nand assessed their alignment with established principles of cortical hierarchy in brain\ndevelopment. To further characterize these subtypes, we analyzed their neurocognitive\nprofiles using cognitive measures and meta-analytic functional decoding. We subsequently\nexplored the biological underpinnings of these FC development patterns by examining their\nassociations with gene expression. Finally, we evaluated the generalizability and\nrobustness of the identified subtypes of FC development using data from an independent\ncohort.”\n8. Figure 1 shows the RBD analysis pipeline, but does not clearly label model inputs and\noutputs, the formula for RBD, or how RBD feeds into downstream clustering and association\nanalyses. Please consider annotating workflow nodes, core equations, and analytic steps to\nimprove interpretability, especially for interdisciplinary readers.\nResponse: Thank you for the constructive suggestion. We have updated Fig. 1 accordingly\nto better illustrate the analytical workflow of this study.\nFig. 1. Overview of the analysis workflow. (A) The regional brain development (RBD) index was\nderived from functional connectivity (FC). Region-specific brain age prediction models were built\nupon regional FC profiles comprising FC measures of each region with all other cortical regions, to\nassess region-wise FC development. These prediction models were trained using a 5-fold cross-\nvalidation approach to estimate brain age for each individual and cortical region. Spatial variability in\nprediction performance and inter-individual heterogeneity in the RBD index were then analyzed to\ngain insights into region-wise FC development. (B) Spatially covarying patterns of the RBD index\nwere identified through clustering, which stratified individuals into subtypes with distinct FC\ndevelopment profiles. These subtypes were further characterized by examining their associations\nwith cortical hierarchy, neurocognition measures, and gene expression.\n9. Currently, subtyping is based only on brain imaging. There is increasing interest in behavioral\n(e.g., LCA) or joint brain-behavior subtyping approaches (see Drysdale et al., 2017, Am J\nPsychiatry; Smith et al., 2022, Trends Cogn Sci). I encourage the authors to discuss prospects\nfor integrating behavioral subtyping with brain-based subtyping, and perhaps even attempt a\npreliminary behavioral subtyping analysis with available data to explore overlap and\ncomplementarity, which could inform future stratified intervention strategies.\nResponse: Thank you for this thoughtful comment. We have further discussed this as a\npromising future direction in the revised manuscript.\n“Lastly, while this study focuses on FC development in youth, a promising direction for\nfuture research involves integrating multi-modal data, including brain imaging, behavioral\nmeasures, and other relevant measures, to investigate the heterogeneity of brain\ndevelopment from a multi-dimensional perspective. Such approaches could offer a more\ncomprehensive understanding of the complex mechanisms underlying brain development\nand maturation and provide valuable insights into the convergence and divergence of\ndevelopmental heterogeneity revealed across different data modalities.”\n10. All major analyses are correlational, and there is no direct causal chain established between\nFC subtype, cognition, or gene expression. Please make a clear distinction between correlation\nand causality in the discussion and avoid over-interpretation. Where possible, supplement with\nsensitivity analyses—such as subgroup analysis by age, sex, or SES; mediation analysis or\nSEM to explore whether FC subtype mediates links between molecular/structural factors and\ncognitive outcomes. Also, highlight that future work leveraging longitudinal data, intervention, or\nmulti-omic integration will be needed to clarify causal mechanisms. Such enhancements would\ngreatly increase the scientific impact and robustness of findings.\nResponse: Thank you for this thoughtful comment. Discussion of causal inference,\nlongitudinal data analysis, multimodal data integration has been clearly acknowledged as\nlimitations and identified as directions for future work in the revised manuscript.\nReviewer #2 (Remarks to the Author):\nIn this manuscript, Li and colleagues examined spatial patterns of functional brain development\nacross two normative samples and characterized data-driven subgroups by their cognitive\nprofiles. The study is well motivated: deriving regional brain age estimates addresses an\nimportant limitation of traditional brain age studies that rely on a single global metric.\nMy primary concern relates to the authors' approach to replicating the data-driven analysis. The\nauthors applied k-means clustering to identify subgroups based on patterns of regional brain\ndevelopment (RBD) indices in their primary cohort (PNC) and compared cognitive scores across\nsubgroups. Because clustering will always produce a solution, independent replication is critical\nto demonstrate that subgroups are meaningful rather than artifacts. To address this, the authors\nclassified individuals from their replication cohort (HCP-D) using PNC-derived centroids and\nconcluded that subgroup replicability was supported by similar cognitive differences across\ndatasets. However, this approach does not validate replicability.\nTrue replication would require rerunning clustering independently in both cohorts and\ndemonstrating that (1) subgroup patterns of RBD indices and (2) associated cognitive profiles\nare similar. Without establishing (1), it remains unclear whether comparable subgroup structures\nexist. Classifying HCP-D participants based on PNC centroids assumes, rather than tests,\nsubgroup generalizability and risks overestimating replicability. For example, it is theoretically\npossible that multiple, drastically different RBD patterns could yield similar distributions of\ncognitive scores without reflecting the same underlying neurodevelopmental patterns. To\nconvincingly demonstrate replicability, the authors should rerun clustering independently in the\nHCP-D cohort and compare both subgroup structures and cognitive profiles to those observed\nin PNC.\nSeparately, the cognitive measures differ substantially between the two datasets (executive\nfunctioning, social cognition, and memory in PNC, versus fluid cognition, crystallized cognition,\nand total cognition in HCP-D). As a result, even if subgroup differences were consistent in\ndirection or magnitude, it would be difficult to interpret this as evidence of replication, given that\nthe constructs are not directly comparable. To strengthen their conclusions, the authors should\nharmonize cognitive domains across cohorts or more explicitly acknowledge this limitation.\nResponse: We appreciate the reviewer’s insightful comments and constructive\nsuggestions. We agree that stratifying individuals in the HCP-D individuals using the PNC-\nderived centroids primarily serves as a generalization analysis. As a large sample of\nparticipants would have to complete both sets of cognitive tests to allow the mapping\nbetween cohorts, harmonizing cognitive measures might be challenging and would require\nnew data collection. To better align cognitive domains between the two cohorts, we\nexamined additional cognitive measures available in the HCP-D cohort. For example, the\nFlanker test is most relevant to the executive function, and the List Sort task is most\nrelevant to the memory function. As shown in the updated Fig. 5, individuals in subtype P2\ndemonstrated higher cognitive performance than those in P1 and P3, particularly on the\nFlanker test (p=0.004) and the List sort task (p=0.032), further supporting the subtype-\nrelated differences in cognitive profiles. We have added a discussion noting that future\nvalidation using harmonized cognitive measures or additional datasets would be beneficial.\nTo further assess the robustness of the identified FC development patterns, we replicated\nthe clustering procedure within the HCP-D cohort. The resulting subtypes showed\nconsistent spatial patterns with the PNC-derived centroids (Fig. S10A), and the cognitive\nprofiles (Fig. S10B) across subtypes were similar to those observed in the PNC-based\nanalysis (Fig. 5D). Individual-level subtype assignments also demonstrated good\nconsistency (Fig. S10C), reinforcing the reproducibility of the finding. These results have\nbeen included in the revised manuscript and supplementary material.\nFig. 5. Generalization analysis in the HCP-Development cohort. (A) Regional prediction accuracy for\nbrain age, measured by correlation coefficient (r). (B) Spatial correspondence between prediction\naccuracy in the PNC and HCP-D cohorts, revealing consistent pattens of regional FC development\n(p<0.0001, Spin-test). (C) Assignment of HCP-D individuals to the FC development subtypes defined\nin the PNC cohort. For each HCP-D individual, the distances between their RBD index map and the\ncentroid RBD index maps of the PNC-defined subtypes were computed, and the individual was\nassigned to the subtype with the smallest distance. For example, individuals i and i+1 were assigned\nto subtypes P2 and P1, respectively. (D) Statistical comparisons of cognition scores across the three\nsubtypes in HCP-D, highlighting subtype-related differences in cognitive profiles. CardSort,\nDimensional Change Card Sort Task; Flanker, Flanker Test; ListSort, List Sort Working Memory\nTask; Reading, Oral Reading Test; Pattern, Pattern Comparison Processing Speed Task; PicSeq,\nPicture Sequence Memory Task; PicVocab, Picture Vocabulary Test.\nFig. S10. Subtypes of functional connectivity development identified based on RBD index in the\nHCP-D cohort through clustering. (A) Centroid RBD index maps of the three subgroups identified. (B)\nStatistical comparisons of cognition scores across the three subgroups. (C) Crosstab of subtype\nassignments based on matching to the PNC-derived subtype centroids (subtyping solution 1) and the\nHCP-D subtype centroids (subtyping solution 2), demonstrating subtype consistency between\ncohorts.\nAdditional comments for improving the manuscript:\n• The use of cross-sectional data to study development is a major limitation that is not discussed\nin the manuscript. This should be acknowledged.\nResponse: Thank you for this thoughtful comment. We agree that longitudinal data is\nindispensable to capture the associations of intra-individual FC and cognition changes\nalong development. We have acknowledged this as a limitation in the revised Discussion\nsection:\n“In addition, the current analysis was based on cross-sectional data to characterize the\ndevelopment differences across individuals. Future research incorporating longitudinal data\nis needed to investigate intra-individual changes in FC and cognition throughout\ndevelopment.”\n• Normative modeling, which similarly derives individualized regional deviations from expected\ntrajectories (and can incorporate non-linearity and covariates, stated limitations of the chosen\napproach), is not mentioned. While I do not expect a change in approach, a brief rationale for\nchoosing RBD over normative modeling would be helpful.\nResponse: Thank you for this thoughtful comment. While normative modeling is\nincreasingly used to derive deviation scores from normative trajectories for individuals, it is\ntypically based on element-wise regression, which poses challenges for modeling FC due\nto its high dimensionality. Prior studies applying normative modeling to FC have either used\nsummary measures such as FC strength [ref r12] or surrogate scores derived from\ndimension reduction techniques [ref r13], which may not fully capture the richness of FC\nprofiles. Although multivariate normative modeling methods exist, their application to FC\ndata remains relatively limited.\nIn contrast, the predictive modeling of brain age provides an interpretable, composite score\nthat reflects the development status of each cortical region based on its full FC profiles,\nwithout relying on FC strength or dimension reduction. Moreover, the calibration of\npredicted brain age against chronological age in the pipeline is an analogy to normative\nmodeling, enabling the estimation of individual deviations in brain development.\nWe have added a discussion regarding this in the revised manuscript.\n[ref r12] Sun, Xiaoyi, et al. \"Mapping neurophysiological subtypes of major depressive disorder using\nnormative models of the functional connectome.\" Biological Psychiatry 94.12 (2023): 936-947.\n[ref r13] Kessler, Daniel, Michael Angstadt, and Chandra Sripada. \"Growth charting of brain\nconnectivity networks and the identification of attention impairment in youth.\" JAMA psychiatry 73.5\n(2016): 481-489.\n• The choice of the Schaefer 400 parcellation should be justified, particularly given evidence\nsuggesting that the 200-parcel resolution produces more representative connectivity networks\n(Luppi et al., 2021). Additionally, the omission of the typically paired Melbourne subcortical atlas\nis a missed opportunity to capture subcortical contributions.\nResponse: Thank you for this thoughtful comment. While the optimal resolution of\nparcellation may vary depending on the application, a recent study [ref r14], which\nsystematically evaluated the fMRI data-processing pipelines, identified nine optimal\npipelines, four of which used a 400-parcel resolution and three used a 200-parcel\nresolution. To assess the robustness of our findings, we conducted a sensitivity analysis\nusing the Schaefer-300 and Schaefer-200 atlases. As shown in Fig. S7, the subtypes\nidentified using the 300-parcel and 200-parcel atlases exhibited spatial patterns of FC\ndevelopment consistent with those derived from the 400-parcel atlas. Additionally,\nindividuals in subgroup P2 consistently demonstrated better cognitive performance across\natlases. These results suggest that our analysis is robust across commonly used cortical\nparcellations ranging from 200 to 400 parcels.\nWe agree that extending the analysis to include subcortical regions is an important future\ndirection, as it would enable investigation of the subcortical-cortical FC development and its\nrelationship to neurocognition. This has been acknowledged as a promising direction for\nfuture work in the revised Discussion section.\n[ref r14] Luppi, Andrea I., et al. \"Systematic evaluation of fMRI data-processing pipelines for\nconsistent functional connectomics.\" Nature Communications 15.1 (2024): 4745.\nFig. S7. Subtypes of functional connectivity development identified based on RBD index using\ndifferent cortical atlases. (A) Centroid RBD index maps of the three subgroups identified using the\nSchaefer-300 atlas. (B) Statistical comparisons of cognition scores across the three subgroups\nidentified using the Schaefer-300 atlas. (C) Centroid RBD index maps of the three subgroups\nidentified using the Schaefer-200 atlas. (D) Statistical comparisons of cognition scores across the\nthree subgroups identified using the Schaefer-200 atlas.\n• Raw correlation values were used as inputs to the ridge regression model. Typically, Fisher r-\nto-z transformation is applied to stabilize variance and normalize correlations; the rationale for\ndeviating from this standard practice should be provided.\nResponse: Thank you for this thoughtful comment. While both raw correlation [refs r15-r17]\nand Fisher’s r-to-z-transformed correlation values are commonly used as FC measures in\npredictive modeling, we conducted an additional sensitivity analysis using the Fisher’s r-to-\nz-transformed correlation values as model input. As illustrated in Fig. S8, the subtyping\nresults obtained using the Fisher’s r-to-z-transformed correlation values were largely\nconsistent with those derived using raw correlation values, in terms of spatial patterns (Fig.\nS8A), cognitive profiles (Fig. S8B), and individual-level subtype assignments (Fig. S8C).\nThese results further support the robustness of the identified FC development patterns.\nThese results have been included in the supplementary material of the revised manuscript.\n[ref r15] Chen, Jianzhong, et al. \"Shared and unique brain network features predict cognitive,\npersonality, and mental health scores in the ABCD study.\" Nature communications 13.1 (2022):\n2217.\n[ref r16] Dhamala, Elvisha, et al. \"Functional brain networks are associated with both sex and gender\nin children.\" Science Advances 10.28 (2024): eadn4202.\n[ref r17] Ooi, Leon Qi Rong, et al. \"Longer scans boost prediction and cut costs in brain-wide\nassociation studies.\" Nature (2025): 1-10.\nFig. S8. Subtypes of functional connectivity development identified based on RBD index maps\nderived using Fisher’s r-to-z-transformed FC measures. (A) Centroid RBD index maps of the three\nsubgroups identified. (B) Statistical comparisons of cognition scores across the three subgroups. (C)\nCrosstab of subtype assignments using raw correlation values (subtyping solution 1) and Fisher’s r-\nto-z-transformed correlation values as FC measures (subtyping solution 2), demonstrating subtype\nconsistency across FC measures.\n• I appreciate the authors' approach to determining the optimal number of clusters. However, it\nwould strengthen the manuscript to include more standard clustering validity metrics (e.g.,\nSilhouette coefficient) in the supplemental material for comparison.\nResponse: Thank you for the constructive suggestion. In addition to stability based (ARI)\nand data-fitting based (within-cluster distance) measures (Fig. S4A), we further justified the\nchoice of number of clusters using the NbClust package [ref r5] and according to the data\ndistribution characteristics. Particularly, NbClust uses a wide variety of indices to determine\nthe optimal number of clusters using Hierarchical clustering or K-means methods. As\nshown in Fig. S4B, the most frequently selected cluster numbers across clustering methods\nwere 2 and 3, suggesting them as optimal choices for our dataset. To further verify this\nfinding, we examined the distribution of individual RBD patterns in a 2D embedding space\nusing t-SNE projection (Fig. S4C). It was observed that clusters P1 and P3 in the 3-cluster\nsolution aligned closely with clusters P1 and P2 in the 2-cluster solution, while cluster P2 in\nthe 3-cluster solution captured fine-grained inter-individual variability between P1 and P3.\nBased on all these quantitative and qualitative assessments, the number of clusters was set\nto 3 in our analysis.\nMoreover, we replicated the clustering procedure in the HCP-D cohort and identified\nconsistent RBD subtypes (Fig. S10), further supporting the robustness and reproducibility of\nthe 3-cluster solution. All these results have been added to the revised manuscript and\nsupplemental materials.\n[ref r5] Charrad, Malika, et al. \"NbClust: an R package for determining the relevant number of clusters in\na data set.\" Journal of statistical software 61 (2014): 1-36.\nFig. S4. Determination of the number of subtypes based on clustering performance and data\ndistribution. (A) Reproducibility and data-fitting measures for FC development subtyping across\ndifferent numbers of subtypes (ranging from 2 to 7), evaluated over 100 split-half replications.\nReproducibility was assessed using the Adjusted Rand index, and data-fitting was quantified as\n1−𝑑𝑎𝑣𝑔, where 𝑑𝑎𝑣𝑔 represents the average of all data points' distances to their corresponding\ncluster centroids. (B) Optimal number of clusters suggested by various metrics in the NbClust\npackage, using both hierarchical clustering and K-means methods. Blue-colored blocks indicate the\nnumber of clusters recommended by each metric. (C) Distribution of individual FC development\npatterns projected into a 2D space using t-Distributed Stochastic Neighbor Embedding (t-SNE), with\noverlaid subtype assignments from both 2-cluster and 3-cluster solutions. Each dot represents one\nindividual. A crosstab of the 2-cluster and 3-cluster assignments is shown on the right.\n• More detail is needed on the clustering methodology: Why was k-means chosen over\nalternative algorithms (e.g., spectral clustering)? Why was the cosine distance metric used?\nResponse: Thank you for the constructive suggestion. Imaging data used for subtyping are\ngenerally unlabeled, making it nontrivial to identify the optimal clustering algorithm solely\nbased on clustering performance, as standard clustering metrics are not necessarily aligned\nwith cognitive or behavior differences across subtypes. A prior study [ref r4] on MRI based\nsubtyping suggested that K-means and Gaussian Mixture Models (GMM) perform well in\nterms of stability and generalizability. Building on this study, we also adopted GMM to\nsubtype FC development patterns in our study. As shown in Fig. S6, the subtyping results\nobtained using both K-means and GMM methods are largely consistent in terms of spatial\npatterns (Fig. S6A), cognitive profiles (Fig. S6B), and individual-level subtype assignments\n(Fig. S6C), demonstrating the robustness of the identified FC development patterns. These\nresults have been included in the supplementary material of the revised manuscript\n[r4] Zhang, Shirui, et al. \"A Comparative Analysis of Subtyping Methodologies on Cross-sectional\nsMRI Data.\" Neuroscience Bulletin (2025): 1-7.\nFig. S6. Subtypes of functional connectivity development identified based on RBD index maps using\nGMM. (A) Centroid RBD index maps of the three subgroups identified in the PNC cohort. (B)\nStatistical comparisons of cognition scores across the three subgroups. (C) Crosstab of subtype\nassignments by K-means (subtyping solution 1) and GMM (subtyping solution 2), demonstrating\nsubtype consistency between clustering methods.\n• Complete statistical reporting is needed throughout the manuscript, including test statistics,\nexact p-values, and effect sizes. Currently, many results are only partially reported (e.g.,\nindicating p < 0.05 without further detail).\nResponse: Thank you for this suggestion. We have included complete test statistics in the\nrevised manuscript.\n• The manuscript should specify that Python and scikit-learn were used for the ridge regression\nmodeling, including software versions.\nResponse: Thank you for this suggestion. We have included more detailed information\nabout the code used for the analysis.\n“The prediction models were implemented using the scikit-learn library (v1.4.1) in Python\n(v3.9.18).”\n• A table summarizing cohort characteristics (age, age range, sex, race/ethnicity, scanner\ninformation, motion parameters, cognitive performance, etc.) should be included.\nResponse: Thank you for the helpful suggestion. We have summarized these\ncharacteristics in Table 1 and Table 2 in the revised manuscript.\nReviewer #2 (Remarks on code availability):\nOther than the code being pretty sparsely commented, the README file was clear and I see no\nissues with the pipeline.\nResponse: Thank you for your observation. We have updated the README file and are\nworking to provide more detailed instructions for data preparation and code execution at\nhttps://github.com/MLDataAnalytics/FC-SpatialHeterogeneity.\nWe sincerely thank the reviewers for their thoughtful feedback and valuable time."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed most of my concerns. However, there are still some minor points needed to be addressed before publication to improve the quality of the manuscript:1. Figure2A: the fancy 3-D visualization of the dimension space actually deliver no useful information. I would personally recommend three scatter plots (Dim1-on-Dim2, Dim2-on-Dim3, Dim1-on-Dim3) to elaborate the distribution of subtypes in the space. It should be very easy to do some web search and learn about the 3-d visualization’s shortfall in scientific paper. 2. The Y-axis of Figure3B is missing.", "reviewer_2": "I thank the authors for their thoughtful revisions and the additional analyses conducted in response to prior feedback. These additions have substantially strengthened the manuscript. The consistent results across multiple rigorous sensitivity analyses further support the robustness of the findings. I have no further concerns."}, "author_response": "Thank you for these suggestions.\n1. We think this comment refers to Fig. 3A, as it mentions ‘the distribution of subtypes’. To\nclarify, the 3D visualizations depict the centroid RBD (regional brain development) index\nmaps for each subgroup, illustrating the typical spatial pattern of functional connectivity\ndevelopment for each subgroup, rather than a three-dimension distribution across\ncortical regions. These visualizations are intended to clearly show subgroup-specific FC\ndevelopment patterns and subgroup differences from both lateral and medial views.\n2. We have reviewed the figures and confirm that all relevant information is included.\nReviewer #2 (Remarks to the Author):\nI thank the authors for their thoughtful revisions and the additional analyses conducted in\nresponse to prior feedback. These additions have substantially strengthened the manuscript.\nThe consistent results across multiple rigorous sensitivity analyses further support the\nrobustness of the findings. I have no further concerns.\nResponse: We greatly appreciate your thoughtful comments and constructive suggestions,\nwhich have substantially strengthened our manuscript."}]

Decision letter