Developmental topography of cortical thickness during infancy.

During the first 2 postnatal years, cortical thickness of the human brain develops dynamically and spatially heterogeneously and likely peaks between 1 and 2 y of age. The striking development renders this period critical for later cognitive outcomes and vulnerable to early neurodevelopmental disorders. However, due to the difficulties in longitudinal infant brain MRI acquisition and processing, our knowledge still remains limited on the dynamic changes, peak age, and spatial heterogeneities of cortical thickness during infancy. To fill this knowledge gap, in this study, we discover the developmental regionalization of cortical thickness, i.e., developmentally distinct regions, each of which is composed of a set of codeveloping cortical vertices, for better understanding of the spatiotemporal heterogeneities of cortical thickness development. We leverage an infant-dedicated computational pipeline, an advanced multivariate analysis method (i.e., nonnegative matrix factorization), and a densely sampled longitudinal dataset with 210 serial MRI scans from 43 healthy infants, with each infant being scheduled to have 7 longitudinal scans at around 1, 3, 6, 9, 12, 18, and 24 mo of age. Our results suggest that, during the first 2 y, the whole-brain average cortical thickness increases rapidly and reaches a plateau at about 14 mo of age and then decreases at a slow pace thereafter. More importantly, each discovered region is structurally and functionally meaningful and exhibits a distinctive developmental pattern, with several regions peaking at varied ages while others keep increasing in the first 2 postnatal years. Our findings provide valuable references and insights for early brain development.

Early postnatal brain development, especially during the first 2 postnatal years, features dynamic growth in brain structures and rapid expansions in behavioral and cognitive abilities (1–3). For example, cortical thickness (CT) develops spatially heterogeneously, with an average increase of 40% during the first 2 y of age (4, 5). It has been found that the early development of CT positively links to later intelligence quotient, as well as the general factors of intelligence in healthy individuals in the late childhood, adolescence, and adulthood (6). Development of CT in the frontal regions has been shown to be predictive of cognitive outcomes in typically developing children and adolescents (7). Meanwhile, increasing prevalence of neurodevelopmental disorders emphasizes the importance and vulnerability of this critical period compared with the entire lifespan (3). For example, abnormalities of CT have been observed in many neurodevelopmental disorders, such as schizophrenia, autism, attention deficit/hyperactivity disorder, Down syndrome, fragile X syndrome, and Angelman syndrome (5, 8), which are likely rooted in early postnatal brain development. Consequently, charting developmental patterns of CT in healthy infants is of great importance in understanding the underpinnings of neurodevelopmental disorders, identifying early biomarkers, and planning early targeted interventions.

Various efforts have been made to study the development of CT, which, however, have shown inconsistent results. Several early studies (9) suggested that CT increases from age 6 to 10 y before starting to thin at 10 to 12 y of age. While more recent studies (10–12) reported monotonic decrease of CT since 3 y of age, which is the youngest age in their datasets. Li et al. (4) reported that CT presents dynamic increase in the first postnatal year, followed by region-specific increase or decrease in the second postnatal year, based on a longitudinal dataset acquired at age 0, 1, and 2 y. Taken together, it is suggested that CT likely peaks between 1 and 2 y of age (3). However, due to the difficulties in acquisition and processing of the infant brain MRI, which typically exhibits extremely low tissue contrast and dynamic imaging appearance, our knowledge still remains very limited regarding the dynamic growth, peak age, and regional heterogeneities of CT during the first 2 y after birth.

Therefore, in this paper, we investigated detailed spatiotemporal developmental patterns of CT during infancy, by leveraging 210 longitudinal MRI scans from 43 typically developing infants, with each having up to 7 densely sampled time points during the first 2 postnatal years (i.e., around 1, 3, 6, 9, 12, 18, and 24 mo of age). Instead of using predefined cortical regions based on prior knowledge from adult populations, we discovered the first topographical regionalization of CT naturally formed during early development by using a data-driven nonnegative matrix factorization (NMF) method (13), which is able to decompose the cerebral cortex into interpretable and meaningful regions. As we applied the NMF method onto all of the longitudinal scans at different ages from all subjects, the discovered cortical topography indicates a set of developmentally distinct regions across subjects and ages. Each of these regions is composed of a group of cortical vertices with codeveloping CT, thus revealing the spatial heterogeneity of the CT development. For each discovered region, we further modeled the region-specific developmental trajectory and peak ages of CT.

Results

We studied a population of 43 infants with a total of 210 longitudinal MRI scans in the first 2 postnatal years. The structural MR images were processed using an infant-dedicated pipeline (14) to reconstruct the cortical surfaces and compute the CT maps. The NMF method (13) was used to discover the developmental regionalization map of CT.

Topography of Developmental Regionalization.

The appropriate region number K was determined according to 3 criteria shown in Fig. 1 . As we mainly focus on relatively large-scale patterns of the developmental regionalization of CT, we set the largest K as 30. We observed that both the reconstruction error (Fig. 1) and instability (Fig. 1) reached local minima at K = 6 and 17, suggesting stable regionalization patterns. This is further confirmed by the silhouette coefficient (Fig. 1), which reaches a local peak at K = 6 and a plateau starting from K = 15 to K = 17. Therefore, we finally chose both K = 6 and K = 17 to present our discovered developmental regionalization of CT, with a primary focus on K = 17 to provide a finer pattern of developmental regionalization. Meanwhile, to inspect the coarsest pattern of developmental regionalization, we also show the result with K = 2.

Fig. 1.

Three criteria for determining the region number K. The reconstruction error (A), instability (B), and silhouette coefficient (C) with respect to different region numbers. Yellow dots indicate the local minimums of reconstruction error and instability, and the local maximums of silhouette coefficient.

As shown in Fig. 2, all of the regions identified with K = (2, 6, 17) are bilaterally located in a quasisymmetrical fashion. Of the 2 regions with K = 2 (Fig. 2), the first region encompasses the lateral and medial prefrontal, lateral temporal, and inferior parietal areas, and the second region consists of the sensorimotor, visual, sensory association (superior parietal cortex), insular, and visual cortices. In the case of K = 6, the discovered regions still exhibit bilaterally relatively symmetric patterns (Fig. 2 and ). These regions approximately encompass: 1) perisylvian areas, inferior parietal lobules, and posterior cingulate cortex; 2) medial occipital and dorsal sensorimotor areas; 3) insula and orbitofrontal areas; 4) medial prefrontal areas and superior parietal lobules; 5) middle, inferior, and medial temporal cortices and fusiform; and 6) dorsal frontal cortex and temporal pole, respectively. When increasing K to 17, the discovered regions still exhibit bilaterally symmetric patterns but only contain spatially continuous areas (Fig. 2 and ). These regions approximately cover: 1) perisylvian areas, 2) medial occipital cortex, 3) medial orbitofrontal cortex, 4) medial prefrontal cortex, 5) medial temporal areas and fusiform, 6) temporal pole, 7) precuneus, 8) inferior parietal lobules, 9) middle insula and anterior superior temporal cortices, 10) lateral orbitofrontal cortices and anterior insula, 11) middle and posterior cingulate cortices, 12) dorsal somatosensory area, 13) inferior frontal triangularis and opercularis, 14) superior parietal lobules, 15) posterior temporal and lateral occipital cortices, 16) sensorimotor areas, and 17) paracentral and superior frontal areas, respectively. Many large or disjoint regions at K = 6 are separated into smaller ones at K = 17. For example, the region 6 at K = 6 is divided into subregions 6, 13, and 17 at K = 17, and the region 1 at K = 6 is divided into subregions 1, 8, and 11 at K = 17.

Fig. 2.

Discovered regions with different region numbers: K = 2 (A), K = 6 (B), and K = 17 (C), where warmer colors correspond to higher values. In B and C, each small brain represents a region, placed according to the center of its corresponding location on the big brain.

Developmental Pattern and Peak Age of Each Discovered Region.

The longitudinal developmental trajectory of the average CT in each hemisphere of each subject is shown as a red line (female) or blue line (male) in Fig. 3 (the top left box). According to the general cross validation (GCV) (15) errors, the general additive mixed models (GAMM) (16) better fitted the population’s developmental trajectories of the average CT in both hemispheres than other 3 parametric models (detailed GCVs are reported in ), and the estimated population’s developmental trajectory is delineated using green dashed curves in Fig. 3 (the top left box). This result suggests that the average CT first increases rapidly after birth and then reaches peak at around 426 d (14 mo), as indicated by the yellow hexagons and arrows, before it decreases at a slow pace. No gender difference was found on the developmental trajectory of the average CT in either hemisphere (P = 0.91 for the left hemisphere and P = 0.75 for the right hemisphere).

Fig. 3.

Developmental trajectories of the average CT of each hemisphere (in the top left box) and each discovered region (i.e., component) in the left hemisphere, as well as the region-specific peak ages. The corresponding result of each region in the right hemisphere is shown in . The y axis stands for CT (1.5 to 4.5 mm for all components), and the x axis represents the age in days. Red lines and blue lines represent females and males, respectively. The dashed green curves illustrate the fitted model of the population’s trajectories. The peak point of each fitted curve is signified using a yellow hexagon and an arrow.

For each discovered region in the regionalization map with K = 17, its longitudinal trajectory of the region-averaged CT of each subject in the left hemisphere is shown in Fig. 3 (i.e., components 1 to 17; herein, a “component” indicates a region). For these regions, the population’s developmental trajectories of the average CT were also best fitted by using the GAMM (detailed GCVs are reported in ), shown as green dash curves in Fig. 3. The estimated peak age of CT in each region is indicated by a yellow hexagon and an arrow. Again, no significant gender difference on the developmental trajectory of CT was found in any region (the P values of the gender differences are reported in ). The peak age and developmental trajectory of each region in the right hemisphere largely resemble those of the corresponding regions in the left hemisphere, which can be found in . The parameters of the fitted curves are detailed in .

As we can see, CT exhibits region-specific and bilaterally quasisymmetric developmental patterns in terms of the starting point, peak age, and developmental velocities. At birth, regional differences are already present (Figs. 3 and 4b), with the CT decreasing along the anterior to posterior (A P) direction and the ventral to dorsal (V D) direction. Such a spatial pattern is largely maintained during the first 2 y. Specifically, thin cortices are located in the sensorimotor cortex (components 12 and 16), parietal lobules (components 8 and 14), middle temporal areas and fusiform (component 5), posterior temporal areas, lateral and medial occipital areas (components 2 and 15).

Fig. 4.

(A) Regional peak days of CT development. (B) Regional CT development at different time points. (C) Regional velocity of CT development.

During the first 2 postnatal years, several regions have a limited amount of CT change (Fig. 3), such as the sensorimotor, paracentral, superior parietal, and primary visual cortices. Several regions show more dramatic changes in CT, including the medial orbitofrontal cortices (component 3), medial prefrontal cortices (component 4), lateral orbitofrontal areas and anterior insula (component 10), as well as the posterior temporal and lateral occipital cortices (component 15). Meanwhile, the development rate, or velocity, of CT is distributed not only spatially heterogeneously but also heterochronously, as shown in Fig. 4. Across all cortical regions, the velocity of CT initially presents the highest positive values at the first month (increasing as much as 0.1 mm per month) during the first 2 y. During this process, the prefrontal cortex, especially its medial part, has the highest velocity in the whole period between 0 and 12 mo, followed by the lateral frontal, posterior temporal, and lateral occipital cortices. The sensorimotor and primary visual cortices increase with low velocities and their developmental trajectories are rather flat (Fig. 3). After 12 mo, for the regions that reach the peak during the second year, their velocities gradually decrease to zero (reaching peak), and change to negative values thereafter. For other regions that do not reach the peak during the first 2 y, their velocities decrease to very small yet still positive values after the first year, suggesting very slow but continuous increase during the second year. Regions exhibiting faster decreases of CT include the lateral temporal, insula and anterior frontal, medial orbitofrontal, and cingulate cortices. Regions showing slow but continuous increase during 12 and 24 mo include the sensorimotor, superior frontal, and superior parietal areas. Overall, CT appears to decrease relatively faster at 18 than 24 mo (Fig. 4).

The peak days of CT of different regions are pointed out by the yellow arrows in Fig. 3 and mapped on the brain in Fig. 4. The regions that do not peak between 0 to 24 mo are colored in white. As can be observed, most regions peak between 350 to 457 d (11.5 to 15 mo), and the regional peak ages on the left hemisphere largely mirror the results on the right hemisphere. The orbitofrontal cortex (component 3) arrives the very first to peak at 350 d (11.5 mo), followed by several regions that peak between 360 and 370 d (12 mo), including the insula (component 9), perisylvian areas (component 1), and posterior temporal and lateral occipital cortices (component 15). Regions that peak relatively late include the inferior parietal lobule, middle and posterior cingulate cortices, and inferior frontal and medial occipital areas (components 8, 11, 13, and 2), which all peak after 14 mo. Other regions show continuous increase without reaching a peak, such as the sensorimotor, superior frontal, superior parietal and medial prefrontal areas (components 4, 7, 12, 14, 16, and 17). To visualize the development of CT in detail, we additionally fitted the developmental trajectory for each cortical vertex and mapped the vertex-wise fitted CT values onto the cortical surface at different ages, as shown in Fig. 5.

Fig. 5.

Vertex-wise maps of CT development at different ages.

Discussion

In this paper, we revealed the cortical topography of developmental regionalization of CT during infancy by discovering vertices that are codeveloping across ages and subjects using a data-driven method. Our discovered regions based on infant CT development are anatomically and functionally meaningful, yet differ from previous cortical areas defined by sulcal-gyral patterns, cytoarchitecture, or functional connectivity. Several distinct characteristics of these discovered regions are discussed below.

Most regions are located rather symmetrically in both hemispheres from K = 2 to K = 17 (Fig. 2), even though the method was conducted on both hemispheres simultaneously, without imposing any constraint on the left–right corresponding (homotopic) vertices. Consequently, our results indicate a predominantly symmetric pattern of CT development during early infancy. Such symmetric patterns in development of CT and genetic organization of CT have also been found in adolescents (17) and adults (18) when using comparable region numbers. This is also consistent with the functional parcellation of the cerebral cortex from nearly all of the cohorts with variable ages, where most of the functional parcels have been found to be symmetric (19). A possible mechanism could be the dense interhemispheric structural connections, which may exert certain constraints to the cortices of both sides, making their developmental patterns similar.

At the coarsest level of regionalization with K = 2, as shown in Fig. 2, the 2 discovered regions are strikingly well aligned with the spatial pattern of evolutionary cortical surface expansion between adult macaques and humans (20). The region 1 in Fig. 2 mainly includes the lateral and medial prefrontal, lateral temporal, and inferior parietal areas, which largely correspond to high-expanding areas; while the region 2 contains the sensorimotor, sensory association (superior parietal cortex), insula, and visual cortices, which approximately correspond to low-expanding areas. Moreover, the region 1 consistently shows low CT but high velocity of CT development during the first 2 y, which has also been found to exhibit high-growth rate of cortical folding during infancy (21) and high intersubject variability of functional connectivity in both infants and adults (22, 23). In contrast, the region 2 develops in a rather low velocity during infancy, which corresponds to regions with low-growth rate of cortical folding (21) and low intersubject variability of functional connectivity (22, 23). These results suggest that these 2 discovered regions indeed exhibit distinct characterizes in brain development, evolution, and intersubject variability, which are likely associated with their underlying cellular, genetic, and functional differences.

At the finer level of regionalization with K = 17, several regions in high-order association cortices show consistency between our developmental regionalization of CT in infants and the genetic topography of CT in adult twins (18). For example, the component 13 (inferior frontal, triangularis, and opercularis) resembles the cluster 11 in ref. 18, the component 5 (medial temporal and fusiform) corresponds to the cluster 9 in ref. 18, the component 8 (inferior parietal) corresponds to the cluster 3 in ref. 18, and the components 6 and 9 form the cluster 8 in ref. 18. Such relatedness may suggest that the early development of CT in these high-order association cortices is more affected by genetic factors than other regions, e.g., sensorimotor and visual regions. However, it would be better to chart the infant-specific genetic topography of CT to provide a more meaningful comparison.

It is interesting that the developmental regionalization of CT during infancy in the present study manifested quite distinct patterns with that of adolescents (17). For example, at K = 2, the cortical regionalization in infants shows a rough separation between the unimodal regions and the high-order association regions, while an anterior–posterior division is shown in that of adolescents. At K = 17, most cortical regions show very different boundaries between infants and adolescents, except the medial occipital cortex (i.e., primary visual area). Such distinction could stem from the fact that these 2 studies investigated different developmental periods with different underlying microstructural mechanisms.

It should be noted that the sensorimotor area in infant CT regionalization was distinguished as one region at both K = 2 and K = 6, but was split along the dorsal–ventral direction at K = 17. This separation is distinct from anterior–posterior separations that were commonly observed in adolescents (17). According to the associated functions, the component 12 (the dorsal part) mainly involves the sensing and controlling of the lower body, trunk, and arms, while the component 16 (the ventral part) is mainly in charge of hands, fingers, head, and other small muscle groups. According to the fitted CT developmental trajectories shown in Fig. 3, the dorsal part continuously shows lower CT compared with the ventral sensorimotor area. Based on the CT development velocity shown in Fig. 4, the ventral part initially shows higher development velocity in the first 12 mo and later shows lower development velocity in the second year, compared with the dorsal part of the sensorimotor area. These observations suggest that the ventral sensorimotor develops earlier than the dorsal part, which is consistent with the cephalocaudal direction, i.e., from head to toe, in the motor skill development of infants (24).

Previous studies of CT during infancy either have very sparse acquisition time points or only use MRI scans after 2 y of age. The main challenge lies in processing infant brain MR images, which typically exhibit low tissue contrast and dynamic appearance (2), thus rendering it difficult to precisely chart the developmental trajectories of CT during this period. We addressed these issues by leveraging a densely sampled longitudinal dataset scanned starting from birth and a set of infant-dedicated computational tools that have been extensively validated in many infant studies (4, 5, 21, 25). We revealed that the development of the average CT of the entire cortex follows an inverted U-shaped trajectory during the first 2 postnatal years, and the overall peak of CT is reached at around 14 mo of age. Our finding is also supported by previous reports that have suggested that CT increases dynamically in the first year while changing subtly in the second year, according to the longitudinal MRI scans at 0, 1, and 2 y of age (4, 5). Moreover, our finding is also in line with the recent hypothesis that CT likely peaks between 1 and 2 y of age (3), which was inferred from latest brain development studies in both in vivo MRI and postmortem histological research.

Although the underlying mechanisms that drive CT changes remain unclear, the initial dynamic increase of CT after birth is likely associated with a set of complicated microstructural changes, e.g., the increases in dendritic arborization, axonal elongation and thickening, synaptogenesis, and glial proliferation (26). The following cortical thinning might be driven by complex yet age-related factors that occur at different ages. For example, during infancy, it is assumed that the decrease of CT may be related to synaptic pruning that follows the synapse overproduction (5). For old children and adolescents, the cortical thinning is also likely related to the myelination and decrease of brain volume (27). Contrary to cortical thinning related to brain atrophy for elder adults, cortical thinning at young ages is consistently believed as a part of normal brain maturation, and earlier cortical thinning is likely a positive sign indicating higher cognitive outcome for adolescents (6).

We found that the major patterns of the spatial distribution of thin and thick cortices have been largely established at birth and even earlier, which is well kept during infancy, consistent with the findings in refs. 4, 28, and 29. Specifically, the prefrontal areas, including the components 4, 10, and 13, as well as the superior temporal cortex, including the component 15, are observed with relatively thick cortex throughout the first 2 y, consistent with the discoveries in fetal cortical development (28). Meanwhile, the fetal prefrontal cortex shows significant change in CT (28), while our study shows that CT of this region increases continuously at high rate throughout the first postnatal year. The components 2, 12, and 16, corresponding to the visual and sensorimotor cortices, continuously exhibit low CT and low CT development velocity during the first 2 y; a similar pattern has also been shown in fetuses of 40 postmenstrual weeks (29).

Meanwhile, all cortical regions show distinct developmental trajectories—although the whole-brain average CT and most regions reach the peak during the first 2 y, several regions (e.g., the sensorimotor, sensory association, motor association areas) tend to continuously increase in CT and do not reach a peak during this period. Among all of the regions, the medial orbitofrontal cortex peaks the very first, consistent with previous findings, where CT in this region has been found to decrease most remarkably from year 1 to year 2 (4). Several primary cortical regions develop relatively earlier, such as the auditory cortex and visual cortices, which are found to peak during the first 2 y. While other primary cortices, e.g., the sensorimotor, as well as motor association and sensory association areas, develop later and do not reach a peak during the first 2 y. This finding is in concordance with ref. 4, which revealed that, in the second postnatal year, CT in the visual and auditory cortices tends to decrease, while CT in the motor cortex likely increases. Our finding is also consistent with ref. 30, in which the posterior midbody of the human corpus callosum shows slow increase in area during infancy and arrives a plateau at around 70 mo of age, compared with that most parts of the corpus callosum stabilize at around 40 mo of age. In our study, the sensorimotor cortex, to which the posterior midbody of the corpus callosum connects, also shows slow development in CT compared with other regions, suggesting a potential relatedness between their early development. However, whether the development of CT in this region after 24 mo would follow a similar growth trajectory as the posterior midbody of the corpus callosum needs to be verified.

Limitations.

Certain limitations in this study should be noted. First, the brain MRI in this study was collected at 1-mm resolution, which likely hampers the accuracy in measuring infant CT that is usually between 1.5 to 4 mm. High-resolution images are highly desired for infant studies. Second, our infant study did not directly compare with adult studies. According to the average CT of 2.63 mm in the Human Connectome Project (31), we find that infant CT initially presents around 80% of the adult level at birth. At the peak level, the infant average CT is around 103% of the adult level (i.e., 2.7 mm) and decreases to a near adult level (101%) at the end of the second year. However, a rigorous comparison between infants and adults should be performed by using infant and adult data with the similar imaging resolution from the similar scanner, following the similar imaging protocol and data-processing pipeline. This is an important topic of our future work. Third, although our dataset has densely sampled time points, the scan time of subjects is relatively concentrated at several ages. Therefore, the modeled CT at in-between time points may be less accurate than those at scanned time points. These limitations can be resolved in the near future by the Baby Connectome Project (32), which is an accelerated longitudinal study acquiring high-resolution, multimodal MR images, thus allowing more detailed investigations of baby brain development.

Conclusions

In summary, this study has 2 main contributions. First, it has unprecedentedly charted the topography of developmental regionalization of CT during infancy, thus delineating developmentally distinct regions. Second, this study has revealed the region-specific trajectory and peak age of CT during the first 2 y. These discoveries confirm the inferences drawn from the recent related imaging and histological studies and provide an important reference for future studies regarding both normal and abnormal brain development.

Materials and Methods

Subjects.

This study was approved by the Institutional Review Board at the University of North Carolina (UNC) at Chapel Hill, School of Medicine. Pregnant mothers were recruited during the second trimester of pregnancy from the UNC hospitals and informed consents were obtained from all of the parents. All infants in this study were free of congenital anomalies, metabolic diseases, and focal lesion. All infants were scanned during natural sleep without sedation, with their heads secured in a vacuum-fixation device. A total of 210 longitudinal brain MRI scans at around 1, 3, 6, 9, 12, 18, and 24 mo of age were acquired from 43 term-born infants (with gestational ages 261 to 294 d), including 21 males and 22 females. For further details, see .

Discovering Spatiotemporal Heterogeneity of CT Development.

A data-driven approach was used to discover the cortical regions formed during the CT growth during infancy. Our aim is to discover the infant-specific cortical topography of developmental regionalization of CT by grouping codeveloping cortical vertices into the same region. To this end, we adopted the NMF method (13) by including the whole longitudinal course of CT of all subjects into a large data matrix . By decomposing using NMF as , the resulting nonnegative elements in each column of the base matrix naturally point out a group of cortical vertices jointly developing across subjects and ages, thus indicating a distinct region during the developmental regionalization. To find the appropriate region number K for the NMF method, we jointly considered 3 widely adopted criteria, i.e., reconstruction error, instability (33), and silhouette coefficient. A region number was chosen based on all 3 criteria, i.e., a low reconstruction error and instability, as well as a high silhouette coefficient. Further details are presented in .

Charting Longitudinal Developmental Trajectories.

We adopted 3 parametric models, i.e., a linear, a quadratic, and a sigmoid model, and one nonparametric model, i.e., the GAMM (16), to fit the complex developmental trajectory of CT in each discovered region, given that CT increases dynamically during the first year and then exhibits region-specific increase or decrease during the second year (4). The GCV error (15) was used to determine the best-fitted model on our data and was estimated according to ref. 15, as it has been commonly used as a metric for comparing nonparametric with parametric models (15). More information regarding the model fitting and model statistics is presented in .