1. Neuroscience
Download icon

An expanding manifold in transmodal regions characterizes adolescent reconfiguration of structural connectome organization

  1. Bo-yong Park  Is a corresponding author
  2. Richard AI Bethlehem
  3. Casey Paquola
  4. Sara Larivière
  5. Raul Rodríguez-Cruces
  6. Reinder Vos de Wael
  7. Neuroscience in Psychiatry Network (NSPN) Consortium
  8. Edward T Bullmore
  9. Boris C Bernhardt  Is a corresponding author
  1. McConnell Brain Imaging Centre, Montreal Neurological Institute and Hospital, McGill University, Canada
  2. Department of Data Science, Inha University, Republic of Korea
  3. Autism Research Centre, Department of Psychiatry, University of Cambridge, United Kingdom
  4. Brain Mapping Unit, Department of Psychiatry, University of Cambridge, United Kingdom
  5. Institute of Neuroscience and Medicine (INM-1), Forschungszentrum Jülich, Germany
Research Article
  • Cited 0
  • Views 833
  • Annotations
Cite this article as: eLife 2021;10:e64694 doi: 10.7554/eLife.64694

Abstract

Adolescence is a critical time for the continued maturation of brain networks. Here, we assessed structural connectome development in a large longitudinal sample ranging from childhood to young adulthood. By projecting high-dimensional connectomes into compact manifold spaces, we identified a marked expansion of structural connectomes, with strongest effects in transmodal regions during adolescence. Findings reflected increased within-module connectivity together with increased segregation, indicating increasing differentiation of higher-order association networks from the rest of the brain. Projection of subcortico-cortical connectivity patterns into these manifolds showed parallel alterations in pathways centered on the caudate and thalamus. Connectome findings were contextualized via spatial transcriptome association analysis, highlighting genes enriched in cortex, thalamus, and striatum. Statistical learning of cortical and subcortical manifold features at baseline and their maturational change predicted measures of intelligence at follow-up. Our findings demonstrate that connectome manifold learning can bridge the conceptual and empirical gaps between macroscale network reconfigurations, microscale processes, and cognitive outcomes in adolescent development.

Introduction

Adolescence is a time of profound and genetically mediated changes in whole-brain network organization (Larsen and Luna, 2018; Menon, 2013). Adolescent development is important for the maturation in cognitive and educational functions and brain health more generally, a notion reinforced by the overlapping onset of several neurodevelopmental and psychiatric disorders (Hong et al., 2019; Khundrakpam et al., 2017; Paus et al., 2008). With increased capacity to carry out longitudinal studies in large samples, it is now possible to track changes in brain network organization within subjects, providing insights into maturational processes, their biological underpinnings, and their effects on behavior and cognition.

By offering an in vivo window into brain organization, neuroimaging techniques, such as magnetic resonance imaging (MRI), offer the ability to track adolescent brain development over time. Several cross-sectional and longitudinal studies in neurodevelopmental cohorts have focused on the analysis of morphological changes (Gogtay et al., 2004; Shaw et al., 2006; Tamnes et al., 2017), including MRI-based cortical thickness (Shaw et al., 2006; Tamnes et al., 2017) and volumetric measures (Gogtay et al., 2004; Tamnes et al., 2017). Studies robustly show initial gray matter increases until mid-late childhood followed by a decline for the rest of the lifespan. During adolescence, cortical thickness decreases in widespread brain regions (Khundrakpam et al., 2013; Shaw et al., 2006; Sotiras et al., 2017; Tamnes et al., 2017). Thus, contextualizing connectome alterations relative to established patterns of cortical thickness findings may establish whether inter-regional network changes occur above and beyond these diffuse effects of regional morphological maturation. More recent work explored changes in intracortical microstructure, capitalizing on myelin-sensitive contrasts such as magnetization transfer ratio (MT) mapping, which generally suggest overall increases in adolescence (Paquola et al., 2019a; Whitaker et al., 2016) together with depth-dependent shifts in intracortical myelin profiles (Paquola et al., 2019a). Besides the increasingly recognized changes in cortico-cortical connectivity, studying subcortical regions offer additional insights for understanding brain maturation during adolescence. Indeed, an increasing body of connectome-level studies emphasizes that subcortical structures contribute significantly to patterns of whole-brain organization, dynamics, and cognition (Hwang et al., 2017; Müller et al., 2020; Shine et al., 2019). In prior neurodevelopmental studies, it has been shown that the volumes of the striatum and thalamus decrease between adolescence and adulthood, potentially paralleling processes resulting in cortical gray matter reduction during this time window (Herting et al., 2018). A close inter-relationship between cortical and subcortical development is also suggested by recent functional connectivity work suggesting that cortico-subcortical pathways are intrinsically remodeled during adolescence (Váša et al., 2020), and these changes affect cognitive functioning. Collectively, these prior findings suggest measurable trajectories of cortical and subcortical structural organization and support associations to cognitive development (Baum et al., 2020; Shaw et al., 2006).

Recent conceptual and methodological advances enable the study of brain organization, development, and substrates underlying cognitive trajectories in humans. One key modality to track developmental changes in structural connectivity is diffusion MRI (dMRI), a technique sensitive to the displacement of water in tissue that allows for the non-invasive approximation of inter-regional white matter tracts. Prior cross-sectional and longitudinal studies in children and adolescents outlined changes in the microstructure of major white matter tracts during development based on the analysis of dMRI-derived tissue parameters (Lebel and Beaulieu, 2011; Schmithorst and Yuan, 2010). These findings have been complemented by assessments of brain network topology using graph-theoretical analysis (Baker et al., 2015; Hagmann et al., 2010; Lebel and Beaulieu, 2011; Oldham and Fornito, 2019), which reported a relatively preserved spatial layout of structural hubs across adolescent development on the one hand (Hagmann et al., 2010), yet with a continued strengthening of their connectivity profiles, likely underpinned by the ongoing maturation of long-range association fibers (Baker et al., 2015; Lebel and Beaulieu, 2011; Oldham and Fornito, 2019).

One emerging approach to address connectome organization and development comes from the application of manifold learning techniques to connectivity datasets. By decomposing whole-brain structural and functional connectomes into a series of lower dimensional axes capturing spatial gradients of connectivity variations, these techniques provide a compact perspective on large-scale connectome organization (Margulies et al., 2016; Paquola et al., 2019b; Vos de Wael et al., 2020a). In addition, these techniques capture multiple, potentially overlapping gradients in connectivity along cortical mantle, which can represent both subregional heterogeneity and multiplicity within a brain region (Haak and Beckmann, 2020). In prior work, we showed that multiple dMRI gradients can illustrate structural underpinnings of dynamic functional communication in the adult human connectome (Park et al., 2021b). In line with prior conceptual accounts, the low-dimensional eigenvectors (i.e., gradients) derived from these techniques provide continuous dimensions of cortical organization, and thus the eigenvectors can jointly generate intrinsic coordinate systems of the brain based on connectivity (Bijsterbosch et al., 2020; Haak et al., 2018; Huntenburg et al., 2018; Margulies et al., 2016; Mars et al., 2018). Beyond these methodological considerations, prior work has shown that the principal gradients estimated from resting-state functional (Margulies et al., 2016), microstructural (Paquola et al., 2019b), and diffusion MRI (Park et al., 2021a) all converge broadly along an established model of sensory-fugal hierarchy and laminar differentiation (Mesulam, 1998), allowing gradient mapping techniques to make conceptual contact to theories of cortical organization, development, and evolution (Buckner and Krienen, 2013; Goulas et al., 2018; Huntenburg et al., 2018; Sanides, 1969; Sanides, 1962). An emerging literature has indeed shown utility of the gradient framework to study primate evolution and cross-species alignment (Blazquez Freches et al., 2020; Valk et al., 2020; Xu et al., 2020), neurodevelopment (Hong et al., 2019; Paquola et al., 2019a), as well as plasticity and structure-function coupling (Park et al., 2021b; Valk Sofie et al., 2020; Vázquez-Rodríguez et al., 2019). In a recent assessment by our team, manifold learning techniques have been applied to myelin sensitive intracortical MT data, showing an increasing myeloarchitectural differentiation of association cortex throughout adolescence (Paquola et al., 2019a). Still, the longitudinal maturation of dMRI connectomes in children and adolescents using manifold techniques has not been tracked.

Imaging-transcriptomics approaches allow for the identification of cellular and molecular factors that co-vary with imaging-based findings (Arnatkeviciute et al., 2019; Fornito et al., 2019; Gorgolewski et al., 2014; Hawrylycz et al., 2015; Thompson et al., 2013). Recently established resources, such as the Allen Human Brain Atlas (Arnatkeviciute et al., 2019; Hawrylycz et al., 2015), can be utilized to spatially associate macroscale imaging/connectome data with the expression patterns of thousands of genes. These findings have already been applied in the study of healthy adults (Hawrylycz et al., 2015; Park et al., 2020) and typically developing adolescents (Mascarell Maričić et al., 2020; Padmanabhan and Luna, 2014; Paquola et al., 2019a; Vértes et al., 2016; Whitaker et al., 2016), as well as individuals suffering from prevalent brain disorders (Altmann et al., 2018; Hashimoto et al., 2015; Klein et al., 2017; Park et al., 2021a; Patel et al., 2021; Romero-Garcia et al., 2019). The gene sets that co-vary with in vivo findings can furthermore be subjected to gene set enrichment analyses to discover potentially implicated molecular, cellular, and pathological processes (Ashburner et al., 2000; Carbon et al., 2019; Chen et al., 2013; Dougherty et al., 2010; Kuleshov et al., 2016; Morgan et al., 2019; Romero-Garcia et al., 2018; Subramanian et al., 2005). For example, studies in newborns have shown that cortical morphology reflects spatiotemporal patterns of gene expression in fetuses, linking molecular mechanisms to in vivo measures of cortical development in early life (Ball et al., 2020). Work in adolescents has furthermore shown that developmental changes in regional cortical thickness measures and myelin proxies spatially co-localize with the expression patterns of genes involved in synaptic and oligodendroglial function (Paquola et al., 2019a; Whitaker et al., 2016). Building on these prior investigations, the current study aimed at exploring whether adolescent structural connectome reconfigurations, assessed using manifold learning techniques, reflect the expression patterns of specific genes in order to identify potential molecular signatures of macroscale structural network development.

Here, we charted developmental changes in structural connectome organization, based on an accelerated longitudinal neuroimaging study involving 208 participants investigated between 14 and 26 years of age (Kiddle et al., 2018; Whitaker et al., 2016). Compared to cross-sectional designs, longitudinal studies track within-subject change, separating developmental effects from between-subject variability (Louis et al., 1986). We first estimated longitudinal changes in structural connectome manifolds across age. This compact and lower dimensional space furthermore allowed for the integration of connectome-level findings with changes in MRI-based measures of cortical morphology and intracortical myelin. We furthermore projected subcortico-cortical connectivity patterns into the manifold space to assess parallel developmental shifts of these pathways in the studied time window. Connectome manifold changes were contextualized at the molecular level via transcriptomic association and developmental enrichment analyses based on post-mortem datasets, which furthermore allowed for data-driven exploration of time windows of spatially co-localized gene sets. To also assess behavioral associations of connectome manifold changes, we utilized supervised machine learning to predict future measures of cognitive function quantified via the intelligence quotient (IQ). IQ is a widely used marker of general cognitive abilities, which shows good test–retest reliability (Brown and May, 1979; Watkins and Smith, 2013; Catron, 1978; G.-Matarazzo et al., 1973; Snow et al., 1989; Wagner and Caldwell, 1979) and has previously been applied to index overall cognitive function during development (Crespi, 2016; Garde et al., 2005; Garde et al., 2000; Koenis et al., 2018; Park et al., 2016; Ramsden et al., 2011; Shaw et al., 2006; Suprano et al., 2020). In the study of neurodevelopment, neuroimaging reports have previously assessed associations between IQ and large-scale network measures in children to adolescents (Koenis et al., 2018; Ramsden et al., 2011; NSPN Consortium et al., 2018; Shaw et al., 2006; Suprano et al., 2020). Multiple sensitivity analyses were conducted at several steps to verify the robustness of our findings, and analytical code is made fully accessible to allow for independent replication of our findings.

Results

These findings were based on the Neuroscience in Psychiatry Network (NSPN) cohort (Kiddle et al., 2018; Whitaker et al., 2016). In brief, we studied 208 healthy individuals enrolled in an accelerated longitudinal study, with approximately equal numbers of males and females in each of five age-related strata that collectively spanned the time period from 14 to 25 years coinciding with transition from adolescence to young adulthood. Participants (48% female) had a mean age of 18.82 years (range = 14–25 years) at baseline and 19.95 years (15–26 years) at follow-up. The average interval between baseline and follow-up scan was 11.28 months (range = 6–12 months). See Materials and methods for details on participant selection, image processing, and analysis.

Macroscale structural connectome manifold

For every participant, we built cortex-wide structural connectome manifolds formed by the eigenvectors displaying spatial gradients in structural connectome organization using non-linear dimensionality reduction techniques (Vos de Wael et al., 2020aVos de Wael et al., 2020bhttps://github.com/MICA-MNI/BrainSpace). Individual manifolds were aligned to a template manifold estimated from a hold-out dataset (see Materials and methods) (Langs et al., 2015; Vos de Wael et al., 2020a). Three eigenvectors (E1, E2, and E3) explained approximately 50% of information in the template affinity matrix (i.e., 20.7/15.8/13.5% for E1/E2/E3, respectively), with each eigenvector showing a different axis of spatial variation across the cortical mantle (Figure 1A). Eigenvectors depicted a continuous differentiation between medial and lateral cortices (E1), between inferior and superior cortices (E2), and between anterior and posterior areas (E3). For each participant and time point, we calculated manifold eccentricity, which depicts how far each node is located from the center of the template manifold (see Materials and methods). It thus quantifies the changes in eigenvectors between the time points in terms of expansion and contraction instead of comparing multidimensional connectome manifolds (Bethlehem et al., 2020). The manifold eccentricity showed high values in frontal and somatomotor regions, while temporoparietal, visual, and limbic regions showed low values (Figure 1B).

Figure 1 with 15 supplements see all
Structural connectome manifolds.

(A) Systematic fiber tracking based on diffusion magnetic resonance imaging generated a cortex-wide structural connectome, which was subjected to diffusion map embedding. As shown in the scree plot, three eigenvectors (E1, E2, E3) accounted for approximately 50% information of connectome data, and each depicted a different gradual transition across the cortical mantle. (B) Manifold eccentricity measured by Euclidean distance between the template center and each data point. Arrows depict average positional change in connectivity space from baseline to follow-up. The color of each arrow represents each brain region mapped on the surface on the bottom. (C) The histogram represents age distribution of all subjects at baseline and follow-up. The colors on brain surfaces indicate t-statistics of regions showing significant longitudinal changes in manifold eccentricity across age, following multiple comparisons correction with a false discovery rate < 0.05. Datapoint colors in the scatter plot represent t-statistics. Identified regions are represented with arrows that originate from baseline to follow-up. (D) Stratification of age-related changes in manifold eccentricity according to prior models of cortical hierarchy (Mesulam, 1998) and functional magnetic resonance imaging communities (Yeo et al., 2011).

Figure 1—source data 1

Source files for connectome manifolds and age-related changes in manifold eccentricity.

https://cdn.elifesciences.org/articles/64694/elife-64694-fig1-data1-v2.xlsx

Changes in manifold eccentricity across age

Leveraging linear mixed effect models that additionally controlled for effects of sex, site, head motion, and subject-specific random intercepts (Worsley et al., 2009), we assessed changes in manifold eccentricity across age (see Materials and methods). Manifold eccentricity expanded as age increased, especially in bilateral prefrontal and temporal areas, as well as left early visual and right lateral parietal cortices (false discovery rate [FDR] < 0.05; Benjamini and Hochberg, 1995; Figure 1C). Stratifying these effects along four cortical hierarchical levels, defined using an established taxonomy based on patterns of laminar differentiation and tract-tracing data in non-human primates (Mesulam, 1998), we identified peak effects in heteromodal association and paralimbic areas (Figure 1D). Convergent findings were observed when analyzing the effects with respect to intrinsic functional communities (Yeo et al., 2011), showing highest effects in default mode and limbic areas followed by visual and frontoparietal cortices. No significant contraction of manifold eccentricity was observed. In addition, we could not find any significant effects when we fitted the model with a quadratic form of age (i.e., age2), indicating the manifold eccentricity linearly increases across age.

To conceptualize the findings derived from manifold eccentricity with respect to conventional network topologies, we correlated manifold eccentricity changes with several graph-theoretical measures of structural connectome (Figure 1—figure supplement 1; Rubinov and Sporns, 2010). We first defined six spatially contiguous clusters within the regions that showed significant age-related changes in manifold eccentricity (see Figure 1C) and correlated within-subject changes in manifold eccentricity with developmental changes in degree centrality, connectivity distance, and modular parameters (i.e., within-module degree and participation coefficient based on modules defined via Louvain’s community detection algorithm [Blondel et al., 2008]; see Materials and methods; Figure 1—figure supplement 2). We found significant positive associations for degree centrality and within-module degree, suggesting that connectome manifold expansion reflects a concurrent increase of overall connectivity, particularly within modules. Stratifying changes in manifold eccentricity, as well as connectome topology measures, according to the discretized age bins confirmed these age-related trends (Figure 1—figure supplement 3). Indeed, except for participation coefficient, values in general increased from childhood to young adulthood.

Effects of cortical morphology and microstructure

Previous studies demonstrated significant changes in cortical morphology and microstructure during adolescence, showing co-occurring reductions in cortical thickness and MT skewness, the latter being an index of depth-dependent intracortical myelin changes in multiple lobes (Gogtay et al., 2004; Khundrakpam et al., 2017; Paquola et al., 2019a; Shaw et al., 2006). We replicated these findings by showing cortical thinning in almost all brain regions across the studied age window as well as reductions in depth-dependent MT skewness, suggestive of supragranular enrichment of myelin (Figure 2A). To evaluate whether the age-related changes in manifold eccentricity were robust above and beyond these regional changes in cortical thickness and MT, we implemented linear mixed effect models including cortical thickness and MT as covariates in the analysis of developmental change in manifold eccentricity (Figure 2B). While we observed virtually identical spatial patterns of manifold eccentricity changes in models that controlled for thickness, MT skewness, and both, age-related effects in regions of significant manifold eccentricity findings (see Figure 1C) were reduced in models that additionally controlled for these covariates (average reduction of t-value in models controlling for thickness/MT skewness/both = 42/18/68%).

Age-related effects on macro- and microstructural metrics of cortical anatomy.

(A) The t-statistics of identified regions that showed significant age-related changes in cortical thickness (upper row) and magnetization transfer ratio MT (bottom row), and stratification of t-statistics according to cortical hierarchy (Mesulam, 1998) and functional community (Yeo et al., 2011). (B) Age-related changes in manifold eccentricity after controlling for cortical thickness and MT.

Figure 2—source data 1

Source files for age-related changes in cortical thickness and magnetization transfer ratio.

https://cdn.elifesciences.org/articles/64694/elife-64694-fig2-data1-v2.xlsx

Age-related changes in subcortico-cortical connectivity

Besides visualizing these changes in cortico-cortical connectivity, we also capitalized on the manifold representation to assess adolescent changes in the connectivity of subcortical regions, to obtain a more holistic insight into whole-brain connectome reconfigurations during this time period, and to examine whether subcortical connectivity patterns undergo parallel developmental trajectories (Hwang et al., 2017; Shine et al., 2019). Specifically, we assessed changes in subcortical-weighted manifolds across age, defined by projecting the streamline strength of subcortical regions to cortical targets to the manifold space (see Materials and methods). Such an analysis situates changes in subcortico-cortical pathways in the macroscale context of cortico-cortical connectivity identified in the previous analyses. After multiple comparisons correction, the caudate and thalamus showed significant age-related effects on subcortical-weighted manifolds (FDR < 0.05; Figure 3), and marginal effects were observed in the putamen, pallidum, and hippocampus (FDR < 0.1).

Longitudinal changes in subcortical-weighted manifolds.

The t-statistics of age-related changes in subcortical-weighted manifolds. The effects of each subcortical region are reported on the radar plot. FDR: false discovery rate.

Figure 3—source data 1

Source files for age-related changes in subcortical-weighted manifolds.

https://cdn.elifesciences.org/articles/64694/elife-64694-fig3-data1-v2.xlsx

Transcriptomic association analysis

Connectome organization, in general, and macroscale gradients, in particular, have been argued to reflect genetic expression profiles, underscoring the close link between the physical layout of the brain and innate transcriptional patterning (Buckner and Krienen, 2013; Fornito et al., 2019). Here, we carried out a transcriptomic association analysis and developmental enrichment analyses to contextualize the age-related manifold eccentricity changes with respect to patterns of post-mortem gene expression from a sample of independent adults (Figure 4A). Specifically, leveraging mixed effect models, we associated the spatial patterns of manifold change across age in the NSPN sample (controlling for covariation of cortical thickness and MT) with cortical maps of post-mortem gene expression data from the Allen Institute for Brain Sciences (Arnatkeviciute et al., 2019; Gorgolewski et al., 2015; Gorgolewski et al., 2014; Hawrylycz et al., 2012; Markello et al., 2020). Among the list of most strongly associated genes (FDR < 0.05), we selected only genes that were consistently expressed across different donors (r > 0.5) (Arnatkeviciute et al., 2019; Hawrylycz et al., 2012; Markello et al., 2020; Supplementary file 1). We performed developmental gene set enrichment analysis using the cell-type-specific expression analysis (CSEA) tool, which compares the selected gene list with developmental enrichment profiles (see Materials and methods) (Dougherty et al., 2010; Xu et al., 2014). This analysis highlights developmental time windows across macroscopic brain regions in which genes are strongly expressed. We found marked expression of the genes enriched from childhood onward in the cortex, thalamus, and cerebellum (FDR < 0.001; Figure 4B). Although signal was reduced, genes were also enriched for expression in the striatum at the transition from childhood to adolescence (FDR < 0.05). On the other hand, identified genes were not found to be expressed in the hippocampus and amygdala.

Transcriptomic analysis.

(A) Gene decoding process by associating t-statistics from the linear mixed effect model with post-mortem gene expression maps. (B) We identified genes that were spatially correlated with the input t-statistic map (false discovery rate [FDR] < 0.05) and selected only those that were furthermore consistently expressed across different donors (r > 0.5). These genes were input to a developmental enrichment analysis, showing strong associations with cortex, thalamus, striatum, and cerebellum during the childhood-to-adulthood time window. The degree of gene expression for developmental windows is reported on the bottom. The curve represents log transformed FDR-corrected p-values, averaged across the brain regions for each of the time windows reported on the bottom. ***FDR < 0.001, ** FDR < 0.01, *FDR < 0.05.

Figure 4—source data 1

Source files for developmental enrichment profiles.

https://cdn.elifesciences.org/articles/64694/elife-64694-fig4-data1-v2.xlsx

Association between connectome manifold and cognitive function

Finally, to establish associations between connectome reconfigurations and cognitive functioning, we utilized supervised machine learning to predict full IQ at follow-up using manifold eccentricity features. Independent variables were combinations of cortical and subcortical manifold features at baseline and their age-related trajectory data. We used elastic net regularization with nested ten-fold cross-validation (Cawley and Talbot, 2010; Parvandeh et al., 2020; Tenenbaum et al., 2000; Varma and Simon, 2006; Zou and Hastie, 2005) (see Materials and methods), and repeated the prediction 100 times with different training and test dataset compositions to mitigate subject selection bias. Across cross-validation and iterations, 6.24 ± 5.74 (mean ± SD) features were selected to predict IQ using manifold eccentricity of cortical regions at baseline, 6.20 ± 5.14 cortical features at baseline and maturational change, 5.45 ± 5.99 cortical and subcortical features at baseline, and 5.16 ± 5.43 at baseline and maturational change, suggesting that adding more independent variables may not per se lead to improvement in prediction accuracy. The manifold eccentricity of cortical regions at baseline significantly predicted future IQ score (mean ± SD r = 0.14 ± 0.04; mean absolute error [MAE] = 8.93 ± 0.16, p=0.09). Prediction performance was slightly improved when we combined the manifold eccentricity both at baseline and differences between follow-up and baseline (r = 0.18 ± 0.04; MAE = 9.10 ± 0.19, p=0.04) (Figure 5A). Notably, prediction accuracy was improved if we additionally considered subcortical manifold features (baseline: r = 0.17 ± 0.03; MAE = 8.74 ± 0.11, p=0.04; baseline and maturational change: r = 0.21 ± 0.02; MAE = 8.86 ± 0.14, p=0.01) (Figure 5B). The regions showing strongest predictive validity for IQ were prefrontal, parietal, and temporal cortices, as well as the caudate and thalamus. The probability map of the selected brain regions (bottom right of Figure 5B) was further decoded using Neurosynth (Yarkoni et al., 2011), revealing strong associations with higher-order cognitive and social terms (Figure 5—figure supplement 1). We compared the prediction performance of our model with a baseline model, where IQ of the test set was simple average of training set (r = −0.15 ± 0.06, MAE = 8.98 ± 0.04, p=0.12; see Materials and methods). We found that our model outperformed this baseline model (Meng’s z-test p < 0.001) (Meng et al., 1992). We also predicted the change of IQ between the baseline and follow-up, instead of IQ at follow-up, using the imaging features. However, we could not find significant results.

Figure 5 with 2 supplements see all
Intelligence quotient (IQ) prediction by baseline and follow-up measures of cortical and subcortical manifolds.

(A) Probability of selected brain regions across ten-fold cross-validation and 100 repetitions for predicting future IQ using only baseline manifold eccentricity (left), and both baseline and maturational change in the feature (right). Correlations between actual and predicted IQ are reported. Black lines indicate mean correlation, and gray lines represent 95% confidence interval for 100 iterations with different training/test dataset. (B) The prediction performance when both cortical and subcortical features were considered. MAE: mean absolute error.

Figure 5—source data 1

Source files for selected probability as well as actual and predicted intelligence quotient.

https://cdn.elifesciences.org/articles/64694/elife-64694-fig5-data1-v2.xlsx

Sensitivity analysis

Spatial scale

Repeating the longitudinal modeling with a different spatial scale (i.e., 300 parcels), findings were highly consistent (Figure 1—figure supplement 4).

Site and sex effects

Furthermore, manifold eccentricity of the identified cortical regions and age consistently correlated positively across different sites and within both biological sexes, yielding non-significant interaction effects (Figure 1—figure supplement 5).

Different parameters for diffusion map embedding

When we changed parameters of diffusion map embedding for generating connectome manifolds (see Materials and methods), t-statistic maps of age-related changes in manifold eccentricity were largely consistent (mean ± SD linear correlation r = 0.92 ± 0.10).

Gradient alignment fidelity

When calculating linear correlations between template and individual manifolds before and after alignment, we found significant increases after alignment (r = 0.92 ± 0.03/0.93 ± 0.03/0.94 ± 0.03) compared to before alignment (−0.02 ± 0.03/–0.001 ± 0.37/0.003 ± 0.12) for E1/E2/E3, respectively, supporting effectiveness of alignment. After excluding 10% of subjects with poor alignment (cutoff r = 0.83; the new set was correlated with the template manifold, r = 0.94 ± 0.01), we found consistent age-related changes in manifold eccentricity (Figure 1—figure supplement 6), with the t-statistic map showing strong correlation to the map derived in the whole sample (r = 0.97, p<0.001).

Connectome manifold generation using principal component analysis

In a separate analysis, we generated eigenvectors using principal component analysis (Wold et al., 1987), instead of diffusion map embedding (Coifman and Lafon, 2006), and found consistent spatial maps (linear correlation = 0.998 ± 0.001 across E1/E2/E3; Figure 1—figure supplement 7A) and longitudinal findings (Figure 1—figure supplement 7B).

Longitudinal changes in graph-theoretical measures

Repeating the longitudinal modeling using graph-theoretical centrality measures, we found significant age-related longitudinal changes in degree and eigenvector centrality, while betweenness centrality did not reveal significant effects, in similar regions to those that had significant age-related changes in manifold eccentricity (Figure 1—figure supplement 8). Correlating the effect size maps for manifold eccentricity and each graph measure, we found a significant yet variable spatial similarity of the effect maps (betweenness centrality: r = 0.18, spin-test p = 0.02; degree centrality: r = 0.57, p < 0.001; eigenvector centrality: r = 0.47, p < 0.001).

Manifold eccentricity based on all eigenvectors

Repeating manifold eccentricity calculation and age modeling using all eigenvectors, instead of using only the first three, we observed relatively consistent results with our original findings (linear correlation of manifold eccentricity r = 0.54, p<0.001; t-statistic map r = 0.68, p<0.001), also pointing to manifold expansion in transmodal cortices (Figure 1—figure supplement 9).

Robustness of group representative structural connectome

We compared gradients derived from the group representative structural connectome, based on (i) distance-dependent thresholding (Betzel et al., 2019) and (ii) consistency thresholding (Wang et al., 2019; Figure 1—figure supplement 10). We found high similarity in spatial maps of the estimated manifolds (r = 0.89 ± 0.01 for E1; 0.93 ± 0.004 for E2; 0.85 ± 0.01 for E3 across six different thresholds), indicating robustness.

Connectome manifolds based on structural parcellation

We repeated our analyses with a structural parcellation, defined using a sub-parcellation of folding based on the Desikan–Killiany atlas (Desikan et al., 2006; Vos de Wael et al., 2020a; Figure 1—figure supplement 11). Despite slight differences in the topography of manifold eccentricity in lateral prefrontal, temporal, and occipital cortices, we could replicate strong age-related effects in heteromodal association areas, together with effects in caudate and hippocampus (FDR < 0.05), and marginally in thalamus (FDR < 0.1).

Longitudinal modeling using edge weights

Repeating the longitudinal modeling across age using connectome edge weights, we found significant increases in edge weights within frontoparietal and default mode networks, as well as in attention and sensory networks (FDR < 0.05; Figure 1—figure supplement 12), consistent with findings based on manifold eccentricity.

Manifold eccentricity and pubertal stages

We repeated the longitudinal modeling within a subset of participants who completed the Tanner scale (n = 73) (Marshall and Tanner, 1970; Marshall and Tanner, 1969) and found relatively consistent albeit weaker age-related changes in manifold eccentricity as for the overall sample (Figure 1—figure supplement 13A). Notably, manifold eccentricity within the identified regions derived from overall sample and Tanner scale revealed a significant interaction effect (t = 2.36, p=0.01; Figure 1—figure supplement 13B), suggesting that participants in early pubertal stages show more marked changes in manifold eccentricity across age compared to those in later stages.

IQ prediction using nonlinear model

We predicted IQ at follow-up using a regression tree method (Breiman et al., 1984), instead of linear regression model, but we could not find improved prediction performance (Figure 5—figure supplement 2).

Discussion

The current study tracked whole-brain structural connectome maturation from adolescence to young adulthood in an accelerated longitudinal imaging cohort (Kiddle et al., 2018; Whitaker et al., 2016). Capitalizing on advanced manifold learning techniques applied to dMRI-derived connectomes, we established that higher-order association cortices in prefrontal, medial and superior temporal areas, as well as parieto-occipital regions, expanded in their connectome manifold representation indicative of an increased differentiation of these systems from the rest of the brain in adolescence. Parallel topological analysis based on graph theory indicated that these changes anatomically coincided with increases in the within-module connectivity of transmodal cortices. Findings were consistent across the different acquisition sites and biological sexes, and similar albeit slightly weaker when correcting connectivity manifolds for MRI-based measures of macrostructure (cortical thickness) and microstructure (skewness of MT depth profile). In addition to the cortical manifold expansion, we found parallel reconfigurations of subcortical connectivity patterns for the caudate and thalamus. Decoding our findings with post-mortem gene expression maps implicated genes enriched in adolescence and young adulthood, again pointing to both cortical as well as subcortical targets. Finally, the combination of both cortical and subcortical manifold measures predicted behavioral measures of intelligence at follow-up, with higher performance than cortical or subcortical data alone. Collectively, our findings provide new insights into adolescent structural connectome maturation and indicate how multiple scales of cortical and subcortical organization can interact in typical neurodevelopment.

Leveraging advanced manifold learning, we depicted macroscale connectome organization along continuous cortical axes. Similar approaches have previously been harnessed to decompose microstructural (Paquola et al., 2019b; Paquola et al., 2019a) and functional MRI (Bethlehem et al., 2020; Hong et al., 2019; Margulies et al., 2016; Murphy et al., 2019; Vos de Wael et al., 2020a). These techniques are appealing as they offer a low-dimensional perspective on connectome reconfigurations in a data-driven and spatially unconstrained manner. In our longitudinal study, we could identify marked connectome expansion during adolescence, mainly encompassing transmodal and heteromodal association cortex in prefrontal, temporal, and posterior regions, the territories known to mature later in development (Gogtay et al., 2004; Shaw et al., 2006). Findings remained consistent when we considered a linear dimensionality reduction technique, suggesting robustness to methodological details of this analysis. Connectome expansion can be understood as an overall greater differentiation of the connectivity of these areas from the rest of the brain as they would then cover wider portions of the corresponding manifold space. Manifold expansion in higher-order areas correlated with an increase in their within-module connectivity, but not with participation coefficient and connectivity distance measures that would be more reflective of their between-module connectivity. In light of potential limitations of dMRI tractography in detecting long-distance fiber tracts (Betzel et al., 2019; Maier-Hein et al., 2017), we cannot rule out a reduced sensitivity of our approach for the study of long-range inter-regional connections. Nevertheless, our diffusion modeling was based on constrained spherical-deconvolution approaches together with SIFT2-based tractogram filtering, in addition to using a distance-dependent thresholding approach that may have partially mitigated these limitations (Betzel et al., 2019). Thus, our findings do overall confirm and extend prior dMRI studies that have focused on specific tracts and that have indicated considerable developmental shifts in diffusion parameters, such as increases in fractional anisotropy and decreases in mean diffusivity in early and late adolescence (Olson et al., 2009). Other studies have furthermore reported increased streamline count estimates (Genc et al., 2020). In this context, our macroscale manifold findings likely reflect an ongoing consolidation of transmodal cortical communities. These findings align with prior graph-theoretical studies, which have pointed to concurrent increases in network integration and consolidation of network hubs from late childhood to early adulthood (Baker et al., 2015; Lebel and Beaulieu, 2011; Oldham and Fornito, 2019). Considering their distributed regional substrate, these network effects are likely driven by the ongoing maturation of fiber bundles that interconnect these higher-order cortices, including superior longitudinal fascicules, but also thalamic and basal ganglia pathways (Tamnes et al., 2010), throughout adolescence.

Projecting manifold solutions back onto cortical surfaces allowed us to integrate our connectome manifold results with morphometric and intracortical intensity indices obtained via structural and quantitative MRI contrasts in the same participants. We were thus able to balance the network-level effects against trajectories of intracortical remodeling. Longitudinal changes in these cortical features were overall in agreement with prior work, suggesting marked reductions in cortical thickness in adolescence (Khundrakpam et al., 2017; Shaw et al., 2006), possibly reflecting synaptic pruning processes (Petanjek et al., 2011) together with decreases in the skewness of intracortical MT profiles, a feature sensitive to preferential myelination of supragranular layers (Paquola et al., 2019a). Although we still observed significant age-related changes in manifold eccentricity after controlling for these intracortical and morphological measures, the effect sizes of our findings were reduced. This effect was particularly evident when running a parallel analysis that additionally controlled for depth-dependent shifts in cortical microstructure, a finding in line with more generally demonstrated links between cortical microstructural depth profiles and inter-regional connectivity (Paquola et al., 2019b). In the context of adolescence and prior findings in the NSPN dataset (Paquola et al., 2019a), these results thus suggest a coupled process that affects depth-dependent shifts in cortical myeloarchitecture, on the one hand, and adolescent shifts in macroscale connectome organization, on the other hand, as shown by our longitudinal manifold analyses.

In addition to emphasizing a distributed set of association cortices and their cortico-cortical connections, analysis of subcortico-cortical connectivity patterns highlighted parallel developmental processes in several subcortical structures and their connections, particularly the caudate and thalamus. These findings were independently supported by transcriptomic association studies and developmental enrichment analyses, which implicated genes expressed in cortical regions and these subcortical structures during late childhood, adolescence, and young adulthood. The caudate nucleus of the striatum has long been recognized to play an important role in mediating large-scale cortical network organization (Aglioti, 1997; Alexander and Crutcher, 1990; Graybiel, 1995), a finding also increasingly recognized in the connectome literature (Hwang et al., 2017; Müller et al., 2020; Shine et al., 2019). It is known to modulate activity in prefrontal association areas during memory-driven internal thought processes (Aglioti, 1997), and higher-order cognitive functions, notably motivational processes, decision making, as well as cognitive control and executive functions more generally (Aglioti, 1997; Graybiel, 1995). Regions of the striatum participate in dense cortico-subcortical feedback loops and exchange neural activity through dense connections with adjacent basal ganglia structures as well as the thalamus (Aglioti, 1997; Alexander and Crutcher, 1990; Shine, 2021). Associating macroscopic changes in manifold eccentricity with post-mortem microarray data provided by the Allen Human Brain Atlas (Arnatkeviciute et al., 2019; Fornito et al., 2019; Gorgolewski et al., 2014; Hawrylycz et al., 2015; Thompson et al., 2013), we identified gene sets expressed in cortical regions and subcortical structures of the thalamus and striatum during late childhood, adolescence, and young adulthood. Despite these findings being associative and based on separate datasets, they overall support our results that brain network maturation from late childhood until early adulthood implicates micro- and macroscale factors in both subcortical and cortical networks. Coupled network and molecular changes may ultimately change subcortical and cortical circuit properties, including the balance of excitation and inhibition (E/I). Human brain development involves spatio-temporal waves of gene expression changes across different brain regions and developmental time windows (Ip et al., 2010; Kang et al., 2011; Shin et al., 2018). In the study of adolescent development, prior studies have suggested shifts in E/I balance, evolving from a dominant inhibitory bias in early developmental stages towards stronger excitatory drive in later stages, and suggested that these may underlie the maturation of cognitive functions such as working memory and executive control (Dorrn et al., 2010; Lander et al., 2017; Liu et al., 2007). In common neurodevelopmental disorders, including autism, schizophrenia, and attention-deficit hyperactivity disorder, imbalances in cortical E/I and cortico-subcortical network function have been demonstrated (Cellot and Cherubini, 2014; Gandal et al., 2018; Lee et al., 2017; Lewis et al., 2005; Nelson and Valakh, 2015; Park et al., 2021a; Sohal and Rubenstein, 2019; Trakoshis et al., 2020), potentially downstream to perturbations of different neurotransmitter systems, such as interneuron-mediated GABA transmission (Bonaventura et al., 2017; Kilb, 2012; Liu et al., 2007; Park et al., 2021a; Silveri et al., 2013; Trakoshis et al., 2020; Tziortzi et al., 2014).

Higher-order cognitive function implicates functionally relevant whole-brain network mechanisms, and its prediction may thus leverage structurally governed principles of network integration and segregation. Application of a supervised machine learning framework with cross-validation and regularization to our cohort demonstrated that it is possible to predict inter-individual variations in future IQ from structural connectome manifold data. These findings complement conceptual accounts linking brain organization to cognitive function (Margulies et al., 2016; Mesulam, 1998) and earlier efforts to predict IQ measures from inter-regional measures of connectivity and graph-theoretical indices of network topology (Greene et al., 2018). Notably, evaluations of several feature combinations highlighted that predictive performance was highest when including both baseline and trajectory data, and when complementing cortical and subcortical manifold features. These findings re-emphasize the benefits of incorporating subcortical nodes in the characterization of large-scale cortical network organization and overall cognitive function (Alves et al., 2019; Müller et al., 2020; Shine, 2021; Shine et al., 2019). Of note, although our model significantly outperformed a baseline model, the relationship between the actual and predicted IQ scores did not locate on the equality line and the strength of the association was rather weak. Further improvements in brain-based IQ prediction in adolescence, for example, through combinations of structural and functional imaging features, will be a focus of future work.

Adolescence is a time characterized by ongoing brain changes (Baum et al., 2020; Gogtay et al., 2004; Larsen and Luna, 2018; Menon, 2013; Shaw et al., 2006), gradually increasing independence from caregivers, accompanied by strong increments in knowledge and our ability to think more abstractly and to cooperate with our peers to achieve common goals. On the other hand, adolescence is also a sensitive time window for risk taking, the development of addictions, and is associated with high rates of onset of several psychiatric disorders (Hong et al., 2019; Khundrakpam et al., 2017). Our study has shown that structural brain network organization continues to mature significantly during this time period, with higher-order association cortices in prefrontal and posterior regions especially showing an expansion of their corresponding connectome manifold signature. Findings were related to an increased strengthening of intra-community connectivity as well as cortico-subcortical connectivity to thalamo-striatal regions. Although the current work was restricted to a longitudinal sample of typically developing adolescents, our framework may be useful to explore multiscale network perturbations in cohorts with a psychiatric diagnosis or those at risk for addiction or atypical neurodevelopment.

Materials and methods

Participants

We obtained imaging and phenotypic data from the NSPN 2400 cohort, which contains questionnaire data on 2402 individuals (with MRI data on a subset of ~300) from adolescence to young adulthood in a longitudinal setting (Kiddle et al., 2018; Whitaker et al., 2016). The NSPN study was ethically approved by the National Research Ethics Service and conducted in accordance with NHS research governance standards. All participants provided informed consent in writing, with additional parental consent for participants aged less than 16 years at enrollment. Included participants completed quality-controlled (see Data preprocessing section) multimodal MRI scans consisting of T1-weighted, MT, and dMRI for at least two time points. Our final sample consisted of a total of 208 participants (48% female; mean [range] age = 18.82 [14–25] years at baseline and 19.95 [15–26] at follow-up; inter-scan interval of 11.28 [6–12] months), collected from three different UK sites: Wolfson Brain Imaging Centre and MRC Cognition and Brain Sciences Unit in Cambridge; and University College London. We divided the participants into template and non-template cohorts with matched age, sex, and site ratio. The template dataset (n = 30; 50% female; mean [range] age = 18.69 [15–24] years at baseline and 19.84 ± 2.66 [16–25] at follow-up) was used for constructing the group mean template manifold and the non-template dataset (n = 178; 47% female; mean [range] age = 18.84 [14–25] years at baseline and 19.97 [15–26] at follow-up) was used for conducting main analyses. Of note, changing the template dataset composition did not markedly affect main findings (Figure 1—figure supplement 14).

MRI acquisition

Request a detailed protocol

Imaging data were obtained using Siemens Magnetom TIM Trio 3T scanners. T1-weighted and MT sequences were acquired using a quantitative multiparameter mapping sequence (repetition time [TR]/flip angle = 18.7 ms/20° for T1-weighted and 23.7 ms/6° for MT; six equidistance echo times [TE] = 2.2–14.7 ms; voxel size = 1 mm3; 176 slices; field of view [FOV] = 256 × 240 mm; matrix size = 256 × 240 × 176) (Weiskopf et al., 2013). The dMRI data were acquired using a spin-echo echo-planar imaging sequence (TR = 8700 ms; TE = 90 ms; flip angle = 90°; voxel size = 2 mm3; 70 slices; FOV = 192 × 192 mm2; matrix size = 96 × 96 × 70; b-value = 1000 s/mm2; 63 diffusion directions; and 6 b0 images).

Data preprocessing

Request a detailed protocol

T1-weighted data were processed using the fusion of neuroimaging preprocessing (FuNP) pipeline integrating AFNI, FSL, FreeSurfer, ANTs, and Workbench (Avants et al., 2011; Cox, 1996; Fischl, 2012; Glasser et al., 2013; Jenkinson et al., 2012; Park et al., 2019), which is similar to the minimal preprocessing pipeline for the Human Connectome Project (Glasser et al., 2013). Gradient nonlinearity and b0 distortion correction, non-brain tissue removal, and intensity normalization were performed. The white matter and pial surfaces were generated by following the boundaries between different tissues (Dale et al., 1999), and they were averaged to generate the midthickness contour, which was used to generate the inflated surface. The spherical surface was registered to the Conte69 template with 164k vertices (Van Essen et al., 2012) and downsampled to a 32k vertex mesh. Quality control involved visual inspection of surface reconstruction of T1-weighted data, and cases with faulty cortical segmentation were excluded. Surface-based co-registration between T1-weighted and MT weighted scans was performed. We generated 14 equivolumetric cortical surfaces within the cortex and sampled MT intensity along these surfaces (Paquola et al., 2019a). The vertex-wise MT profiles for each surface depth were averaged based on the Schaefer atlas with 200 parcels (Schaefer et al., 2018). The dMRI data were processed using MRtrix3 (Tournier et al., 2019), including correction for susceptibility distortions, head motion, and eddy currents. We visually inspected the quality of co-registration between the adolescence data and adult-driven surface template as well as parcellation atlas, and all data showed reasonable registration results.

Structural connectome manifold identification

Request a detailed protocol

Structural connectomes were generated from preprocessed dMRI data (Tournier et al., 2019). Anatomically constrained tractography was performed using different tissue types derived from the T1-weighted image, including cortical and subcortical gray matter, white matter, and cerebrospinal fluid (Smith et al., 2012). Multi-shell and multi-tissue response functions were estimated (Christiaens et al., 2015), and constrained spherical-deconvolution and intensity normalization were performed (Jeurissen et al., 2014). The tractogram was generated with 40 million streamlines, with a maximum tract length of 250 and a fractional anisotropy cutoff of 0.06. Subsequently, spherical-deconvolution informed filtering of tractograms (SIFT2) was applied to reconstruct whole-brain streamlines weighted by the cross-section multipliers, which considers the fiber bundle’s total intra-axonal space across its full cross-sectional extent (Smith et al., 2015). The structural connectome was built by mapping the reconstructed cross-section streamlines onto the Schaefer 7-network based atlas with 200 parcels (Schaefer et al., 2018) then log-transformed to adjust for the scale (Fornito et al., 2016). We opted for this atlas as it (i) allows contextualization of our findings within macroscale intrinsic functional communities (Yeo et al., 2011), (ii) incorporates the option to assess results across different granularities, and (iii) aligns the current study with previous work from our group (Benkarim et al., 2020; Paquola et al., 2020; Park et al., 2021b; Park et al., 2021a; Rodríguez-Cruces et al., 2020) and others (Baum et al., 2020; Betzel et al., 2019; Osmanlıoğlu et al., 2019).

Cortex-wide structural connectome manifolds were identified using BrainSpace (https://github.com/MICA-MNI/BrainSpace; Vos de Wael et al., 2020a). First, a template manifold was estimated using a group representative structural connectome of the template dataset. The group representative structural connectome was defined using a distance-dependent thresholding that preserves long-range connections (Betzel et al., 2019). An affinity matrix was constructed with a normalized angle kernel, and eigenvectors were estimated via diffusion map embedding (Figure 1A), a nonlinear dimensionality reduction technique (Coifman and Lafon, 2006) that projects connectome features into low-dimensional manifolds (Margulies et al., 2016). This technique is only controlled by a few parameters, computationally efficient, and relatively robust to noise compared to other nonlinear techniques (Errity and McKenna, 2007; Gallos et al., 2020; Hong et al., 2020; Tenenbaum et al., 2000), and has been extensively used in the previous gradient mapping literature (Hong et al., 2019; Hong et al., 2020; Huntenburg et al., 2017; Larivière et al., 2020a; Margulies et al., 2016; Müller et al., 2020; Paquola et al., 2019a; Park et al., 2021b; Valk et al., 2020; Vos de Wael et al., 2020a). It is controlled by two parameters α and t, where α controls the influence of the density of sampling points on the manifold (α = 0, maximal influence; α = 1, no influence) and t controls the scale of eigenvalues of the diffusion operator. We set α = 0.5 and t = 0 to retain the global relations between data points in the embedded space, following prior applications (Hong et al., 2019; Margulies et al., 2016; Paquola et al., 2019a; Paquola et al., 2019b; Vos de Wael et al., 2020a). Briefly, the eigenvectors estimated from the decomposition technique generate a connectivity coordinate system (Bijsterbosch et al., 2020; Haak et al., 2018; Huntenburg et al., 2018; Margulies et al., 2016; Mars et al., 2018) – the diffusion map, where Euclidean distances in the manifold correspond to diffusion times between the nodes of the network (Coifman and Lafon, 2006). In this manifold space, interconnected brain regions with similar connectivity patterns are closely located, and regions with weak similarity in connectivity patterns are located farther apart. After generating the template manifold, individual-level manifolds were estimated from the non-template dataset and aligned to the template manifold via Procrustes alignment (Langs et al., 2015; Vos de Wael et al., 2020a). To analyze change in the low-dimensional manifold space, we simplified the multivariate eigenvectors into a single scalar value that is., manifold eccentricity (Figure 1B). Manifold eccentricity was calculated as the Euclidean distance between the manifold origin and all data points (i.e., brain regions) in manifold space. The template center was defined as the centroid of the first three eigenvectors, which explained 50% variance. Specifically, manifold eccentricity was defined as follows:

(1) CT=1Ni=1NTE1i, i=1NTE2i, i=1NTE3i
(2) ME=e=13IEe-CTe2

CT is the template manifold origin, N the number of brain regions, T the template manifold, ME the manifold eccentricity, I the individual manifold, and CTe the origin of e-th template manifold. Simply, as shown in Figure 1—figure supplement 15, each brain region (i.e., each dot in the scatter plot) is described as a vector from the manifold origin (i.e., triangular mark in the scatter plot), and manifold eccentricity is simply a length (i.e., Euclidean distance) of that vector. Shifts in connectivity patterns of a given region thus will lead to shifts in the vectors, which in turn changes the manifold eccentricity. Thus, manifold eccentricity quantifies global brain organization based in the connectivity space.

Age-related changes in structural manifolds

Request a detailed protocol

We assessed changes in manifold eccentricity across age using a linear mixed effect model (Worsley et al., 2009), controlling for effects of sex, site, head motion, and subject-specific random intercept to improve model fit in accelerated longitudinal designs. The t-statistics of each brain region were computed, and we corrected for multiple comparisons by using an FDR threshold of q < 0.05 (Figure 1C; Benjamini and Hochberg, 1995). We stratified age-related effects based on a seminal model of neural organization and laminar differentiation that contains four cortical hierarchical levels (Mesulam, 1998), as well as seven intrinsic functional communities (Yeo et al., 2011; Figure 1D). To assess the effects with respect to age2, we repeated implementing a linear mixed effect model by adding a quadratic term of age to the model.

To provide the underlying structure of manifold eccentricity, we compared the changes in manifold eccentricity with those in connectome topology measures. We first defined clusters within the identified regions based on their spatial boundaries (Figure 1—figure supplement 1A). Then, we calculated degree centrality, as well as modular measures of within-module degree and participation coefficient using the Brain Connectivity Toolbox (https://sites.google.com/site/bctnet/) (Rubinov and Sporns, 2010) and connectivity distance using a recently published approach (Larivière et al., 2020b). Degree centrality is defined as the row-wise sum of the weighted connectivity matrix, representing the connection strength of a given node (Rubinov and Sporns, 2010). Connectivity distance is a given brain region’s geodesic distance to its structurally connected brain areas within the cortex (Oligschläger et al., 2019), and it is defined as the multiplication between the geodesic distance and the binarized structural connectome (Hong et al., 2019; Oligschläger et al., 2019). Within-module degree and participation coefficient are nodal measures reflecting different facets of community organization (Rubinov and Sporns, 2010). For each individual subject, community structure was defined using Louvain’s algorithm (Blondel et al., 2008) and a consistency matrix was constructed, where each element of the matrix represents whether the two different nodes are involved in the same community (i.e., 1) or not (i.e., 0) (Figure 1—figure supplement 2A). We constructed the group-wise consistency matrix by averaging the consistency matrix of all subjects and applied k-means clustering (Figure 1—figure supplement 2B). The optimal number of clusters was determined using the silhouette coefficient, that is, the k that maximized the silhouette coefficient (Kannan et al., 2010). We calculated within-module degree and participation coefficient based on these modules. Within-module degree is the degree centrality within a community, indicating intra-community connection strength, while participation coefficient represents inter-community connectivity (Rubinov and Sporns, 2010). We calculated linear correlations between changes in manifold eccentricity and those in each graph-theoretical measure for each cluster (Figure 1—figure supplement 1B). The significance of the correlation was corrected using 1000 permutation tests by randomly shuffling subject indices in one of the data vectors, and we corrected for multiple comparisons across clusters using an FDR procedure (Benjamini and Hochberg, 1995). To visualize age-related changes in these parameters, we stratified each measure according to discretized age bins (<17, 17–19, 19–21, 21–23, ≥23; Figure 1—figure supplement 3).

Cortical morphology and microstructure

Request a detailed protocol

It has been shown that macroscale cortical morphology and microstructure significantly change during development (Gogtay et al., 2004; Khundrakpam et al., 2017; Paquola et al., 2019a; Shaw et al., 2006). Here, we confirmed these changes by assessing age-related changes in MRI-based cortical thickness measures and intracortical measures of MT, an index sensitive to myelin content (Weiskopf et al., 2013), using linear mixed effect models (Figure 2A; Worsley et al., 2009). We further regressed out cortical thickness and MT from the connectome manifold eccentricity metric. We then implemented linear mixed effect models using the residuals of manifold measures to assess whether age-related connectome manifold effects exist above and beyond age-related effects on cortical morphology and microstructure (Figure 2B).

Subcortico-cortical connectivity

Request a detailed protocol

To assess age-related changes in subcortical manifold organizations in addition to cortical manifold structures, we first parcellated the accumbens, amygdala, caudate, hippocampus, pallidum, putamen, and thalamus for each individual (Patenaude et al., 2011), and approximated cross-sectionl streamlines connect each subcortical region to the rest of the brain. For each individual and each subcortical region, we projected the streamline strength to cortical manifold space by weighting the cortical manifolds with the streamline strength of the connection between each subcortical region and cortical parcels, yielding a matrix with the form of (number of brain regions × number of cortical manifolds). We averaged the matrix across the axis of cortical manifolds to construct subcortical-weighted manifold vector. We assessed age-related changes in the subcortical-weighted manifold using a linear mixed effect model (Worsley et al., 2009), controlling for sex, site, head motion, and subject-specific random intercept, and FDR corrected for multiple comparisons (Figure 3Benjamini and Hochberg, 1995).

Transcriptomic analysis

Request a detailed protocol

We performed spatial correlation analysis to post-mortem gene expression data and carried out a developmental enrichment analysis (Figure 4). In brief, we first correlated the t-statistics map, which represents age-related changes in manifold eccentricity that controlled for cortical morphology and microstructure, with the post-mortem gene expression maps provided by the Allen Institute using the Neurovault gene decoding tool (Gorgolewski et al., 2015; Gorgolewski et al., 2014; Hawrylycz et al., 2012). Leveraging mixed effect models to associate the input t-statistic map with the genes of six donor brains, Neurovault yields the gene symbols associated with the input spatial map. Gene symbols that passed for a significance level of FDR < 0.05 were further tested whether they are consistently expressed across different donors using abagen toolbox (Markello et al., 2020; copy archived at https://github.com/rmarkello/abagenArnatkeviciute et al., 2019; Hawrylycz et al., 2012). For each gene, we estimated whole-brain gene expression map and correlated it between all pairs of donors. Leveraging CSEA developmental expression tool (http://genetics.wustl.edu/jdlab/csea-tool-2; Dougherty et al., 2010; Xu et al., 2014), we evaluated the significance of overlap between the genes showing consistent whole-brain expression pattern across donors (FDR < 0.05) with RNAseq data obtained from BrainSpan dataset (http://www.brainspan.org). The significance was calculated based on Fisher’s exact test (Fisher, 1922) with FDR correction (Benjamini and Hochberg, 1995). The CSEA tool provides simplified results of gene enrichment profiles along six major brain regions (i.e., cortex, thalamus, striatum, cerebellum, hippocampus, amygdala) across 10 developmental periods (from early fetal to young adulthood) approximated from mouse data, yielding a total of 60 combinations of developmental enrichment profiles (Xu et al., 2014). We repeated developmental enrichment analysis using the genes identified from the rotated maps of the age-related changes in manifold eccentricity (100 spherical rotations). For each iteration, we obtained developmental expression profiles using the identified genes, where the FDR-corrected p-values built a null distribution. For each brain division and developmental period, if the actual p-value is placed outside 95% of the null distribution, it was deemed significant. As the Allen Brain Institute repository is composed of adult post-mortem datasets, it should be noted that the associated gene symbols represent indirect associations with the input t-statistic map derived from the developmental data.

Association with the development of cognitive function

Request a detailed protocol

Leveraging a supervised machine learning with ten-fold cross-validation, we predicted full IQ score measured by the Wechsler Abbreviated Scale of Intelligence (Wechsler, 1999) at follow-up using cortical and subcortical features. Four different feature sets were evaluated: (i) manifold eccentricity of the identified cortical regions at baseline and (ii) manifold eccentricity at baseline and its longitudinal change (i.e., differences between follow-up and baseline), and (iii) cortical manifold eccentricity and subcortical-weighted manifold of the identified regions at baseline and (iv) manifold eccentricity and subcortical-weighted manifold at baseline and their longitudinal changes. For each evaluation, a subset of features that could predict future IQ was identified using elastic net regularization (ρ=0.5) with optimized regularization parameters (L1 and L2 penalty terms) via nested ten-fold cross-validation (Cawley and Talbot, 2010; Parvandeh et al., 2020; Tenenbaum et al., 2000; Varma and Simon, 2006; Zou and Hastie, 2005). We split the dataset into training (9/10) and test (1/10) partitions, and each training partition was further split into inner training and testing folds using another ten-fold cross-validation. Within the inner fold, elastic net regularization finds a set of non-redundant features to explain the dependent variable. Using a linear regression, we predicted the IQ scores of inner fold test data using the features of the selected brain regions by controlling for age, sex, site, and head motion. The model with minimum MAE across the inner folds was applied to the test partition of the outer fold, and the IQ scores of outer fold test data were predicted. The prediction procedure was repeated 100 times with different training and test sets to reduce subject selection bias. Prediction accuracy was indexed by computing linear correlations between the actual and predicted IQ scores as well as MAE. A 95% confidence interval of the accuracy measures was also reported. Permutation-based correlations across 1000 tests were conducted by randomly shuffling subject indices to check whether the prediction performance exceeded chance levels. To assess whether our model outperforms baseline model, we predicted IQ of test data using average of IQ of training data (i.e., predicted IQ = mean(training set IQ)). The improvement of prediction performance was assessed using Meng’s z-test (Meng et al., 1992). In addition to predicting future IQ, we performed the same prediction analysis to predict the change of IQ between the baseline and follow-up.

Sensitivity analysis

Spatial scale

Request a detailed protocol

To assess the consistency of our findings across spatial scales, we additionally performed the linear mixed effect modeling using a finer parcellation scheme of 300 parcels (Figure 1—figure supplement 4; Schaefer et al., 2018).

Site and sex effect

Request a detailed protocol

Participants were recruited from three different sites. To assess whether the longitudinal changes in manifold eccentricity across age are consistent across different sites, we calculated interaction effects of the relationship between age and manifold eccentricity of the identified regions across sites (Figure 1—figure supplement 5B). In addition, we computed interaction effect of the relationship between age and manifold eccentricity across male and female subjects to assess whether the age-related changes are affected by biological sexes (Figure 1—figure supplement 5C).

Different parameters for diffusion map embedding

Request a detailed protocol

To assess the sensitivity of our findings, we generated connectome manifolds with different parameters for diffusion map embedding (α = 0.25, 0.5, 0.75; t = 0, 1, 2, 3). We assessed age-related changes of the newly defined manifold eccentricity and calculated linear correlation with t-statistic map of the default setting (α = 0.5; t = 0; Figure 1C).

Gradient alignment fidelity

Request a detailed protocol

To assess robustness of individual alignment, we computed linear correlations between the template and individual manifolds before and after alignment. We also repeated the linear mixed effect modeling after excluding 10% of subjects with the lowest alignment to the template manifold (Figure 1—figure supplement 6).

Connectome manifold generation using principal component analysis

Request a detailed protocol

To explore consistency of our results when using different dimensionality reduction techniques, we generated connectome manifolds using principal component analysis (Wold et al., 1987), instead of relying on diffusion map embedding (Coifman and Lafon, 2006), and performed longitudinal modeling (Figure 1—figure supplement 7). We compared the eigenvectors estimated from diffusion map embedding and principal component analysis using linear correlations.

Longitudinal changes in graph-theoretical measures

Request a detailed protocol

To compare longitudinal changes in manifold eccentricity with those in graph-theoretical centrality measures, we calculated betweenness, degree, and eigenvector centrality of the structural connectomes and built similar linear mixed effects models to assess longitudinal change (Figure 1—figure supplement 8). Betweenness centrality is the number of weighted shortest paths between any combinations of nodes that run through that node, degree centrality is the sum of edge weights connected to a given node, and eigenvector centrality measures the influence of a node in the whole network (Lohmann et al., 2010; Rubinov and Sporns, 2010; Zuo et al., 2012). Spatial similarity between t-statistics of centrality and manifold measures was assessed with 1000 spin tests that account for spatial autocorrelation (Alexander-Bloch et al., 2018).

Manifold eccentricity analysis based on all eigenvectors

Request a detailed protocol

We repeated our analysis by calculating manifold eccentricity from all eigenvectors to assess consistency of the findings (Figure 1—figure supplement 9).

Robustness of group representative structural connectome

Request a detailed protocol

We compared the distance-dependent thresholding (Betzel et al., 2019) that was adopted for the main analysis with a consistency thresholding approach (Wang et al., 2019). The latter averages subject-specific matrices, in addition to performing a 50, 40, 30, 20, and 10% thresholding, as well as simple averaging (i.e., 0% thresholding) (Figure 1—figure supplement 10).

Connectome manifolds based on structural parcellation

Request a detailed protocol

To confirm whether functional and structural parcellation schemes yield consistent results, we repeated our main analyses using 200 cortical nodes structural parcellation scheme, which preserves the macroscopic boundaries of the Desikan–Killiany atlas (Desikan et al., 2006; Vos de Wael et al., 2020aFigure 1—figure supplement 11).

Longitudinal modeling using edge weights

Request a detailed protocol

In addition to the analyses based on manifold eccentricity, linear mixed effect modeling using connectome edge weights assessed age-related longitudinal changes in streamline strength (Figure 1—figure supplement 12).

Manifold eccentricity and pubertal stages

Request a detailed protocol

To assess the relationship between manifold eccentricity and pubertal stages, we selected a subset of participants who completed Tanner scale (Marshall and Tanner, 1970; Marshall and Tanner, 1969), which quantifies pubertal stages from 1 (pre-puberty) to 5 (final phase of physical maturation). However, the score was collected at baseline and for 73/208 participants only. To confirm robustness, we performed linear mixed effect modeling using this subset (Figure 1—figure supplement 13A). In addition, we assessed interaction effects of Tanner scale and manifold eccentricity restricted to the regions identified from the overall sample (Figure 1—figure supplement 13B).

IQ prediction using nonlinear model

Request a detailed protocol

We additionally predicted future IQ score using decision tree learning, a nonlinear approach that builds a regression tree model a root node and split leaf nodes, where the leaf nodes contain the response variables (Breiman et al., 1984; Figure 5—figure supplement 2).

Data and code availability

Request a detailed protocol

The imaging and phenotypic data were provided by the NSPN 2400 cohort. As stated in https://doi.org/10.1093/ije/dyx117, the NSPN project is committed to make the anonymised dataset fully available to the research community, and participants have consented to their de-identified data being made available to other researchers. A data request can be made to openNSPN@medschl.cam.ac.uk. Codes for connectome manifold generation are available at https://doi.org/10.1038/s42003-020-0794-7; https://github.com/MICA-MNI/BrainSpace (copy archived at swh:1:rev:1fb001f4961d3c0b05b7715f42bcc362b31b96a5; Vos de Wael et al., 2020b), and those for calculating manifold eccentricity and subcortical-weighted manifold, as well as performing linear mixed effect modeling to assess age-effects on these features, at out GitHub (https://github.com/MICA-MNI/micaopen/tree/master/manifold_features; copy archived at swh:1:rev:d3988d51e01940007595761dab6b846ce2506433; Park, 2021). Source data are provided with this paper.

Data availability

The imaging and phenotypic data were provided by the NSPN 2400 cohort. As stated in https://doi.org/10.1093/ije/dyx117, the NSPN project is committed to make the anonymised dataset fully available to the research community, and participants have consented to their de-identified data being made available to other researchers. A data request can be made to openNSPN@medschl.cam.ac.uk. Codes for connectome manifold generation are available at https://doi.org/10.1038/s42003-020-0794-7; https://github.com/MICA-MNI/BrainSpace (copy archived at https://archive.softwareheritage.org/swh:1:rev:1fb001f4961d3c0b05b7715f42bcc362b31b96a5/), and those for calculating manifold eccentricity and subcortical-weighted manifold, as well as performing linear mixed effect modeling to assess age-effects on these features at our GitHub (https://github.com/MICA-MNI/micaopen/tree/master/manifold_features; copy archived at https://archive.softwareheritage.org/swh:1:rev:d3988d51e01940007595761dab6b846ce2506433/).

References

  1. Book
    1. Breiman L
    2. Friedman J
    3. Stone CJ
    4. Olshen RA
    (1984)
    Classification and Regression Trees
    Wadsworth, Inc.
    1. Carbon S
    2. Douglass E
    3. Dunn N
    4. Good B
    5. Harris NL
    6. Lewis SE
    7. Mungall CJ
    8. Basu S
    9. Chisholm RL
    10. Dodson RJ
    11. Hartline E
    12. Fey P
    13. Thomas PD
    14. Albou LP
    15. Ebert D
    16. Kesling MJ
    17. Mi H
    18. Muruganujan A
    19. Huang X
    20. Poudel S
    21. Mushayahama T
    22. Jc H
    23. LaBonte SA
    24. Siegele DA
    25. Antonazzo G
    26. Attrill H
    27. Brown NH
    28. Fexova S
    29. Garapati P
    30. Jones TEM
    31. Marygold SJ
    32. Millburn GH
    33. Rey AJ
    34. Trovisco V
    35. Dos Santos G
    36. Emmert DB
    37. Falls K
    38. Zhou P
    39. Goodman JL
    40. Strelets VB
    41. Thurmond J
    42. Courtot M
    43. Osumi DS
    44. Parkinson H
    45. Roncaglia P
    46. Acencio ML
    47. Kuiper M
    48. Lreid A
    49. Logie C
    50. Lovering RC
    51. Huntley RP
    52. Denny P
    53. Campbell NH
    54. Kramarz B
    55. Acquaah V
    56. Ahmad SH
    57. Chen H
    58. Rawson JH
    59. Chibucos MC
    60. Giglio M
    61. Nadendla S
    62. Tauber R
    63. Duesbury MJ
    64. Del NT
    65. Meldal BHM
    66. Perfetto L
    67. Porras P
    68. Orchard S
    69. Shrivastava A
    70. Xie Z
    71. Chang HY
    72. Finn RD
    73. Mitchell AL
    74. Rawlings ND
    75. Richardson L
    76. Sangrador-Vegas A
    77. Blake JA
    78. Christie KR
    79. Dolan ME
    80. Drabkin HJ
    81. Hill DP
    82. Ni L
    83. Sitnikov D
    84. Harris MA
    85. Oliver SG
    86. Rutherford K
    87. Wood V
    88. Hayles J
    89. Bahler J
    90. Lock A
    91. Bolton ER
    92. De Pons J
    93. Dwinell M
    94. Hayman GT
    95. Laulederkind SJF
    96. Shimoyama M
    97. Tutaj M
    98. Wang SJ
    99. D’Eustachio P
    100. Matthews L
    101. Balhoff JP
    102. Aleksander SA
    103. Binkley G
    104. Dunn BL
    105. Cherry JM
    106. Engel SR
    107. Gondwe F
    108. Karra K
    109. MacPherson KA
    110. Miyasato SR
    111. Nash RS
    112. Pc N
    113. Sheppard TK
    114. Shrivatsav Vp A
    115. Simison M
    116. Skrzypek MS
    117. Weng S
    118. Wong ED
    119. Feuermann M
    120. Gaudet P
    121. Bakker E
    122. Berardini TZ
    123. Reiser L
    124. Subramaniam S
    125. Huala E
    126. Arighi C
    127. Auchincloss A
    128. Axelsen K
    129. Argoud GP
    130. Bateman A
    131. Bely B
    132. Blatter MC
    133. Boutet E
    134. Breuza L
    135. Bridge A
    136. Britto R
    137. Bye-A-Jee H
    138. Casals-Casas C
    139. Coudert E
    140. Estreicher A
    141. Famiglietti L
    142. Garmiri P
    143. Georghiou G
    144. Gos A
    145. Gruaz-Gumowski N
    146. Hatton-Ellis E
    147. Hinz U
    148. Hulo C
    149. Ignatchenko A
    150. Jungo F
    151. Keller G
    152. Laiho K
    153. Lemercier P
    154. Lieberherr D
    155. Lussi Y
    156. Mac-Dougall A
    157. Magrane M
    158. Martin MJ
    159. Masson P
    160. Natale DA
    161. Hyka NN
    162. Pedruzzi I
    163. Pichler K
    164. Poux S
    165. Rivoire C
    166. Rodriguez-Lopez M
    167. Sawford T
    168. Speretta E
    169. Shypitsyna A
    170. Stutz A
    171. Sundaram S
    172. Tognolli M
    173. Tyagi N
    174. Warner K
    175. Zaru R
    176. Wu C
    177. Chan J
    178. Cho J
    179. Gao S
    180. Grove C
    181. Harrison MC
    182. Howe K
    183. Lee R
    184. Mendel J
    185. Muller HM
    186. Raciti D
    187. Van Auken K
    188. Berriman M
    189. Stein L
    190. Sternberg PW
    191. Howe D
    192. Toro S
    193. Westerfield M
    194. The Gene Ontology Consortium
    (2019) The gene ontology resource: 20 years and still GOing strong
    Nucleic Acids Research 47:D330–D338.
    https://doi.org/10.1093/nar/gky1055
    1. Cawley GC
    2. Talbot NLC
    (2010)
    On over-fitting in model selection and subsequent selection Bias in performance evaluation
    Journal of Machine Learning Research : JMLR 11:2079–2107.
    1. Coifman RR
    2. Lafon S
    (2006) Diffusion maps
    Applied and Computational Harmonic Analysis 21:5–30.
    https://doi.org/10.1016/j.acha.2006.04.006
    1. Patel Y
    2. Parker N
    3. Shin J
    4. Howard D
    5. French L
    6. Thomopoulos SI
    7. Pozzi E
    8. Abe Y
    9. Abé C
    10. Anticevic A
    11. Alda M
    12. Aleman A
    13. Alloza C
    14. Alonso-Lana S
    15. Ameis SH
    16. Anagnostou E
    17. McIntosh AA
    18. Arango C
    19. Arnold PD
    20. Asherson P
    21. Assogna F
    22. Auzias G
    23. Ayesa-Arriola R
    24. Bakker G
    25. Banaj N
    26. Banaschewski T
    27. Bandeira CE
    28. Baranov A
    29. Bargalló N
    30. Bau CHD
    31. Baumeister S
    32. Baune BT
    33. Bellgrove MA
    34. Benedetti F
    35. Bertolino A
    36. Boedhoe PSW
    37. Boks M
    38. Bollettini I
    39. Del Mar Bonnin C
    40. Borgers T
    41. Borgwardt S
    42. Brandeis D
    43. Brennan BP
    44. Bruggemann JM
    45. Bülow R
    46. Busatto GF
    47. Calderoni S
    48. Calhoun VD
    49. Calvo R
    50. Canales-Rodríguez EJ
    51. Cannon DM
    52. Carr VJ
    53. Cascella N
    54. Cercignani M
    55. Chaim-Avancini TM
    56. Christakou A
    57. Coghill D
    58. Conzelmann A
    59. Crespo-Facorro B
    60. Cubillo AI
    61. Cullen KR
    62. Cupertino RB
    63. Daly E
    64. Dannlowski U
    65. Davey CG
    66. Denys D
    67. Deruelle C
    68. Di Giorgio A
    69. Dickie EW
    70. Dima D
    71. Dohm K
    72. Ehrlich S
    73. Ely BA
    74. Erwin-Grabner T
    75. Ethofer T
    76. Fair DA
    77. Fallgatter AJ
    78. Faraone SV
    79. Fatjó-Vilas M
    80. Fedor JM
    81. Fitzgerald KD
    82. Ford JM
    83. Frodl T
    84. Fu CHY
    85. Fullerton JM
    86. Gabel MC
    87. Glahn DC
    88. Roberts G
    89. Gogberashvili T
    90. Goikolea JM
    91. Gotlib IH
    92. Goya-Maldonado R
    93. Grabe HJ
    94. Green MJ
    95. Grevet EH
    96. Groenewold NA
    97. Grotegerd D
    98. Gruber O
    99. Gruner P
    100. Guerrero-Pedraza A
    101. Gur RE
    102. Gur RC
    103. Haar S
    104. Haarman BCM
    105. Haavik J
    106. Hahn T
    107. Hajek T
    108. Harrison BJ
    109. Harrison NA
    110. Hartman CA
    111. Whalley HC
    112. Heslenfeld DJ
    113. Hibar DP
    114. Hilland E
    115. Hirano Y
    116. Ho TC
    117. Hoekstra PJ
    118. Hoekstra L
    119. Hohmann S
    120. Hong LE
    121. Höschl C
    122. Høvik MF
    123. Howells FM
    124. Nenadic I
    125. Jalbrzikowski M
    126. James AC
    127. Janssen J
    128. Jaspers-Fayer F
    129. Xu J
    130. Jonassen R
    131. Karkashadze G
    132. King JA
    133. Kircher T
    134. Kirschner M
    135. Koch K
    136. Kochunov P
    137. Kohls G
    138. Konrad K
    139. Krämer B
    140. Krug A
    141. Kuntsi J
    142. Kwon JS
    143. Landén M
    144. Landrø NI
    145. Lazaro L
    146. Lebedeva IS
    147. Leehr EJ
    148. Lera-Miguel S
    149. Lesch KP
    150. Lochner C
    151. Louza MR
    152. Luna B
    153. Lundervold AJ
    154. MacMaster FP
    155. Maglanoc LA
    156. Malpas CB
    157. Portella MJ
    158. Marsh R
    159. Martyn FM
    160. Mataix-Cols D
    161. Mathalon DH
    162. McCarthy H
    163. McDonald C
    164. McPhilemy G
    165. Meinert S
    166. Menchón JM
    167. Minuzzi L
    168. Mitchell PB
    169. Moreno C
    170. Morgado P
    171. Muratori F
    172. Murphy CM
    173. Murphy D
    174. Mwangi B
    175. Nabulsi L
    176. Nakagawa A
    177. Nakamae T
    178. Namazova L
    179. Narayanaswamy J
    180. Jahanshad N
    181. Nguyen DD
    182. Nicolau R
    183. O'Gorman Tuura RL
    184. O'Hearn K
    185. Oosterlaan J
    186. Opel N
    187. Ophoff RA
    188. Oranje B
    189. García de la Foz VO
    190. Overs BJ
    191. Paloyelis Y
    192. Pantelis C
    193. Parellada M
    194. Pauli P
    195. Picó-Pérez M
    196. Picon FA
    197. Piras F
    198. Piras F
    199. Plessen KJ
    200. Pomarol-Clotet E
    201. Preda A
    202. Puig O
    203. Quidé Y
    204. Radua J
    205. Ramos-Quiroga JA
    206. Rasser PE
    207. Rauer L
    208. Reddy J
    209. Redlich R
    210. Reif A
    211. Reneman L
    212. Repple J
    213. Retico A
    214. Richarte V
    215. Richter A
    216. Rosa PGP
    217. Rubia KK
    218. Hashimoto R
    219. Sacchet MD
    220. Salvador R
    221. Santonja J
    222. Sarink K
    223. Sarró S
    224. Satterthwaite TD
    225. Sawa A
    226. Schall U
    227. Schofield PR
    228. Schrantee A
    229. Seitz J
    230. Serpa MH
    231. Setién-Suero E
    232. Shaw P
    233. Shook D
    234. Silk TJ
    235. Sim K
    236. Simon S
    237. Simpson HB
    238. Singh A
    239. Skoch A
    240. Skokauskas N
    241. Soares JC
    242. Soreni N
    243. Soriano-Mas C
    244. Spalletta G
    245. Spaniel F
    246. Lawrie SM
    247. Stern ER
    248. Stewart SE
    249. Takayanagi Y
    250. Temmingh HS
    251. Tolin DF
    252. Tomecek D
    253. Tordesillas-Gutiérrez D
    254. Tosetti M
    255. Uhlmann A
    256. van Amelsvoort T
    257. van der Wee NJA
    258. van der Werff SJA
    259. van Haren NEM
    260. van Wingen GA
    261. Vance A
    262. Vázquez-Bourgon J
    263. Vecchio D
    264. Venkatasubramanian G
    265. Vieta E
    266. Vilarroya O
    267. Vives-Gilabert Y
    268. Voineskos AN
    269. Völzke H
    270. von Polier GG
    271. Walton E
    272. Weickert TW
    273. Weickert CS
    274. Weideman AS
    275. Wittfeld K
    276. Wolf DH
    277. Wu MJ
    278. Yang TT
    279. Yang K
    280. Yoncheva Y
    281. Yun JY
    282. Cheng Y
    283. Zanetti MV
    284. Ziegler GC
    285. Franke B
    286. Hoogman M
    287. Buitelaar JK
    288. van Rooij D
    289. Andreassen OA
    290. Ching CRK
    291. Veltman DJ
    292. Schmaal L
    293. Stein DJ
    294. van den Heuvel OA
    295. Turner JA
    296. van Erp TGM
    297. Pausova Z
    298. Thompson PM
    299. Paus T
    300. Writing Committee for the Attention-Deficit/Hyperactivity Disorder
    301. Autism Spectrum Disorder
    302. Bipolar Disorder
    303. Major Depressive Disorder
    304. Obsessive-Compulsive Disorder
    305. and Schizophrenia ENIGMA Working Groups
    (2021) Virtual histology of cortical thickness and shared neurobiology in 6 psychiatric disorders
    JAMA Psychiatry 78:47.
    https://doi.org/10.1001/jamapsychiatry.2020.2694
  2. Book
    1. Sanides F
    (1962)
    Die Architektonik des Menschlichen Stirnhirns, Monographien aus dem Gesamtgebiete der Neurologie und Psychiatrie
    Berlin, Heidelberg: Springer Berlin Heidelberg.
  3. Book
    1. Wechsler D
    (1999)
    Wechsler Abbreviated Scales of Intelligence (WASI
    San Antonio, TX: Psychological Corporation.

Decision letter

  1. Lucina Q Uddin
    Reviewing Editor; University of Miami University, United States
  2. Timothy E Behrens
    Senior Editor; University of Oxford, United Kingdom
  3. Ben D Fulcher
    Reviewer; University of Sydney, Australia

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This manuscript describes a longitudinal study of the adolescent structural connectome. The authors find strong effects of expansion of structural connectomes in transmodal brain regions during adolescence. They also report findings centered on the caudate and thalamus, and supplement the structural connectivity analyses with transcriptome association analyses revealing genes enriched in specific brain regions. Finally, intelligence measures are predicted from baseline structural measures. These findings help to bridge large-scale brain network configurations, microscale processes, and cognition in adolescent development.

Decision letter after peer review:

Thank you for submitting your article "An expanding manifold characterizes adolescent reconfigurations of structural connectome organization" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Timothy Behrens as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Ben Fulcher (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

Summary:

This manuscript describes a longitudinal study of the adolescent structural connectome. Park et al. report on an analysis of existing semi-longitudinal NSPN 2400 data to learn how the projections of high-dimensional structural connectivity patterns onto a three-dimensional subspace change with age during adolescence. They employ a non-linear manifold learning algorithm (diffusion embedding), thereby linking the maturation of global structural connectivity patterns to an emerging approach in understanding brain organization through spatial gradient representations. The authors find strong effects of expansion of structural connectomes in transmodal brain regions during adolescence. They also report findings centered on the caudate and thalamus, and supplement the structural connectivity analyses with transcriptome association analyses revealing gens enriched in specific brain regions. Finally, intelligence measures are predicted from baseline structural measures.

This is an interesting and comprehensive set of analyses on an important topic. Overall, the figures are lovely. The sensitivity analyses are particularly commendable.

The paper is well written, the data are fantastic, and the analyses are interesting. Some suggestions and points for clarification (both theoretical and methodological) are below.

Essential revisions:

– There is not much in the introduction about why co-localized gene sets are of interest to explore. What is already known about brain development using this approach, and how does the current work fill a gap in our knowledge?

– Similarly, the introduction states that the study aims to "predict future measures of cognitive function". What cognitive functions specifically were of interest in this study, and why? No rationale or background is provided for conducting these analyses.

– The authors claim that their study examines "the entire adolescent time period", however some would argue that age 14 does not represent the earliest age at which adolescence onsets. I think it would be more accurate to say the study covers the mid to late adolescent period.

– In the results it is stated that three eigenvector explained approximately 50% of the variance in the template affinity matrix. Here it would be helpful to report exactly how much of the variance was explained by each (E1, E2, E3).

– Pubertal development occurs across the age range investigated, and affects brain structure and function. Was information on pubertal stage of participants available? Did some participants undergo changes in pubertal status from timepoint 1 to timepoint 2?

– The introduction does not mention cortical thickness much, therefore these analyses come as a bit of surprise in the results.

– As in the introduction, there is not much interpretation of the transcriptome findings in the discussion.

– For constructing the structural connectome, the Schaefer 7-network atlas was utilized. Can the authors comment on why a functional atlas (rather than a structural atlas) was used here?

– While this work touches on an important topic, ties nicely with the increasing body of papers on global brain gradients, and its overall conclusions are warranted, I am not (yet) convinced that it offers fundamentally new insights that could not have been gleaned from previous work (after all, manifold learning simply displays a shadow of the underlying patterns; if the patterns change, so does their shadow). I am also not convinced by the rationale for employing diffusion embedding: the authors state that the ensuing gradients are heritable, conserved across species, capture functional activation patterns during task states, and provide a coordinate system to interrogate brain structure and function, but that would be true for any method that adequately captures biologically meaningful variance in the structural connectivity patterns.

– The authors show that the maturational change of the manifold features predict intelligence at follow-up but did not show that intelligence itself exhibited changes that exceeded the error bounds of the regression line. Why not predict IQ change?

– The slight improvements in prediction accuracy observed after adding maturational change and subcortical features to the features at baseline will necessarily happen by adding more regression parameters and may not be meaningful.

– Interpretability and Lack of Comparison

The authors claim repeatedly that they are "capitalizing on advanced manifold learning techniques".

One could imagine an infinite number of papers that take a dataset, use a technique to extract a metric, X (e.g., eccentricity), and then write about the changes in X with some property of interest, Y (e.g., age). Given this set of papers (and the non-independence between the set of possible Xs), the reader ought to be most interested in those Xs that provide the best performance and simplest interpretation, with other papers being redundant.

Thus, a nuanced approach to presenting a paper like this is to demonstrate that the metric used represents an advance over alternative, simpler-to-compute, or clearer-to-interpret metrics that already exist.

In this paper, however, the authors do not demonstrate the benefits of their particular choice of applying a specific nonlinear dimensionality reduction method using 3 dimensions alignment to a template manifold and then computing an eccentricity metric. For example:

i. Is the nonlinearity required (e.g., does it outperform PCA or MDS)?

ii. Is there something special about picking 3 dimensions to do the eccentricity calculation? Is dimensionality reduction required at all (e.g., would you get similar results by computing eccentricity in the full-dimensional space?)

iii. Does it outperform basic connectome measures (e.g., the simple ones the authors compute)?

There is a clear down-side of how opaque the approach is (and thus difficult to interpret relative to, say, connectivity degree), so one would hope for a correspondingly strong boost in performance. The authors could also do more to develop some intuition for the idea of a low-dimensional connection-pattern-similarity-space, and how to interpret taking Euclidean distances within such a space.

– Developmental Enrichment Analysis

Both in the main text and in the Materials and methods, this is described as "genes were fed into a developmental enrichment analysis". Can some explanation be provided as to what happens between the "feeding in" and what comes out? Without clearly described methods, it is impossible to interpret or critique this component of the paper. If the methodological details are opaque, then the significance of the results could be tested numerically relative to some randomized null inputs being 'fed in' to demonstrate specificity of the tested phenotype.

– IQ prediction

The predictions seem to be very poor (equality lines, y = x, should be drawn in Figure 5, to show what perfect predictions would look like; linear regressions are not helpful for a prediction task, and are deceptive of the appropriate MAE computation). The authors do not perform any comparisons in this section (even to a real baseline model like `predicted_IQ = mean(training_set_IQ)`). They also do not perform statistical tests (or quote p-values), but nevertheless make a range of claims, including of "significant prediction" or "prediction accuracy was improved", "reemphasize the benefits of incorporating subcortical nodes", etc. All of these claims should be tested relative to rigorous statistics, and comparisons to appropriate baseline/benchmark approaches.

– Group Connectome

Given how much the paper relies on estimating a group structural connectome, it should be visualized and characterized. For example, a basic analysis of the distribution of edge weights and degree, especially as edge weights can vary over orders of magnitude and high weights (more likely to be short distances) may therefore unduly dominate some of the low-dimensional components. The authors may also consider testing robustness performed to alternative ways of estimating the connectome [e.g., Oldham et al. NeuroImage 222, 117252 (2020)] and its group-level summary [e.g., Roberts et al. NeuroImage 145, 1-42 (2016)].

– Individual Alignment

The paper relies on individuals being successfully aligned to the template manifold. Accordingly, some analysis should be performed quantifying how well individuals could be mapped. Presumably some subjects fit very well onto the template, whereas others do not. Is there something interesting about the poorly aligned subjects? Do your results improve when excluding them?

https://doi.org/10.7554/eLife.64694.sa1

Author response

Essential revisions:

– There is not much in the introduction about why co-localized gene sets are of interest to explore. What is already known about brain development using this approach, and how does the current work fill a gap in our knowledge?

We updated Introduction by adding a new paragraph motivating the transcriptomic analysis and outlining the novelty of our work:

“Imaging-transcriptomics approaches allow for the identification of cellular and molecular factors that co-vary with imaging-based findings (Arnatkeviciute et al., 2019; Fornito et al., 2019; Gorgolewski et al., 2014; Hawrylycz et al., 2015; Thompson et al., 2013). Recently established resources, such as the Allen Human Brain Atlas (Arnatkeviciute et al., 2019; Hawrylycz et al., 2015), can be utilized to spatially associate macroscale imaging/connectome data with the expression patterns of thousands of genes. These findings have already been applied in the study of healthy adults (Hawrylycz et al., 2015; Park et al., 2020b; Q. Xu et al., 2020) and typically developing adolescents (Mascarell Maričić et al., 2020; Padmanabhan and Luna, 2014; Paquola et al., 2019a; Vértes et al., 2016; Whitaker et al., 2016), as well as individuals suffering from prevalent brain disorders (Altmann et al., 2018; Hashimoto et al., 2015; Klein et al., 2017; Park et al., 2020a; Patel et al., 2021; Romero-Garcia et al., 2019). The gene sets that co-vary with in vivo findings can furthermore be subjected to gene set enrichment analyses to discover potentially implicated molecular, cellular, and pathological processes (Ashburner et al., 2000; Carbon et al., 2019; Chen et al., 2013; Dougherty et al., 2010; Kuleshov et al., 2016; Morgan et al., 2019; Romero-Garcia et al., 2018; Subramanian et al., 2005). For example, studies in newborns have shown that cortical morphology reflects spatiotemporal patterns of gene expression in foetuses, linking molecular mechanisms to in vivo measures of cortical development in early life (Ball et al., 2020). Work in adolescents has furthermore shown that developmental changes in regional cortical thickness measures and myelin proxies spatially co-localize with the expression patterns of genes involved in synaptic and oligodendroglial function (Paquola et al., 2019a; Whitaker et al., 2016). Building on these prior investigations, the current study aimed at exploring whether adolescent structural connectome reconfigurations, assessed using manifold learning techniques, reflect the expression patterns of specific genes, in order to identify potential molecular signatures of macroscale structural network development.”

– Similarly, the introduction states that the study aims to "predict future measures of cognitive function". What cognitive functions specifically were of interest in this study, and why? No rationale or background is provided for conducting these analyses.

We thank the Reviewers for this comment.

We incorporated additional justification and details on the cognitive prediction in the Introduction:

“To also assess behavioral associations of connectome manifold changes, we utilized supervised machine learning to predict future measures of cognitive function quantified via the intelligence quotient (IQ). IQ is a widely used marker of general cognitive abilities, which shows good test-retest reliability (Brown and May, 1979; Canivez and Watkins, 1998; Catron, 1978; Matarazzo et al., 1973; Snow et al., 1989; Wagner and Caldwell, 1979) and has previously been applied to index overall cognitive function during development (Crespi, 2016; Garde et al., 2005, 2000; Koenis et al., 2018; Park et al., 2016; Ramsden et al., 2011; Shaw et al., 2006; Suprano et al., 2020). In the study of neurodevelopment, neuroimaging reports have previously assessed associations between IQ and large-scale network measures in children to adolescents (Koenis et al., 2018; Ramsden et al., 2011; Seidlitz et al., 2018; Shaw et al., 2006; Suprano et al., 2020).”

– The authors claim that their study examines "the entire adolescent time period", however some would argue that age 14 does not represent the earliest age at which adolescence onsets. I think it would be more accurate to say the study covers the mid to late adolescent period.

We appreciate the Reviewer’s comment.

We clarified the sample in the Introduction:

“Here, we charted developmental changes in structural connectome organization, based on an accelerated longitudinal neuroimaging study involving 208 participants investigated between 14 to 26 years of age (Kiddle et al., 2018; Whitaker et al., 2016)”

– In the results it is stated that three eigenvector explained approximately 50% of the variance in the template affinity matrix. Here it would be helpful to report exactly how much of the variance was explained by each (E1, E2, E3).

As suggested, we provided further details in the Results:

“Three eigenvectors (E1, E2, and E3) explained approximately 50% of information in the template affinity matrix (i.e., 20.7/15.8/13.5% for E1/E2/E3, respectively), with each eigenvector showing a different axis of spatial variation across the cortical mantle (Figure 1A).”

– Pubertal development occurs across the age range investigated, and affects brain structure and function. Was information on pubertal stage of participants available? Did some participants undergo changes in pubertal status from timepoint 1 to timepoint 2?

The NSPN consortium collected Tanner scale (Marshall and Tanner, 1970, 1969), which quantifies pubertal stages from 1 (pre-puberty) to 5 (final phase of physical maturation). However, the score was collected at baseline and for 73/208 participants only. We could thus not assess longitudinal changes in manifold features according to changes in pubertal stages indexed by Tanner staging. Nevertheless, linear mixed effect modeling in the subset of these 73 participants confirmed our overall findings. While the statistical power was reduced in this smaller subset, spatial patterns of age-effects were consistent with the larger sample (Figure 1—figure supplement 13A). Notably, manifold eccentricity within the identified regions derived from overall sample and Tanner scale revealed significant interaction effect (t = 2.36, p = 0.01; Figure 1—figure supplement 13B), suggesting participants in early pubertal stage show stronger changes in manifold eccentricity across age compared to those in later stages.

We updated the Results:

“k) Manifold eccentricity and pubertal stages. We repeated the longitudinal modeling within a subset of participants who completed the Tanner scale (n = 73) (Marshall and Tanner, 1970, 1969), and found relatively consistent albeit weaker age-related changes in manifold eccentricity as for the overall sample (Figure 1—figure supplement 13A). Notably, manifold eccentricity within the identified regions derived from overall sample and Tanner scale revealed a significant interaction effect (t = 2.36, p = 0.01; Figure 1—figure supplement 13B), suggesting participants in early pubertal stages show more marked changes in manifold eccentricity across age compared to those in later stages.”

as well as Materials and methods:

“k) Manifold eccentricity and pubertal stages. To assess the relationship between manifold eccentricity and pubertal stages, we selected a subset of participants who completed Tanner scale (Marshall and Tanner, 1970, 1969), which quantifies pubertal stages from 1 (pre-puberty) to 5 (final phase of physical maturation). However, the score was collected at baseline and for 73/208 participants only. To confirm robustness, we performed linear mixed effect modeling using this subset (Figure 1—figure supplement 13A). In addition, we assessed interaction effects of Tanner scale and manifold eccentricity restricted to the regions identified from the overall sample (Figure 1—figure supplement 13B).”

– The introduction does not mention cortical thickness much, therefore these analyses come as a bit of surprise in the results.

As suggested, we now contextualize adolescent cortical thickness changes in the Introduction:

“Early cross-sectional and longitudinal studies in neurodevelopmental cohorts focused on the analysis of morphological changes (Gogtay et al., 2004; Shaw et al., 2006; Tamnes et al., 2017), including MRI-based cortical thickness (Shaw et al., 2006; Tamnes et al., 2017) and volumetric measures (Gogtay et al., 2004; Tamnes et al., 2017). Studies robustly show initial grey matter increases until mid-late childhood followed by a decline for the rest of the lifespan. During adolescence, cortical thickness decreases in widespread brain regions (Khundrakpam et al., 2013; Shaw et al., 2006; Sotiras et al., 2017; Tamnes et al., 2017). Thus, contextualizing connectome alterations relative to established patterns of cortical thickness findings may establish whether inter-regional network changes occur above and beyond these diffuse effects of regional morphological maturation.”

– As in the introduction, there is not much interpretation of the transcriptome findings in the discussion.

As suggested, we expanded our interpretation of the transcriptome findings in the Discussion:

“Associating macroscopic changes in manifold eccentricity with post-mortem microarray data provided by the Allen Human Brain Atlas (Arnatkeviciute et al., 2019; Fornito et al., 2019; Gorgolewski et al., 2014; Hawrylycz et al., 2015; Thompson et al., 2013), we identified gene sets expressed in cortical regions and subcortical structures of the thalamus and striatum during late childhood, adolescence, and young adulthood. Despite these findings being associative and based on separate datasets, they overall support our results that brain network maturation from late childhood until early adulthood implicates micro- and macroscale factors in both subcortical and cortical networks. Coupled network and molecular changes may ultimately change subcortical and cortical circuit properties, including the balance of excitation and inhibition (E/I). Human brain development involves spatio-temporal waves of gene expression changes across different brain regions and developmental time windows (Ip et al., 2010; Kang et al., 2011; Shin et al., 2018). In the study of adolescent development, prior studies have suggested shifts in E/I balance, evolving from a dominant inhibitory bias in early developmental stages towards stronger excitatory drive in later stages, and suggested that these may underlie the maturation of cognitive functions such as working memory and executive control (Dorrn et al., 2010; Lander et al., 2017; Liu et al., 2007). In common neurodevelopmental disorders, including autism, schizophrenia, and attention deficit hyperactivity disorder, imbalances in cortical E/I and cortico-subcortical network function have been demonstrated (Cellot and Cherubini, 2014; Gandal et al., 2018; Lee et al., 2017; Lewis et al., 2005; Nelson and Valakh, 2015; Park et al., 2020a; Sohal and Rubenstein, 2019; Trakoshis et al., 2020), potentially downstream to perturbations of different neurotransmitter systems, such as interneuron-mediated GABA transmission (Bonaventura et al., 2017; Kilb, 2012; Liu et al., 2007; Park et al., 2020a; Silveri et al., 2013; Trakoshis et al., 2020; Tziortzi et al., 2014).”

– For constructing the structural connectome, the Schaefer 7-network atlas was utilized. Can the authors comment on why a functional atlas (rather than a structural atlas) was used here?

We opted for the widely used Schaefer parcellation (Schaefer et al., 2018) as it (i) allows contextualization of our findings within the context of intrinsic functional communities (Yeo et al., 2011), (ii) incorporates the option to assess results across different parcel numbers, and (iii) aligns the current study with prior work from our group (Benkarim et al., 2020; Paquola et al., 2020; Park et al., 2021, 2020a; Rodríguez-Cruces et al., 2020) and others (Baum et al., 2020; Betzel et al., 2019; Osmanlıoğlu et al., 2019).

This justification is now included in the Materials and methods:

“We opted for this atlas as it (i) allows contextualization of our findings within macroscale intrinsic functional communities (Yeo et al., 2011), (ii) incorporates the option to assess results across different granularities, and (iii) aligns the current study with previous work from our group (Benkarim et al., 2020; Paquola et al., 2020; Park et al., 2021, 2020a; Rodríguez-Cruces et al., 2020) and others (Baum et al., 2020; Betzel et al., 2019; Osmanlıoğlu et al., 2019).”

We agree with the reviewer that diffusion MRI studies have often been conducted based on structural atlases. We thus repeated our analyses using a structural parcellation scheme with 200 nodes that preserves the macroscopic boundaries of the Desikan-Killiany atlas (Desikan et al., 2006; Vos de Wael et al., 2020). We found largely consistent age-effects on manifold eccentricity and subcortical-weighted manifold in heteromodal association areas, as well as in caudate, hippocampus, and thalamus (Figure 1—figure supplement 11).

We updated Results:

“i) Connectome manifolds based on structural parcellation. We repeated our analyses with a structural parcellation, defined using a sub-parcellation of folding based on the Desikan-Killiany atlas (Desikan et al., 2006; Vos de Wael et al., 2020) (Figure 1—figure supplement 11). Despite slight differences in the topography of manifold eccentricity in lateral prefrontal, temporal, and occipital cortices, we could replicate strong age-related effects in heteromodal association areas, together with effects in caudate and hippocampus (FDR < 0.05), and marginally in thalamus (FDR < 0.1).”

as well as Materials and methods:

“i) Connectome manifolds based on structural parcellation. To confirm whether functional and structural parcellation schemes yield consistent results, we repeated our main analyses using 200 cortical nodes structural parcellation scheme, which preserves the macroscopic boundaries of the Desikan-Killiany atlas (Desikan et al., 2006; Vos de Wael et al., 2020) (Figure 1—figure supplement 11).”

– While this work touches on an important topic, ties nicely with the increasing body of papers on global brain gradients, and its overall conclusions are warranted, I am not (yet) convinced that it offers fundamentally new insights that could not have been gleaned from previous work (after all, manifold learning simply displays a shadow of the underlying patterns; if the patterns change, so does their shadow). I am also not convinced by the rationale for employing diffusion embedding: the authors state that the ensuing gradients are heritable, conserved across species, capture functional activation patterns during task states, and provide a coordinate system to interrogate brain structure and function, but that would be true for any method that adequately captures biologically meaningful variance in the structural connectivity patterns.

We thank the reviewer for the positive evaluation of our work and happy to provide further rationale for the use of gradients. Being derived from connectomes, we agree that gradient mapping techniques do, in part, recapitulate findings that have been previously shown with more established techniques, for example, graph-theoretical analyses that consider the brain as a graph with nodes and edges and estimate connectome patterns by calculating graph parameters, such as centrality and modular measures (Bullmore and Sporns, 2009; Rubinov and Sporns, 2010). However, we believe that gradient mapping approaches offer a promising combination of methodological and conceptual properties: (i) Gradient mapping based on manifold learning approaches reduce high dimensional (i.e., n x n) connectomes into a series of k spatial maps (i.e., gradients, each of the form n x 1) in a data-driven way and thus provide a simplified view on large-scale connectome organization. Of note, the dimensionality reduction procedure we opted for (i.e., diffusion map embedding) is computationally efficient and robust to noise when computing a globally optimal solution compared to linear approaches (Tenenbaum et al., 2000). (ii) As in the current work, the solution of these techniques is not necessarily a single gradient, but it can be multiple, potentially overlapping gradients. Representing multiple gradients has been suggested to offer the ability to characterize both subregional heterogeneity and functional multiplicity of brain regions (Haak and Beckmann, 2020). In prior work, we showed that the use of multiple diffusion MRI gradients can help to better understand dynamic functional communication patterns in the adult human brain connectome (Park et al., 2021). (iii) As the gradients derived from these techniques can serve as continuous axes of cortical organization, the use of several gradients jointly allows to generate intrinsic coordinate systems in connectivity space, in line with prior conceptual accounts (Bijsterbosch et al., 2020; Haak et al., 2018; Huntenburg et al., 2018; Margulies et al., 2016; Mars et al., 2018). (iv) In addition to the aforementioned methodological advances, prior work has shown that principal gradients, e.g., those estimated from task-free functional MRI (Margulies et al., 2016), follow established models of cortical hierarchy and laminar differentiation (Mesulam, 1998). These observations have allowed gradient mapping approaches to make conceptual contact with seminal work on cortical organization, evolution, and development (Buckner and Krienen, 2013; Goulas et al., 2018; Huntenburg et al., 2018; Sanides, 1969, 1962). What’s more, an emerging literature has shown utility of the gradient framework to study primate evolution and cross-species alignment (Blazquez Freches et al., 2020; Sofie L. Valk et al., 2020; T. Xu et al., 2020), neurodevelopment (Hong et al., 2019; Paquola et al., 2019a), as well as plasticity and structure-function coupling (Park et al., 2021; Vázquez-Rodríguez et al., 2019).

We updated Introduction:

“One emerging approach to address connectome organization and development comes from the application of manifold learning techniques to connectivity datasets. […] An emerging literature has indeed shown utility of the gradient framework to study primate evolution and cross-species alignment (Blazquez Freches et al., 2020; Sofie L. Valk et al., 2020; T. Xu et al., 2020), neurodevelopment (Hong et al., 2019; Paquola et al., 2019a), as well as plasticity and structure-function coupling (Park et al., 2021; Sofie L Valk et al., 2020; Vázquez-Rodríguez et al., 2019).”

To build cortex-wide structural connectome manifolds, our work opted for diffusion map embedding, a non-linear dimensionality reduction technique (Coifman and Lafon, 2006), following a seminal connectivity gradient study (Margulies et al., 2016). Follow-up work from our group and others adopted the same approach (Hong et al., 2019, 2020; Huntenburg et al., 2017; Larivière et al., 2020; Margulies et al., 2016; Müller et al., 2020; Paquola et al., 2019a; Park et al., 2021; Sofie L. Valk et al., 2020; Vos de Wael et al., 2020). In addition to ensuring continuity with these prior studies, non-linear techniques do in general explain more information with the same number of components than linear techniques, which may increase performance on group classification, cluster identification, and phenotypic score prediction (Errity and McKenna, 2007; Gallos et al., 2020; Hong et al., 2020). Diffusion map embedding, in particular, is only controlled by few parameters, and thus computationally efficient and rather robust to noise (Errity and McKenna, 2007; Gallos et al., 2020; Hong et al., 2020; Tenenbaum et al., 2000). Of note, the non-linearity is simply introduced by running the singular value decomposition, which is also at the heart of PCA, on the Markov matrix of the normalized connectome, which also lends a direct interpretation to the reduced manifold as being interpretable as a ‘diffusion map’ between the nodes.

We nevertheless agree with the Reviewer’s point that non-linear and linear techniques often converge, and that it is not fully established which approach should be prioritized (Hong et al., 2020; Vos de Wael et al., 2020). To assess consistency, we thus repeated our analysis based on principal component analysis (PCA) and observed consistent manifold dimensions (linear correlation = 0.998 ± 0.001 across E1/E2/E3; Figure 1—figure supplement 7A), as well as similar age-effects as in the main analysis, notably a similar expansion of the manifold space (Figure 1—figure supplement 7B).

We updated Results:

“e) Connectome manifold generation using principal component analysis. In a separate analysis, we generated eigenvectors using principal component analysis (Wold et al., 1987), instead of diffusion map embedding (Coifman and Lafon, 2006), and found consistent spatial maps (linear correlation = 0.998 ± 0.001 across E1/E2/E3; Figure 1—figure supplement 7A) and longitudinal findings (Figure 1—figure supplement 7B).”

as well as Materials and methods:

“An affinity matrix was constructed with a normalized angle kernel, and eigenvectors were estimated via diffusion map embedding (Figure 1A), a non-linear dimensionality reduction technique (Coifman and Lafon, 2006) that projects connectome features into low dimensional manifolds (Margulies et al., 2016). This technique is only controlled by few parameters, computationally efficient, and relatively robust to noise compared to other non-linear techniques (Errity and McKenna, 2007; Gallos et al., 2020; Hong et al., 2020; Tenenbaum et al., 2000), and has been extensively used in the previous gradient mapping literature (Hong et al., 2019, 2020; Huntenburg et al., 2017; Larivière et al., 2020a; Margulies et al., 2016; Müller et al., 2020; Paquola et al., 2019a; Park et al., 2021; Sofie L. Valk et al., 2020; Vos de Wael et al., 2020).”

“e) Connectome manifold generation using principal component analysis. To explore consistency of our results when using different dimensionality reduction techniques, we generated connectome manifolds using principal component analysis (Wold et al., 1987), instead of relying on diffusion map embedding (Coifman and Lafon, 2006), and performed longitudinal modeling (Figure 1—figure supplement 7). We compared the eigenvectors estimated from diffusion map embedding and principal component analysis using linear correlations.”

– The authors show that the maturational change of the manifold features predict intelligence at follow-up, but did not show that intelligence itself exhibited changes that exceeded the error bounds of the regression line. Why not predict IQ change?

As suggested, we additionally predicted the change of IQ between the baseline and follow-up, instead of IQ at follow-up, using the imaging features. However, we could not find significant results for predicting ∆IQ using cortical features at baseline (mean ± SD r = -0.10 ± 0.04, MAE = 6.62 ± 0.06, p = 0.27) nor at both baseline and maturational change (r = -0.08 ± 0.04, MAE = 6.60 ± 0.06, p = 0.34). Adding subcortical regions did not improve the performance (baseline: r = -0.01 ± 0.03, MAE = 7.02 ± 0.12, p = 0.73; baseline and maturational change: r = 0.03 ± 0.03, MAE = 6.53 ± 0.06, p = 0.67). The low performance might be due to small variations in ∆IQ, where many participants (42%) showed IQ changes less than 5 and 76% less than 10.

We updated Results:

“We also predicted the change of IQ between the baseline and follow-up, instead of IQ at follow-up, using the imaging features. However, we could not find significant results.”

as well as Materials and methods:

“In addition to predicting future IQ, we performed the same prediction analysis to predict the change of IQ between the baseline and follow-up.”

– The slight improvements in prediction accuracy observed after adding maturational change and subcortical features to the features at baseline will necessarily happen by adding more regression parameters and may not be meaningful.

To improve robustness of prediction analysis, as well as avoid overfitting, the revised work employed a nested ten-fold cross-validation framework (Cawley and Talbot, 2010; Parvandeh et al., 2020; Tenenbaum et al., 2000; Varma and Simon, 2006) with elastic net regularization (Zou and Hastie, 2005). Specifically, we split the dataset into training (9/10) and test (1/10) partitions, and each training partition was further split into inner training and testing folds using another ten-fold cross-validation. Within the inner fold, elastic net regularization finds a set of non-redundant features to explain the dependent variable. Using a linear regression, we predicted the IQ scores of inner fold test data using the features of the selected brain regions by controlling for age, sex, site, and head motion. The model with minimum mean absolute error (MAE) across the inner folds was applied to the test partition of the outer fold, and the IQ scores of outer fold test data were predicted. The prediction procedure was repeated 100 times with different training and test sets to reduce subject selection bias. Across cross-validation and iterations, 6.24 ± 5.74 (mean ± SD) features were selected to predict IQ using manifold eccentricity of cortical regions at baseline, 6.20 ± 5.14 cortical features at baseline and maturational change, 5.45 ± 5.99 cortical and subcortical features at baseline, and 5.16 ± 5.43 at baseline and maturational change. In this scenario, adding more independent variables may not per se lead to improvement in prediction accuracy.

We updated Results:

“We used elastic net regularization with nested ten-fold cross-validation (Cawley and Talbot, 2010; Parvandeh et al., 2020; Tenenbaum et al., 2000; Varma and Simon, 2006; Zou and Hastie, 2005) (see Methods), and repeated the prediction 100 times with different training and test dataset compositions to mitigate subject selection bias. Across cross-validation and iterations, 6.24 ± 5.74 (mean ± SD) features were selected to predict IQ using manifold eccentricity of cortical regions at baseline, 6.20 ± 5.14 cortical features at baseline and maturational change, 5.45 ± 5.99 cortical and subcortical features at baseline, and 5.16 ± 5.43 at baseline and maturational change, suggesting that adding more independent variables may not per se lead to improvement in prediction accuracy.”

as well as Materials and methods:

“For each evaluation, a subset of features that could predict future IQ was identified using elastic net regularization (𝜌 = 0.5) with optimized regularization parameters (L1 and L2 penalty terms) via nested ten-fold cross-validation (Cawley and Talbot, 2010; Parvandeh et al., 2020; Tenenbaum et al., 2000; Varma and Simon, 2006; Zou and Hastie, 2005). We split the dataset into training (9/10) and test (1/10) partitions, and each training partition was further split into inner training and testing folds using another ten-fold cross-validation. Within the inner fold, elastic net regularization finds a set of non-redundant features to explain the dependent variable. Using a linear regression, we predicted the IQ scores of inner fold test data using the features of the selected brain regions by controlling for age, sex, site, and head motion. The model with minimum mean absolute error (MAE) across the inner folds was applied to the test partition of the outer fold, and the IQ scores of outer fold test data were predicted. The prediction procedure was repeated 100 times with different training and test sets to reduce subject selection bias.”

– Interpretability and Lack of Comparison

The authors claim repeatedly that they are "capitalizing on advanced manifold learning techniques".

One could imagine an infinite number of papers that take a dataset, use a technique to extract a metric, X (e.g., eccentricity), and then write about the changes in X with some property of interest, Y (e.g., age). Given this set of papers (and the non-independence between the set of possible Xs), the reader ought to be most interested in those Xs that provide the best performance and simplest interpretation, with other papers being redundant.

Thus, a nuanced approach to presenting a paper like this is to demonstrate that the metric used represents an advance over alternative, simpler-to-compute, or clearer-to-interpret metrics that already exist.

In this paper, however, the authors do not demonstrate the benefits of their particular choice of applying a specific nonlinear dimensionality reduction method using 3 dimensions alignment to a template manifold and then computing an eccentricity metric.

We thank the reviewers for this suggestion, please see the points below. In brief, both linear and non-linear dimensionality reduction techniques compress high dimensional connectome data (n x n) into a series of lower-dimensional eigenvectors (i.e., gradients, each of the form n x 1), offering a synoptic view on connectome development. In our case, we used these manifold learning techniques to derive a three-dimensional connectivity space to which we were able to sensitively track developmental changes in adolescence. Beyond these methodological considerations, prior work has shown that principal gradients estimated from task-free functional (Margulies et al., 2016), microstructural (Paquola et al., 2019b), and diffusion MRI (Park et al., 2020a) broadly converge along an established model of sensory-fugal neural hierarchy and laminar differentiation (Mesulam, 1998), allowing the gradient mapping approach to align with theories of cortical organization and evolution (Buckner and Krienen, 2013; Goulas et al., 2018; Huntenburg et al., 2018; Sanides, 1969, 1962). We believe that the relative analytical simplicity and ability to contextualize work in classic theory represent an attractive justification for gradient mapping techniques. Below, we furthermore addressed the specific comments:

For example:

i. Is the nonlinearity required (e.g., does it outperform PCA or MDS)?

Non-linear dimensionality techniques are theoretically appealing, as they can explain more information of the original data with the same number of components than linear techniques, which may boost performance on group classification, cluster identification, and phenotypic score prediction (Errity and McKenna, 2007; Gallos et al., 2020; Hong et al., 2020). Of note, diffusion map embedding is analytically straightforward, as it simply performs a singular value decomposition on the Markov matrix, which can be readily interpreted as a diffusion process throughout the connectome. Researchers have previously adopted both linear (Hong et al., 2020; Murphy et al., 2019; Shine et al., 2019; Tian et al., 2020) and non-linear techniques (Guell et al., 2018; Haak and Beckmann, 2020; Hong et al., 2019; Huntenburg et al., 2017; Larivière et al., 2020; Margulies et al., 2016; Müller et al., 2020; Paquola et al., 2019b, 2019a; Sofie L. Valk et al., 2020), and it thus remains an open question as to which technique should be preferred (Hong et al., 2020; Vos de Wael et al., 2020). We opted to stick with the diffusion maps in our main analysis to stay consistent with prior work from our group and others (Hong et al., 2019, 2020; Huntenburg et al., 2017; Larivière et al., 2020; Margulies et al., 2016; Müller et al., 2020; Paquola et al., 2019a; Park et al., 2021; Sofie L. Valk et al., 2020; Vos de Wael et al., 2020). In the revised sensitivity analyses, we also derived gradients via PCA, and showed consistent manifold dimensions using this technique (linear correlation = 0.998 ± 0.001 across E1/E2/E3; Figure 1—figure supplement 7A), as well as consistent age-effects on manifold eccentricity (Figure 1—figure supplement 7B).

We updated Results:

“e) Connectome manifold generation using principal component analysis. In a separate analysis, we generated eigenvectors using principal component analysis (Wold et al., 1987), instead of diffusion map embedding (Coifman and Lafon, 2006), and found consistent spatial maps (linear correlation = 0.998 ± 0.001 across E1/E2/E3; Figure 1—figure supplement 7A) and longitudinal findings (Figure 1—figure supplement 7B).”

Materials and methods:

“An affinity matrix was constructed with a normalized angle kernel, and eigenvectors were estimated via diffusion map embedding (Figure 1A), a non-linear dimensionality reduction technique (Coifman and Lafon, 2006) that projects connectome features into low dimensional manifolds (Margulies et al., 2016). This technique is only controlled by few parameters, computationally efficient, and relatively robust to noise compared to other non-linear techniques (Errity and McKenna, 2007; Gallos et al., 2020; Hong et al., 2020; Tenenbaum et al., 2000), and has been extensively used in the previous gradient mapping literature (Hong et al., 2019, 2020; Huntenburg et al., 2017; Larivière et al., 2020a; Margulies et al., 2016; Müller et al., 2020; Paquola et al., 2019a; Park et al., 2021; Sofie L. Valk et al., 2020; Vos de Wael et al., 2020).”

“e) Connectome manifold generation using principal component analysis. To explore consistency of our results when using different dimensionality reduction techniques, we generated connectome manifolds using principal component analysis (Wold et al., 1987), instead of relying on diffusion map embedding (Coifman and Lafon, 2006), and performed longitudinal modeling (Figure 1—figure supplement 7). We compared the eigenvectors estimated from diffusion map embedding and principal component analysis using linear correlations.”

as well as Discussion:

“In our longitudinal study, we could identify marked connectome expansion during adolescence, mainly encompassing transmodal and heteromodal association cortex in prefrontal, temporal, and posterior regions, the territories known to mature later in development (Gogtay et al., 2004; Shaw et al., 2006). Findings remained consistent when we considered a linear dimensionality reduction technique, suggesting robustness to methodological details of this analysis.”

ii. Is there something special about picking 3 dimensions to do the eccentricity calculation? Is dimensionality reduction required at all (e.g., would you get similar results by computing eccentricity in the full-dimensional space?)

Based on the screen plot in Figure 1A, we selected the first three eigenvectors as (i) they explained approximately 50% of variance in the affinity matrix, (ii) each explained >10% (i.e., 20.7/15.8/13.5% for E1/E2/E3, respectively), and (iii) the choice of three eigenvectors captured a clearly visible eigengap. Prior work has chosen a similar approach, as higher order components may sometimes not show clear spatial patterns and/or higher noise contributions (Haak et al., 2018; Hong et al., 2019, 2020; Margulies et al., 2016; Paquola et al., 2019b, 2019a; Park et al., 2021, 2020a; Vos de Wael et al., 2018). In response to the reviewers’ suggestion, we also assessed manifold eccentricity from the full dimensional space (Figure 1—figure supplement 9A-B). The spatial pattern of manifold eccentricity correlated with the original solution based on three eigenvectors (r = 0.54, p < 0.001). Moreover, longitudinal changes in manifold eccentricity were similar to the original findings (linear correlation of t-statistic map = 0.68, p < 0.001; Figure 1—figure supplement 9C), confirming manifold expansion in transmodal cortices.

We updated Results:

“g) Manifold eccentricity based on all eigenvectors. Repeating manifold eccentricity calculation and age modeling using all eigenvectors, instead of using only the first three, we observed relatively consistent results with our original findings (linear correlation of manifold eccentricity r = 0.54, p < 0.001; t-statistic map r = 0.68, p < 0.001), also pointing to manifold expansion in transmodal cortices (Figure 1—figure supplement 9).”

as well as Materials and methods:

“g) Manifold eccentricity analysis based on all eigenvectors. We repeated our analysis by calculating manifold eccentricity from all eigenvectors to assess consistency of the findings (Figure 1—figure supplement 9).”

iii. Does it outperform basic connectome measures (e.g., the simple ones the authors compute)?

To compare the age-effects on manifold eccentricity with graph-theoretical effects, we calculated betweenness, degree, and eigenvector centrality measures (Rubinov and Sporns, 2010). While betweenness centrality did not reveal significant effects, degree and eigenvector centrality showed effects in similar regions as manifold eccentricity (Figure 1—figure supplement 8). Calculating linear correlation between the effect size maps, we found a graded pattern of similarity, which was low but significant for betweenness centrality (r = 0.18, p = 0.02; significance determined by 1,000 spin tests that account for spatial autocorrelation (Alexander-Bloch et al., 2018)) and moderate for degree and eigenvector centrality (r = 0.57, p < 0.001; r = 0.47, p < 0.001).

We updated Results:

“f) Longitudinal changes in graph-theoretical measures. Repeating the longitudinal modeling using graph-theoretical centrality measures, we found significant age-related longitudinal changes in degree and eigenvector centrality, while betweenness centrality did not reveal significant effects, in similar regions to those that had significant age-related changes in manifold eccentricity (Figure 1—figure supplement 8). Correlating the effect size maps for manifold eccentricity and each graph measure, we found a significant yet variable spatial similarity of the effect maps (betweenness centrality: r = 0.18, spin-test p = 0.02; degree centrality: r = 0.57, p < 0.001; eigenvector centrality: r = 0.47, p < 0.001).”

as well as Materials and methods:

“f) Longitudinal changes in graph-theoretical measures. To compare longitudinal changes in manifold eccentricity with those in graph-theoretical centrality measures, we calculated betweenness, degree, and eigenvector centrality of the structural connectomes and build similar linear mixed effects models to assess longitudinal change (Figure 1—figure supplement 8). Betweenness centrality is the number of weighted shortest paths between any combinations of nodes that run through that node, degree centrality is the sum of edge weights connected to a given node, and eigenvector centrality measures the influence of a node in the whole network (Lohmann et al., 2010; Rubinov and Sporns, 2010; Zuo et al., 2012). Spatial similarity between t-statistics of centrality and manifold measures were assessed with 1,000 spin tests that account for spatial autocorrelation (Alexander-Bloch et al., 2018).”

There is a clear down-side of how opaque the approach is (and thus difficult to interpret relative to, say, connectivity degree), so one would hope for a correspondingly strong boost in performance. The authors could also do more to develop some intuition for the idea of a low-dimensional connection-pattern-similarity-space, and how to interpret taking Euclidean distances within such a space.

CT=1N[i=1NT(E1)i,i=1NT(E2)i,i=1NT(E3)i]
ME=e=13{I(Ee)CT(e)}2

Here, we used eigenvectors estimated from the decomposition technique to generate a new coordinate system in connectivity space (Bijsterbosch et al., 2020; Haak et al., 2018; Huntenburg et al., 2018; Margulies et al., 2016; Mars et al., 2018). Notably, in this manifold space, interconnected brain regions with similar connectivity patterns are closer together, while regions that do not have significant connectivity nor similarity in connectivity patterns are located farther apart. In our work, the specific space is a diffusion map, where Euclidean distances in the manifold correspond to diffusion times between the nodes of the network. To analyze how the multi-dimensional manifold structures change in the low dimensional manifold space, we simplified the multivariate eigenvectors into a single scalar value i.e., manifold eccentricity. Manifold eccentricity was calculated as the Euclidean distance between the manifold origin and all data points (i.e., brain regions) in manifold space. For example, using three eigenvectors, manifold eccentricity was defined as follows:CT is the template manifold origin; N number of brain regions; T() template manifold; ME manifold eccentricity; I() individual manifold; CT(e) origin of eth template manifold. This concept is visualized in Figure 1—figure supplement 15. Each brain region (i.e., each dot in the scatter plot) is described as a vector from the manifold origin (i.e., triangular mark in the scatter plot), and manifold eccentricity is simply a length (i.e., Euclidean distance) of that vector. Shifts in connectivity patterns of a given region thus will lead to shifts in the vectors, which in turn changes the manifold eccentricity. Thus, manifold eccentricity quantifies global brain organization in connectivity space.

We updated Materials and methods:

“Briefly, the eigenvectors estimated from the decomposition technique generate a connectivity coordinate system (Bijsterbosch et al., 2020; Haak et al., 2018; Huntenburg et al., 2018; Margulies et al., 2016; Mars et al., 2018) – the diffusion map, where […] Shifts in connectivity patterns of a given region thus will lead to shifts in the vectors, which in turn changes the manifold eccentricity. Thus, manifold eccentricity quantifies global brain organization based in the connectivity space.”

– Developmental Enrichment Analysis

Both in the main text and in the Materials and methods, this is described as "genes were fed into a developmental enrichment analysis". Can some explanation be provided as to what happens between the "feeding in" and what comes out? Without clearly described methods, it is impossible to interpret or critique this component of the paper. If the methodological details are opaque, then the significance of the results could be tested numerically relative to some randomized null inputs being 'fed in' to demonstrate specificity of the tested phenotype.

We implemented developmental enrichment analysis using the cell-type specific expression analysis (CSEA) developmental expression tool (http://genetics.wustl.edu/jdlab/csea-tool-2) (Dougherty et al., 2010; Xu et al., 2014). This technique evaluates the significance of the overlap between the identified gene list and RNAseq data obtained from BrainSpan dataset (http://www.brainspan.org) across six brain regions (i.e., cortex, thalamus, striatum, cerebellum, hippocampus, amygdala) and ten developmental periods (from early fetal to young adulthood) approximated from mouse data, yielding a total of 60 combinations of developmental enrichment profiles (Xu et al., 2014). The significance is calculated based on Fisher’s exact tests (Fisher, 1922) with FDR correction (Benjamini and Hochberg, 1995).

We updated Results:

“We performed developmental gene set enrichment analysis using the cell-type specific expression analysis (CSEA) tool, which compares the selected gene list with developmental enrichment profiles (see Methods) (Dougherty et al., 2010; Xu et al., 2014). This analysis highlights developmental time windows across macroscopic brain regions in which genes are strongly expressed. We found marked expression of the genes enriched from childhood onwards in the cortex, thalamus, and cerebellum (FDR < 0.001; Figure 4B).”

as well as Materials and methods:

“Leveraging cell-type specific expression analysis (CSEA) developmental expression tool (http://genetics.wustl.edu/jdlab/csea-tool-2) (Dougherty et al., 2010; Xu et al., 2014), we evaluated the significance of overlap between the genes showing consistent whole-brain expression pattern across donors (FDR < 0.05) with RNAseq data obtained from BrainSpan dataset (http://www.brainspan.org). The significance was calculated based on Fisher’s exact test (Fisher, 1922) with FDR correction (Benjamini and Hochberg, 1995). The CSEA tool provides simplified results of gene enrichment profiles along six major brain regions (i.e., cortex, thalamus, striatum, cerebellum, hippocampus, amygdala) across ten developmental periods (from early fetal to young adulthood) approximated from mouse data, yielding a total of 60 combinations of developmental enrichment profiles (Xu et al., 2014).”

As suggested, we additionally assessed spatial correlations between the post-mortem gene expression maps and rotated maps of the age-related changes in manifold eccentricity (100 spherical, random rotations). For each iteration, the selected genes were used for developmental enrichment analysis using CSEA tool (Xu et al., 2014), and we obtained FDR-corrected p-values for each brain division and developmental period. We, thus, built a null distribution using the rotated p-values into which the actual p-value was placed. If the actual p-value did not belong to the 95% of the null distribution, it was deemed significant. The results remained consistent.

We updated Materials and methods:

“We repeated developmental enrichment analysis using the genes identified from the rotated maps of the age-related changes in manifold eccentricity (100 spherical rotations). For each iteration, we obtained developmental expression profiles using the identified genes, where the FDR-corrected p-values built a null distribution. For each brain division and developmental period, if the actual p-value placed outside 95% of the null distribution, it was deemed significant.”

– IQ prediction

The predictions seem to be very poor (equality lines, y = x, should be drawn in Figure 5, to show what perfect predictions would look like; linear regressions are not helpful for a prediction task, and are deceptive of the appropriate MAE computation). The authors do not perform any comparisons in this section (even to a real baseline model like `predicted_IQ = mean(training_set_IQ)`). They also do not perform statistical tests (or quote p-values), but nevertheless make a range of claims, including of "significant prediction" or "prediction accuracy was improved", "reemphasize the benefits of incorporating subcortical nodes", etc. All of these claims should be tested relative to rigorous statistics, and comparisons to appropriate baseline/benchmark approaches.

We thank the reviewer for the comment. To assess the significance of our prediction model, we compared the prediction performance to the suggested baseline model (i.e., predicted IQ = mean(training set IQ)). The prediction performance for the baseline model was not good showing a low negative correlation between actual and predicted IQ scores (r = -0.15 ± 0.06, MAE = 8.98 ± 0.04, p = 0.12). Our model based on cortical features (baseline: r = 0.14 ± 0.04, MAE = 8.93 ± 0.16, p = 0.09; baseline and maturational change: r = 0.18 ± 0.04, MAE = 9.10 ± 0.19, p = 0.04), and based on both cortical and subcortical features (baseline: r = 0.17 ± 0.03, MAE = 8.74 ± 0.11, p = 0.04; baseline and maturational change: r = 0.21 ± 0.02, MAE = 8.86 ± 0.14, p = 0.01) outperformed the baseline model (Meng’s z-test p < 0.001 for all cases) (Meng et al., 1992). The results suggest that while manifold features are only weakly related to IQ, they still provide information on future IQ above and beyond the baseline model.

We updated the Results:

“We compared the prediction performance of our model with a baseline model, where IQ of the test set was simple average of training set (r = -0.15 ± 0.06, MAE = 8.98 ± 0.04, p = 0.12; see Methods). We found that our model outperformed this baseline model (Meng’s z-test p < 0.001) (Meng et al., 1992).”

Materials and methods:

“To assess whether our model outperforms baseline model, we predicted IQ of test data using average of IQ of training data (i.e., predicted IQ = mean(training set IQ)). The improvement of prediction performance was assessed using Meng’s z-test (Meng et al., 1992).”

as well as Discussion:

“Of note, although our model significantly outperformed a baseline model, the relationship between the actual and predicted IQ scores did not locate on the equality line and the strength of the association was rather weak. Further improvements in brain-based IQ prediction in adolescence, for example through combinations of structural and functional imaging features, will be a focus of future work.”

We additionally performed the prediction analysis using a regression tree method (i.e., decision tree learning), a non-linear approach based on tree model constituting root node and split leaf nodes, where the leaf nodes contain the response variables (Breiman et al., 1984). However, the prediction results were not improved compared to linear model (Figure 5—figure supplement 2).

We updated Results:

“l) IQ prediction using non-linear model. We predicted IQ at follow-up using a regression tree method (Breiman et al., 1984), instead of linear regression model, but we could not find improved prediction performance (Figure 5—figure supplement 2).”

as well as Materials and methods:

“l) IQ prediction using non-linear model. We additionally performed prediction analysis for predicting future IQ score using regression tree method (i.e., decision tree learning), a non-linear approach based on tree model constituting root node and split leaf nodes, where the leaf nodes contain the response variables (Breiman et al., 1984) (Figure 5—figure supplement 2).”

– Group Connectome

Given how much the paper relies on estimating a group structural connectome, it should be visualized and characterized. For example, a basic analysis of the distribution of edge weights and degree, especially as edge weights can vary over orders of magnitude and high weights (more likely to be short distances) may therefore unduly dominate some of the low-dimensional components). The authors may also consider testing robustness performed to alternative ways of estimating the connectome [e.g., Oldham et al. NeuroImage 222, 117252 (2020)] and its group-level summary [e.g., Roberts et al. NeuroImage 145, 1-42 (2016)].

We thank the Reviewers for these suggestions. We calculated structural connectome in individual level, and estimated group representative connectome using a distance-dependent thresholding (Betzel et al., 2019), which is visualized in Figure 1A. As suggested, we performed the longitudinal modeling using the connectome edge weights and found significant increases in weights within frontoparietal and default mode networks, as well as attention and sensory networks (FDR < 0.05; Figure 1—figure supplement 12). These results are consistent with our main findings based on manifold eccentricity that pointed to mainly transmodal effects.

We updated Results:

“j) Longitudinal modeling using edge weights. Repeating the longitudinal modeling across age using connectome edge weights, we found significant increases in edge weights within frontoparietal and default mode networks, as well as in attention and sensory networks (FDR < 0.05; Figure 1—figure supplement 12), consistent with findings based on manifold eccentricity.”

as well as Materials and methods:

“j) Longitudinal modeling using edge weights. In addition to the analyses based on manifold eccentricity, linear mixed effect modeling using connectome edge weights assessed age-related longitudinal changes in streamline strength (Figure 1—figure supplement 12).”

To assess robustness of group representative structural connectomes, we built the connectomes with two different approaches: (i) Distance-dependent thresholding (Betzel et al., 2019) was adopted for the main analysis, and (ii) consistency thresholding was assessed in a supplementary analysis. Consistency thresholding approach averages subject specific matrices, in addition to performing a 50, 40, 30, 20, 10% thresholding, as well as simple averaging (i.e., 0% thresholding) (Wang et al., 2019) (Figure 1—figure supplement 10). We calculated linear correlations between spatial maps of eigenvectors derived from distance-dependent thresholding and the group consistency method and found high similarity (r = 0.89 ± 0.01 for E1; 0.93 ± 0.004 for E2; 0.85 ± 0.01 for E3 across consistency thresholds).

We updated Results:

“h) Robustness of group representative structural connectome. We compared gradients derived from the group representative structural connectome, based on (i) distance-dependent thresholding (Betzel et al., 2019) and (ii) consistency thresholding (Wang et al., 2019) (Figure 1—figure supplement 10). We found high similarity in spatial maps of the estimated manifolds (r = 0.89 ± 0.01 for E1; 0.93 ± 0.004 for E2; 0.85 ± 0.01 for E3 across six different thresholds), indicating robustness.”

as well as Materials and methods:

“h) Robustness of group representative structural connectome. We compared the distance-dependent thresholding (Betzel et al., 2019) that was adopted for the main analysis with a consistency thresholding approach (Wang et al., 2019). The latter averages subject specific matrices, in addition to performing a 50, 40, 30, 20, 10% thresholding, as well as simple averaging (i.e., 0% thresholding) (Figure 1—figure supplement 10).”

Individual Alignment

The paper relies on individuals being successfully aligned to the template manifold. Accordingly, some analysis should be performed quantifying how well individuals could be mapped. Presumably some subjects fit very well onto the template, whereas others do not. Is there something interesting about the poorly aligned subjects? Do your results improve when excluding them?

Individual manifolds were aligned to template manifold with Procrustes alignment, which makes eigenvectors from different individuals more comparable e.g., by flipping eigenvector signs (Langs et al., 2015; Vos de Wael et al., 2020). To evaluate the procedure, we calculated linear correlations between template and individual manifolds before/after alignment. We found significant improvement in correlations after aligning individual manifolds to the template (0.92±0.03/0.93±0.03/0.94±0.03 after alignment; -0.02±0.03/-0.001±0.37/0.003±0.12 before alignment for E1/E2/E3).

As suggested, we also performed the longitudinal modeling after excluding 10% of subjects with poor alignment (cutoff r = 0.83; the new set had a linear correlation with template manifold of mean ± SD r = 0.94 ± 0.01). We still could find consistent results (Figure 1—figure supplement 6), where the t-statistics showed strong correlation with those derived using the whole subjects (r = 0.97, p < 0.001), suggesting robustness.

We updated Results:

“d) Gradient alignment fidelity. When calculating linear correlations between template and individual manifolds before and after alignment, we found significant increases after alignment (r = 0.92±0.03/0.93±0.03/0.94±0.03) compared to before alignment (-0.02±0.03/-0.001±0.37/0.003±0.12), for E1/E2/E3, respectively, supporting effectiveness of alignment. After excluding 10% of subjects with poor alignment (cutoff r = 0.83; the new set was correlated with the template manifold, r = 0.94 ± 0.01), we found consistent age-related changes in manifold eccentricity (Figure 1—figure supplement 6), with the t-statistic map showing strong correlation to the map derived in the whole sample (r = 0.97, p < 0.001).”

as well as Materials and methods:

“d) Gradient alignment fidelity. To assess robustness of individual alignment, we computed linear correlations between the template and individual manifolds before and after alignment. We also repeated the linear mixed effect modeling after excluding 10% of subjects with the lowest alignment to the template manifold (Figure 1—figure supplement 6).”

https://doi.org/10.7554/eLife.64694.sa2

Article and author information

Author details

  1. Bo-yong Park

    1. McConnell Brain Imaging Centre, Montreal Neurological Institute and Hospital, McGill University, Montreal, Canada
    2. Department of Data Science, Inha University, Incheon, Republic of Korea
    Contribution
    Conceptualization, Software, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    by9433@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7096-337X
  2. Richard AI Bethlehem

    1. Autism Research Centre, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    2. Brain Mapping Unit, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    Contribution
    Conceptualization, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0714-0685
  3. Casey Paquola

    1. McConnell Brain Imaging Centre, Montreal Neurological Institute and Hospital, McGill University, Montreal, Canada
    2. Institute of Neuroscience and Medicine (INM-1), Forschungszentrum Jülich, Jülich, Germany
    Contribution
    Conceptualization, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0190-4103
  4. Sara Larivière

    McConnell Brain Imaging Centre, Montreal Neurological Institute and Hospital, McGill University, Montreal, Canada
    Contribution
    Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  5. Raul Rodríguez-Cruces

    McConnell Brain Imaging Centre, Montreal Neurological Institute and Hospital, McGill University, Montreal, Canada
    Contribution
    Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  6. Reinder Vos de Wael

    McConnell Brain Imaging Centre, Montreal Neurological Institute and Hospital, McGill University, Montreal, Canada
    Contribution
    Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
  7. Neuroscience in Psychiatry Network (NSPN) Consortium

    Contribution
    Data curation
    Additional information
    A complete list of investigators from the Neuroscience in Psychiatry Network (NSPN) Consortium can be found in the Supplementary Information.
    1. Edward Bullmore, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    2. Raymond Dolan, Max Planck University College London Centre for Computational Psychiatry and Ageing Research, University College London, London, United Kingdom
    3. Ian Goodyer, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    4. Peter Fonagy, Research Department of Clinical, Educational and Health Psychology, University College London, London, United Kingdom
    5. Peter Jones, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    6. Michael Moutoussis, Max Planck University College London Centre for Computational Psychiatry and Ageing Research, University College London, London, United Kingdom
    7. Tobias Hauser, Max Planck University College London Centre for Computational Psychiatry and Ageing Research, University College London, London, United Kingdom
    8. Sharon Neufeld, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    9. Rafael Romero-Garcia, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    10. Michelle St Clair, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    11. Petra Vértes, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    12. Kirstie Whitaker, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    13. Becky Inkster, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    14. Gita Prabhu, Max Planck University College London Centre for Computational Psychiatry and Ageing Research, University College London, London, United Kingdom
    15. Cinly Ooi, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    16. Umar Toseeb, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    17. Barry Widmer, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    18. Junaid Bhatti, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    19. Laura Villis, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    20. Ayesha Alrumaithi, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    21. Sarah Birt, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    22. Aislinn Bowler, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    23. Kalia Cleridou, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    24. Hina Dadabhoy, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    25. Emma Davies, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    26. Ashlyn Firkins, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    27. Sian Granville, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    28. Elizabeth Harding, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    29. Alexandra Hopkins, Max Planck University College London Centre for Computational Psychiatry and Ageing Research, University College London, London, United Kingdom
    30. Daniel Isaacs, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    31. Janchai King, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    32. Danae Kokorikou, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    33. Christina Maurice, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    34. Cleo McIntosh, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    35. Jessica Memarzia, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    36. Harriet Mills, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    37. Ciara O’Donnell, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    38. Sara Pantaleone, Wellcome Centre for Human Neuroimaging, University College London, London, United Kingdom
    39. Jenny Scott, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    40. Beatrice Kiddle, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    41. Ela Polek, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    42. Pasco Fearon, Research Department of Clinical, Educational and Health Psychology, University College London, London, United Kingdom
    43. John Suckling, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    44. Anne-Laura van Harmelen, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    45. Rogier Kievit, Max Planck University College London Centre for Computational Psychiatry and Ageing Research, University College London, London, United Kingdom
    46. Sam Chamberlain, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
  8. Edward T Bullmore

    Brain Mapping Unit, Department of Psychiatry, University of Cambridge, Cambridge, United Kingdom
    Contribution
    Conceptualization, Resources, Methodology, Writing - review and editing
    Competing interests
    ETB. serves on the scientific advisory board of Sosei Heptares and as a consultant for GlaxoSmithKline.
  9. Boris C Bernhardt

    McConnell Brain Imaging Centre, Montreal Neurological Institute and Hospital, McGill University, Montreal, Canada
    Contribution
    Conceptualization, Supervision, Writing - original draft, Writing - review and editing
    For correspondence
    boris.bernhardt@mcgill.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9256-6041

Funding

Canada Research Chairs

  • Boris C Bernhardt

National Research Foundation of Korea (NRF2020R1A6A3A03037088)

  • Bo-yong Park

Fonds de la Recherche du Quebec – Santé

  • Bo-yong Park
  • Casey Paquola
  • Raul Rodríguez-Cruces
  • Boris C Bernhardt

Montreal Neurological Institute and Hospital (MNI) (Molson Neuro-Engineering fellowship)

  • Bo-yong Park

British Academy (Post-Doctoral Fellowship)

  • Richard AI Bethlehem

Autism Research Trust

  • Richard AI Bethlehem

Canadian Institutes of Health Research

  • Sara Larivière

National Institute for Health Research (Senior Investigator award)

  • Edward T Bullmore

Natural Sciences and Engineering Research Council of Canada (NSERC Discovery-1304413)

  • Boris C Bernhardt

Canadian Institutes of Health Research (FDN-154298)

  • Boris C Bernhardt

SickKids Foundation (NI17-039)

  • Boris C Bernhardt

Azrieli Center for Autism Research (ACAR-TACC)

  • Boris C Bernhardt

BrainCanada

  • Boris C Bernhardt

MNI-Cambridge collaborative award

  • Bo-yong Park
  • Richard AI Bethlehem
  • Casey Paquola
  • Boris C Bernhardt

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

Dr. Bo-yong Park was funded by the National Research Foundation of Korea (NRF-2020R1A6A3A03037088), Molson Neuro-Engineering fellowship by Montreal Neurological Institute and Hospital (MNI), and a postdoctoral fellowship of the Fonds de la Recherche du Quebec – Santé (FRQ-S). Dr. Richard AI Bethlehem was funded by a British Academy Post-Doctoral Fellowship and the Autism Research Trust. Dr. Casey Paquola and Dr Raul R Cruces were funded through postdoctoral fellowships of the FRQ-S. Ms. Sara Larivière acknowledges funding from the Canadian Institutes of Health Research (CIHR). Dr. Edward T Bullmore was supported by a Senior Investigator award from the National Institute of Health Research (NIHR). Dr. Boris C Bernhardt acknowledges research support from the National Science and Engineering Research Council of Canada (NSERC Discovery-1304413), the CIHR (FDN-154298), SickKids Foundation (NI17-039), Azrieli Center for Autism Research (ACAR-TACC), BrainCanada, FRQ-S, and the Tier-2 Canada Research Chairs program. Drs. Bo-yong Park, Richard A I Bethlehem, Casey Paquola, and Boris C Bernhardt are jointly funded through an MNI-Cambridge collaborative award. The Neuroscience and Psychiatry Network (NSPN) study was funded by a Wellcome Trust award to the University of Cambridge and University College London. The data were curated and analyzed using a computational facility funded by an MRC research infrastructure award (MR/M009041/1) and supported by the NIHR Cambridge Biomedical Research Centre. The views expressed are those of the authors and not necessarily those of the NHS, the NIHR, or the Department of Health and Social Care.

Ethics

Human subjects: Participants provided informed written consent for each aspect of the study, and parental consent was obtained for those aged 14-15 years old. Ethical approval was granted for this study by the NHS NRES Committee East of England-Cambridge Central (project ID 97546). The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008.

Senior Editor

  1. Timothy E Behrens, University of Oxford, United Kingdom

Reviewing Editor

  1. Lucina Q Uddin, University of Miami University, United States

Reviewer

  1. Ben D Fulcher, University of Sydney, Australia

Publication history

  1. Received: November 7, 2020
  2. Accepted: March 30, 2021
  3. Accepted Manuscript published: March 31, 2021 (version 1)
  4. Version of Record published: April 30, 2021 (version 2)

Copyright

© 2021, Park et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 833
    Page views
  • 150
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Neuroscience
    Behrad Noudoost et al.
    Research Article

    Visually guided behavior relies on the integration of sensory input and information held in working memory (WM). Yet it remains unclear how this is accomplished at the level of neural circuits. We studied the direct visual cortical inputs to neurons within a visuomotor area of prefrontal cortex in behaving monkeys. We show that the efficacy of visual input to prefrontal cortex is gated by information held in WM. Surprisingly, visual input to prefrontal neurons was found to target those with both visual and motor properties, rather than preferentially targeting other visual neurons. Furthermore, activity evoked from visual cortex was larger in magnitude, more synchronous, and more rapid, when monkeys remembered locations that matched the location of visual input. These results indicate that WM directly influences the circuitry that transforms visual input into visually guided behavior.

    1. Developmental Biology
    2. Neuroscience
    Jessica Douthit et al.
    Research Article Updated

    As neural circuits form, growing processes select the correct synaptic partners through interactions between cell surface proteins. The presence of such proteins on two neuronal processes may lead to either adhesion or repulsion; however, the consequences of mismatched expression have rarely been explored. Here, we show that the Drosophila CUB-LDL protein Lost and found (Loaf) is required in the UV-sensitive R7 photoreceptor for normal axon targeting only when Loaf is also present in its synaptic partners. Although targeting occurs normally in loaf mutant animals, removing loaf from photoreceptors or expressing it in their postsynaptic neurons Tm5a/b or Dm9 in a loaf mutant causes mistargeting of R7 axons. Loaf localizes primarily to intracellular vesicles including endosomes. We propose that Loaf regulates the trafficking or function of one or more cell surface proteins, and an excess of these proteins on the synaptic partners of R7 prevents the formation of stable connections.