Individual identification and individual variability analysis based on cortical folding features in developing infant singletons and twins.

Studying the early dynamic development of cortical folding with remarkable individual variability is critical for understanding normal early brain development and related neurodevelopmental disorders. This study focuses on the fingerprinting capability and the individual variability of cortical folding during early brain development. Specifically, we aim to explore (a) whether the developing neonatal cortical folding is unique enough to be considered as a "fingerprint" that can reliably identify an individual within a cohort of infants; (b) which cortical regions manifest more individual variability and thus contribute more for infant identification; (c) whether the infant twins can be distinguished by cortical folding. Hence, for the first time, we conduct infant individual identification and individual variability analysis involving twins based on the developing cortical folding features (mean curvature, average convexity, and sulcal depth) in 472 neonates with 1,141 longitudinal MRI scans. Experimental results show that the infant individual identification achieves 100% accuracy when using the neonatal cortical folding features to predict the identities of 1- and 2-year-olds. Besides, we observe high identification capability in the high-order association cortices (i.e., prefrontal, lateral temporal, and inferior parietal regions) and two unimodal cortices (i.e., precentral gyrus and lateral occipital cortex), which largely overlap with the regions encoding remarkable individual variability in cortical folding during the first 2 years. For twins study, we show that even for monozygotic twins with identical genes and similar developmental environments, their cortical folding features are unique enough for accurate individual identification; and in some high-order association cortices, the differences between monozygotic twin pairs are significantly lower than those between dizygotic twins. This study thus provides important insights into individual identification and individual variability based on cortical folding during infancy.

INTRODUCTION

The human cerebral cortex is a highly convoluted and complex structure, which presents remarkable individual variability in both cortical folding pattern and function (Duan, Xia, et al., 2017; Duan et al., 2018; Finn et al., 2017; MacDonald, Nyberg, & Bäckman, 2006; Mangin et al., 2004). Although it is still not well understood, several studies found that the individual variability of the cerebral cortex might be able to link human behavior and cognition to brain anatomy (Frost & Goebel, 2012; Kanai & Rees, 2011; J. Liu, Liao, Xia, & He, 2018; MacDonald et al., 2006; Mueller et al., 2013; Nordahl et al., 2007; Tavor et al., 2016), and the cerebral cortex is altered in structure and morphology in various neurodevelopmental disorders, for example, autism, cerebral palsy, attention deficit hyperactivity disorder and schizophrenia (Bigler et al., 2007; Keshavan et al., 1998; Lee et al., 2007; Takahashi et al., 2009). Thus, studying the patterns and characteristics of early brain development, for example, the developmental individual variability pattern and fingerprinting capability of infant cortical folding, would be of great importance for understanding many neurodevelopmental disorders. Current evidences indicate that the individual variability is nonuniform across the cortex, with high‐order association regions being more variable than the unimodal cortex in both cortical folding pattern and functional connectivity (Hill et al., 2010; Finn et al., 2015; Dubois & Adolphs, 2016). Accordingly, there has been a growing interest in the fingerprinting capability of the human brain, that is, whether the features of the human brain are unique and distinguishable for individual identification. Several studies performed individual identification based on structural or functional characteristics of the human brain, for example, morphological shape description of cortical and subcortical structures, named BrainPrint (Wachinger et al., 2015; Wachinger, Golland, & Reuter, 2014), shape description of white matter fibers, named FiberPrint (Kumar, Desrosiers, Siddiqi, Colliot, & Toews, 2017; Kumar, Toews, Chauvin, Colliot, & Desrosiers, 2018), and functional connectivity networks, named functional connectome fingerprint (Biazoli Jr et al., 2017; Finn et al., 2015; Horien, Shen, Scheinost, & Constable, 2019; Kaufmann et al., 2017; Liu et al., 2018; Tavor et al., 2016).

Of note, most of the above studies focus on adult individual identification, in which the brain structure and function have been relatively stable in the duration across different scans. Only a few studies based on functional connectome involve the developing brain from adolescent datasets (Horien et al., 2019; Kaufmann et al., 2017). To date, no study has investigated the fingerprinting property of cortical folding of the dynamic developing infant brain. In fact, studies on infant identification based on cortical folding are of great neuroscientific significance as they will allow the examination of: (a) whether the neonatal cortical folding, which undergoes dynamic development during the first postnatal years, is reliable and unique enough to be considered as a “fingerprint” of an infant; (b) which cortical regions particularly manifest more individual variability of cortical folding during the early brain development and are more capable of infant identification; (c) whether infant twins can be correctly distinguished based solely on the cortical folding. Addressing these questions remains a potential challenge since a plethora of evidences suggest that the first 2 years is the most dynamic phase of the postnatal brain development (Gilmore et al., 2007), during which both cortical volume and surface area expand rapidly, and new tertiary cortical folding emerges (Duan, Rekik, et al., 2017; Gilmore et al., 2011; Li, Nie, Wang, et al., 2013a; Li, Nie, Wang, Shi, Lyall et al., 2013b; Nie et al., 2011) as shown in Figure 1. Hence, individual identification involving neonatal brains should be more challenging than adult individual identification. Furthermore, we note that the large‐scale longitudinal infant datasets, which can be used for infant identification, are scarce and invaluable, due to the challenges in recruitment of eligible pregnant mothers and the difficulties of collecting quality motion‐free infant MR images and their follow‐up scans. Besides, it is also difficult to properly process the infant MR images, which typically exhibit dynamic appearance and extremely low tissue contrast (Li et al., 2019).

Figure 1

Longitudinal infant cortical surfaces at 0, 1, and 2 years of age of a pair of monozygotic twins. From (a) to (c), the cortical surfaces are color‐coded by mean curvature, average convexity, and sulcal depth, respectively

This study presents the first attempt to leverage the morphological features of cortical folding for infant identification and also explore the cortical regions charactering more individual variability for identification. Our particular choice of investigating the fingerprinting potential of cortical folding is first motivated by the fact that the major cortical folding patterns and their individual variability in the adult brain are already present at birth and largely preserved during postnatal brain development (Duan et al., 2018; Gilmore, Knickmeyer, & Gao, 2018; Hill et al., 2010; Rekik, Li, Lin, & Shen, 2016). Second, the morphology of cortical folding varies remarkably across individuals, since its formation process is influenced by many complex genetic and environmental factors. Even for the monozygotic (MZ) twins with identical genes and similar developmental environments, their cortical folding morphologies are different as shown in Figure 1. Thus, we hypothesize that the cortical folding is unique and reliable enough for infant identification and cortical regions exhibiting large individual variability have high fingerprinting potential.

In this study, we perform the infant individual identification and individual variability study in a relatively large‐scale longitudinal dataset, including 1,141 longitudinal MRI scans from 472 infants scanned at birth, 1 and 2 years of age. In particular, first, we propose an intuitive ROI‐based voting method based on the combination of multiple cortical folding features, that is, mean curvature, average convexity, and sulcal depth, to identify infant individuals. We apply the proposed method in the whole dataset to identify the 1‐year‐olds and 2‐year‐olds based on their neonatal cortical folding features. Second, we establish the maps of the regional identification capability and also the longitudinal maps of individual variability of cortical folding, thus inspecting the relationship between the regions that contribute most to the individual identification and the regions that manifest high individual variability. Third, we perform individual identification and analyze the individual variability of infant cortical folding for MZ and dizygotic (DZ) twins.

MATERIALS AND METHODS

Subjects and MR image acquisition

The Institutional Review Board of the University of North Carolina (UNC) School of Medicine approved this study. Pregnant mothers during their second trimester of pregnancy were recruited in the UNC Hospitals. The informed consent forms were signed by both parents of each infant. The subjects with abnormalities on fetal ultrasound, or major medical diseases or psychiatric illness on mother were excluded. All the infants in this study were free of congenital anomalies, metabolic disease, and focal lesions. All the infants were fed, swaddled, and fitted with ear protection before MRI scanning. And none of them was sedated for MRI scanning. More information on this cohort can be found in (Gilmore et al., 2011).

In this study, 1,141 longitudinal MRI scans were acquired from 472 healthy neonates (Year 0) with their incomplete follow‐up scans in the first two postnatal years. Specifically, 387 neonates were longitudinally scanned at 1 year of age (Year 1 for short), 282 neonates were longitudinally scanned at 2 years of age (Year 2 for short), and 197 neonates were longitudinally scanned at both 1 and 2 years of age. The mean gestational age at birth of the 472 neonates is 37.2 ± 2.8 weeks, and the percentages of males (45.3%) and females (54.7%) are relatively comparable. This infant cohort includes 40 pairs of MZ twins and 59 pairs of DZ twins. More demographic information of the infant cohort is shown in Table 1.

Table 1

Demographic information of the infant cohort

	Year 0	Year 1	Year 2	Total scans
Subject	472	387	282	1,141
Singleton	274	225	166	665
Twin	198	162	116	476
MZ twin	80	56	54	190
DZ twin	118	106	62	286
Age at MRI scan (months)	1.1 ± 0.6	13.3 ± 1.7	24.8 ± 1.0	—

Abbreviations: DZ: dizygotic; MZ: monozygotic.

T1‐weighted and T2‐weighted MR images were acquired on a Siemens head‐only 3T scanner (Allegra, Siemens Medical System, Erlangen, Germany) with a circularly polarized head coil. For T1‐weighted images, 160 sagittal slices were acquired with a 3D magnetization‐prepared rapid gradient echo (MPRAGE) sequence: TR = 1900 ms, TE = 4.38 ms, inversion time = 1,100 ms, flip angle = 7°, and resolution = 1 × 1 × 1 mm3. For T2‐weighted images, 70 transverse slices were acquired by using a turbo spin‐echo (TSE) sequences: TR = 7,380 ms, TE = 119 ms, flip angle = 150°, and resolution = 1.25 × 1.25 × 1.95 mm3 (Gilmore et al., 2011).

Image processing and surface mapping

Image preprocessing

All the MR images were processed using an infant dedicated computational pipeline (Li et al., 2015), which has been extensively validated on more than 2000 infant brain MRI scans. Specifically, the pipeline using the following main procedures for image preprocessing: (a) stripping the brain skull through a learning‐based method (Shi et al., 2012); (b) removing the cerebellum and brain stem using a registration method (Shen & Davatzikos, 2002); (c) correcting the intensity inhomogeneity using N3 algorithm (Sled, Zijdenbos, & Evans, 1998); (d) rigidly aligning all images at each age to the age‐matched infant brain atlas (Shi et al., 2011); (e) segmenting the brain MR images into gray matter (GM), white matter (WM) and cerebrospinal fluid (CSF) with an infant‐specific patch‐driven level‐set‐based method (Wang et al., 2013; Wang et al., 2014); (f) masking and filling the noncortical structures (i.e., lateral ventricles and subcortical structures), and then further separating each brain into left and right hemispheres.

Cortical surface reconstruction

Based on the tissue segmentation results, topologically and geometrically correct cortical surfaces for each hemisphere of each scan were reconstructed by a deformable surface method (Li, Nie, Wang, Shi, Gilmore, et al., 2014; Li, Nie, Wu, et al., 2012). Specifically, first, the topological defects in the WM were corrected to ensure a spherical topology for each hemisphere (Hao, Li, Wang, Meng, & Shen, 2016; Sun et al., 2019). Then, the corrected WM was tessellated as a triangular mesh to form the initial inner cortical surface. Next, the inner cortical surface was further deformed by preserving its initial topology to refine itself and reconstruct the outer cortical surface, using forces derived from Laplace's equation (Li, Nie, Wang, Shi, Gilmore, et al., 2014; Li, Nie, Wu, et al., 2012). The inner cortical surface of each hemisphere was further smoothed, inflated and mapped onto a standard sphere (Fischl, Sereno, & Dale, 1999).

Cortical surface registration and parcellation

For the longitudinal study of cortical folding, we need to establish both longitudinal cortical correspondences across different ages and cross‐sectional cortical correspondences across different individuals. To this end, for each hemisphere of each scan, we registered its spherical cortical surface onto the age‐matched instance in the 4D infant cortical surface atlases (Li et al., 2015; Wu et al., 2017) using Spherical Demons (Yeo et al., 2010). The inter‐subject and intra‐subject vertex‐to‐vertex correspondences of cortical surfaces were thus established across all the scans at three ages in our infant dataset. Consequently, each cortical surface was accordingly resampled to a standard mesh tessellation with 40,962 vertices (subdivision level 6 of an icosahedron). Meanwhile, the Desikan‐Killiany parcellation (Desikan et al., 2006) maps with 34 regions of interest (ROIs) in the 4D surface atlases were propagated to each cortical surface. Of note, the Desikan‐Killiany parcellation is generally adopted for adults, not dedicated for infants. However, all the primary and secondary cortical folds of the adult brain have already been established at term birth (Chi, Dooling, & Gilles, 1977; Duan et al., 2018; Li et al., 2016; Li, Wang, Shi, Lin, & Shen, 2014b), thus it is also reliable to properly propagate the Desikan‐Killiany parcellation based on major cortical folds to infant cortical surface.

Cortical folding features

To comprehensively characterize the complex morphology of the cortical folding, we chose three complementary morphological features of the inner cortical surface, that is, mean curvature, average convexity, and sulcal depth. The reason of choosing these features is that they are widely and jointly used in surface registration and morphometry, and are reliable and informative for characterizing the inter‐subject and intra‐subject variability in cortical folding. Specifically, the mean curvature is an intuitive indication of cortical gyri and sulci and reflects the finer‐scale local geometric information of cortical folding, including primary, secondary, and tertiary folds (Cachia et al., 2003; Dubois et al., 2007; Habas et al., 2011; Li, Nie, Wang, Shi, Lyall, et al., 2013b; Nie et al., 2010). The average convexity indicates the integrated normal movement of a vertex during cortical surface inflation and mainly reflects the larger‐scale morphological information of primary cortical folds (Fischl et al., 1999). As for the sulcal depth, which was computed as the distance from each vertex on the inner cortical surface to the closest point on the cerebral hull surface, encodes the continuously varying morphology of both coarse and fine cortical folding (Li et al., 2015; Li, Nie, Wang, Shi, Lyall, et al., 2013b; Meng, Li, Lin, Gilmore, & Shen, 2014; Van Essen, 2005). Thus, these features provide comprehensive and complementary information about cortical folding. Other features, for example, cortical thickness, which is very sensitive to noise, and gyrification, which is a relative global description of cortical folding, are unsuitable for such individual identification and individual variability study. Figure 1 shows the longitudinal infant cortical surfaces at 0, 1, and 2 years of age of a pair of MZ twins. From (a) to (c), the cortical surfaces are color‐coded by mean curvature, average convexity, and sulcal depth, respectively. Herein, the positive and negative mean curvatures/average convexity on the cortical surface correspond to sulcal bottoms and gyral crests, respectively. Note that the sulcal depth increases remarkably with the brain development during this stage.

Folding‐based infant identification

To investigate whether cortical folding morphology is able to fingerprint each infant, we proposed to identify infant individuals at 1 or 2 years of age using their neonatal cortices based on cortical folding features, that is, mean curvature, average convexity, and sulcal depth. Specifically, we proposed an intuitive ROI‐based voting framework for infant identification. Compared to the global‐based method, this strategy would largely improve the identification accuracy, since the negative effect of the regions with less identification capability would be neglected and suppressed during the voting process, while the positive effect of the regions with higher identification capability would be preserved and strengthened. Herein, the global‐based method denotes that the identification was conducted by comparing the similarities of whole‐brain cortical feature maps between the to‐be‐identified scan and all other scans in the database, and the identity of the most similar scan in the database would be the final identification result. Besides, the feature combinations in the proposed framework would outperform the single feature in infant identification, since multi‐view features jointly and comprehensively characterize the cortical folding. In this study, for a new coming scan at a later age, we performed the individual identification by finding its identity based on its most similar scan in a database at an early age. Herein, we carried out three sets of experimental tasks to evaluate the identification accuracy of the proposed method. Specifically, for each feature or feature combination, the settings about the early age and later age of the datasets in three tasks are: (a) Year 0 and Year 1; (b) Year 0 and Year 2; (c) Year 1 and Year 2, respectively. Herein, in each task, each scan at a later age has its corresponding scan in the database at an early age. The identification accuracy was measured as the percentage of infants whose identity was correctly predicted out of the total number of infants to be identified in each task. Since the cerebral cortex develops much faster during the first year than during the second year (Gilmore et al., 2018; Li et al., 2014), the first two tasks are much more challenging than the third one. Thus, the first two tasks are our primary focus.

Let denote cortical feature maps, given a to‐be‐identified infant at a later age . First, for each feature type in each cortical ROI , we computed the pairwise similarity of feature maps between and each of the subjects { at an early age in the database, where denotes the total number of scans. To account for the dynamic brain size changes during development, we specifically used the Pearson correlation coefficient to measure the similarity between two feature maps, denoted as . Herein, =1, 2, and 3, correspond to features: mean curvature, average convexity, and sulcal depth, respectively; , correspond to 68 ROIs from the Desikan‐Killiany parcellation (Desikan et al., 2006) for left and right hemispheres; { correspond to the age groups in our three tasks, respectively. Second, for each ROI, all similarities calculated between the to‐be‐identified subject and all subjects in the database at an early age were sorted in decreasing order, and the identity of the subject with the largest similarity was regarded as a potential identity for the to‐be‐identified subject. In total, 68 potential identities were obtained from all the ROIs based on each specific feature type. For the feature combinations, the potential identities from all feature types were concatenated together, thus combining the contributions of multi‐view morphological features. Third, through a simple majority voting across all ROIs, the potential identity presenting the highest frequency among all the potential identities was considered as the final predicted identity for the subject to be identified. Of note, a threshold of the ratio between the frequencies of the first ranked potential identity and the second ranked potential identity was set to reject the to‐be‐identified subject for which there is no corresponding scan in the dataset. The threshold is empirically set as 2. Based on this threshold, if the frequency of the first ranked potential identity is not much higher than that of the second one, the to‐be‐identified subject would be rejected and would not obtain the final predicted identity. Finally, we evaluated our identification framework using each morphological feature type independently, as well as their different combinations.

Patterns of cortical identification capability and individual variability

Different cortical regions contribute differentially to individual identification. Their capacity to fingerprint the brain might be associated with how much they vary across subjects. To investigate this potential relationship, that is, understanding whether and how the regions with high identification capability are in accordance with the regions encoding high individual variability, for each cortical folding feature, we computed its ROI‐specific identification capability maps in the first two tasks and also computed its individual variability maps in all age groups. Thus, we can intuitively inspect and compare these two patterns. Of note, to explore their relationship, an alternative choice is to first find the regions with the largest individual variability, then use these regions for infant identification, and further determine their identification capability. However, this approach might miss out on key fingerprinting regions that do not necessarily largely vary across subjects. Hence, in this study, we opted for exploring the fingerprinting regions and the variable regions across subjects separately.

Mapping the ROI‐specific identification capability

To investigate which cortical ROIs most likely encode the individual fingerprint during dynamic postnatal brain development and thus contribute more to the individual identification, we charted the ROI‐specific identification capability maps. For each ROI, we calculated its identification accuracy as the percentage of the number of correctly identified infants by this ROI out of the total number of infants to be identified in each task. Then, we mapped the resulting ROI identification accuracies based on three types of cortical folding features onto the cortical surfaces, respectively. Herein, since the first two tasks are more challenging than the third one, we only show the representative maps of the first two tasks.

Mapping the longitudinal vertex‐wise individual variability

To explore the relationship between the regions with high individual variability and regions with high identification capability, we further computed vertex‐wise individual variability map of each cortical folding feature in each age group. Of note, the reason for choosing to establish the vertex‐wise rather than ROI‐wise map is that the former map would contain more fine‐scale information thus would more comprehensively display the individual variability distribution. To be consistent with the similarity measure in infant identification, instead of calculating the standard deviation of features at each vertex across all subjects, we computed the individual variability of each vertex based on Pearson correlation of its local feature maps of N‐ring neighborhood among individuals at the same age. Specifically, for each vertex , the variability across all subjects in the ‐th feature at age was quantified by averaging the feature map differences between 's N‐ring neighborhoods on cortical surfaces across all pairs of subjects (:where is the subject number; denotes the Pearson correlation coefficient. Of note, the ‐th ROI herein denotes the vertex 's N‐ring neighborhood, and N was empirically set as 15. Then, for each feature, the resulted variability values of all vertices at each age were mapped onto the atlas surface with the corresponding age, thereby producing the longitudinal maps of the vertex‐wise individual variability of cortical folding.

Twins study: Individual identification and cortical variability

Twins, including MZ twins and DZ twins, are ideal candidates for investigating the potential genetic and environmental influences on cortical folding, since they share similar genes and developmental environments. Recent studies demonstrate that MZ twins share 100% of their genes, while DZ twins share on average 50% of their segregating genes (Peper, Brouwer, Boomsma, Kahn, & Hulshoff Pol, 2007; Wachinger et al., 2015). In this work, we compared the differences between MZ twins and DZ twins in two aspects. On one hand, we calculated the identification accuracies of MZ twins and DZ twins in the above‐mentioned three tasks, respectively, and examined whether the difficulties for identifying the individuals from MZ and DZ twins are significantly different using two‐sample t‐test. Herein, similar to the experiments in Section 2.4, we identified the 1‐ and 2‐year‐old brains of each infant, who is one of the MZ or DZ twins, based on its own neonatal brain. Of note, the experiment on identifying one of the twins based on its twin brother/sister is not in the scope of this study. On the other hand, we calculated the difference between MZ twin pairs based on each specific cortical folding feature in each ROI at each age, and compared it with the corresponding difference between DZ twin pairs, thus investigating differences between different kinds of twins during early brain development. Since MZ twins share more genes compared to DZ twins, we infer that the difference of cortical regions between MZ twin pairs would be lower than that of between DZ twin pairs. To examine this inference, specifically, we defined the cortical folding difference across ‐th pair of MZ (or DZ) twins ( in ‐th ROI with ‐th feature as:where denotes the Pearson correlation coefficient. Herein, the regional cortical folding difference between MZ twin pairs and the difference between DZ twin pairs were calculated respectively. Then, the difference between the average difference of MZ twin pairs and the average difference of DZ twin pairs was defined as , where () is the total number of pairs of MZ (DZ) twins at age . For intuitive comparison of the ROI‐wised difference degree, we mapped the difference onto the longitudinal cortical surfaces. Besides, we further performed one‐tailed two‐sample t‐tests for testing in which ROIs the differences between MZ twin pairs are statistically significantly lower than those between DZ twin pairs. We also performed the one‐tailed test in the other direction to check whether there is any ROI in which the differences between MZ twin pairs are unexpectedly significantly higher than those between DZ twin pairs. We have controlled the false discovery rate (FDR) with the Benjamini‐Hochberg procedure for correction of multiple comparisons. We also mapped the ROIs with significant differences onto the longitudinal cortical surfaces.

RESULTS

Infant identification based on cortical folding features

Table 2 displays the infant identification accuracies in three tasks by using the proposed ROI‐based voting framework based on different combinations of cortical folding features. The corresponding identification accuracies of the global‐based method are also displayed in this table. We can clearly see that all three cortical folding features, that is, mean curvature, average convexity and sulcal depth, are capable of identifying infants in the three tasks using the ROI‐based voting framework. Specifically, most of accuracies based on the three kinds of features in the first two tasks are nearly 100%, and the accuracies in the third task are all 100%. We notice that the identification accuracies based on average convexity and sulcal depth are slightly higher than those based on mean curvature, and the identification accuracies in the first two tasks are slightly lower than those in the third task, but no statistical significant difference was found based on one‐tailed two‐sample t‐test with a significance level of p < .05. When using the combinations of multiple cortical folding features, the identification accuracies would further improve. For instance, the accuracies in the three target tasks are all 100%, when using the combination of all three folding features together. As for the global‐based method, the identification accuracies in the first two challenging tasks are statistically significantly lower than those (all 100%) in the third task, and are also statistically significantly lower than those in the first two tasks based on the ROI‐based voting framework, according to the results of one‐tailed t‐tests (p < .05).

Table 2

Infant identification accuracies of the proposed ROI‐based voting framework as well as the global‐based method using different combinations of cortical folding features

Feature combinations	Identification accuracy (%)
	Year 0 → Year 1 (387)		Year 0 → Year 2 (282)		Year 1 → Year 2 (197)
	ROI	Global	ROI	Global	ROI	Global
Curv	95.09	83.46	97.52	87.59	100.00	100.00
Conv	99.74	94.83	99.65	95.74	100.00	100.00
Depth	99.22	95.61	98.94	95.04	100.00	100.00
Curv & Conv	99.74	91.99	100.00	93.26	100.00	100.00
Curv & Depth	99.48	93.54	100.00	94.33	100.00	100.00
Conv & Depth	100.00	95.87	99.65	96.10	100.00	100.00
Curv & Conv & Depth	100.00	94.83	100.00	95.74	100.00	100.00

Note: The number in each bracket denotes the number of subjects to be identified in each task. Year 0 → Year 1: using the dataset with scans at Year 0 to predict the identities of scans at Year 1; Year 0 → Year 2: using the dataset with scans at Year 0 to predict the identities of scans at Year 2; Year 1 → Year 2: using the dataset with scans at Year 1 to predict the identities of scans at Year 2.

Abbreviations: Curv, mean curvature; Conv, average convexity; Depth, sulcal depth.

Besides, we also evaluated the false positive rate (FPR) and false negative rate (FNR) of the proposed method. To calculate the FPR and FNR, we performed the experiments in the inverse direction of the three tasks, that is, Year 1 → Year 0, Year 2 → Year 0, and Year 2 → Year 1, in which some to‐be‐identified subjects in Year 0/Year1 have no corresponding scans in Year 1/Year 2 dataset. When the threshold of the ratios is 2.0, the identification accuracies of three inverse tasks are 99.15%, 99.78% and 99.22%; the FPRs are 0.00%, 0.53% and 1.58%; the FNRs are 1.03%, 0.00% and 0.00%, respectively, as displayed in Table S1.

ROI‐specific identification capability patterns

For each cortical folding feature, Figure 2 displays corresponding maps of the identification capabilities of 68 ROIs in the first two tasks, respectively. Of note, we did not provide maps of the third task since in which the identification capabilities of most ROIs are very high due to the lower identification difficulty and thus the corresponding colors in the maps are almost red and indistinguishable. Table 3 displays the Pearson correlation coefficients of the 68 ROIs' identification accuracies across different features and different tasks, respectively, which quantitatively denote that these accuracy distributions are significantly positively correlated (the corresponding values of the correlations are all less than .05). According to Figure 2 and Table 3, we can conclude that the distributions of the regions with higher or lower identification accuracy are largely consistent across different tasks and different cortical folding features: mean curvature (Figure 2a), average convexity (Figure 2b), and sulcal depth (Figure 2c). Specifically, the precentral gyrus, rostral middle frontal gyrus, superior frontal gyrus, lateral occipital cortex, caudal middle frontal gyrus, and parietal cortex in both hemispheres and the middle and inferior temporal gyri in the left hemisphere show higher identification accuracies (described in roughly decreasing order of the identification accuracies). Figure 2d displays the average identification accuracies of three cortical folding features in each task. To numerically inspect these ROI accuracies, we also list the first 20 ROIs with the highest average identification accuracies in each task and the ROIs with higher average accuracies across both tasks in Table S2. From Figure 2d and Table S2, we can see that the rostral middle frontal gyrus in both hemispheres present the highest average identification capability in most cases. Then in decreasing order of the average accuracy of both hemispheres, the lateral occipital cortex, precentral gyrus, superior frontal gyrus, inferior parietal cortex, superior parietal cortex show higher identification capabilities. Besides, interestingly, Figure 2 and Table S2 both show that the identification accuracies of the precentral gyrus, inferior parietal gyrus, middle and inferior temporal gyri in the left hemispheres are higher than those in the right hemisphere, while the identification accuracies of the superior frontal gyrus in the left hemisphere are lower than those in the right hemisphere.

Figure 2

Maps of the identification accuracies of 68 ROIs in first two tasks based on cortical folding features: (a) mean curvature, (b) average convexity, and (c) sulcal depth, respectively. (d) displays the average identification accuracies of the above three features in each task. The first two columns show the lateral view, and the last two columns show the medial view. (Two tasks: Year 0 → Year 1, that is, using the dataset with scans at Year 0 to predict the identities of scans at Year 1; Year 0 → Year 2, that is, using the dataset with scans at Year 0 to predict the identities of scans at Year 2. Acc.: accuracy)

Table 3

Pearson correlation coefficients of 68 ROIs' identification accuracies across different features and different tasks

		Across features				Across tasks
		Curv vs Conv	Curv vs Depth	Conv vs Depth		Curv	Conv	Depth
r	Task 1	0.9	0.83	0.92	Task 1 vs 2	0.85	0.90	0.97
	Task 2	0.83	0.79	0.94	—	—	—	—

Note: Task 1: Year 0 → Year 1, that is, using the dataset with scans at Year 0 to predict the identities of scans at Year 1; Task 2: Year 0 → Year 2, that is, using the dataset with scans at Year 0 to predict the identities of scans at Year 2.

Abbreviations: Curv, mean curvature; Conv, average convexity; Depth, sulcal depth.

Longitudinal vertex‐wise individual variability patterns

Longitudinal maps of vertex‐wise individual variability based on mean curvature, average convexity, and sulcal depth, are respectively displayed in Figure 3. We can thus visually inspect the distribution patterns of the regions with high individual variability, and also intuitively compare it with that of the regions with high identification capability. It can be clearly seen that the regions with higher individual variability are largely consistent at 0, 1, and 2 years of age within each type of cortical folding features, and are also largely consistent across different cortical folding features. The above consistency of the individual variability distributions across different features and different ages is quantitatively measured by the Pearson correlation coefficients displayed in Table 4, which denote that all the distributions are consistent and significantly correlated with p < .05. Specifically, the middle frontal gyrus (including both rostral and caudal middle frontal gyrus), the bank of the superior temporal sulcus (STS), supramarginal gyrus, parts of the inferior and superior parietal cortices, the inferior temporal gyrus, and parts of the lateral occipital cortex exhibit higher individual variability in Figure 3. In particular, the middle and inferior temporal gyri and the bank of the STS in the left hemispheres show higher individual variability than those in the right hemispheres. These are largely in accordance with the discovered regions with high identification capability as shown in Figure 2 and Table S2. Interestingly, the middle temporal gyri in both hemispheres show lower individual variability than the inferior temporal gyrus; and the precentral gyrus, which presents high capability of infant identification, does not exhibit high individual variability based on cortical folding features. In contrast, the supramarginal gyrus, which does not show apparent high identification capability in Figure 2, shows high individual variability here, and meanwhile its variability is lower in the left hemisphere than in the right hemisphere. Additionally, in the medial view, most regions show lower individual variability, particularly in sulcal depth maps.

Figure 3

Longitudinal vertex‐wise individual variability maps at 0, 1, and 2 years of age with three types of cortical folding features: (a) mean curvature; (b) average convexity; (c) sulcal depth, based on 15‐ring neighborhoods of each vertex across all pairs of subjects at each age. The first two columns show the lateral view, and the last two columns display the medial view. (Var.: variability)

Table 4

Pearson correlation coefficients of vertex‐wise maps of individual variabilities across different features and different ages

		Across features				Across ages
		Curv vs Conv	Curv vs Depth	Conv vs Depth		Year 0 vs Year 1	Year 0 vs Year 2	Year 1 vs Year 2
r	Year 0	0.78	0.69	0.87	Curv	0.85	0.90	0.97
	Year 1	0.83	0.70	0.87	Conv	0.94	0.94	0.99
	Year 2	0.85	0.68	0.87	Depth	0.92	0.92	0.99

Abbreviations: Curv, mean curvature; Conv, average convexity; Depth, sulcal depth.

Furthermore, we also charted ROI‐wise individual variability maps in Figure 4 based on the correlation coefficients defined similarly for infant identification. Specifically, the individual variability value of each ROI was computed as the average difference (1‐correlation) of ROI‐wise feature maps across all pairs of subjects at each age. As we can see, the overall distribution patterns of ROI‐wise individual variability are largely consistent with that of vertex‐wise individual variability in Figure 3, while the latter ones display much finer‐scale information. Besides, to intuitively compare the proposed vertex‐wise neighborhood‐based measure in Equation (1) with the more common measure, we charted vertex‐wise individual variability maps based on standard deviation (SD) across subjects, as shown in Figure S1. As can be seen, both measures obtain similar individual variability patterns in all tasks. However, the vertex‐wise SD measure is more sensitive to noise compared to our proposed neighborhood‐based measure, especially in the maps of mean curvature feature. In addition, we also charted longitudinal vertex‐wise individual variability maps based on the correlation coefficients calculated with other different neighborhood sizes, for example, 10‐ring and 20‐ring, as shown in Figures S2‐S3. As we can see, the resulted individual variability patterns are largely consistent with the patterns revealed using 15‐ring neighborhood in Figure 3. Hence, we finally chose the middle one, that is, 15‐ring, to explore the individual variability patterns.

Figure 4

Longitudinal ROI‐wise individual variability maps at 0, 1, and 2 years of age with three types of cortical folding features: (a) mean curvature; (b) average convexity; (c) sulcal depth, based on the average difference (1‐correlation) of ROI‐wise feature maps across all pairs of subjects at each age. The first two columns show the lateral view, and the last two columns display the medial view. (Var.: variability)

Table 5 displays the identification accuracies for MZ and DZ twins using different combinations of cortical folding features, respectively. According to the results of two‐sample t‐tests with a significance level of p < .05, no statistical significant difference of the identification accuracy was found across different kinds of twins using both ROI‐ or global‐based methods. Moreover, based on the sulcal depth, we mapped the difference in each ROI between the average difference (defined in Equation ) across all pairs of MZ twins and the average difference across all pairs of DZ twins onto the cortical surface at the corresponding age, as shown in Figure 5. As we can see, most ROIs, including prefrontal gyrus, parietal gyrus, temporal gyrus et al., are in blue, which denote that the average differences between MZ twin pairs are lower than those between DZ twin pairs. Besides, these patterns are largely preserved during brain development in the first 2 years. Furthermore, as shown in Figures S4‐S5, we also displayed the numerical comparisons of the regional average differences of cortical folding between MZ twin pairs and that of between DZ twin pairs at each age in 20 cortical ROIs using sulcal depth and average convexity, respectively. Herein, these ROIs are the first 20 ROIs with the highest average identification accuracies as shown in Table S2. Notably, we can see from these figures that in the majority of cortical regions, the difference between MZ twin pairs is lower than that between DZ twin pairs, independent of both age and cortical folding feature. In addition, longitudinal cortical maps in Figure 6 show the results of the one‐tailed two‐sample t‐test (p < .05, FDR‐corrected) with the alternative hypothesis that the differences between MZ twin pairs are statistically significantly lower than those between DZ twin pairs. In these magenta ROIs in Figure 6, the above alternative hypothesis was accepted, and most of these ROIs overlapped with those ROIs in blue in Figure 5. As for the results of the one‐tailed test with an alternative hypothesis in the other direction, none of the ROIs has been found where the differences between MZ twin pairs are statistically significantly higher than those between DZ twin pairs.

Table 5

Identification accuracies for MZ twins and DZ twins based on both ROI‐based and global‐based methods with different combinations of cortical folding features, respectively

Features combination	Identification accuracy (%)
	Year 0 → year 1				Year 0 → year 2
	MZ twins (56)		DZ twins (106)		MZ twins (54)		DZ twins (62)
	ROI	Global	ROI	Global	ROI	Global	ROI	Global
Curv	96.43	83.93	96.23	85.85	98.15	88.89	98.39	88.71
Conv	100.00	98.21	100.00	95.28	100.00	96.30	98.39	95.16
Depth	100.00	96.43	99.06	96.23	100.00	94.44	98.39	95.16
Curv & Conv	100.00	94.64	100.00	92.45	100.00	94.44	100.00	93.55
Curv & Depth	100.00	94.64	99.06	95.28	100.00	96.30	100.00	93.55
Conv & Depth	100.00	98.21	100.00	96.23	100.00	96.30	98.39	95.16
Curv & Conv & Depth	100.00	98.21	100.00	96.23	100.00	96.30	100.00	95.16

Note: Number in bracket denotes the number of subjects to be identified in each task. Year 0 → Year 1: using the dataset with scans at Year 0 to predict the identities of scans at Year 1; Year 0 → Year 2: using the dataset with scans at Year 0 to predict the identities of scans at Year 2. (The identification accuracies in the third tasks [Year 1 → Year 2] are all 100% as shown in Table 2, thus they are not displayed here).

Abbreviations: Curv, mean curvature; Conv, average convexity; Depth, sulcal depth; MZ, monozygotic; DZ, dizygotic.

Figure 5

Longitudinal cortical maps of the ROI‐wise difference between the average sulcal depth difference across all pairs of monozygotic (MZ) twins and the average difference across all pairs of dizygotic (DZ) twins. Herein, the blue color indicates that the average difference of MZ twins is smaller than that of DZ twins. (Diff.: difference)

Figure 6

Longitudinal cortical maps of the ROIs (magenta) in which the differences between monozygotic (MZ) twin pairs are significantly lower than those between dizygotic (DZ) twin pairs. Herein, one‐tailed two‐sample t‐test is adopted with a significance level of p < .05, and the false discovery rate (FDR) has been controlled with the Benjamini‐Hochberg method for correction of multiple comparisons

DISCUSSION

The first main contribution of this study resides in the finding that the cortical folding morphology fingerprints the dynamic developing infant brain and is reliable for individual identification during the first postnatal years. Despite the dramatic global change in cortical size, shape and folding features during birth and 2 years of age (Li, Nie, Wang, Shi, Lin, et al., 2013a; Li, Wang, Shi, Lyall et al., 2014; Lyall et al., 2014; Meng et al., 2014), as also shown in Figure 1, we achieved promising accuracies in identifying 1‐ and 2‐year‐old brains from neonatal cortices using the combinations of cortical folding features (Table 2 ). More importantly, we can thus anticipate that the evidenced fingerprinting power of the neonatal brain of a specific subject can be carried out across the whole lifespan. The reasons for this assumption are in two aspects. First, all major cortical folds and individual variability patterns of the human brain are established at term birth (Duan et al., 2018; Hill et al., 2010; Li, Wang, Shi, Lin, & Shen, 2014a). Second, the most dynamic phase of postnatal brain development is the first 2 years of life, and the folding patterns only undergo minor changes during later childhood and adulthood, thus the 2‐year‐olds' brains largely resemble the adult brains in cortical folding (Gilmore et al., 2007; Li, Nie, Wang, Shi, Lin, et al., 2013a; Li, Nie, Wang, Shi, Gilmore, et al., 2014). Hence, once the 2‐year‐olds can be correctly identified, the possibility of identifying the adult brains based on their neonatal cortical folding patterns would be very high. However, further investigations are required to validate this assumption using a larger longitudinal dataset covering both developing and adult periods.

Table 2 provides us further insights into the infant identification tasks from neonatal cortical folding. Specifically, first, the combination of three kinds of cortical folding features can slightly improve identification accuracies compared to any single feature. Though the improvement is not significant, we prefer to adopt the combinations of three features into the proposed individual identification framework because of two reasons: (a) the mean curvature, average convexity, and sulcal depth provide complementary morphological information of cortical folding from different aspects, as described at the beginning of Section 2.3; (b) the accuracies are all 100% in all tasks, outperforming any single feature and other feature combinations. Second, the identification accuracies in the first two tasks using neonatal brain to identify 1‐ and 2‐year‐olds are lower than that in the third task using 1‐year‐olds to identify 2‐year‐olds. Compared to the first two tasks, the third task is more similar to the adult individual identification, due to the moderate change of cortical folding from Year 1 to Year 2. Thus, these results indirectly validate that the infant individual identification involving neonates is much more challenging than the adult individual identification. Furthermore, the identification accuracies in the first task (i.e., using the dataset with scans at Year 0 to predict the identities of scans at Year 1) are sometimes lower than that of the second task (i.e., using the dataset with scans at Year 0 to predict the identities of scans at Year 2). It might seem less reasonable, since the first task should be easier than the second one, because of the smaller brain development in the first year, in comparison to the first 2 years. To analyze whether this result is caused by the imbalanced datasets in the first two tasks, we repeated experiments with balanced datasets based on both ROI‐based and global‐based methods, as shown in Table S3. Here, to obtain the balanced Year 1 and Year 2 datasets, we randomly selected 200 subjects for 10 times from their original datasets, respectively. Table S3 shows the averaged accuracies of 10 times experiments. As we can see, the accuracies in Task 1 are still lower than Task 2 in some experiments. Excluding the reason of imbalanced datasets, we speculate that the different surface quality in Year 1 and Year 2 datasets might be responsible for this unexpected result in Table 2. Specifically, in the image processing pipeline, the surfaces in Year 0 dataset are reconstructed based on the segmentation results obtained from T2‐weighted images, which show better tissue contrast than the T1‐weighted images of neonatal brains, while the surfaces in Year 1 and Year 2 datasets are reconstructed from T1‐weighted images. Due to the poorer contrast of T1‐weighted images at Year 1 compared with those at Year 2, the surface quality of images in Year 1 dataset in the first task is poorer than that in Year 2 dataset in the second task, which thus might lead to the unexpected slightly lower identification accuracies in the first task.

To handle the case where the subject to be identified has no corresponding scan in the dataset, we set a threshold of the ratio between the frequencies of the first ranked potential identity and the second ranked potential identity. We recorded the ratios during all subjects' identification procedures, and found that the minimum ratios in the first two tasks are 2.0 and 2.2, respectively. The distributions of the ratios are displayed in histograms in Figure 7. Of note, choosing a proper threshold of the ratios is important for the individual identification method. If the threshold is too large, the FPR would be 0, but the FNR would be large; if it is too small, the FPR would be large, and the FNR would be 0. In both situations with improper thresholds, the identification accuracies would be low. The accuracies, FPRs and FNRs based on different thresholds in inverse tasks are displayed in Table S1. Here, we set the threshold to 2.0 according to the above minimum ratio in the first two tasks. Based on thresholding, if a new coming scan has no corresponding scan in the early dataset, we would reject it, thus controlling the false discovery rate.

Figure 7

Histogram distributions of the ratios between the frequencies of the first ranked potential identity and the second ranked potential identity during the individual identification procedure in the first two tasks: (a) using the dataset with scans at Year 0 to predict the identities of scans at Year 1; (b) using the dataset with scans at Year 0 to predict the identities of scans at Year 2

Relation of cortical identification capability and individual variability patterns

The second main contribution of our study is that we found that the infant cortex shows regionally variable identification accuracy; and the regions with high identification capability are largely overlapped with the regions exhibiting high individual variability. Interestingly, these structural regions are also highly in line with the functional regions showing high individual identification capability and high individual variability (Finn et al., 2015; Finn et al., 2017; Liu et al., 2018; Mueller et al., 2013). This validates that the individual variability might be able to link brain structure to brain function (Frost & Goebel, 2012; MacDonald et al., 2006).

Cortical identification capability patterns

Based on the identification accuracies of ROIs from the Desikan‐Killiany parcellation (Desikan et al., 2006), we found that the regions with higher identification capability on all cortical folding features are mostly located in high‐order association cortices (including the rostral middle frontal gyrus, superior frontal gyrus, caudal middle frontal gyrus, inferior parietal cortex, middle and inferior temporal gyrus) and two unimodal cortices (including precentral gyrus and lateral occipital cortex). Furthermore, this discovery is consistent to some extent with related studies on functional connectome fingerprinting. For example, Miranda‐Dominguez et al. found that the unique connections in high‐order association regions in frontal and parietal cortices mainly contribute to the individual identification (Miranda‐Dominguez et al., 2014), and they further demonstrated that the high‐order systems, including the frontoparietal, dorsal attention, ventral attention, cingulo‐opercular and default systems, present substantial heritability variability (Miranda‐Dominguez et al., 2017). Horien et al. found that both medial frontal and frontoparietal networks are unique and stable for individual identification during infancy (Horien et al., 2019). Finn et al. found that the high‐order frontoparietal functional network, comprised of frontal, parietal and temporal lobes, is of the highest identification power (Finn et al., 2015). However, in their study, the primary visual network and the motor network are highly consistent across subjects, thus do not contribute substantially to the individual identification when using functional MRI, which contrasts our observation on cortical folding morphology. Interestingly, there is another functional MRI study largely supporting our finding. Specifically, Liu et al. found that not only the frontoparietal network, but also other two high‐order cognitive systems (i.e., default mode and dorsal attention networks) and two primary systems (i.e., visual and sensorimotor networks) exhibited high identification power (Liu et al., 2018).

Cortical individual variability patterns

Based on the longitudinal vertex‐wise individual variability maps, we found that the regions with higher individual variability on cortical folding features are mostly distributed in high‐order association cortices (including the middle frontal gyrus, the bank of the STS, middle and inferior temporal gyrus, part of the inferior and superior parietal cortices and part of the lateral occipital cortex); and the regions with lower individual variability are mainly located in unimodal cortices (including sensorimotor cortex and auditory cortex). Several studies on individual variability support this finding. Specifically, for individual variability of cortical folding, Hill et al. found that it is higher in association areas and lower in the motor cortex (Hill et al., 2010). As for the individual variability of functional connectivity, it is significantly higher in the high‐order association cortex and lower in unimodal cortices (Finn et al., 2017; Gao et al., 2014; Miranda‐Dominguez et al., 2017; Mueller et al., 2013). We also found that the individual variability patterns of infant cortical folding are relatively consistent during the first 2 years of life, and are also relatively consistent across different features. This discovery is also largely in line with previous studies (Cao, Huang, & He, 2017; Gao et al., 2014).

Notably, the regions with high identification capability and those with high individual variability largely overlap in high‐order association cortices. This is reasonable since the regions characterizing high individual variability should be more distinctive across different individuals. Interestingly, the precentral gyrus (motor cortex) presents high identification capability (Figure 2), despite its low vertex‐wise individual variability (Figure 3). This might be explained by the junction structure of the middle frontal gyrus and the precentral gyrus, which provides essential distinct morphological information with high priority for individual identification. This also supports our primary choice to separately identify the most fingerprinting regions and establish the individual variability maps.

The third main contribution of our study is the discovery that both the MZ and DZ twins' brains can be correctly identified using the cortical folding features despite similar genetic and environmental influences. Table 5 demonstrates that the cortical folding features are reliable for identifying both infant MZ and DZ twins, and there is no statistical significant difference in the difficulty degree between their identification. Besides, the accuracies show slightly higher values than the corresponding identification accuracies in Table 2, but no significant improvement was found. The slight difference might be caused by the largely imbalanced datasets in the tasks of individual and twin identification. Moreover, Figure 5 and Figures S4‐S5 show that the discordance between MZ twin pairs in most cortical regions, especially the high‐order association cortices, is generally lower than that between DZ twin pairs. Figure 6 further validates that in most of these high‐order association cortices, the degree of discordance between MZ twin pairs is significantly lower than that of between DZ twin pairs, which is in line with the universal biological principles (Kaminsky et al., 2009). In a few regions in Figure 5, the differences between MZ twin pairs are slightly larger than the difference between DZ twin pairs, but no significance was found according to the results of the one‐tailed test. According to the results of these statistical tests, we can have an interesting conclusion that only in some high‐order association cortices, the differences between MZ twin pairs are significantly lower than those between DZ twin pairs, while in other ROIs, there is no statistical significant difference between the discordance of MZ and DZ twins.

It is interesting that the MZ twins, who share the identical genetic makeup (DNA) from a single fertilized egg (Jain, Prabhakar, & Pankanti, 2002; Patwari & Lee, 2008), present distinctive cortical folding patterns in infants. Though the underlying reason is still unclear, recent studies found that the genetically‐identical cells and organisms are not an entirely genetic characteristic, but influenced by both genetic and environmental factors in a dynamic and complex manner (Jha et al., 2018; Raser & O'shea, 2005). Specifically, first, the variation in gene expression may contribute to the phenotypic variability (Patwari & Lee, 2008; Raser & O'shea, 2005) of cortical folding patterns. Second, the prenatal environmental factors (Patwari & Lee, 2008), including the umbilical cord length, access to the nutrition, blood pressure, and position in the womb, also play import roles during the prenatal dynamic development of cortical folding. To this end, we can conclude that (a) both genetic and environmental factors could influence the early development of cortical folding morphology; (b) individual identification based on cortical folding is valid and promising, because there are no two identical brains, even for MZ twins.

Additional considerations

There are several issues that require further considerations as listed below.

Cortical parcellation choice

From a methodological perspective, the scale and definition of the ROI might influence the patterns of regions with high identification accuracy to some extent. With a cortical region showing high identification accuracy, it might be hard to know which specific part of this region is more critical, given the large sizes of some ROIs from Desikan‐Killiany parcellation (Desikan et al., 2006). Future work may explore the patterns of identification accuracy with a finer‐scale ROI parcellation to better understand which sub‐regions are more or less contributive for individual identification and to better inspect the relation of cortical identification capability and individual variability patterns.

Longitudinal individual cortical folding study across lifespan

It remains unclear to which extent the individual cortical folding is consistent across the whole lifespan. This would be ideally explored using a longitudinal brain dataset with follow‐up MR images from birth to adulthood—which is currently nonexistent. Future studies should include further collecting follow‐up images, and exploring the consistent aspects and developmental aspects of cortical folding during brain development. Additionally, since our experiments with promising identification accuracies were carried out on healthy subjects, it is still unclear that whether specific neurodevelopmental disorder would influence the individual variability and fingerprinting power of the cortical folding. However, it is promising since some studies found that the descriptor of brain morphology can be used to effectively identify adult individuals with Alzheimer's Disease (Peper et al., 2007; Wachinger et al., 2015). Further studies would be required to validate this assumption based on datasets including both healthy subjects and subjects with neurodevelopmental disorders. It would constitute a formidable step forward to demonstrate this in brain development and maturation research.

Infant identification

Though our proposed method based on cortical folding features achieves promising identification accuracies (100%), we do realize that it is not a realistic way for infant identification at present, since MRI is a relatively slow and expensive imaging examination until now. Of note, our main innovation of this study is not the real application of infant identification but rather the three neuroscientific discoveries we reported. Thus, we concisely review the background of current infant identification methods as follows. To our knowledge, this is the first study to leverage the developing cortical folding as the biometric trait for infant identification. There are a few infant identification studies based on other conventional biometric traits, for example, fingerprint (Jain et al., 2016), footprint (E. Liu, 2017), face (Bharadwaj, Bhatt, Singh, Vatsa, & Singh, 2010), or iris (Corby et al., 2006). Compared to the cortical folding features, these biometric traits are more convenient to acquire. However, their performance is less promising, especially when involving neonates (typically < 90% in accuracy) due to the rapidly changing biometric traits during infancy. Besides, these exterior biometric traits are typically unstable and easy to be artificially changed or imitated on purpose in the real application. In future, once brain MRI becomes fast, convenient and cheap to acquire, cortical folding could potentially be a reliable biometric trait for infant identification.

CONCLUSION

In this study, we investigated infant identification and individual variability using a large‐scale longitudinal dataset with 1,141 scans. The contribution of this work is threefold. First, for the first time, we reported that the cortical folding during dynamic early brain development is unique and effective for individual identification. The identification accuracies of using neonates to identify 1‐year‐olds and 2‐year‐olds all reach 100% with the combinations of cortical folding features, that is, mean curvature, average convexity, and sulcal depth. We can also infer that the neonatal brain could be used for individual identification during the whole lifespan. Second, we discovered that the regions with high identification capability mainly distributed in high‐order association cortices (including prefrontal, lateral temporal, and inferior parietal regions) and two unimodal cortices (including precentral gyrus and lateral occipital cortex). Furthermore, the patterns of individual variability are age‐consistent and region‐specific, which largely overlap with the distribution patterns of identification capability. More interestingly, these discovered regions with large structural individual variability are also largely in line with the regions encoding functional individual variability. Third, we also found that even for MZ twins, who share identical genes and similar developmental environment, their cortical folding patterns are unique fingerprints, as the cortical folding might be shaped by complex gene expression and prenatal environmental factors. As for the cortical variability of twins, we found that only in some high‐order association cortices, the differences between MZ twin pairs are significantly lower than those between DZ twin pairs, while in other cortices, no significant difference was found. This study thus provides important insights into individual identification and individual variability based on cortical folding during early brain development.

CONFLICT OF INTERESTS

The authors report no conflict of interest.

Data S1: Supporting Information

Click here for additional data file.