Generation of Infant-dedicated Fine-grained Functional Parcellation Maps of Cerebral Cortex

Infancy is a dynamic and immensely important period in human brain development. Studies of infant functional development using resting-state fMRI rely on precisely defined cortical parcellation maps. However, available adult-based functional parcellation maps are not applicable for infants due to their substantial differences in functional organizations. Fine-grained infant-dedicated cortical parcellation maps are highly desired but remain scarce, due to difficulties ranging from acquiring to processing of infant brain MRIs. In this study, leveraging 1,064 high-resolution longitudinal rs-fMRIs from 197 infants from birth to 24 months and advanced infant-dedicated processing tools, we create the first set of infant-specific, fine-grained cortical functional parcellation maps. Besides the conventional folding-based cortical registration, we specifically establish the functional correspondences across individuals using functional gradient densities and generate both age-specific and age-common fine-grained parcellation maps. The first set of comprehensive brain functional developmental maps are accordingly derived, and reveals a complex, hitherto unseen multi-peak fluctuation development pattern in temporal variations of gradient density, network sizes, and local efficiency, with more dynamic changes during the first 9 months than other ages. Our proposed method is applicable in generating fine-grained parcellations for the whole lifespan, and our parcellation maps will be available online to advance the neuroimaging field.


Introduction 1
The dynamic brain functional development during the first two postnatal years is important for 2 establishing cognitive abilities and behaviors that could last a lifetime (1-3). As a prerequisite for 3 understanding how the brain works and develops, cortical parcellation maps provide a repository 4 that helps cortical area localization, network node definition, inter-subject comparison, inter-study 5 communication, and comparison, as well as reducing data complexity while improving statistical 6 sensitivity and power (4). In the functional aspect, researchers used to reveal and understand the 7 cortical network topography by clustering cortical vertices into parcels that are different from each 8 other in functional architecture using adult resting-state fMRI (rs-fMRI) data (5, 6). Although these 9 clustering-based methods can produce convincing results given a limited number of clusters, they 10 are not suitable for fine-grained parcellations (e.g., ≫ 100 parcels), as they usually result in 11 considerable disjointed fragments that are hardly explainable. To this point, recent adult 12 parcellations (7-10) started to use gradient-based methods, i.e., the functional gradient density, to 13 delineate sharp changes of resting-state functional connectivity (RSFC) patterns to promote the 14 meaningfulness and accuracy of parcel boundaries. 15 All the abovementioned studies derived functional parcellation maps using adult data, which 16 are not suitable for infant studies featuring dynamic brain structural and functional development, 17 due to enormous differences in brain functional organization between infants and adults (11,12). 18 Therefore, infant-specific cortical functional parcellation maps are highly desired, but remain 19 scarce, due to difficulties in both acquiring high-resolution infant brain multi-modal MR images 20 and challenges in processing infant MR images that typically have prominently dynamic imaging 21 appearance and extremely poor tissue contrast (1,2,13). Of note, another critical issue of using the 22 abovementioned methods for generating infant functional parcellations is that these methods 23 typically computed the functional gradient density map for a cohort directly based on cortical 24 folding-based registration and extensive spatial smoothing of functional connectivity, which, 25 4 however, cannot lead to accurate functional alignment across individuals, due to large variation 26 between folding and functional areas. Thus many vital details of the functional architecture are 27 blurred and inherently missed in the resulting functional parcellation maps. 28 In this paper, we aim to generate the first set of infant-specific, high-resolution, fine-grained 29 functional parcellation maps on the cortical surface to significantly accelerate early brain 30 development studies. To this end, this study leverages a large high-resolution dataset with 1,064 rs-31 fMRI scans and 394 T1-weighted and T2-weighted structural MRI scans from birth to 2 years of 32 age, as part of the UNC/UMN Baby Connectome Project (14). To ensure accuracy, all MR images 33 are processed using an extensively validated, advanced infant-dedicated cortical surface-based 34 pipeline (15). To establish accurate cortical functional alignment across individuals, we propose a 35 novel method to first compute the functional gradient density map of each infant scan, rather than 36 for the whole cohort in the traditional way, to capture fine-grained functional patterns, and then co-37 register all functional gradient density maps across individuals based on both cortical folding and 38 functional gradient information. Following steps detailed in Fig. 1, our derived group-average 39 functional gradient density maps capture much more details of cortical functional architecture than 40 the conventional method, thus enabling us to generate fine-grained age-specific cortical parcellation 41 maps of infants at multiple ages, i.e., 3,6,9,12,18, and 24 months of age. To facilitate infant 42 studies requiring parcel-to-parcel correspondences across ages, we also generate age-common 43 parcellation maps that are suitable for all ages during the first two years. Our infant functional 44 parcellation maps will soon be released to the public to greatly contribute to the pediatric 45 neuroimaging research community. 46 5 2. Results 47 Fig. 1. The procedure of infant parcellation using functional gradient density. Our major steps include structural and functional MRI processing, computing gradient density map for each scan, function-based cortical registration, and generating fine-grained functional parcellation.

48
We unprecedentedly investigated the fine-grained cortical surface-based functional parcellation 49 maps of the infant cerebral cortex using 1,064 high-resolution (2 × 2 × 2 mm 3 ) resting-state fMRI 50 scans from 197 healthy infants, with subject demographics shown in Table 1 and Fig. 8. To capture 51 detailed patterns of sharp transition between cortical areas, after the conventional cortical folding-52 based inter-individual cortical registration, the gradient density map of cortical functional 53 connectivity was computed on each scan of each individual and further used as a reliable functional 54 6 feature for function-based registration for establishing functionally more meaningful cortical 55 correspondences across individuals. This resulted in considerably detailed visualization of 56 functional boundaries on the cerebral cortex, which was used to create the infant-dedicated fine-57 grained cortical functional parcellation maps. Detailed steps of the proposed method are illustrated 58 in Fig. 1. 59

103
To evaluate the consistency of gradient density across different age groups, we thresholded and 104 binarized the age-specific functional gradient density maps to their top 50% and 25% gradient 105 density. These binary maps were summed up, resulting in a gradient density overlap map indicating 106 its age consistency shown in Fig. 4 (b). In these maps, "one" stands for high gradient densities that 107 appeared in only one age group, and "six" represents high gradient densities that appeared in all six 108 10 age groups. It is worth noting that most high gradient densities are repeatedly detected in all six age 109 groups, suggesting the high consistency of majorities of high gradient densities in all age groups. 110 Further, to better illustrate the functional architecture development, we computed the temporal 111 variability of gradient density maps between neighboring age groups, as shown in Fig. 4 (a). In 112 general, the temporal variabilities of functional gradient density are at a relatively low level (<=0.05) 113 at all age intervals. Across all ages, high temporal variabilities are mainly presented in high-order 114 (b) Consistency of high gradient density across ages. (c) Variabilities between each age-specific functional gradient density map and that of the age-common map. association areas, including the left middle and inferior frontal, middle temporal, right superior 115 frontal, precuneus medial prefrontal, and bilateral supramarginal, posterior superior temporal, and 116 medial frontal areas. Other regions mostly exhibit low temporal variabilities, especially in the 117 sensorimotor and medial occipital regions. Keeping this spatial distribution, the temporal variability 118 shows a multi-peak fluctuation, where the gradient density decreases from 3-6 to 6-9 months, 119 followed by an increase during 9-12 and 12-18 months, and drops again during 18-24 months. 120 121

Age-common Functional Gradient Density and Parcellation Maps 122
Since infant functional MRI studies typically involve multiple age groups, it is highly desired to 123 have an age-common functional parcellation map that features parcel-to-parcel correspondences 124 across ages, so that it can be conveniently employed for all ages during infancy. Therefore, we also 125 computed the age-common gradient density map ( Fig. 5 (a)) as the average of the functional 126 gradient density maps of all six age groups. The variabilities between the age-common gradient 127 density map and each age-specific gradient density map are illustrated in Fig. 4 (c). Compared to 128 the temporal variability between neighboring age groups ( Fig. 4 (a)), the age-common gradient 129 density map shows small variability to all age-specific maps. The spatial distributions of high and 130 low variabilities remain mostly similar to that of the temporal variabilities between neighboring 131 age groups, with high variability presented in some high-order association cortices and low 132 variability in unimodal cortices. Consequently, it can be speculated that the age-common gradient 133 density map can be used to generate an age-common parcellation map that is suitable for all subjects 134 from birth to 2 years of age. 135 136 12 The resulted age-common functional parcellation map based on the age-common gradient 137 density map is shown in Fig. 5 (b), which has 864 parcels in total (L: 432, R: 432) excluding the 138 medial wall. It should be noted that we manually removed some apparently over-segmented regions 139 13 after using the watershed algorithm, and prior to that, we had 903 parcels in total (L: 448, R: 455). 140 The parcel boundaries of the age-common parcellation map are well aligned with high gradient 141 density regions and show largely bilaterally symmetric patterns of the areal organization. In the 142 following development-related analyses in this study, we mainly employed the age-common 143 parcellation map to facilitate comparisons of infants across ages. 144 Compared to existing fine-grained parcellation maps, such as the multi-modal adult parcellation 145 (4), the age-common infant parcellation map has comparatively smaller and more evenly 146 distributed parcel sizes and shapes. Also, as shown in Fig. 5 (f), some areas of our parcels show 147 substantial overlap with the known cortical areas of adults, such as the visual areas V1, MT, MST, 148 sensorimotor areas 2, 3, 4, auditory areas A1, LBelt, and language areas 44, 45. To further examine 149 the validity of our parcellation map, we compared it with 1,000 null parcellation maps in terms of 150 variance and homogeneity, with the results shown in Fig. 5 (c) and (d). It can be observed that our 151 parcellation map shows significantly higher homogeneity (p=2e-10) and lower variance (p=4e-05), 152 indicating the meaningfulness of the resulting parcellation map. 153

Network Organization and Development 154
We performed network clustering of the generated parcels in each age group to reveal the early 155 development of functional network organization. The number of networks for each age group is 156 determined separately according to the random split-half stability analysis. Empirically, the 157 network number is set as 2 to 30, and the stability plots are shown in Fig. 6 (a). Higher stability 158 suggests a better clustering result, hence a more meaningful network organization. Overall, when 159 the network number surpasses 15, the stability does not show a substantial raise or decrease, 160 indicating the network numbers likely hold less than 15. Hence, we look for the cluster number on 161 a 'peak' or prior to a 'descending cliff', which guarantees high stability or more significant network 162 numbers. As a result, we find that the most suitable cluster numbers for different age groups are 7 163 networks for 3 months, 9 networks for 6 months, 10 networks for 9, 12, 18, and 24 months. Of note, 164 14 we choose10 networks for 18 months so as to be consistent during development, even though it is 165 neither a peak nor a cliff. 166 The spatiotemporal patterns of the discovered functional network organization are shown in Fig.  167 6 (b). Overall, changes in network structure from 3 to 9 months are more extensive than those from 168 9 to 24 months. Specifically, the sensorimotor network splits into two subnetworks from 3 to 6 169 months, and the boundary between them moved toward the ventral direction from 6 to 9 months. 170 The hand sensorimotor expands, while the mouth sensorimotor shrinks and both stabilize after 9 171 months. The auditory network is distinguished at 3 months and merges into the hand sensorimotor 172 at 6 months. The visual network splits into peripheral and central visual subnetworks from 6 to 9 173 months and remains stable until a slight shrinkage at 24 months. 174 15 Other networks exhibit more complex development with multi-peak fluctuation of the size in 175 certain networks. Specifically, the anterior default mode network expands from 3 to 6 months, and 176 shrinks from 6 to 9 months and from 12 to 18 months, and expands thereafter. The lateral posterior 177 default mode network that emerged at 6 months shrinks from 6 to 9 months and then expands from 178 9 to 18 months; while the medial posterior default mode network emerged at 9 months only lightly 179 shrinks thereafter. The anterior and posterior default mode networks develop to the adult-like 180 pattern at 18 months, while till 24 months, they are still detected as two separate networks. The 181 superior temporal network shrinks from 3 to 6 months, expands from 6 to 9 months, and then 182 shrinks again from 9 to 24 months. The anterior frontoparietal network gradually shrinks from 3 to 183 24 months, except for an expansion from 12 to 18 months. On the medial surface, the posterior 184 frontoparietal network expands to include the parahippocampal gyrus from 3 to 6 months and then 185 disappears by 9 months. On the lateral surface, the posterior frontoparietal network expands from 186 3 to 9 months to include the inferior temporal part and becomes stable thereafter. The dorsal 187 attention network is seen at 6 months and evolves to the adult-like pattern at 9 months and keeps 188 stable thereafter. 189

Parcel-wise Development 190
Homogeneity of functional connectivity can be used as a criterion for characterizing functional 191 development. Fig. 7 (a) shows the parcel-wise homogeneity development during infancy. Our 192 results suggest that the overall parcel-wise homogeneity shows a monotonic decrease trend during 193 the first two years by maintaining similar relative spatial distribution. Higher homogeneities are 194 located in the sensorimotor, paracentral, posterior insula, inferior parietal, posterior superior 195 temporal, lateral occipital, and occipital pole. Low homogeneities are presented in lateral prefrontal, 196 medial frontal, anterior insula, inferior temporal, and temporal pole.

Discussion 203
In this study, we created the first set of both age-specific and age-common, infant-dedicated, fine-204 grained, and cortical surface-based functional parcellation maps using functional gradient density 205 maps. We analyzed the spatiotemporal patterns of age-specific functional gradient density maps 206 and found that age-common functional gradient density maps are suitable for creating fine-grained 207 functional parcellation maps for all ages in the infant cohort. We validated the meaningfulness of 208 the parcellation and showed that its boundaries substantially reproduced known areal boundaries, 209 and its parcels featured high homogeneity and low variance. Finally, we illustrated the infantile 210 development in network structure, parcel homogeneity, and parcel local efficiency. 211 This study used the functional gradient density as a feature for improving functional alignment 212 across individuals, in addition to the conventional cortical folding features used by previous adult 213 functional parcellations (7,9). As a result, our method not only captured important coarse gradient 214 17 patterns discovered by previous methods, but also revealed much more detailed areal boundaries at 215 a remarkable resolution, as compared in Fig. 2. The main reason is that previous studies solely 216 relied on cortical folding-based registration, thus inevitably suffered from significant inter-subject 217 variability in the relation between cortical folding and functional areas, leading to less accurate 218 inter-subject functional correspondences. For infant-dedicated functional parcellations, the only 219 one available is the volumetric-based parcellation based on image registration and clustering 220 generated by Shi et al. (16), without any advanced surface-based processing and registration. In 221 contrast, our parcellation maps are generated based on the cortical surface, which well respects the 222 topology of the convoluted cerebral cortex, and avoids mixing signals from opposite sulcal banks 223 and different tissues, leading to more accurate functional signals resampling, smoothing, 224 computation and registration. Moreover, our parcellations leveraged high-quality 2 mm isotropic 225 fMRI data that densely covers the first two years, instead of data with a coarse resolution of 4 mm 226 isotropic centering at birth, 1 and 2 years of age (16). 227 Standing upon the detailed gradient density patterns by the proposed method, we generated age-228 specific fine-grained parcellation maps for 3, 6, 9, 12, 18, and 24 months of age. We found that the 229 temporal variability (temporal changes) of the functional gradient density generally decreases 230 during most age intervals, except for a slight increase from 9 to 12 months. This may suggest that 231 the development of functional architecture gradually slows down during the first two years. We 232 also found that high temporal variabilities mostly presented in high-order association cortices, 233 implying that they are developing at a more considerable pace compared to unimodal cortices. 234 Besides the age-specific parcellations, we also generated an age-common parcellation that suites 235 infants at all ages to help brain development-related studies under two considerations. First, infant 236 studies typically involve subjects of different ages, and it is not convenient to use different 237 parcellations for comparison between different ages. Second, our age-common gradient density 238 map shows low variability to all age-specific gradient density maps and therefore can generate the 239 18 representative fine-grained functional parcellation map that is suitable for all ages during infancy. 240 However, we will make both the age-specific and age-common functional parcellation maps 241 accessible to the public in the case that some researchers may still prefer age-specific parcellation 242 maps. 243 This study successfully augmented the resolution of the existing cortical parcellations from ~300 244 to ~900 areas, which represents a finer architecture of brain functional organizations compared to 245 previous ones. This fine-grained cortical organization is also in line with Eickhoff et al. (17), where 246 they believe that 200-300 areas are not the ultimate resolution for cortical parcellations due to the 247 multi-hierarchical formation of the brain. Glasser et al. (4) also consider 360 as a lower bound for 248 cortical parcellations since each parcel can be represented as a combination of several smaller 249 regions. Consequently, our fine-grained infant cortical parcellation maps provide a great platform 250 for analyzing pediatric neuroimaging data with a greatly boosted resolution, thus leading to more 251 meaningful discoveries on the fine-scaled functional architecture of infant brains. 252 Besides, it is worth noting that our parcellation increased the resolution in a meaningful way. 253 First, our functional gradient density maps are highly reproducible. By separating subjects into non-254 overlapping parts, their gradient density patterns are repeated with a dice ratio of ~0.93. Second, 255 our age-common infant parcellation shows high accordance in some specific cortical areas defined 256 by Glasser et al. (4), which is recognized as the state-of-the-art adult parcellation map. As illustrated 257 in Fig. 5 (d), our gradient density map-derived parcellation contains parcels that have substantial 258 overlap with the known adult area V1 defined by Glasser et al. (4). Other known cortical areas of 259 adults, such as sensorimotor areas 2, 3, and 4, were also overlapped with a combination of several 260 parcels in our parcellation. These observations that parcel borders conform to some adult cortical 261 areas lend substantial visual validity to the parcellation. 262 When applying parcellation as a tool to explore infant brain functional development, our results 263 reveal complex multi-peak fluctuations in several aspects, including parcel number, temporal 264 19 variation of gradient density, network organization, and local efficiency. To the best of our 265 knowledge, this complex fluctuation development trend is not reported in previous literature and 266 should fill an important knowledge gap for infantile brain functional development. These functional 267 developmental patterns are very different from early brain structure development, where the 268 cortical thickness follows an inverted-"U" shaped trajectory, while the surface area and cortical 269 volume monotonically increase following a logistic curve. The multi-peak fluctuations potentially 270 mirror different milestones of behavioral/cognitive abilities, which likely emerge at different ages 271 during infancy (30). However, the underlying mechanisms of such developmental patterns remain 272 to be further investigated. 273 For network organization ( Fig. 6 (b)), at 3 months, networks likely groups vertices with close 274 spatial locations, resulting in networks being more dependent on the local anatomy. After 9 months, 275 the primary functional systems reach steady and present adult-like patterns, while high-order 276 functional networks still show substantial differences compared to the adult-like pattern. Our results 277 suggest that a primitive form of brain functional networks is present at 3 months, which is largely 278 consistent with recent studies suggesting that most resting-state networks are already in place at 279 term birth (18-20). Besides, our results also suggest that, compared to high-order functional 280 networks, the primary functional system is more developed in infants. This confirms previous 281 findings in infant cortical thickness development (21), suggesting that the primary functional 282 systems develop earlier than high-order systems. 283 At the network level, the sensorimotor system splits into two sub-networks, i.e., the mouth-and 284 hand-sensorimotor at 6 months, which were also observed in infants and toddlers (22). The visual 285 network is split into central (primary) visual and peripheral (high-order) visual cortices at 9 months 286 and well maintains this pattern until 24 months. This subdivision of mouth-and hand-sensorimotor 287 networks is also found in adults (5). The high-order functional systems, including the default mode, 288 frontoparietal, and dorsal attention network, exhibit considerable development during 3 to 9 months, 289 20 followed by some minor adjustments from 12 to 24 months. A previous infant study (19) also  290 demonstrated that functional network development shows more considerable change in the first 291 year compared to the second year. At 24 months, both default mode and frontoparietal networks 292 show a lack of strong cross-lobe connections. Though several studies identified some prototypes of 293 cross-lobe connection (19, 23), their links seem not as strong as to be stably distinguished (24-26). 294 Our results suggest that the high-order functional networks are far more from established at 24 295 months of age. It is worth noting that, the size changes of networks can be quite subtle between a 296 short time interval, which emphasizes the importance of using a fine-grained parcellation map. 297 The parcel homogeneity measures the development within parcels. Our result shows ( Fig. 7 (a)) 298 that unimodal cortices, including the sensorimotor, auditory, and visual areas, show high 299 homogeneity, which is largely consistent with adults (7). However, the inferior parietal and 300 posterior superior temporal cortices, which show high homogeneity in infants, are observed low 301 homogeneity in adults (7). Besides, the prefrontal area, which shows relatively low homogeneity 302 in infants, seems to develop to a medium-to-high homogeneity in adults. Almost all parcels are 303 observed decreased homogeneity with age. This is likely related to the development of brain 304 function, especially in high-order cortices, which show increased heterogeneity, and consequently 305 decreased homogeneity. Among the high-order association cortices, the prefrontal area has the 306 lowest homogeneity, followed by temporal and then parietal regions, suggesting different levels of 307 functional development. 308 Local efficiency measures a different developmental aspect -it represents the connection of 309 parcels to neighbors. Higher local efficiency is usually related to higher functional segregation (11). 310 Our results (Fig. 7 (b)) suggest that local efficiency shares certain similar spatial distribution with 311 homogeneity -they both increase in anterior to posterior and ventral to dorsal directions. During 312 development, the local efficiency shows a complex developmental trend: although 24 months 313 shows a strong increase compared to 1 month, there is a dip from 12 to 21 months that should be 314 21 noted. The age-related increase of local efficiency was previously found from 18 months to 18 315 years (27), 5 to 18 years (28), and 12 to 30 years (29), and is likely explained by progressive white 316 matter maturation (27). This trajectory of local efficiency is not contradictory to the previous 317 studies (31, 32), since they only measured the averaged local efficiency of all nodes to reflect 318 network characterization, thus missing important characteristics of parcel-level local efficiency. 319 This further stresses the importance of performing parcel-wise analyses and the significance of 320 fine-grained infant cortical parcellations. 321

Conclusion 322
In summary, for the first time, this study constructed a comprehensive set of cortical surface-based 323 infant-dedicated fine-grained functional parcellation maps. To this end, we developed a novel 324 method for establishing functionally more accurate inter-subject cortical correspondences. We 325 delineated age-specific parcellation maps at 3, 6, 9, 12, 18, and 24 months of age as well as an age-326 common parcellation map to facilitate studies involving infants at different ages. Our parcellation 327 maps were demonstrated meaningful by comparing with known areal boundaries and through 328 quantitative evaluation of homogeneity and variance of functional connectivity. Leveraging our 329 infant parcellation, we provide the first comprehensive visualizations of the infant brain functional 330 developmental maps on the cortex and reveal a complex multi-peak fluctuation functional 331 development trend, which will serves as valuable references for future early brain developmental 332 studies. Our generated fine-grained infant cortical functional parcellation maps will be released to 333 the public soon to greatly advance pediatric neuroimaging studies. 334

Fig. 8. Longitudinal distribution of scans. Each point represents a scan with its scanned age (in months)
shown in the x-axis, with males in blue and females in red, and each horizontal line represents one subject, with males in blue and females in red.

337
Subjects in this study are from the UNC/UMN Baby Connectome Project (BCP) data90set (14). 338 The BCP focuses on normal early brain development, where all infants were born at the gestational 339 age of 37-42 weeks and free of any major pregnancy and delivery complications. In this study, 394 340 high-resolution longitudinal structural MRI scans were acquired from 197 (90 males and 107 341 females) typically developing infants, as demonstrated in Fig. 8. Images were acquired on 3T 342 Siemens Prisma MRI scanners using a 32-channel head coil during natural sleeping. T1-weighted 343 images (208 sagittal slices) were obtained by using the three-dimensional magnetization-prepared

Resting-State fMRI Processing 368
Infant rs-fMRI processing was conducted according to an infant-specific functional pipeline (31, 369 45, 46). The head motion was corrected using FSL, as well as the spatial distortions due to gradient 370 24 non-linearity. The rs-fMRI scans were then registered to the T1-weighted structural MRI of the 371 same subject using a boundary-based registration approach (47). All of the transformations and 372 deformation fields were combined and used to resample the rs-fMRI data in the native space 373 through a one-time resampling strategy. After conservative high-pass filtering with a sigma of 374 1,000 s to remove linear trends in the data, individual independent component analysis was 375 conducted to decompose each of the preprocessed rs-fMRI data into 150 components using 376 MELODIC in FSL. An automatic deep learning-based noise-related component identification 377 algorithm was used to identify and remove non-signal components to clean the rs-fMRI data (48). 378

Cortical Surface Reconstruction and Mapping 379
Based on the tissue segmentation results, inner, middle and outer cortical surfaces of each 380 hemisphere of each MRI scan were reconstructed and represented by triangular meshes with correct 381 topology and accurate geometry, by using a topology-preserving deformable surface method (33, 382 49). Before cortical surface reconstruction, topology correction on the whiter matter surface was 383 performed to ensure the spherical topology of each surface (50). After surface reconstruction, the 384 inner cortical surface, which has vertex-to-vertex correspondences with the middle and outer 385 cortical surfaces, was further smoothed, inflated, and mapped onto a standard sphere (51). 386 To ensure the accuracy in longitudinal analysis during infancy, it is necessary to perform 387 longitudinally-consistent cortical surface registration (15). Specifically, 1) for each subject, we first 388 co-registered the longitudinal cortical surfaces using Spherical Demons (52) based on cortical 389 folding-based features, i.e., average convexity and mean curvature. 2) Longitudinal cortical 390 attribute maps were then averaged to obtain the intra-subject mean surface maps. 3) For each 391 hemisphere, all intra-subject mean surface maps were then co-registered and averaged to get the 392 population-mean surface maps. 4) The population-mean surface maps were aligned to the HCP 393 32k_LR space through registration to the "fsaverage" space as in (53). By concatenating the three 394 deformation fields of steps 1, 3 and 4, we directly warped all cortical surfaces from individual scan 395 25 spaces to the HCP 32k_LR space. These surfaces were further resampled as surface meshes with 396 32,492 vertices, thus establishing vertex-to-vertex correspondences across individuals and ages. 397 All results were visually inspected to ensure sufficient quality for subsequent analysis. The inner 398 and outer cortical surfaces were used as a constraint to resample the rs-fMRI time courses onto the 399 middle cortical surface with 32,492 vertices using the HCP workbench (54), and the time courses 400 were further spatially smoothed on the middle cortical surface with a small Gaussian kernel ( = 401

Generation of Fine-grained Cortical Functional Parcellation Maps 403
In this section, we describe detailed steps for generating fine-scaled infant cortical functional 404 parcellation maps (see Fig. 1). Specifically, we first describe the computation of the gradient 405 density map of functional connectivity for each scan, followed by a function-based registration step 406 based on gradient density maps. Then, we detail the computation of both "age-specific" and "age-407 common" parcellation maps based on the functional registration results and our evaluation scheme. 408 At last, we introduce how we use the parcellation maps to discover the functional network 409 organization development, as well as parcel homogeneity and local efficiency during infancy. 410

Computation of Individual Functional Gradient Density Map 411
The gradient density of functional connectivity (7) identifies sharp changes of RSFC, thus 412 intrinsically representing the transition from one functional parcel to another, and is widely used in 413 generating meaningful fMRI-based cortical parcellations in adult studies (7-9). For each fMRI scan 414 of each infant subject, the computation of the gradient density of functional connectivity on the 415 cortical surface is summarized in the following steps. 1) For each fMRI scan, the functional 416 connectivity matrix is built by pair-wise correlating each vertex with all other cortical vertices in 417 the CIFTI file to create a 32k×64k RSFC matrix for each hemisphere. 2) Each RSFC matrix is 418 transformed to z scores using Fisher's r-to-z transformation. 3) For each fMRI scan, the z-419 transformed RSFC of each vertex is correlated with all cortical vertices within the same hemisphere, 420 26 creating a 2 nd order correlation matrix (RSFC-2 nd ) sized 32k×32k for each hemisphere. 4) For some 421 scan visits consisting of both AP and PA scans, all RSFC-2 nd matrices of the same visit from the 422 same subject are averaged, so that all subjects contribute equally, even they may have different 423 numbers of scans in one visit. 5) The gradient of functional connectivity is computed on the RSFC-424 2 nd as in (4), resulting in a 32k×32k gradient matrix per hemisphere. 6) By performing the 425 watershed-based boundary detection (7) on the gradient matrix, we obtain 32k binary boundary 426 maps per hemisphere. 7) The functional gradient density map is defined as the average of 32k 427 binary boundary maps. 428

Cortical Surface Registration based on Functional Gradient Density 429
Previous studies mostly computed population-based functional gradient density map, where 430 cortical surfaces were usually co-registered to a common space using only cortical folding-based 431 features. However, due to the highly-variable relationship between cortical folds and functions, 432 especially in high-order association regions, researchers recently are getting more aware of the 433 necessity of functional features-based registration (55,56). To this end, in addition to cortical 434 folding-based co-registration, we further use the gradient density of functional connectivity as a 435 meaningful functional feature to perform a second-round of co-registration of cortical surfaces for 436 the purpose of more accurate functional alignment. 437 Specifically, based on cortical folding-based surface co-registration, 1) the functional gradient 438 density maps of all scans are averaged to generate the population-mean functional gradient density 439 map. 2) To improve inter-individual cortical functional correspondences, the functional gradient 440 density map of each scan is then aligned onto the current population-mean functional gradient 441 density map using Spherical Demons (52) by incorporating functional gradient density as a feature. 442 3) All warped functional gradient maps are then resampled and averaged to obtain the newly 443 improved population-mean functional gradient density map with sharper and more detailed 444 functional architecture. 4) Steps 2 and 3 are repeated iteratively until no visually observed changes 445 27 in the population-mean functional gradient density map (4 iterations in our experiment). After this 446 procedure, all individual functional gradient density maps are co-registrated, thus establishing 447 functionally more accurate cortical correspondences across individuals. 448

Generation of Parcellation Maps based on Functional Gradient Density 449
Age-specific Parcellation Maps: To capture the spatiotemporal changes of fine-grained cortical 450 functional maps during infancy, we group all scans into 6 representative age groups, i.e., 3,6,9,451 12, 18, and 24 months of age based on the distribution of scan ages. For each age group, we compute 452 the age-specific group-average functional gradient density maps by averaging the gradient density 453 maps of all scans within the group, without any smoothing. Detailed information of each age group 454 is reported in Table 1. A watershed method is then applied on each age-specific group-average 455 functional gradient density map to generate the corresponding functional parcellation maps (7). 456 This watershed segmentation algorithm starts by detecting local minima in 3-ring neighborhoods, 457 and iteratively grows the region until reaching ambiguous locations, where vertices can be assigned 458 to multiple regions. These locations appear to be borders that separate parcels and reflect putative 459 boundaries of functional connectivity according to the functional gradient density maps. 460 Age-common Parcellation Maps: Ideally, the age-specific parcellation maps are the more 461 appropriate representation of the cortical functional architecture at the concerned age. However, 462 many neuroimaging studies involve infants across multiple ages, thus the age-specific parcellation 463 maps may not be proper choices due to different parcel numbers and variation in parcel boundaries 464 across ages, thus inducing difficulties in across-age comparisons. To facilitate infant studies 465 involving multiple age groups, we also compute an age-common gradient density map, which is 466 the average of all 6 age-specific functional gradient density maps without any smoothing, so that 467 each age group contributes equally to the age-common map. According to the age-common 468 functional gradient density map, we generate the age-common functional parcellation map using 469 the watershed segmentation method as well. The subsequent parcellation evaluation, functional 470 28 network architecture and longitudinal development analyses are performed using the age-common 471 parcellation maps. 472

Evaluation of Parcellation Maps 475
Reproducibility: Ideally, a functional gradient density-based parcellation map should extract 476 robust common gradient information that shows the transition between parcels. We thus test if the 477 gradient density map is reproducible on different subjects. Therefore, randomly divided "generating" 478 and "repeating" groups (7, 9) are used to calculate mean gradient density map, separately. These 479 two maps are then binarized by keeping only 25% highest gradient density as in (7,8), and the dice 480 ratio overlapping index between the two binarized maps is calculated to evaluate the reproducibility 481 of the functional gradient map. This process is repeated multiple times (1,000 times in this study) 482 to get a reliable estimation. 483 Homogeneity: The functional gradient density-based parcellation identifies large gradients, 484 representing sharp transition in functional connectivity pattern and avoiding large gradients inside 485 parcels as much as possible. Meanwhile, a parcel that accurately represents a cortical area should 486 not only be distinct from its neighbors in functional connectivity pattern, but also has a homogenous 487 29 functional connectivity pattern across all vertices inside. Therefore, we estimate the homogeneity 488 of each parcel as in (7). Specifically, we first compute the mean correlation profile of each vertex 489 across all subjects. Next, the correlation patterns of all vertices within one parcel is entered into a 490 principal component analysis; the percentage of the variance that is explained by the largest 491 principal component is used to represent the homogeneity of this parcel. 492 Variance: As the functional connectivity pattern within a parcel should be relatively uniform, we 493 also measure the variability of the connectivity pattern within each parcel, with smaller variability 494 indicating greater uniformity and hence higher parcellation quality. Specifically, for each parcel, 495 we first obtain a matrix with each column representing subject-average z-score of functional 496 connectivity profile of one vertex in the parcel. Then we compute the sum of standard deviation of 497 each row to represent the variability of this parcel. The average variability of all parcels is used to 498 represent the variability of the parcellation map. 499 As parcellation maps usually have different numbers, sizes and shapes in parcels, to have fair 500 comparison and be consistent with (7, 8), we compare our parcellation maps with 'null 501 parcellations'. The null parcellations are generated by rotating by a random amount along x, y and 502 z axes on the 32,492 spherical surfaces, which relocate each parcel while keeping the same number 503 and size of parcels. We compare both variability and homogeneity of our parcellation and that of 504 the random rotated null parcellations. Notably, in any random rotation, some parcels will inevitably 505 be rotated into the medial wall, where no functional data exist. The homogeneity/variance of a 506 parcel rotated into the medial wall is not calculated; instead, we assign this parcel the average 507 homogeneity/variance of all random versions of the parcel that were rotated into non-medial-wall 508 cortical regions. 509 Variability Between Functional Gradient Density Maps: A variability map visualizes the 510 variability or dissimilarity between two functional gradient density maps, and is estimated as 511 follows. For a vertex , a surface patch centering at is extracted (10-ring neighborhood in this 512 30 study), and two vectors 1 and 2 within this patch are then extracted from two functional gradient 513 density maps. Their variability at is computed as 0.5 × (1 − ( 1 , 2 )), where (⋅,⋅) 514 stands for Pearson's correlation. As a result, the variability/dissimilarity is within the range of [0, 515 1], where high value stands for high variability/dissimilarity and vice versa. In this study, we mainly 516 measure the variability between functional gradient density maps in two aspects: 1) the temporal 517 variability, which computes the variability of functional gradient density maps between two 518 consecutive age groups to reflect the developmental changes of the gradient density maps; 2) the 519 variability between the age-common functional gradient density map and each age-specific 520 functional gradient density map, for quantitatively evaluating whether it is appropriate to use the 521 age-common parcellation maps for all 6 age groups. 522

Functional Development Analysis 523
Functional Network Detection: To discover the developmental evolution of large-scale cortical 524 functional networks, we employ a network discovery method (5) to each of the 6 age groups. 525 Specifically, for each subject in each age group, given parcels, we first compute the average time 526 course of each parcel (excluding the medial wall), and compute the correlation of the average time 527 courses between any two parcels. This results in a × matrix, which is further binarized by 528 setting the top 10% of the correlations to one and the rest to zero. For each age group, all × 529 matrices are averaged across individuals independently. A clustering algorithm (57) is then applied 530 to estimate networks of parcels with similar connectivity profiles. 531 To determine the optimal cluster number for each age group, we employ the random split-half 532 test to compute the stability for each , with higher stability corresponding to more meaningful 533 clustering results. Specifically, for each age group, we randomly split all subjects into two folds 534 and run the clustering algorithm separately to obtain two independent clustering results 1 and 2 , 535 and the similarity between 1 and 2 is evaluated using the Amari-type distance (58). This 536 experiment is repeated 200 times for each age group, and the resulted similarities are averaged to 537 31 represent the stability for . During this process, the range of is set to [2,30] according to the 538 existing literature of functional network discovery (5, 6). 539

Parcel-wise Development:
We computed the homogeneity and local efficiency of each parcel in 540 the age-common parcellation to characterize infantile parcel-wise developmental patterns regarding 541 functional homogeneity and functional segregation, respectively. The homogeneity is computed as 542 described in Section 3.4 for each subject, where higher parcel homogeneity indicates more unified 543 connectivity pattern within the parcel. The local efficiency is computed using the GRETNA Toolkit 544 (59) for each subject. Herein, multiple thresholds are used, keeping 50% to 5% connections with 545 1% as a step, and the area under curve (AUC) is calculated to represent the local efficiency to avoid 546 the influence of connectivity densities. The local efficiency corresponds to the mean information 547 transfer efficiency between a particular parcel and all its connected nodes, which is proportional to 548 the clustering coefficient. Parcels with higher local efficiency can more effectively share 549 information to its connected parcels, and thus help build effective segregated networks. To have 550 intuitive and spatiotemporally detailed views of their development, we use the sliding window 551 technique to compute homogeneity and local efficiency in each age window by averaging all scans 552 within the same age window. The windows are centered at each month, with a window width of 90 553 days (±45 days) at 2 months of age, increasing 4 days in width for each following month and 554 reaching 182 days (±91 days) at 2 years of age. 555

Competing Interests 560
The authors declare no competing interests.