Introduction

What are the core principles of human brain organisation? As we develop, when are those principles established? And how do those principles vary from person to person, according to differences in cognition and behaviour? These are central questions in developmental and systems neuroscience. At a molecular level, early developmental patterning of the cortex arises through interacting gradients of morphogens and transcription factors (see Cadwell et al., 2019). The resultant areal and progenitor specialisation produces a diverse pool of neurones, glia, and astrocytes (Hawrylycz et al., 2015). Across childhood, an initial burst in neuronal proliferation is met with later protracted synaptic pruning (Bethlehem et al., 2022), the dynamics of which are governed by an interplay between experience-dependent synaptic plasticity and genomic control (Gottlieb, 2007). An emerging approach is to conceptualise brain organisation through gradients, revealing latent dimensions (see Huntenburg et al., 2018). This provides an empirical extension to seminal work that demonstrated a sensory-fugal axis of structural and functional brain organisation (Mesulam, 1998). However, most studies have focussed on gradients of organisation in adulthood and at a group-level. Our aim was to better understand how those organisational principles emerge throughout childhood and adolescence, and to map variability in that organisation to key individual differences in cognition and psychopathology.

With the advent and rapid proliferation of neuroimaging technologies, it is now possible to image the developing brain at scale. White matter axonal inter-regional projections form a structural connectome, which reconfigures during childhood and adolescence towards a small-world topology (Khundrakpam et al., 2013), and optimises a wiring cost-trade-off (Betzel et al., 2016; Bullmore & Sporns, 2009). Highly intercorrelated structural cores, primarily within the posterior medial cortex, provide a structural backbone (Hagmann et al., 2008) for functional connectivity (FC). Like SC, FC also reconfigures during childhood and adolescence, characterised by increasing segregation between default-mode, frontoparietal and cingulo-opercular networks, alongside integration between dorsal anterior cingulate cortices and the cingulo-opercular network (Fair et al., 2007, 2008), producing a more distributed topology (Fair et al., 2009). In parallel, anatomical constraints on functional networks weaken across development, with strengthening of long-range functional connections between anatomically-distant regions, and maintenance of short-range connectivity (Fair et al., 2009).

The trends described above reflect group-level developmental trends, but how do we capture these broad anatomical and functional organisational principles at the level of an individual? Gradients of organisation may be derived through non-linear dimensionality reduction techniques, such as diffusion map embedding (DME, Coifman & Lafon, 2006), providing a summary metric of organisational variability. The result is a set of eigenvectors, within a low-dimensional coordinate system, representing organisational axes upon which inter-regional structural and functional connections sit. Several studies have now used DME to derive organisational gradients of structural and functional connectivity (reviewed in Royer et al., 2024). In typical development, structural gradients derived from microstructural profile covariance (Paquola et al., 2019; Valk et al., 2022), histology (Paquola et al., 2019), and probabilistic tractography (Park et al., 2021) extend from primary sensorimotor regions towards limbic and paralimbic regions. Functional gradients in childhood are centred in unimodal cortices, but are then anchored in transmodal cortices at the onset of puberty (Dong et al., 2021), along which the default-mode network is organised (Margulies et al., 2016). Functional constraints on structural manifolds follow an anterior-posterior axis, strongest in highly-differentiated unimodal regions and weakest in undifferentiated transmodal regions (Paquola et al., 2019; Valk et al., 2022). Structure-function correspondence develops hierarchically across this axis, such that regions with the strongest structural differentiation have the weakest correspondence to FC across childhood and adolescence (Park et al., 2022).

The increasing availability of large-scale multi-modal neuroimaging databases, in both neurotypical and neurodivergent samples, creates the opportunity to comprehensively quantify canonical trajectories and sensitive periods of variability in cortex-wide coordinate systems of structural and functional gradients (see Huntenburg et al., 2018), with sufficient statistical power to make conclusions across development and phenotypes. Applying diffusion-map embedding as an unsupervised machine-learning technique, we derived gradients of structural and functional brain organisation in children and adolescents, spanning the full neurotypical-neurodivergent range. Crucially, we employ a mixed cross-sectional and longitudinal design, leveraging 887 structural and 727 functional scans, across a broad developmental window (6-19 years old), allowing us to robustly chart developmental trajectories. First, we derived structural and functional axes of organisation for each child. We anticipated a primary unimodal-transmodal structural gradient dominant across development, accompanied by a primary unimodal-centred functional gradient in early childhood which is replaced by a unimodal-transmodal functional gradient at the onset of adolescence (Dong et al., 2021; Xia et al., 2022). Next, as an exploratory analysis, we charted the development of individual variability in these gradients, and statistically quantified sensitive periods. Third, we examined the patterning of the intersection of structural and functional gradients, termed structure-function coupling, across the cortex, anticipating a unimodal-transmodal organisational axis (Baum et al., 2020; Suárez et al., 2020; Vázquez-Rodríguez et al., 2019). Fourth, we mapped structure-function coupling with multiple dimensions of cognition and neurodevelopmental symptomatology. Based on evidence suggesting that transmodal association cortices exhibit protracted development (Sydnor et al., 2021) and therefore greater developmental sensitivity, with the greatest rate of evolutionary expansion and variability (Buckner & Krienen, 2013; Mueller et al., 2013), we hypothesised that regions of the cortex with the weakest structure-function coupling, namely higher-order transmodal association cortices, would be the largest source of individual differences in psychopathology and cognition.

Results

We leveraged diffusion-weighted imaging (DWI) and resting-state functional magnetic resonance imaging (rsfMRI) scans from two mixed-design neurodevelopmental datasets with overlapping age ranges, common cognitive and psychopathology measures, but different recruitment approaches. The first was the Nathan Kline Institute (NKI) Rockland Sample Longitudinal Discovery of Brain Development Trajectories sub-study. This is a community-ascertained sample of neurotypical children aged between 6 and 17 years old, from Rockland County, USA. We processed 447 DWI scans from 258 participants (56.82% male, 28.35% with one scan, 33.48% with two, and 38.17% with three) and 424 rsfMRI scans from 256 participants (54.35% male, 31.29% with one scan, 37.65% with two, and 31.06% with three). The second dataset was the Centre for Attention, Learning, and Memory (CALM), based in Cambridge (UK). This comprised children and adolescents aged between 6 and 19 years old. Unlike the NKI sample, the CALM sample is intentionally designed to capture the breadth of cognitive profiles present in children who present to children’s professional clinical and educational services. All children in CALM were referred from educational or clinical professional services because of difficulties with cognitive, learning or behaviour. From CALM we processed 440 DWI scans from 352 participants (70% male, 60% with one scan, 40% with two) and 304 rsfMRI scans from 256 participants (66.45% male, 68.42% with one scan, 31.58% with two). See the corresponding participants section in the Methods for full sample details.

Using the Schaefer 200-node 7-network parcellation, we reconstructed individual-level and group-level structural and functional connectomes (Figure 1a), using best practices and several sensitivity analyses to ensure robustness (see MRI and rsfMRI processing sub-sections in the Methods). Structural connectomes were constructed using streamline counts from probabilistic tractography. As diffusion-map embedding requires a smooth fully-connected matrix, and in preparation for further downstream analyses examining coupling of structural and functional gradients, we transformed streamline counts into communicability matrices. These communicability matrices reflect a graph theory measure of information transfer previously shown to maximally predict functional connectivity (Esfahlani et al., 2022; Seguin et al., 2022). Functional connectomes were reconstructed from rsfMRI blood oxygenation-level dependent (BOLD) response amplitude timeseries for each region correlated with that of all other regions. We calculated the normalised angle of each structural and functional connectome to derive symmetric affinity matrices. We applied diffusion-map embedding to these affinity matrices (Figure 1b), as a non-linear dimensionality reduction technique to obtain eigenvectors (‘axes’) along which maximal variance in structural and functional organisation, respectively, is accounted for. See ‘statistical analyses and data processing’ within the Methods section for full details.

Figure 1.

Universal patterning of communicability and functional connectivity gradients across phenotypes

We first examined how group-level gradients of structural communicability and rsfMRI differed between neurotypical (NKI) and neurodivergent (CALM) samples (Figure 1c). Averaged across hemispheres, the first three functional connectivity gradients accounted for 42.11%, 36.16%, and 21.73% variance in NKI, respectively, and a similar proportion in CALM (G1_FC = 42.32%, G2_FC = 32.08%, G3_FC = 25.60%). In both NKI and CALM, we replicated canonical functional gradients (Margulies et al., 2016). That is, the principal functional gradient (G1_FC) extended from bilateral somato-motor regions towards the prefrontal cortex. Spin tests preserving spatial autocorrelation revealed that the first functional gradients in CALM and NKI were strongly and significantly correlated (r_s = .880, p_spin < .001), with regions of maximal divergence between gradients, quantified as the largest absolute difference in eigenvectors, centred within bilateral somatomotor and visual networks. The second functional gradient (G2_FC) extended from bilateral visual regions towards the prefrontal cortex and frontal operculum and insula. It displayed high spatial correspondence between CALM and NKI (r_s = .666, p_spin = .001), with the largest between-group divergence centred within bilateral default mode and control networks. Finally, the third functional gradient (G3_FC) extended from bilateral dorsal attention networks towards somato-motor regions and displayed strong spatial correspondence between groups (r_s = .753, p_spin < .001), with the largest between-group divergences centred within the parietal cortex, precuneus, and posterior cingulate cortex divisions of the right default mode network, alongside the right somatomotor network.

Further, just like the functional gradients, gradients of structural communicability were also organised similarly across datasets. Averaged across hemispheres, the first three communicability gradients accounted for 36.21%, 33.40%, and 30.39% variance in NKI, respectively, and a similar proportion in CALM (G1_SC = 36.44%, G2_SC = 32.84%, G3_SC = 30.72%). In NKI, G1_SC was anchored at one end by visual regions in the right hemisphere, and by lateral and dorsal-medial prefrontal cortices at the other end. In CALM, the same gradient was anchored between bilateral prefrontal cortices and right visual networks and displayed strong significant spatial correspondence to NKI (r_s = .830, p_spin < .001). The largest between-group divergence occurred within the left orbitofrontal cortex and precuneus divisions of the control network, and within the precuneus and prefrontal cortex divisions of the left default-mode network. The second communicability gradient (G2_SC) in NKI extended from visual regions and the precuneus in the left hemisphere towards regions in the dorsal attention network and temporo-occipital-parietal lobes in the right hemisphere. The corresponding gradient in CALM displayed high spatial correspondence (r_s = .848, p_spin < .001), with the largest divergence between gradients centred in the left prefrontal cortex and control networks. Finally, the third communicability gradient (G3_SC) in NKI spanned prefrontal, and somato-motor cortices, and was statistically significantly and strongly correlated with the corresponding CALM gradient (r_s = .719, p_spin< .001), which was anchored at one end by the right temporal cortex of the default mode network, and at the other end by divisions of the left default mode, control, and ventral attention systems.

The largest between-group divergence occurred primarily within the parietal and temporal cortex divisions of the left default-mode network. Together, gradients of functional connectivity and communicability displayed similar organisational principles across datasets and phenotypes.

Despite strong spatial correspondence of gradients between datasets, alignment was not perfect. We investigated the hypothesis that this ‘misalignment’ was driven by technical site variation or recruitment protocol differences. Helpfully, the CALM study also incorporated a smaller comparison sample, comprising age-matched children recruited from the same schools and neighbourhoods as those referred to the main cohort. This extra non-referred sample provides an opportunity to disentangle these explanations, since they were scanned on the same machine and protocol as the main CALM cohort. We derived group-level communicability and functional connectivity gradients across 91 structural and 79 structural scans, respectively, from CALM comparison children. Gradients in the non-referred sample were almost identical to the referred sample, both in terms of communicability (G1_SC r_s = .999; G2_SC r_s = .991; G3_SC r_s = .995, all p_spin < .001) and functional connectivity (G1_FC r_s = .983; G2_FC r_s = .986; G3_FC r_s = .968, all p_spin < .001). In other words, any small differences in alignment between NKI and CALM are unlikely to be the result of the phenotypic differences between the cohorts. Across datasets with varying phenotypes and scanning protocols, at least at a group-level averaged across development during childhood and adolescence, the same broad organisational principles of communicability and functional connectivity emerge.

Stability of individual-level gradients across developmental time

After establishing group-level gradients, we applied diffusion-map embedding to 887 individual structural and 728 functional connectomes, across data sets and timepoints. Within modalities and components, gradients explained a similar proportion of variance across datasets (Figure 1d). The median (range) variance explained by the first 3 structural components were 35.75% (34.55% - 38%), 33.61% (32.44%-35%), and 30.67% (27.8%-32%), respectively, in NKI, and 36.50% (35.22% - 38%), 32.94% (31.83% - 35%), and 30.52% (28.03% - 32%) in CALM. For the functional gradients the median (range) variance explained by the first 3 components were 41.12% (35.61%-53%), 33.61% (32.44%-35%), and 30.67% (27.8%-32%), respectively, in NKI, and 41.60% (36.11% - 51%), 32.77% (25.95%-38%), and 25.40% (20.07% - 31%) in CALM. In other words, across the two datasets there is a remarkable consistency in the ordering and the variance explained at an individual level, for both structural and functional gradients.

But does the ordering of these gradients change with developmental time? Dong and colleagues (2021) demonstrated that the primary (unimodal) functional gradient accounted for 38% variance in FC in childhood but only 11% in adolescence, whilst the primary (association) functional gradient accounted for 38% variance in adolescence but only 12% in childhood. We were unable to replicate this effect (Figure 1e). Instead, across datasets, we observed that the proportion of variance explained by the principal structural and functional gradients, respectively, remained highly stable over time, without any sudden change in variance explained at the onset of adolescence. Put simply, not only is the ordering and the patterning of the gradients highly similar across both datasets, the structural and functional gradients of brain organisation were established early in life and are highly stable across childhood and adolescence.

Manifold eccentricity as a measure of individual variability in gradient organisation

To examine variability in axis organisation across time, we used manifold eccentricity as a dependent variable (Figure 2a). For each modality and dataset, we computed the group manifold origin as the mean of the first three gradients. For each participant’s 3D embedding in manifold space, we calculated the Euclidean distance between each node and the group manifold origin. This measure has previously been reported in structural connectivity (Park et al., 2021), and a similar measure in functional connectivity (Bethlehem et al., 2020). To explore properties of manifold eccentricity as a function of data set and modality, we calculated the nodal coefficient of variation (Figure 2b), averaged across all participants. This was non-normally distributed across data sets and modalities, as evidenced by a non-significant Shapiro-Wilk test (W = .92, p < .001). Therefore, we conducted a Scheirer-Ray-Hare test (Scheirer et al., 1976) as a non-parametric equivalent to a two-way Analysis of Variance. The Scheirer-Ray-Hare test revealed a significant main effect of data set (F(1,176) = 87.31, p < .001) and modality (F(1,196) = 3687.40, p < .001), but a non-significant interaction between dataset and modality (F(1,796) = 2.72, p = .090). A post-hoc Dunn test with Benjamini-Hochberg correction for multiple comparisons revealed that NKI functional manifold eccentricity had the largest coefficient of variation (Mean_CV = .340), and this was significantly greater than that of CALM (Mean_CV = .304, p_FDR = .001). Further, the mean coefficient of variation for NKI structural manifold eccentricity (Mean_CV = .139) was also significantly greater than that of CALM (Mean_CV = .114, p_FDR < .001). Together, this suggests that the children referred by clinical or educational professionals showed a reduced variance in manifold organisation, and this is particularly pronounced in functional connectivity.

Figure 2.

To better understand what manifold eccentricity measures in terms of connectome organisation, we investigated its relationship with nodal graph metrics (Figure 2c). We included two graph theory metrics previously associated with structural manifold eccentricity in a neurotypical sample of children and adolescents: participation coefficient, describing the uniformity of dispersion of edges amongst different modules (Guimerà & Nunes Amaral, 2005), and within-module degree Z-score, a standardised measure of the number of intra-modular connections. Across datasets and modalities, higher manifold eccentricity was significantly positively correlated with within-module degree Z-scores, but negatively correlated with participation coefficient. The consistency across datasets and modalities suggests that higher manifold eccentricity corresponds to greater segregation.

Developmental trajectories of structural manifold contraction and functional manifold expansion

Next we charted age-related changes in organisation of structural and functional manifolds, operationalised as manifold eccentricity, using generalised additive mixed models (GAMMs). We did this globally as well as for 7 intrinsic connectivity networks (ICNs) (Yeo et al., 2011). Within each GAMM, age and an age-dataset interaction were smooth covariates, whilst mean framewise displacement, sex, and dataset were parametric linear covariates. Dataset was a significant predictor of structural manifold eccentricity across all levels of analysis (Figure 3a), with the largest cohort effects globally (t = −74.351, p_FDR < .001) and within the default-mode network (t = −42.476, p_FDR < .001), driven by higher manifold eccentricity in CALM than NKI. Further, remodelling of structural manifolds was consistently and significantly associated with development globally and across all networks, except the somato-motor network. Structural manifold eccentricity decreased with age, indicating greater integration. The strongest developmental effects were observed globally (F = 52.644, p_FDR <.001), alongside limbic (F = 29.486, p_FDR < .001) and dorsal attention (F = 27.148, p_FDR < .001) networks. The non-linear GAMM procedure also produces a measure of the non-linearity of an effect, or its effective degrees of freedom (EDF), whereby an EDF of 1 indicates a linear effect, with larger deviations from 1 indicating greater non-linearity. This metric revealed considerable differences in the nature of the developmental effect on structural manifold eccentricity, ranging from linear globally and within the frontal-parietal network (EDF = 1.000), to non-linear within the visual (EDF = 2.434) and dorsal attention (EDF = 2.056) networks.

Figure 3.

In line with previous work (Sydnor et al., 2023; Tervo-Clemmens et al., 2023; Tooley et al., 2023), but not yet explored in terms of multi-modal developmental neuroimaging data, we precisely quantified sensitive periods of development by calculating the first derivative of the smooth age term for each data set separately, based on the original GAMM specification and data used for model fitting subset by dataset. We then extracted developmental periods in which the simultaneous confidence intervals excluded zero (p < .05, two-sided), computed for each unique age in each modality and data set. Developmental effects on structural manifold eccentricity were statistically significant across the entire developmental age range for both CALM and NKI globally, alongside within dorsal attention, ventral attention, and limbic networks. Within the visual network, developmental effects were statistically significant until mid-adolescence, in both CALM (6.00 – 14.58 years) and NKI (6.68 – 15.27 years). Within the frontoparietal network, age-related change was significant across the entire developmental period for CALM (6.00 – 19.17 years), but not statistically significant for NKI. A similar picture emerged for the default-mode network, where age-related change was statistically significant for the entire developmental range in CALM (6.00 – 19.17 years), but only for a brief adolescent period within NKI (14.04 – 14.59 years). This suggests that within higher-order frontoparietal and default-mode networks, structural manifold remodelling could be reliably predicted by age within children referred for problems with attention, learning, or memory, and this was distinct from the developmental trajectory followed by children within a broader community sample. Finally, we examined interactions between age and dataset. Such interactions were significant predictors of structural manifold eccentricity globally (F = 11.665, p_FDR = .001), alongside default-mode (F = 30.074, p_FDR < .001) and limbic (F = 1.504, p_FDR = .046) networks. Across all three levels of analysis, effects were driven by a stronger age effect within CALM (F_Global = 37.589, p < .001; F_Default = 27.655, p < .001; F_Limbic = 34.867, p < .001) than NKI (F_Global = 18.045, p < .001; F_Default = 1.305, p = .222; F_Limbic = 14.863, p < .001). In summary, the structural manifold is gradually contracting with developmental time, most markedly in the limbic and dorsal attention systems. Manifold eccentricity is higher in neurodivergent than neurotypical participants, but differential age effects suggest that the ongoing contraction is more marked in neurodivergent participants.

Having demonstrated significant expansion of structural manifolds and sensitivity to development and dataset, we repeated the same modelling procedure with functional manifold eccentricity as the outcome (Figure 3b). Dataset was a significant predictor of functional manifold contraction globally (t = −6.975, p_FDR < .001, Mean_CALM = .108, Mean_NKI = .104), and in visual (t = −4.066, p_FDR < .001, Mean_CALM = .127, Mean_NKI = .120), dorsal attention (t = −4.903, p_FDR < .001, Mean_CALM = .103, Mean_NKI = .097), ventral attention (t = −2.906, p_FDR = .015, Mean_CALM = .100, Mean_NKI = .098), fronto-parietal (t = −4.712, p_FDR < .001, Mean_CALM = .100, Mean_NKI = .096), and default-mode (t = −8.885, p_FDR < .001, Mean_CALM =.111, Mean_NKI = .101) networks. Across these partitions, CALM had significantly higher functional manifold eccentricity, albeit to a small degree, with the largest cohort effects globally and within the default-mode network, indicating greater segregation. Whilst structural manifolds were sensitive to development across almost all levels of analysis, we only observed significant age-related functional manifold expansion globally (F = 13.369, p_FDR = .002) and within the fronto-parietal network (F = 14.627, p_FDR= .001). Within both levels of analysis, statistically significant first derivatives of age revealed significant age-related change across the whole developmental period in CALM (6.17 – 19.17 years), but no significant change in NKI. Finally, we observed significant age-dataset interactions globally (F = 7.971, p_FDR = .018) and within the fronto-parietal network (F = 9.816, p_FDR = .008). In both cases, CALM exhibited larger age effects (F_Global = 6.276, p = .013; F_{Fronto-parietal} = 8.567, p = .004) than NKI (F_Global = 2.934, p = .087; F_{Fronto-parietal} = 2.122, p = .202). Together, we demonstrate that structural manifolds contract across childhood and adolescence, indicating greater integration whilst, in parallel, functional manifolds expand, indicating greater segregation. Further, structural manifold expansion shows more widespread effects of development and cohort than functional manifold contraction. Full GAMM statistics are available in Supplementary Table 1.

Structure-function coupling follows a unimodal-transmodal axis, with developmental effects concentrated within higher-order association networks

Thus far, we tracked the organisation of structural and functional manifolds separately. However, extensive evidence points towards structural connectivity providing a back-bone for functional brain organisation (Hagmann et al., 2008; Honey et al., 2007, 2009). Therefore, we turned our attention towards structure-function coupling. Structure-function coupling was operationalised as the Spearman-rank correlation between Euclidean distances between each node and all others within structural and functional manifolds, respectively, repeated for each region and participant (Figure 4a). Thus, a large structure-function coupling value suggests that a node occupies a similar location in relation to all other nodes, within structural and functional manifold space, respectively. Further, a negative coupling value suggests that a node’s position in structural manifold space, relative to all other nodes, is anticorrelated to its relative position in functional manifold space. 199 CALM participants had coupling values for one timepoint, with 84 additional coupling values for a second timepoint. 121 NKI participants had a baseline coupling value, with 126 and 99 additional points for a second and third timepoint, respectively. In line with prior studies (Baum et al., 2020; Suárez et al., 2020; Vázquez-Rodríguez et al., 2019), across data sets, we observed a unimodal-transmodal axis of structure-function coupling magnitude and inter-individual variability (Figure 4b). The smallest and least variable coupling was centred around the unimodal anchor, including bilateral somatomotor and visual networks, whilst the largest and most variable coupling was centred around the transmodal anchor, including bilateral default-mode and limbic networks. Further, stratifying structure-function coupling by network and dataset (Figure 4c) revealed consistent coupling across datasets (see Supplementary Table 2).

Figure 4.

Next, we examined the developmental trajectories of structure-function coupling within the GAMM framework, with the same model specification as before (Figure 4d). Dataset was a statistically significant predictor of structure-function coupling within the fronto-parietal network only (t = - 4.177, p_FDR < .001), with stronger coupling within NKI, whilst age was not a statistically significant predictor globally or across any network partition. However, an interaction between age and dataset was a significant predictor of coupling within dorsal attention (F = 9.579, p_FDR = .008) and default-mode (F = 7.245, p_FDR = .007) networks. Within the dorsal attention network, both datasets exhibited a linear age effect (EDF = 1), with a stronger age effect from NKI (F = 9.948, p = .002) than CALM (F = .484, p = .487), and spanning the entire developmental period examined (6.90 – 17.96 years). Similarly, within the default-mode network, the age effect in NKI (F = 5.110, p = .011) was stronger than CALM (F = 2.267, p = .198). The age effect in NKI was highly non-linear (EDF = 2.341), with a later onset period of statistically significant development (14.39 – 17.96 years) than CALM (11.75 – 14.58 years). Together, these effects demonstrate that the development of structure-function coupling is somewhat dataset-dependent, especially within the higher-order association networks, with developmental change being reduced in the neurodivergent sample.

Structure-function coupling variability predicts dimensions of cognition, but not psychopathology

Thus far, we have charted the canonical trajectories of structural and functional manifold reorganisation across childhood and adolescence and examined structure-function relationships. However, what is the behavioural significance, if any, of differences in structure-function coupling? The final step of our analysis related structure-function coupling to dimensions of psychopathology and cognition. To this end, we took advantage of the full range of behavioural and cognitive heterogeneity offered by combining NKI and CALM. To derive orthogonal dimensions of psychopathology (Figure 5a), we conducted a principal component analysis (PCA) with varimax rotation on 6 age-standardised t-score measures from the Conners III Questionnaire (Conners et al., 2011), measuring inattention, hyperactivity/impulsivity, learning problems, executive functioning, aggression, and peer relations. Component 1 explained 44.73% variance, and captured inattention, learning problems, and executive functioning. Component 2 explained 23.09% variance and was dominated by aggression. Component 3 explained 18.23% variance and was dominated by peer relationships. A series of two-tailed continuity-corrected Mann-Whitney U tests revealed significantly larger loadings across CALM than NKI for psychopathology component 1 (U = 87320, p < .001), component 2 (U = 58213, p < .001), and component 3 (U = 63832, p < .001). Two dimensions of cognition were derived using the same method (Figure 5b), based on 6 z-score age-standardised broad measures of cognition: forward and backward digit recall from the Automated Working Memory Assessment (Alloway, 2007a), alongside Tower achievement score, visual scanning speed, number-letter switching accuracy, and motor speed from the Delis-Kaplan Executive Function System (Delis et al., 2001). The first component explained 32.30% of variance in the cognitive measures. Of the 4 measures which most strongly loaded onto the first component, visual scanning, number-letter switching, and motor speed loaded positively, whilst the Tower task loaded negatively. The second component explained 29.63% of variance, onto which forward and backward digit recall loaded most strongly, suggesting that this component captures working memory. Like psychopathology, the difference between median loadings for participants from CALM and NKI for cognitive components 1 (U = 65086, p < .001) and 2 (U = 96732 p < .001) were statistically significant. Together, this suggests that CALM and NKI differed significantly in terms of psychopathology and cognitive profiles (PCA loadings are reported in Supplementary Table 4).

Figure 5.

We then examined the extent to which dimensions of psychopathology predicted structure-function coupling, using the GAMM framework. Alongside sex, mean framewise displacement, and age as covariates, we additionally included main effects of each psychopathology dimension and interactions with age to parse development-specific effects. Neither the main effects of psychopathology dimension or their interaction with age were statistically significant predictors of structure-function coupling globally, or across any network partition, following multiple comparison correction. However, when we applied the same modelling framework to two dimensions of cognition, we found that higher loadings on the second cognitive dimension, capturing working memory, significantly and linearly predicted stronger structure-function coupling within the default-mode network (EDF = 1.000, F = 13.896, p_FDR= .006). Furthermore, we found several significant interactions between age and this second cognitive factor (working memory). This interaction predicted structure-function coupling globally (F = 9.659, p_FDR = .014), as well as within somato-motor (F = 8.199, p_FDR = .024), dorsal attention (F = 7.175, p_FDR = .039), and default-mode networks (F = 12.425, p_FDR = .006). To parse these interactions, we ran a median split by cognitive ability, and repeated the GAMM procedure within each half (Figure 5c). For global structure-function coupling, we observed a stronger age effect in participants with lower (F = 9.206, p = .003), rather than higher (F = 2.009, p = .157), cognitive dimension scores. On the other hand, for somatomotor structure-function coupling, we observed a stronger developmental effect in participants with higher (F = 4.865, p = .033) working memory dimension scores, rather than lower (F = 3.160, p = .076). The difference in developmental effects on working memory loadings as a predictor of structure-function coupling was perhaps most pronounced within the dorsal attention network, with the age effect on structure-function coupling within the high ability group was over ten-fold (F = 12.138, p = .001) that of the low ability group (F = .130, p = .719). Finally, whilst age tended to be a significant predictor of global, dorsal attention, and somatomotor network structure-function coupling within participants with either high or low scores on the second cognitive dimension, developmental effects were statistically significant for both low (F = 6.548, p = .012) and high (F = 8.084, p = .005) loading groups when predicting default-mode network coupling. In summary, cognitive dimensions, specifically working memory, predict stronger structure-function coupling across several network partitions in a developmentally-specific manner, with larger between-dataset differences earlier in development. Those with the highest cognitive scores exhibited a gradual reduction in coupling, whilst those with the lowest cognitive scores exhibited a gradual increase in coupling. Further, those with higher working memory scores tend to show the strongest developmental changes in coupling, with higher coupling earlier in development, but weaker towards the end of adolescence, and more marked in the dorsal attention network.

Discussion

Using diffusion-map embedding we established highly stable principles of structural and functional brain organisation, with the ordering of key gradients being entirely insensitive to age. Furthermore, group-level structural and functional gradients were highly consistent across datasets – differences between datasets only emerged when taking development into account, through differing rates of structural manifold contraction and functional manifold expansion, respectively. Similarities in group-level gradients and their developmental trajectories suggests that such gradients are established early in life, refined through development, and may be examples of universal organisational principles of cortical development. Finally, we quantified the intersection between these two forms of low-dimensional embedding, generating a metric of the coupling between structural and functional manifolds. This coupling displayed a characteristic unimodal-transmodal organisational axis and was highly consistent across datasets. Interactions between age and the second cognitive dimension, capturing working memory, significantly predicted structure-function coupling globally, alongside within somatomotor, dorsal attention, and default-mode networks, with the largest and most consistent age effects concentrated within higher-order late-developing association cortices.

First, we demonstrated the stable ordering of structural and functional gradients across childhood and adolescence, and across samples making up a neurodivergent-neurotypical continuum. That this consistent pattern emerges despite considerable differences in phenotype suggests that there is a degree of species universality in these organisational principles. The stable functional gradient ordering is in contrast to a small body of developmental work suggesting that the primary functional gradient in childhood differentiates between the somato-motor and auditory cortices from the visual cortex, and is replaced by a gradient anchored in unimodal and heteromodal association cortices during late childhood (Xia et al., 2022) and early adolescence (Dong et al., 2021). Despite the growing evidence in support of different primary gradients for childhood and adolescence, inconsistencies exist in the literature. For example, Dong and colleagues (2021) reported that the principal functional gradient in childhood accounted for 38% variance, with the second accounting for 11%. Such a stark difference between the principal and secondary functional gradient is in contrast to both the current study and prior work (Xia et al., 2022).

The gradual functional manifold expansion we observed, indicating functional segregation, builds upon previous reports of cortex-wide functional segregation across youth (Tooley et al., 2022), but at the level of overarching organisational gradients, rather than individual graph theory measures, providing a more unifying account of cortical reorganisation. In parallel, the gradual structural manifold contraction, indicating structural integration, supports existing generative models of structural brain development predicting network integration across youth (Akarca et al., 2023). However, the developmental changes in structural manifold eccentricity are not uniformly distributed across cortex or participants; the developmental effects on structural manifold contraction were concentrated in higher-order association regions, with significantly greater manifold eccentricity, and thus segregation, in neurodivergent participants, who also exhibited stronger age effects. Tentatively, this raises the possibility that, relative to neurotypical participants, neurodiversity may be associated with greater sensitivity to developmental change in structural gradients. Globally and across intrinsic connectivity networks, a greater proportion of structural manifold contraction effects were sensitive to phenotype and age, compared to functional manifold expansion. This may reflect greater variability in functional, rather than structural, brain development, and therefore reorganisation of the low-dimensional structural manifold space may be a more consistent and reliable marker of development and phenotype.

In line with prior work (Baum et al., 2020; Gu et al., 2021; Park et al., 2022; Preti & Van De Ville, 2019; Suárez et al., 2020), an axis of structure-function coupling emerged across datasets, strongest in unimodal lower-order regions and weakest in heteromodal higher-order regions. This sensorimotor-associative axis of structure-function coupling follows patterns of evolutionary expansion in the cortex, such that regions with the smallest evolutionary expansion and highest cross-species similarity have functionally-specific roles, anchored in unimodal regions, whilst regions with the largest evolutionary expansion and smallest cross-species similarity have functionally-diverse roles, anchored in association cortices (Sydnor et al., 2021; Xu et al., 2020). Cross-species tritiated thymidine experiments (Finlay & Darlington, 1995) have demonstrated that brain regions with disproportionately large volumes, such as association cortices, experience protracted neurogenesis. Speculatively, this may be associated with greater sensitivity to age, which we demonstrated, relative to lower-order networks, as significant interactions between age and data set for structure-function coupling in higher-order networks, particularly those incorporating the frontal lobes. Further, whilst the neurodivergent dataset displayed stronger developmental effects in structural manifold remodelling, we observed stronger developmental effects in structure-function coupling within the neurotypical sample, concentrated in higher-order networks. Speculatively, this may point towards protracted structural manifold development underpinning more flexible structure-function coupling, which is more sensitive to age and associated with neurotypicality.

Additionally, dimensions of cognitive ability, specifically working memory, were robust predictors of stronger structure-function coupling within higher-order association regions. This supports the tethering hypothesis of brain development (Buckner & Krienen, 2013), which posits that association cortices are untethered by, and therefore less susceptible to, the intrinsic morphogenic signalling patterns and extrinsic activity from sensory systems that otherwise constrain development of evolutionarily-ancient regions, such as sensory cortices. This also supports existing structure-function coupling studies which demonstrated that variability in structure-function coupling derived from an N-back working memory task follows a unimodal-transmodal axis, and exhibits developmental effects principally concentrated within higher-order association areas (Baum et al., 2020). The present work expands these findings by testing associations between structure-function coupling and broad dimensions of psychopathology and cognition in neurotypical and neurodivergent participants, across datasets. Across development, we also observed that stronger structure-function coupling within the default-mode network was associated with significantly higher loadings for the second cognitive component, but that this effect was development-dependent. This might suggest that better performance on working memory tasks is supported by stronger alignment between low-dimensional structural and functional axes: this effect is strongest early in development for high-achievers, and strongest in later adolescent development for low-achievers.

Our results have several limitations and avenues for further research. First, longitudinal data in the current study forms a small proportion of the overall sample, resulting in statistical power insufficient to detect smaller effects. This may be particularly evident in the case of functional manifold reorganisation, which has greater variability than structural manifold reorganisation. Second, an important future direction is exploring generative, rather than descriptive, models of structure-function relationships. For example, neuromorphic engineering approaches, endowing neural networks with biological properties of connectomes, have been used to explain how function may emerge from structural constraints (Achterberg et al., 2023; Suárez et al., 2021). For example, in the context of reservoir computing, memory capacity in a working memory task has been shown to be significantly greater in reservoirs constrained by empirical brain networks, compared to rewired connectomes (Suárez et al., 2021). Further, in the context of recurrent neural networks, manipulating task performance, signal communication, and economic wiring constraints has been shown to converge on a narrow parameter window in which both structural and functional brain organisation features emerge (Achterberg et al., 2023). Future research may adapt such techniques to include gradients as major organisational features when evaluating model performance. Note, however, that eigenmodes from cortical geometry, rather than brain connectivity, have been suggested as a simpler and more intuitive explanation of structural constraints on functional connectivity, compared to alternative measures relying on connectomes, such as connectome-based eigenmodes and neural mass models (Pang et al., 2023). Therefore, particularly in relation to structure-function coupling, further research charting changes in geometric manifolds, rather than connectomics, may yield more representative measures of structure-function coupling. Third, the current study lays a foundation for exploring potential biological mechanisms mediating the structure-function relationships. For example, open-access libraries of neurotransmitter receptor densities derived from positron-emission topography (Hansen et al., 2022) may lay the foundation for examining how excitation-inhibition balances along the structural connectome backbone constrains functional connectivity dynamics, both at rest and during tasks.

In summary, by triangulating clinical and community developmental samples, we robustly charted the development of structural and functional manifold reorganisation across childhood and adolescence, utilising 887 structural and 728 functional scans, respectively, by modelling the linear and non-linear functional forms of age. We suggest that childhood and adolescence is characterised by consistently-ordered structural and functional gradients established early in life, refined through development with the greatest change in higher-order association networks, and whose coupling reflects dimensions of cognition, specifically working memory. Our findings also highlight the importance of studying development in examining individual differences in brain organisation across time, particularly during youth, and highlights equifinality. That is, whilst group-level trends were very similar across phenotypes, the trajectories to reach such organisation varied.

Methods

Participants

To evaluate how gradients varied across time and neurodevelopmental profiles, we used two mixed-design samples with overlapping age ranges. The first recruiting neurodivergent children, and the second recruiting neurotypical children (6.4 – 19 years old). The first data set was the Centre for Attention, Learning, and Memory (CALM). Native English-speaking children, referred by educational and clinical professionals, with difficulties in attention, language, and/or memory were included regardless of diagnosis, medication, and health conditions, except with uncorrected sensory impairments or health conditions affecting cognition (see Holmes et al., 2019 for a detailed protocol). Our final sample consisted of participants referred to the service due to difficulties with attention, learning, or memory, and with T1w and rsfMRI or DTI data across at least one time point, excluding repeated scans. Thus, we processed 440 DWI scans from CALM for 352 participants (70% male), where 60% had one scan (Mean_Age = 10.72 ± 1.74 years), and the remainder had two (Mean_Age = 13.76 ± 2.26 years). Further, we processed 304 rsfMRI scans from 256 participants (66.45% male), where 68.42% had one scan (Mean_Age = 10.72 ± 1.81 years) and the remainder had two (Mean_Age = 13.84 ± 2.24 years).

The second sample was the Longitudinal Discovery of Brain Development Trajectories sub-study of the NKI-Rockland Sample (see Tobe et al., 2022 for detailed protocols). Rather than recruiting neurodivergent children, this study contained a community-ascertained neurotypical sample of 339 children (55.16% male) aged between 6 and 17 years old at recruitment from Rockland, Orange, Bergen, and Westchester Counties in the United States. Baseline recruitment started in December 2013, with follow-ups at 12- or 15-month intervals. 136 children had a single session, 76 had one follow-up, and 127 had two follow-ups. In contrast to referral by health and education professionals in CALM, NKI children were referred to the service through word of mouth, flyers, prior engagement with NKI events and education days. Exclusion criteria included a history of autism spectrum disorder, psychiatric hospitalizations, an intelligence quotient below 70, and prior treatment with antidepressants, neuroleptics, and mood stabilisers (see Tobe et al., 2022 for further details). Across two 6-hour days, children completed a broad cognitive battery, provided biological samples for genetic analysis, and multi-modal neuroimaging. After September 2015, this protocol was amended to one 8-hour day. The study was approved by the NKI for Psychiatric Research Institutional Review Board. Written consent was obtained by the caregiver, and then the participant if they returned for follow-up after their 18^th birthday. We processed 448 DWI scans from 259 NKI participants (56.92% male), where 28.35% had one scan (Mean_Age = 11.94 ± 3.04 years), 33.48% had two (Mean_Age = 12.20 ± 2.69 years), and the remainder had three (Mean_Age = 12.46 ± 2.52 years). Further, we processed 425 rsfMRI scans from 257 NKI participants (54.35% male), where 31.29% had one scan (Mean_Age = 12.28 ± 2.97 years), 37.65% had two (Mean_Age = 12.34 ± 2.82 years), and 31.06% had three (Mean_Age = 12.94 ± 2.76 years).

Pre-processing of Behavioural and Cognitive Data

We pre-processed behavioural and cognitive data in R v4.2.2 using RStudio. Using predictive mean matching in the mice v3.16 package (van Buuren & Groothuis-Oudshoorn, 2011), for baseline CALM with neuroimaging data (N = 443), we imputed missing gender (N = 2) and age-standardised scores for object counting (N = 7) from the Phonological Assessment Battery (PhAB; Frederickson et al., 1997), alongside digit recall (N = 4) and backward digit span (N = 5) from the Automated Working Memory Assessment (AWMA; Alloway, 2007a). For longitudinal CALM with neuroimaging data (N = 151), we imputed age standardised PhAB alliteration (N = 2) and object counting (N = 2) scores, alongside AWMA digit recall (N = 3) scores.

MRI Data Acquisition

CALM

Diffusion tensor imaging (DTI) and rsfMRI data were used in the present study, with T1-weighted (T1w) and T2-weighted (T2w) scans for registration. All modalities were collected using a 3T Siemens Prisma scanner with 32-channel head coil, within a 1-hour session at the MRC Cognition and Brain Sciences Unit, Cambridge, United Kingdom (see Holmes et al., 2019 for additional details). Baseline scans had all 4 scan types, whilst T2w scans were not collected in the longitudinal sample. Further, scanner software was upgraded from VD13D to VE11E in March 2021, after longitudinal data collection commenced. Children were first familiarised with a mock scanner and played an interactive game designed to minimise head motion. Across 4 minutes and 32 seconds, 192 slices (1mm isotropic resolution with 256mm field of view (FOV)) of T1w structural data were collected using a 3D magnetization prepared rapid-echo sequence (MP-RAGE), with repetition time (TR) of 2500ms, echo time (TE) of 3.02 seconds, inversion time (TI) of 900ms, 9° flip angle, partial Fourier, and generalized auto-calibrating partial parallel acquisition (GRAPPA) acceleration factor 4. At baseline, T2w data was collected using a 3D sampling perfection with application optimised contrast using different flip angle evolution turbo spin-echo sequence. 29 slices (.6875mm x .6875mm x 5.2mm voxel resolution) were collected for 1 minute and 38 seconds, with TR of 5060ms, TE of 102.9ms, and GRAPPA acceleration factor 4. DTI data of resolution 2mm isotropic voxels and 192mm FOV was acquired using 64 gradient diffusion directions of 1000s/mm² and 4 interleaved reference b-value = 0s/mm² images, with TR of 8500ms, TE of 90ms, GRAPPA acceleration factor 2, 192mm FOV read, and partial Fourier. Whilst participants lay with their eyes closed, 270 volumes with 32 interleaved axial slices of T2* gradient-echo echo planar imaging (EPI) was collected across 9 minutes and 6 seconds, with TR of 2000ms, TE of 30ms, 78° flip angle, 3mm isotropic resolution, and 192mm FOV.

NKI-Rockland Sample Longitudinal Discovery of Brain Development Trajectories

All modalities were collected using a Siemens Magnetom Trio-Tim (syngo MRB17) 3T scanner with 32-channel head coil (see Tobe et al., 2022) at the Nathan Kline Institute for Psychiatric Research, New York, United States. The neuroimaging protocol consisted of a T1w MP-RAGE, DTI, and 6 fMRI sequences of varying TRs, three of which were resting state. Across 4 minutes and 18 seconds, 176 slices of T1w structural data were collected using a single-shot 3D MP-RAGE sequence and GRAPPA acceleration factor of 2 with TR of 1900ms, TE of 2.52ms, TI of 900ms, 9° flip angle, and 250mm FOV. Whilst participants listed to music, 64 slices of interleaved DTI data of 2mm isotropic resolution were collected across 5 minutes and 43 seconds using 128 b = 1500s/mm² diffusion directions and 9 b=0s/mm² images, with TR of 2400ms, TE of 85ms, 90° flip angle, multi-band (MB) acceleration factor of 4, partial Fourier, and 212mm FOV. Whilst participants lay still in silence, across 9 minutes and 45 seconds, 64 slices of an interleaved multi-slice EPI sequence with 2mm isotropic resolution were acquired with TE of 30ms, 65° flip angle, partial Fourier, and MB acceleration factor 4. To maximise spatial, rather than temporal, resolution, we selected the 1400ms TR sequence for our analyses, instead of the 2500ms TR sequence.

MRI Data Processing

Across modalities and data sets, processing pipelines were identical. The following sections contain information from QSIprep (Cieslak et al., 2021) and fMRIprep (Esteban et al., 2019) boilerplates, respectively.

Structural MRI

We processed structural and diffusion-weighted data using QSIprep 0.14.2 (Cieslak et al., 2021), implemented through Nipype 1.6.1 (Gorgolewski et al., 2011), Nilearn 0.8.0 (Abraham et al., 2014) and Dipy 1.4.1 (Garyfallidis et al., 2014). In brief, using Advanced Normalization Tools 2.3.1 (ANTs; Avants et al., 2008), the T1w image was corrected for intensity non-uniformity using the N4 algorithm, skull-stripped with a target from the Open Access Series of Imaging Studies, and non-linearly registered to the ICBM 152 Nonlinear Asymmetrical 2009c template (Fonov et al., 2011). Finally, the image was segmented into white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF) using FMRIB’s Automated Segmentation Tool (FAST; Zhang et al., 2001).

DWI

Using MRtrix3 (Tournier et al., 2019), DWI scans were denoised using Marchenko-Pastur Principal Component Analysis, corrected for field inhomogeneity using the N4 algorithm (Tournier et al., 2019), and adjusted to the mean intensity of b = 0s/mm² images. Corrections for head motion and Eddy currents were implemented using FSL 6.0.3 (Andersson & Sotiropoulos, 2016), and motion-associated outliers replaced (Andersson et al., 2016). A rigid ANTs registration created transformations from the b = 0s/mm² to the T1w image. In preparation for constructing affinity matrices, which requires fully connected networks, we conducted whole-brain probabilistic tractography with 10 million streamlines ranging between 30 and 250mm in length, and fibre orientation density (FOD) power of .33. Whilst probabilistic tractograms are considerably denser than deterministic, they suffer from high false-positive rates, and inaccurate reconstruction for crossing and fanning fibres (Maier-Hein et al., 2017). Therefore, we implemented anatomically constrained tractography (Smith et al., 2012), which uses biologically-informed priors to terminate streamlines, such as when entering cortical or subcortical GM, or CSF. To further improve the estimate accuracy of the contribution of WM, GM, and CSF FODs to the DWI signal, single-shell multi-tissue constrained spherical deconvolution was implemented in the MRtrix3Tissue (https://3Tissue.github.io) fork of MRtrix3 (Tournier et al., 2019), and tissue components normalised against inhomogeneity using a log-domain version of an algorithm developed by Raffelt and collleagues (2017). Finally, spherical-deconvolution informed filtering of tractograms (SIFT-2; Smith et al., 2015) ensured that the density of the reconstructed streamlines more accurately mirrored that of the underlying fibre density.

We constructed structural connectomes using the SIFT2-weighted streamline counts, parcellated into 200 cortical regions, each assigned to one of 7 intrinsic functional connectivity networks (Schaefer et al., 2018). We used this parcellation because it matches a previous similar study assessing manifold expansion in adolescents (Park et al., 2021), and provides sufficient spatial resolution to clearly visualise the separation of the extremes in each DME-derived axis. For each participant, we retained the strongest 10% of connections per row, thus creating fully connected networks required for building affinity matrices. To generate a group-representative connectome, we concatenated connectomes across time points, and conducted distance-dependent consensus-thresholding designed to retain the distribution of short-range connections (Betzel et al., 2019). This produced a binary mask, from which we extracted mean streamline counts. Streamlines in the brain provide tracts through which information is transferred between regions through synaptic transmission. One measure of information transfer is communicability, where the communicability of node i is the weighted sum of all possible paths between a pair of nodes passing through it, with shortest paths weighted most strongly (Estrada & Hatano, 2008). Therefore, to better reflect the utility of these tracts in terms of network functioning, and to increase the predictive validity of SC for FC when examining structure-function relationships (Seguin et al., 2020), we converted the thresholded connectomes into weighted communicability matrices (Crofts & Higham, 2009).

Resting-State Functional Image Processing

We processed resting-state functional MRI data using fMRIprep 21.0.1 (Esteban et al., 2019) through Nilearn 0.8.1 (Abraham et al., 2014). In brief, each BOLD run was slice-time corrected using the 3dtshift tool from the Analysis of Functional Images software, head motion parameters estimated using FSL’s mcflirt, resampled to original space, and non-linearly co-registered to the anatomical reference using boundary-based registration with 6 degrees of freedom. The resultant signal was smoothed with an isotropic Gaussian kernel of 6mm full-width half-maximum, and then non-aggressively denoised using independent component analysis-based Automatic Removal of Motion Artifacts (ICA-AROMA; Pruim et al., 2015), during which non-steady states were automatically removed. Total WM, GM, and CSF signals were extracted. Much debate exists around fMRI denoising strategies. However, we selected ICA-AROMA due to prior evidence (Parkes et al., 2018) that it performs well against other denoising strategies, such as regressing head motion parameters, anatomical or temporal derivatives of noise components, in terms of minimising the dependence between distance and FC, the difference in mean signal between high- and low-motion participants, temporal degrees of freedom, and test-retest reliability. Note that whilst these researchers demonstrated that censoring marginally outperformed ICA-AROMA, this decreased temporal degrees of freedom and required at least 4 minutes of uncensored data for each participant. Based on this evidence, we simultaneously regressed out global signal, WM and CSF signals from ICA-AROMA denoised data using Nilearn 0.6.0 (Abraham et al., 2014).

To construct functional connectomes, we averaged across the BOLD time-series, conducted pair-wise Pearson correlations between each pair of ROIs from the Schaefer 200-node 7-network parcellation, and z-transformed the edge weights. Whilst prior DME work retained the top 10% of connections per row (Margulies et al., 2016), this was in vertex-wise data, rather than parcellated. Therefore, to ensure that retained connections were all positive and non-zero, as required for an affinity matrix for diffusion-map embedding, we retained the 10% absolute strongest connections per row. The individual thresholded connectomes were averaged to produce a group-representative connectome.

DWI Sample Construction

Across all modalities, we included participants with low in-scanner motion (described below) and whose connectomes were successfully thresholded across 3 parcellations: the Brainnetome 246-node parcellation (Fan et al., 2016), alongside 100-node and 200-node Schaefer 7-network parcellations (Schaefer et al., 2018).

CALM

At baseline, 411 participants had DWI scans and anatomical data (T1w/T2w) for registration, from which we reconstructed structural connectomes for 408. 1 participant was removed due to high in-scanner motion (mean FD > 3 mm; Power et al., 2012), producing a final sample of 407 participants (Mean_FD = .47 mm, SD_FD = .40 mm). At follow-up, 131 participants had DWI and T1w scans, from which we reconstructed structural connectomes for 129, all with low in-scanner motion (Mean_FD = .33 mm, SD_FD = .20 mm).

NKI

Out of 268 participants with a single scanning session, we reconstructed structural connectomes for 239, and retained 223 across all parcellations (Mean_FD = .74, SD_FD = .34). For those with two scanning sessions, we reconstructed structural connectomes for 141 of 155, all with low in-scanner motion (Mean_FD = .84 mm, SD_FD = .40 mm). For those with three scanning sessions, we reconstructed connectomes for 90 of 104, 2 of which excluded due to motion, resulting in 84 participants (Mean_FD = .98 mm, SD_FD = .42 mm).

rsfMRI Sample Construction

Across both datasets, we retained functional connectomes for participants who exhibited low in-scanner motion (mean FD < .5mm), and less than 20% of spikes exceeding .5mm FD.

CALM

Out of 373 baseline participants with rsfMRI BOLD scans and anatomical data, we reconstructed functional connectomes for 372, and retained 259 after motion quality-control (Mean_FD = .20, SD_FD = .09). At follow-up, from 148 participants with rsfMRI BOLD and anatomical scans, we reconstructed 148 functional connectomes, and retained 127 after motion quality-control (Mean_FD = .19, SD_FD = .08).

NKI

Out of 301 participants with a single scanning session of rsfMRI and T1w data, we reconstructed functional connectomes for 292, and retained 214 following motion quality-control (Mean_FD = .24 mm, SD_FD = .08 mm). From 176 participants with two scanning sessions, we reconstructed functional connectomes for 169 and retained 133 (Mean_FD = .23 mm, SD_FD = .08 mm). From 112 participants with three scanning sessions, we reconstructed functional connectomes for 107 and retained 78 (Mean_FD = .25 mm, SD_FD = .09 mm).

Harmonising CALM Neuroimaging Data

To correct for the scanner software update in March 2021 following baseline collection, which affected 22.83% of functional and 25.58% of DWI scans included in this study, respectively, we harmonized structural and functional connectomes at the edge-level using ComBat with parametric adjustments (Fortin et al., 2017, 2018; Johnson et al., 2007). For each node, ComBat estimates the statistical variance of additive and multiplicative site effects using empirical Bayes. This, alongside a covariance matrix, are regressed out from each node. Our covariates were age at scan, sex, and age-standardised scores for matrix reasoning from the Weschler Abbreviated Scale of Intelligence II (Wechsler, 2011), PhAB (Frederickson et al., 1997) alliteration and object naming, alongside AWMA (Alloway, 2007a) digit recall, backward digit recall, and dot matrix tests.

To assess the effectiveness of harmonisation for thresholded structural and functional connectomes separately, we conducted general linear models of the 8 covariates above and scanner type (2 levels – before and after harmonisation) to predict single global graph theory metrics, computed using the Brain Connectivity Toolbox (Rubinov & Sporns, 2010) in MATLAB R2022a (The MathWorks, Inc., 2022). If harmonisation was successful, we’d expect any significant effects of scanner type before harmonisation to be non-significant after harmonisation. Before harmonisation, scanner type was not a statistically significant predictor for structural efficiency (β < .001, p = .921), assortativity (β < .001, p = .911), density (β < .001, p = .584), transitivity (β = −.003, p = .566), and modularity (β = .001, p = .855). Following harmonisation, scanner type remained a non-significant predictor of structural efficiency (β < .001, p = .935), assortativity (β = .002, p = .803), density (β < .001, p = .139), transitivity (β < −.001, p = .970), and modularity (β < .001, p = .999). For functional connectomes, before harmonisation scanner type was not a statistically significant predictor for efficiency (β = .003, p = .849), assortativity (β < .001, p = .955), density (β < −.001, p = .352), transitivity (β = .004, p = .740), or modularity (β = −.002, p = .869). Following harmonisation, scanner type remained a non-significant predictor of functional efficiency (β = .005, p = .746), assortativity (β = −.005, p = .779), density (β < -.001, p = .841), transitivity (β = .008, p = .486), and modularity (β = .004, p = .727). We therefore concluded that the effects of the scanner upgrade in CALM were minimal and remained so after harmonisation.

Statistical Analyses and Data Processing

Manifold Generation

Using the Brain Space package (Vos de Wael et al., 2020) in Python 3.7.3, we derived group-level and individual-level gradients for structural and functional connectivity, respectively. For each modality and hemisphere, we first constructed an affinity matrix from the group-level connectome, using a normalised angle kernel without additional sparsity, producing a nroi/2 x nroi/2 matrix. This affinity matrix represents statistical similarities between each pair of nodes, normalised between 0.5 (least similarity) and 1 (greatest similarity), and represents the data’s local geometry. We created affinity matrices in each hemisphere to avoid capturing a left-right hemisphere split. We then used diffusion-map embedding as a non-linear dimensionality reduction technique (Coifman & Lafon, 2006) to extract the underlying latent structure of these similarities. The algorithm uses a Markov chain to model transition probabilities between different nodes, using a random walker. The eigenfunctions of this Markov matrix produces a low-dimensional embedding (Coifman et al., 2005), where nodes with similar properties are organised along axes represented by eigenvectors, ordered by decreasing amount of total variance explained. Inter-nodal distances are diffusion distances, where larger distances represent less similarity for a feature, such as connectivity strength. It is controlled by two parameters. The first is diffusion time t, which controls the scale of the embedding. That is, the diffusion operator sums all paths of length t connecting each pair of nodes. The second is anisotropic diffusion α, which controls the contribution of the data distribution to the manifold. The data distribution exerts minimal influence on the manifold when α = 1 (scale-free), and maximum influence when α = 0. In line with prior work (Margulies et al., 2016; Park et al., 2021), we set t = 0, meaning that transitions between nodes i and j along all walk lengths are considered, and α = 0.5. Compared to principal component analysis, diffusion-map embedding does not assume a linear underlying data structure, is robust to noise, and preserves local structure (Coifman & Lafon, 2006; Margulies et al., 2016). Crucially, the pairwise diffusion distances produced are analogous to Euclidean distances, and hence can be mapped back onto the cortex. For each modality, we aligned the right hemisphere manifold to the left using Procrustes rotation. We constructed individual-level manifolds for each modality and hemisphere using individual connectomes and aligned them to the corresponding group template using Procrustes rotation.

Manifold Eccentricity

To examine individual differences in the orientation of structural and functional gradients, we calculated manifold eccentricity, originally conceptualised by Park and colleagues (2021). Consider a 3D structural manifold space, where each participant’s structural organisation can be represented by a nodal 3D coordinate. Manifold eccentricity is the Euclidean distance between the group centroid, defined as the intersection between the 3 group-level templates, and each node’s 3D coordinate. For each participant and modality, we calculated nodal eccentricity, producing a nroi x ngradient vector.

Null Models to Assess Inter-Gradient Relationships

To assess the significance of Spearman-rank correlations between gradients, we used a permutation test preserving spatial autocorrelation, developed for parcellated data by Váša and colleagues (2018). In brief, the parcellated map is projected onto a sphere. FreeSurfer vertices within each parcellated region on the sphere, excluding the medial wall, are averaged to create a centroid, and randomly rotated 10,000 times in each hemisphere. The correlation between the rotated and original coordinates in each hemisphere produces a null distribution.

Measures of Psychopathology and Cognition

Using k-nearest neighbours imputation (k = 5), we imputed missing values for the following measures across datasets: forward (5.56%) and backward (5.56%) digit span tasks from the Automated Working Memory Assessment (Alloway, 2007a); Tower Task total achievement scores (17.94%), visual scanning speed (13.81%), number-letter switching accuracy (16.35%), and motor speed (14.13%) from the Delis-Kaplan assessment suite (Delis et al., 2001); inattention (2.22%), hyperactivity (2.22%), peer relations (2.22%), learning problems (2.06%), executive function (2.06%), and aggression (2.22%) from the Conners III scale (Conners et al., 2011). Only raw cognitive scores were provided, and hence we age-standardised the scores by regressing out age and age².

GAMMs

Using R 4.2.2, we modelled linear and non-linear relationships between age and manifold eccentricity using generalised additive mixed modelling (GAMMs). Like generalised additive models, GAMMs model covariates as a smooth penalised spline, the weighted sum of k basis functions. However, GAMMs additionally include random effects to model correlations between observations, for example in data sets with cross-sectional and longitudinal data. The smooth is decomposed into a combination of fixed unpenalized linear components and random penalized non-linear components (Sørensen et al., 2021). We used the gamm4 R package to specify smooths and parametric coefficients (Lin & Zhang, 1999), and the mgcv R package to specify tensors (Wood, 2017). Mean framewise displacement, sex, and data set were parametric covariates, whilst age and an interaction between age and data set were smooths. Subject-specific random intercepts were included as random effects. To prevent over-fitting, we set a maximum of 3 basis functions. When comparing manifold eccentricity across 7 intrinsic functional connectivity networks (Schaefer et al., 2018), we corrected for multiple comparisons by controlling the false discovery rate.

To visualise the interaction between age and data set as a smooth spline, we predicted manifold eccentricity using n increments of age from youngest to oldest for each data set, where n was the number of participants in that data set. Median mean framewise displacement for each data set and a single factor of sex were used as covariates to generate predictions. Following previous work (Sydnor et al., 2023; Tervo-Clemmens et al., 2023), to generate 95% credible intervals, we took 10,000 draws from a multivariate gaussian distribution with a covariance matrix of the fitted GAMM and generated a critical value by calculating the maximum absolute standardised deviation of the draws from the gaussian distribution to the true GAMM distribution. To assess the significance of the interaction between age and data set, we used the pbkrtest R package to implement a parametric bootstrapped likelihood ratio test of the full model compared to a reduced model (Halekoh & Højsgaard, 2014). We calculated the first derivative of the age smooth for each data set and defined sensitive periods of development as those with non-zero simultaneous confidence intervals (p < .05, two-tailed).

Code/data availability

Analysis scripts and de-identified derivatives from diffusion-map embedding are available at: https://github.com/AlicjaMonaghan/neurodevelopmental_gradients

Acknowledgements

We thank the families of NKI and CALM for their participation, alongside Theodore Satterthwaite and Valerie Sydnor for valuable discussions.