Development of Microstructural and Morphological Cortical Profiles in the Neonatal Brain

Abstract Interruptions to neurodevelopment during the perinatal period may have long-lasting consequences. However, to be able to investigate deviations in the foundation of proper connectivity and functional circuits, we need a measure of how this architecture evolves in the typically developing brain. To this end, in a cohort of 241 term-born infants, we used magnetic resonance imaging to estimate cortical profiles based on morphometry and microstructure over the perinatal period (37–44 weeks postmenstrual age, PMA). Using the covariance of these profiles as a measure of inter-areal network similarity (morphometric similarity networks; MSN), we clustered these networks into distinct modules. The resulting modules were consistent and symmetric, and corresponded to known functional distinctions, including sensory–motor, limbic, and association regions, and were spatially mapped onto known cytoarchitectonic tissue classes. Posterior regions became more morphometrically similar with increasing age, while peri-cingulate and medial temporal regions became more dissimilar. Network strength was associated with age: Within-network similarity increased over age suggesting emerging network distinction. These changes in cortical network architecture over an 8-week period are consistent with, and likely underpin, the highly dynamic processes occurring during this critical period. The resulting cortical profiles might provide normative reference to investigate atypical early brain development.


Introduction
Brain maturation over the perinatal period is rapid and complex. Although the majority of neurons are at their terminal location, synaptogenesis and myelination are ongoing, and (limited) migration of interneurons continues (Petanjek et al. 2011;Paredes et al. 2016). In the cortex, these processes are coordinated with sensory cortex and pathways developing earliest, prefrontal and association areas later (Flechsig of Leipsic 1901;Huttenlocher 1979). This developing architecture and connectivity is critical for efficient functional signaling, supporting and enabling the later development of cognitive and behavioral abilities.
Although alterations in perinatal neurodevelopmental processes have been associated with later cognitive and behavioral difficulties, quantifying myelo-or cytoarchitecture anatomical circuit development in the living human neonate is challenging. A relatively simple hypothesis of cortical connectivity is "similar prefers similar" (Goulas et al. 2016), meaning areas with similar cytoarchitecture preferentially connect. In this context, magnetic resonance imaging (MRI)-based methods such as regional structural covariance (Evans 2013) offer a proxy measure of brain connectivity, with structural similarity of spatially distinct regions of cortex reflecting coordinated maturation (Alexander-Bloch et al. 2013a). This inter-regional similarity is supported by similar genetic or maturational profiles (Alexander-Bloch et al. 2013b) and transcriptomic profiles (Yee et al. 2018) or associated with changes in disease (e.g., Zuo et al. 2018).
Understanding early structural regional similarity can shed light on how the brain develops into an efficient multifunctional system (Cao et al. 2016) and allows the detection of perturbations in normal development (Di Martino et al. 2014;Morgan et al. 2018). Recent years have seen a rise in studies on brain connectivity networks in the perinatal period (Zhao et al. 2019). Structural covariance networks (SCNs) based on measures of gray matter (GM) volume, cortical thickness (CT), cortical folding, and fiber density in the first 2 years of life have been described (Fan et al. 2011;Nie et al. 2014;Geng et al. 2017). Graph measures of these networks at birth show similar organizational properties to brain networks described in adulthood, suggesting that the nonrandom, efficient architecture of the brain is an inherent characteristic of the system, evident from very early on (Fan et al. 2011;Shi et al. 2012). Moreover, measures derived from GM volume covariance networks were found to differ in neonates with familial risk of schizophrenia . As different anatomical measures are capturing distinct developmental processes (Rakic 1988) and are controlled by different genetic mechanisms (Panizzon et al. 2009;Chen et al. 2013), different measures of cortical development show specific spatial and temporal patterns in early childhood Li et al. 2013;Nie et al. 2014;Lyall et al. 2015). Therefore, as might be expected, the resulting single-feature SCNs can be inconsistent.
An alternative approach using multiple anatomical measures to elucidate regional structural similarity has recently been examined in adolescents and adults Seidlitz et al. 2018). Compared with SCNs derived from a single technique, this approach allows for the construction of networks for individuals, rather than over a group, enabling the assessment of individual variability masked by group templates or case-control studies (Seghier and Price 2018). The resulting morphometric similarity networks (MSNs) have superior spatial consistency with cortical cytoarchitecture compared with networks based on diffusion imaging or on one structural measure (CT) ). Furthermore, cortical regions shown to be connected (i.e., similar) based on this approach have complementary expression of specific genes Seidlitz et al. 2018). In adults, these networks are also sensitive to alterations in common ) and rare  neurodevelopmental disorders. They also have neurobiological specificity; cortical MSNs differ between patients with genetic syndromes and are aligned with the regional expression of disorder-related genes . The only work we found to apply MSNs in the perinatal period is a recent study in neonates using a variant of MSNs, where authors were able to successfully predict postmenstrual age (PMA) at scan and differentiated infants born prematurely, with superior performance compared with using single predictive measures (Galdi et al. 2020).
The aim of this work therefore was to examine cortical organization from a multi-morphometric perspective, utilizing MSNs; hypothesizing this approach will provide a more sensitive and informative means to describe the postnatal structure and cytoarchitecture of the brain. We constructed cortical MSNs based on structural/morphological and microstructural diffusion parameters in a sample of 241 infant scanned between 37 and 44 weeks PMA. We hypothesized that, in the period following birth, brain development will be reflected in these networks. Moreover, we were interested in the community structure, or modularity in the neonatal brain, assuming these early clusters will provide the skeleton for developing functionality.

Subjects
This work included a sample of neonates participating in the Developing Human Connectome Project (dHCP) (http://www. developingconnectome.org/), scanned at the Newborn Imaging Centre at Evelina London Children's Hospital, London, UK. Images are available for download and analysis at the project website. This project has received ethical approval (14/LO/1169, IRAS 138070), and written informed consent was obtained from parents. At the beginning of this specific analysis (October 2018), 383 singleton term-born babies had undergone successfully structural and diffusion acquisitions, reconstruction, and early preprocessing. We next removed 87 subjects who failed the structural or diffusion pipelines as described below or who had missing data (e.g., T1 image). Finally, we removed subjects that have not passed diffusion QC: 17 subjects at the lower 5% of motion parameters and 46 subjects whose brains moved out of the field of view during scanning, leading to slices missing on the superior surface of the brain. One subject with major basal ganglia lesion was excluded, resulting in a final sample of 241 subjects (born at between 37 and 42 weeks; 128 males) with MR images acquired at PMA 40.92 ± 1.58 weeks (mean ± sd), range 37.43-44.71 weeks, henceforth "age." No major brain abnormalities were detected in review of the MRI data by a neonatal neuroradiologist in the subjects included.

Image Processing
Cortical surface processing and extraction of individual cortical features followed the pipeline described in Makropoulos et al. (2018)). Briefly, motion-and bias-corrected T2w images were brain extracted and segmented. White, pial, and midthickness surfaces were extracted, inflated, and projected onto a sphere. This was followed by estimation of cortical features including CT, pial surface area (SA), mean curvature (MC), and a gross proxy of myelin content, defined as the ratio between the T1w and T2w images (Glasser and Van Essen 2011), performed after registering the individual T1 and T2 images together, henceforth "myelin index (MI)," as detailed in Makropoulos et al. (2018). All brains were aligned to the 40-week dHCP surface template  using Multimodal Surface Matching (MSM) (Robinson et al. 2013(Robinson et al. , 2014, run with higher-order regularization constraints , to match coarse scale cortical folding (sulcal depth) maps. All other metrics were resampled to the template using this transformation and adaptive barycentric resampling (implemented using Human Connectome Project (HCP) tools, Connectome Workbench (https://www. humanconnectome.org/software/connectome-workbench).
Diffusion images were denoised (Veraart et al. 2016), Gibbsringing suppressed (Kellner et al. 2016), and corrected for subject motion and image distortion with slice-to-volume reconstruction in the multi-shell spherical harmonics and radial decomposition (SHARD) basis (Christiaens et al. 2019a), as described in Christiaens et al. (2019b) and using an image-based field map (Andersson et al. 2003). A tensor model for diffusion data was then fitted using a single shell (b = 1000), and fractional anisotropy (FA) and mean diffusivity (MD) maps were generated using MRtrix3 (Tournier et al. 2019). Neurite density index (NDI) and orientation dispersion index (ODI) maps were calculated using the NODDI toolbox (Zhang et al. 2012) as previously applied in the neonatal brain (Batalle et al. 2019). We acknowledge that the default parameters used to model NODDI here may not be the optimal for our sample (Guerrero et al. 2019). However, there is currently no agreed standard for infant-specific parameters; therefore we opted for the default to at least permit comparisons with previous literature.
Diffusion maps were registered onto individual T2w images using FSL's epi_reg (FLIRT) (https://fsl.fmrib.ox.ac.uk) and then projected onto cortical surface using Connectome Workbench in order to sample imaging features with the same spatial representation. All raw images were visually inspected for motion or image artifact and artifacted data excluded, and processed images were inspected for registration errors.

MSNs Construction
Cortical regions were defined using approximately equal-sized cortical parcellations with Voronoi decomposition at different granularities (n = 50, 100, 150, 200, 250, 300). This option was chosen over parcellation with a pre-existing atlas in order to minimize the effect of variable regional size when calculating the MSNs. For clarity and presentation, all results shown (other than Supplementary Fig. 1) were based on the n = 150 parcellation. Network construction followed steps introduced before Seidlitz et al. 2018): For each subject, each of the eight features was first z-scored across the cortex (to account for feature variation) and then averaged across each ROI in the parcellation, resulting in an eight-feature vector including mean normalized values of CT, MC, MI, SA, FA, MD, NDI, and ODI characterizing each node. Correlation (Pearson's r) between the eight-feature vectors for every two pairs of regions was calculated using Matlab, resulting in an n regions × n regions similarity matrix for each subject. A group structural similarity matrix was produced by averaging the 241 individual similarity matrices ( Fig. 1), and clustering was performed with affinity propagation. This process was also applied for a "leave-one-out" analysis, where similarity was based on a series of seven-feature vectors, leaving one feature out each time. In addition, typical SCNs for the individual features which are much more common in the literature were also generated by correlating each pair of regions across all subjects, resulting in a group similarity matrix that is based on a single morphometric feature. The values used in this analysis were the raw, un-normalized values.

Association Between MSNs with Age at Scan and Sex
To illustrate the effect of age on the individual morphometric features, Spearman's correlation (Spearman's rho, ρ) between the averaged regional single features and age was calculated and results corrected with false discovery rate (FDR) at 5% (Benjamini and Hochberg 1995). Nonparametric correlations were chosen for all age analyses as age at scan was not normally distributed (c) Each pair of regions is correlated using Pearson's r, resulting in a subject-specific similarity matrix with size n regions × n regions, which are then averaged to create a group mean similarity matrix, and the resulting group matrix is clustered using affinity propagation algorithm to examine network modularity.
The effects of age on internodal similarity were examined by correlations between age and the internodal edge-strength. In addition, we investigated the same correlation with mean node strength (the average of a node's edge-strength with every other node). Sex differences in MSNs for individual edge-strength and mean nodal edge-strength were explored using the Mann-Whitney test. In these contrasts, results were corrected with pFDR (Storey 2002;Storey and Tibshirani 2003), which is advantageous over FDR in larger samples.

MSN Clustering Analysis
To investigate network modularity, the group similarity matrix was clustered using affinity propagation (Frey and Dueck 2007, https://psi.toronto.edu/?q=tools). This clustering method has two merits in this current work: First, clustering could be performed without thresholding the group matrix, therefore avoiding information loss due to binarization, and second, it returns the optimal number of clusters (k) for a given preference value. In our case we used the median of the similarity matrix as preference value, making no prior assumptions: While a prespecified k is not mandatory, the algorithm does require an input chosen by the user, termed "preference value." By choosing the median of the similarity matrix as a shared value, we instruct the algorithm to consider all data points or in our case nodes as "exemplars" for the clusters. In addition, we also performed the same clustering approach but restricted the solution to a k = 7, the number of tissue classes in our von Economo atlas (see below), allowing us to compare this (fixed) clustering solution to clustering solutions from other modalities (see below). The rationale behind k = 7 clusters is also based on the maximum rank of the data in the context of a PCA, where the maximum number of components out of eight features from which the MSN was derived (as is the case here) is k − 1.
To measure robustness of the MSN clusters to the type of MRI modalities included, and to rank the importance of input modalities to the eventual solution, agreement between the clustering solution for the eight-feature regional vectors and the clustering solutions for seven-feature regional vectors (leaveone-out analysis) and clustering solutions for the eight single features (i.e., SCNs) was examined using normalized variation of information (Meilȃ 2007;Reichart and Rappoport 2009) using the 'nvi' code from the Pattern Recognition and Machine Learning Toolbox for Matlab (http://prml.github.io/). As in the leave-onefeature-out analysis, we fixed k = 7. In order to rule out the possibility that the clustering solution is derived only from age effects, we also performed clustering on the correlation between the edge-strength and age.
As MSNs in adults have been reported to correspond with cortical cytoarchitecture , we sought to examine the overlap between the resulting clusters with gross von Economo cytoarchitectural classifications (seven classes) (von Economo and Koskinas 1925) using the Dice coefficient to estimate overlap as in Whitaker et al. (2016). We also attempted to examine the spatial overlap between the seven clusters and von Economo classes using a more rigorous method: spatial permutations of the spherical representation of the cortical surface (Alexander-Bloch et al. 2018). We generated 500 null distributions for assessing correspondence between the maps as measured by normalized variation of information, while correcting for multiple comparisons as described in Alexander-Bloch et al. (2018).

Association Between MSN Clusters and Age at Scan
To investigate developing network integration and segregation, for each cluster the averaged edge-strength of connections within a cluster (with higher values indicating cluster distinction within the entire network) and the averaged edge-strength of connections between nodes within the cluster and nodes external to that cluster (indicating integration or segregation of that cluster) were calculated. Spearman's ρ coefficients were calculated for changes in these modular integration and segregation measures against age. Results were FDR corrected at 5% (Benjamini and Hochberg 1995).

Age-Related Associations with Single-Feature Maps
The correlation between parcellated single-feature maps and age at scan revealed metric-specific age associations: SA and MI had a brain-wide positive association with age, followed by CT, ODI, and NDI that also showed a positive association, but to a lesser extent. FA exhibited both positive and negative associations with age, while MD displayed only a negative association. Significant correlations with MC results were sparser, seen in limited perisylvian, frontal, and temporal regions (Fig. 2). Results are overlaid on a 41-week-old neonatal template .

Age-and Sex-Related Associations with MSNs
The correlation between inter-regional edge-strength of MSNs (n = 150 parcels) and age at scan showed both positive and negative associations throughout the brain. They numbered less in frontal and anterior temporal regions (see Fig. 3a for a node-bynode count of significant edges). For mean nodal edge-strength, this spatial gradient was clearer, with anterior cingulate and limbic regions negatively associated with age and lateral and medial parietal positively associated (Fig. 3b). Individual edge-strength did show some sex associations (Fig. 3c), though less extensive, with higher edge-strength in lateral frontal and temporal regions in males and cingulum and medial temporal regions in females (Fig. 3d). Age at birth (U = 6910.5, P = 0.55) and at scan (U = 6817, P = 0.44) did not differ between males and females.

Modularity of Neonatal MSNs
Using affinity propagation to perform clustering resulted in 12 broadly symmetric cortical modules, aligned with sensorymotor, fronto-temporal, anterior frontal, limbic, cingulate, insular, and visual systems (Fig. 4 left). Comparable spatial findings were also found when the initial parcellation included different number of regions (n = 50, 100, 150, 200, 250, 300), but the number of clusters increased linearly with parcellation density. For clarity, results are presented only for the middle density (n = 150).
A similar spatial partition was observed when the number of clusters was fixed at seven (Fig. 4 right) and across parcellation densities ( Supplementary Fig. 1). The main difference between the 7-and 12-cluster solution was a subdivision of the frontotemporal and occipital and parietal clusters in the k = 12 solution. Therefore, as these changes were minor, with no major loss of information on the overall emerging picture, and as a middle ground between complexity and interpretability, further analyses were performed on the k = 7 cluster solution.  To confirm the stability of this solution, clustering was calculated using bootstrapping of 20 random subjects each time, for 500 iterations, and individually for each subject. Node coincidence (how often one node fell in the same node as another) was very high (Supplementary Fig. 2b), and even in MSNs of individual infants, the consistency was very good (Supplementary Fig. 2a) with the whole sample cluster structure clearly evident in individuals and bootstrap samples.
Comparing the resulting MSNs (for k = 7) and von Economo cytoarchitectural classes, similarities were found, more remarkably in limbic and primary sensory areas, as compared with association areas (Fig. 5). Using the spatial permutation test, the normalized variation of information, between the MSN clusters and von Economo classes, was 0.85, ranging z = 0.84-0.93, suggesting relativity modest overlap; however this overlap was more significant than would be expected by chance (P = 0.02).

Feature Contribution
Examination of the clustering solution (for k = 7) for each of the single morphometric measures demonstrated that clustering of FA covariance had the highest agreement with the eightfeature solution (normalized variation of information, z = 0.73), while clustering of MC covariance had the lowest agreement (z = 0.94) (Fig. 6a). The top four metrics in particular are strongly correlated with each other cross-sectionally, derived from the same base diffusion-weighted images, and therefore this may indicate redundancy. This pattern remained when the singlefeature SCNs were calculated with partial correlation, taking into account age, suggesting this is not exclusively derived by age effects in the individual modalities (data not shown).
In the "leave-one-out" seven-feature solutions, the agreement between modules was much higher than for single modalities, as might be expected. Removing SA and MI from MSN estimation had the most profound effect (e.g., resulted in more different clusters), indicating their higher importance to the eventual solution. Removing either MD, NDI, or ODI had less of an effect, with near identical solutions (Fig. 6b).

MSN Clusters and Age at Scan
To demonstrate the age-related changes in cluster differentiation over time, we performed correlations between individual measures of within-module similarity and between-module similarity. Between-module analysis revealed six cluster pairs with significant positive correlation with age, showing betweencluster integration over time-occipital-parietal and anterior frontal, occipital-parietal and fronto-temporal, limbic and somatosensory-auditory, limbic and cingulate, anterior-frontal and cingulate, and insular-medial frontal and cingulate-while eight cluster pairs showed significant negative correlation with age. The limbic and fronto-temporal clusters demonstrated increased within-module similarity with age, suggesting cluster distinction, or increased internal similarity, while no withinmodule negative correlation with age was found (Fig. 7). Exact statistical values appear in Supplementary Table 1.

Age-Related Associations with Clustering Solution
Given the clear association between age at scan and connection strength between edges and clusters, we investigated the resulting MSN modules when calculated on edge-by-edge age trajectories. We performed clustering on the correlation matrix of each node-to-node MSN edge against age. This solution demonstrated both a clearly different result (Supplementary Fig. 3) and also has a different interpretation, showing areas that are similar in their direction and strength of maturation. These results showed more long-range connections and much more heterogeneity in frontal cortex. Primary sensory regions clustered together, and bilateral symmetric fronto-temporal and frontoparietal networks were also evident.

Discussion
Cortical microstructure and morphology go through extensive changes postnatally. Though the majority of neurons are in situ at term age, myelination and synaptogenesis continue to refine the cortical architecture that will persist throughout the lifespan. Using a densely sampled dataset across the perinatal period, we were able for the first time to combine multiple structural imaging features and showed their utility in exploring the dynamic progression of brain maturation. We demonstrated a clear and consistent modular structure, with both local and distributed networks of brain regions having similar cortical profiles. This modular structure showed hemispheric symmetry, with correspondence to known long-range cortical functional networks. This structure partly mirrored von Economo classes, indicative of an underlying cellular composition upon which specialized cortical systems will emerge. To our knowledge, this is the first study in neonates using a multiparametric approach to examine the community structure of the neonatal brain.
MSNs are able to capture the extensive changes occurring throughout a short period of just 8 weeks following birth. Age had a complex effect on edge-strength showing both positive and negative correlations, i.e., regional profiles became more and less similar with age. These age associations were confined to occipital, parietal, and temporal areas and were less evident in frontal areas. Whereas these associations were positive (greater morphometric similarity with the rest of the brain), limbic and anterior cingulate regions showed negative correlations with age (less similarity). These findings suggest that in the neonatal period, sensory and limbic areas and posterior parietal regions have the largest maturational changes, as opposed to regions related to higher executive functions (e.g., prefrontal cortex), qualitatively similar to the pattern seen in Lebenberg et al. (2019) with older infants. Changes may only be evident in higher-order areas at a later stage, in line with the more prolonged synaptogenesis, dendritic arborization, and myelination patterns (Huttenlocher 1990;Huttenlocher and Dabholkar 1997;Teffer and Semendeferi 2012). This may also reflect work in cortical GM volume maturation in neonates, showing a posterioranterior gradient in the first weeks following birth (Gilmore et al. 2007). In all, age-dependent changes observed in MSNs in the newborn are in line with neurobiological findings, as well as other neuroimaging findings.
In a mixed sample of term and preterm infants investigating MSNs (derived from structural and extended diffusion metrics, Figure 7. Inter-and intramodular similarity changes with age at scan: association between clusters' edge-strength and age. Width of line is indicative of the strength of significant association. but not surface measures), Galdi et al. (2020) also found that brain structures became both more similar and dissimilar to each other over the same age range. While that study used different feature vectors, atlas, and populations here, findings are clearly complementary: Occipital connections increased with age, parietal edges were also more positively than negatively associated with age, while temporal edges showed mixed associations with age. In contrast to our findings, frontal regions also showed a strong relationship with age; however, this was more evident in white matter (WM) rather than GM. By choosing to focus only on healthy term-born infants, we attempted to minimize any effects of preterm status confounders on brain maturation (Ball et al. 2017;Batalle et al. 2017).
By exploring the developmental patterns of the individual metrics on which cortical covariance was based, as well as their contribution to the resulting cortical patterning, we were able to dissect early micro and macro properties of the brain. When examining each morphometric measure in isolation, unique developmental trajectories emerged, replicating prior work (Dean et al. 2017;Batalle et al. 2019). However, we also show that no one measure was able to capture the abundance of unique information that is reflected in the combination of morphometric features. Typical measures of cortical anatomy, thickness, and area, as well as myelin content, showed strong increase with age . Age however showed little significant association with curvature, in agreement with previous reports showing relatively little change in MC compared with other surface measures in the months following birth Batalle et al. 2019), suggesting that cortical folding remains relatively static in this period. Though no single individual measure mirrored the clustering solution achieved by combining all measures, the diffusion metrics had the better correspondence ( Fig. 6a) but were apparently less important to the clustering solution (Fig. 6b). This apparent contradiction is due to redundancy. For example, MD and NDI show very similar patterns of correlation with age (Fig. 2), so removal of one or two diffusion measures does not reduce the effective information used to calculate correlations in MSNs.
Cortical diffusion metrics show complex maturational trajectories around term birth. Ball et al. (2013) and Batalle et al. (2019) have shown that MD decreases linearly with age, while FA decreases rapidly until week 38 PMA but then begins to slowly increase. This might explain why at the age range examined here (37-44 weeks) we observed both negative and positive associations between FA and age. Moreover, Batalle et al. found that NDI is positively associated with age, similar to the results here. ODI seemed to stabilize after 38 weeks, whereas we still observe a strong association with age after 38 weeks. We might hypothesize that this is related to the focus on preterm neonates in that work, where microstructural trajectories could differ due to earlier birth and the very different early environmental exposures that entail (both clinical and ex utero related).
Sex differences were less widespread compared with age effects. While sex differences in cortical volumes are evident from birth (Gilmore et al. 2007), they are subtle when looking at single modality covariance networks in early development (Zielinski et al. 2010;Geng et al. 2017). Using the MSNs approach, we were able to detect sex effects in some cortical areas, most pronounced in frontal and temporal regions. This corresponds to results of a neonatal study using tensor-based morphometry (Knickmeyer et al. 2014), and even in adults, multimodal imaging sex differences were mainly localized in the frontal lobe, followed by parietal and temporal lobes (Feis et al. 2013), so this may represent a consistent gross pattern throughout the lifespan.
Generally, any network construction and follow-up clustering analyses are highly reliant on several factors, including population and age range examined Morgan et al. 2019), type of morphometric measures used , parcellation (Arslan et al. 2018), threshold (Bordier et al. 2017), wiring cost (Betzel et al. 2017), and type of modularity estimation (Sporns and Betzel 2016) and are therefore likely to vary accordingly. For module definition, we applied affinity propagation to discover meaningful cortical clusters, without eliminating or changing any connections for network construction. In young adults, utilizing MSNs and the Louvain algorithm, Seidlitz et al. (2018) reported four cortical modules consistent with lobular division. Using affinity propagation to cluster similar developing regions in terms of CT and curvedness between the ages 3 and 20 years, Nie et al. (2013) found eight and five clusters, respectively, also with some alignment to lobular division, and indeed our parcellation based just on CT was qualitatively similar.
While most of the resulting modules presented in these studies tend to be spatially contiguous and local, our clustering solution showed also long-range connections and might be more indicative of possible functional and structural connectivity. The modules were broadly aligned with known functional systems and were relatively stable: These included sensorymotor, fronto-parietal, temporal, limbic, cingulate, and visual regions. Our correspondence between MSNs and functional systems may reflect the underlying architecture for later developed functional networks (van den Heuvel et al. 2015;Geng et al. 2017). Our resulting MSN clusters show some overlap to the von Economo tissue classifications, implying the possible origins of the cortical profiles obtained in our analysis, and provide some reassurance of their validity. Further supporting this, the solution revealed here also shows several analogies to the clustering solution of genetic contribution to SA and CT reported by Chen et al. (2013) and might be especially relevant in early development where external environmental effects on these parameters are still relatively small.
The clustering solutions were spatially robust, by and large not affected by age (confirmed by clustering separately only individuals scanned in the top and bottom PMA-range quartiles), nor was it affected by the neonates' ex utero experience, as the clustering solution for neonates scanned within a week from birth resembled that of the entire sample (data not shown). We also show that the clustering solution is not driven exclusively by age-related changes in the morphometric measures ( Supplementary Fig. 3). So, although age did not alter the location of structure of the MSNs, it was associated with their internal coherence. Age was linked to increased regional similarity within the clusters and both increased and decreased regional similarity between clusters (cluster integration and segregation). Though the exact pattern was complex, the summary was that cingulate showed integration with the limbic and insular clusters while segregating from most of the neocortical clusters, perhaps suggesting an advancement towards a more broad representation of paralimbic structures, while the limbic and fronto-temporal clusters showed also age-dependent increase in intrasimilarity.
While our study's advantages include a large neonatal sample, data acquired and analyzed with optimized protocols for this age group, as well as the use of a number of morphometric measures to characterize the brain, it is not without limitations. In order to fully describe development, longitudinal data is needed. Due to the nature of the surface morphometric measures utilized, we were not able to assess whole-brain connectivity, here excluding subcortical structures, cerebellum, and brain stem. Future work should investigate the structural and functional coupling of neonatal MSNs, as well as their genetic and environmental correlates. The implication of recent work using MSNs is that changes in neurodevelopmental disorders must occur early in development. Using a large sample, our study provides a means to detect meaningful alterations in structural coupling at an early stage in infants with either a known genetic or high likelihood to have a later neurodevelopmental disorder.
To conclude, we present cortical MSNs in the neonatal brain, their association with age, and their community structure, characterizing evolving cortical tissue architecture. Following birth, age is strongly related to edge-strength, in a posterior-anterior gradient reflecting the directionality of axonal and dendritic growth, myelination, and synaptogenesis in this period. Clustering of MSNs into modules demonstrated correspondence to known functional systems and cytoarchitectural atlases and was relatively stable across subjects and over the age range examined. Within these clusters, similarity generally seemed to increase postnatally, implying increased network coherence in this short period. We also find a complex relationship of between-clusters similarity, suggesting some rearrangement in segregation and integration of the network, indicative of cortical maturational processes. Here we provide a template of neonatal structural cortical co-variance profiles. By utilizing both shape and microstructural information, we highlight the complex nature of perinatal cortical development.

Supplementary Material
Supplementary material is available at Cerebral Cortex online.