Heritability of Structural Patterning in the Human Cerebral Cortex

Genetic influences that govern the spatial patterning of the human cortex and its structural variability are still incompletely known. We analyzed structural MR images in twins, siblings, and pairs of unrelated subjects. A comprehensive set of methods was employed to quantify properties of cortical features at different spatial scales. Measures were used to assess the influence of genetic similarity on structural patterning. Results indicated that: (1) Genetic effects significantly influence all structural features assessed here at all spatial resolutions, albeit at different strengths. (2) While strong genetic effects were found at the whole-brain and hemisphere level, effects were weaker at the regional and vertex level, depending on the measure under study. (3) Besides cortical thickness, sulcal (geodesic) depth was found to be under strong genetic control. The local pattern indicated that two axes along (a) the anterior-posterior direction (insula to parieto-occipital sulcus), and (b) superior-inferior direction (central sulcus to callosal sulcus) presumably determine the segregation of four quadrants in each hemisphere early in development. (4) While strong structural asymmetries were found at the regional level, genetic influences on laterality were relatively minor.


Introduction
Billions of human brains develop with astounding similarity, yet intricate individual detail. Decades of research have provided an increasingly detailed understanding of the microscopic development of the human brain (reviews: Bae et al., 2015 ;Geschwind and Rakic, 2013 ), and recent recent large-scale genome-wide association studies enumerated an ever-expanding list of genes involved brain development (e.g., Fiddes et al., 2018 ;Elliot et al., 2018 ). However, determinants of the spatiotemporal gyrification process of the human brain and its regional specialization are largely unknown (reviews: Llinares-Benadero and Borrell, 2019 ;O'Leary et al., 2007 ;Silbereis et al., 2016 ).
Cross-sectional MRI studies involved pairs of monozygotic and dizygotic twins to study the influence of genetic factors on the brain (re-ation study identified 187 loci influencing surface area and 12 influencing cortical thickness ( Grasby and Jahanshad, 2020 ). In a nutshell, genetic influences on cortical thickness and regional surface size were spatially uniform, with minimal regional differences. However, local genetic correlations were often not significant or even negative, indicating that an areal expansion leads to a thinner cortex. Genetic mechanisms that control thickness and size considered as independent Chen et al., 2013;Grasby and Jahanshad, 2020;Winkler et al., 2010 ). Recent advances in the neurobiology of brain development offer some insight into the modulation of cortical thickness and its relation to regional expansion, leading to the hypothesis that a spatially differential cortical expansion drives gyrification ( Llinares-Benadero and Borrell, 2019;Namba and Huttner, 2017;Rash et al., 2019 ). Considerable doubt remains that mechanisms governing regional differences in cortical thickness and areal expansion alone explain the astounding similarity of the brain's folding pattern. Surprisingly, few studies used shape measures to examine developmental and genetic influences on the cortex (e.g., Awate et al., 2010 ;Fish et al., 2017 ;Le Guen et al., 2018a,b ).
In this study, we aimed at providing a comprehensive, quantitative assessment of structural features of the human cortex. We developed a comprehensive battery of 11 methods to assess features across three spatial scales: from the whole brain and its hemispheres to the regional level (10-100 cm 2 ), and finally to the vertex (mm 2 ) level. Besides the wellexamined parameters cortical thickness and size, we studied geodesic depth, surface curvature, and the myelin fraction (T1/T2 ratio), in addition to several composite measures characterizing sulcal and gyral shape. Methods of this battery were designed with some redundancy to allow comparisons across spatial scales and different levels of complexity to better detect and understand method-related influences on the results. In addition to the determination of genetic influences per se, we assessed sex-related differences and hemispheric asymmetries.
To perform a statistical evaluation across individual brains, an approach for spatial normalization and registration is necessary. However, a "perfect " method would remove all inter-subject differences, and current methods for non-linear registration of surface use surface properties that we actually plan to assess. In contrast to previous studies, we did not use an atlas-based scheme to define cortical regions, because non-linear registration is necessary for atlas adaption.
Based on our previous work, we used sulcal basins ( Kruggel, 2018;Yang and Kruggel, 2008 ) that denote a cortical patch, centered around locally deepest points, including the neighboring sulcal walls up to the gyral crowns ( Fig. 1 C). It has been hypothesized that the deepest local points of the basins ( sulcal pits ) are location-invariant and expressed early in development ( Im et al., 2010;Im and Grant, 2019;Le Guen et al., 2018a;Meng et al., 2014;Regis et al., 1995;. Basins are automatically segmented on an individual surface representing the gray/white matter (GM/WM) interface using local properties surface curvature and geodesic depth ( Fig. 1 A, B). While basins are individual structures, we demonstrated that basins cluster into 13 communities ( Fig. 2 ) that are presumably common within a population and provide a registration-free assessment of regional properties ( Kruggel, 2018;Kruggel and Solodkin, 2019 ).
We used this test battery to analyze structural features of the human cortex in twins, siblings, and pairs of unrelated subjects included in the large and well documented data base of the Human Connectome Project (2020) . We aimed at providing a comprehensive assessment of genetic influences on structural patterning in the human cortex highlighting the relationship to brain development.

Subjects and Imaging Data
Based on imaging data and genetic testing of the 1113 subjects in the "1200 Subjects Release " of the Human Connectome Project (2020) , we selected 992 subjects grouped into pairs: Group MZ: 138 pairs of monozygotic twins verified by genetic testing; Group DZ: 78 pairs of dizygotic twins verified by genetic testing; Group SB: 178 pairs of siblings that randomly selected from their known family background; Group UN: 102 pairs of unrelated subjects. The MZ group included 81 female and 57 male pairs. All pairs of the DZ group had the same sex: 48 female and 30 male pairs. Sex was balanced within pairs in groups SB and UN. The mean age was similar between groups (Wilcoxon test).

Segmentation of Cortical Features
As the analytical procedures described in this section were extended from our previous work ( Kruggel, 2018;Kruggel and Solodkin, 2019 ), we included only relevant details.
Generation of hemispheric surfaces: Step 1: Native T1-and T2weighted structural images were co-registered, corrected for intensity inhomogeneities, and the intracranial space was extracted.
Step 2: This space was classified into four compartments, roughly corresponding to GM, WM, cerebro-spinal fluid (CSF) and connective tissue. The inner cavities of the WM segmentation (e.g., inner ventricles, basal ganglia) were filled to form a binary object with genus zero. The cerebellum and brain stem were clipped 15 mm below the AC-PC plane, and split into hemispheres at the mid-sagittal plane.
Step 3: A triangulated surface was computed from this object, and optimally adapted to the GM/WM interface using the intensity-corrected T1-weighted data set. Cortical thickness ( Osechinskiy and Kruggel, 2012 ) and ratio of the voxel-wise T1-weighted/T2-weighted intensity were determined at each vertex location. This ratio was used as a proxy for myelin content ( Glasser and van Essen, 2011 ) Segmentation of cortical basins: Surface curvature was computed from the hemispheric mesh, represented by the shape index ( Fig. 1 A). Geodesic depth was determined in image space, using a constrained distance transform on the sulcal compartment and interpolated at vertex positions ( Fig. 1 B). Finally, basins were segmented by a watershedregion growing process guided by surface curvature and geodesic depth. Starting from seeds at locally deepest vertices in convex regions, this process grew regions by adding neighboring vertices, prioritized by decreasing depth and curvature until regions match at cortical rims. Each basin received a unique label stored at each vertex ( Fig. 1 C). Regions that were clipped at the mid-sagittal plane and the brain stem were excluded from the process. This procedure yielded 100-140 basins per hemisphere.
Inter-subject alignment: Step 4: To compare data across individuals, each hemispheric mesh was unfolded to a unit sphere while minimizing angle and area distortion ( Kruggel, 2008 ). The overall correspondence between individual spherical meshes was maximized by finding a rotation that optimizes the normalized mutual information (NMI) of the vertex-wise basin labels with an arbitrarily chosen reference ( Fig. 1 D). Basin labels, geodesic depth, curvature, cortical thickness, and myelin ratio were re-sampled on a common icosahedral mesh while (approximately) retaining the spatial resolution. Note that this step introduces a spatial normalization. Thus, data for each hemisphere were stored as a 3D matrix of the 992 subjects by 163842 vertices by 5 features, along with the spherical mesh that represented vertex positions and neighborhood relationships. To maximize pair-wise similarity, data of member B was linearly registered to member A, using NMI as similarity metric.
Community segmentation: Step 5: Neighboring basins were clustered into communities based on their mutual overlap across individuals using a graph-based algorithm ( Campigotto et al., 2014 ), as previously described ( Kruggel, 2018 ). Using each subject as a reference, 992 individual segmentations into communities were obtained. Individual segmentations were combined into a cohort-level map that expressed the vertex-wise variability of community labels. Thus, regions of low variability are commonly found at similar locations with similar features across the cohort.

Table 1
Overview of methods M1-M11 included in the test battery, in relation to the structural level (column 1), the data type assessed (column 2), and measure type (columns 3 and 4). Glossary of abbreviations: AR = hemispheric surface area; BGS = pair-wise similarity of basin graphs; BLD = pair-wise vertex-vertex distance of sulcal bottom lines; BLS = pair-wise similarity of basin labels; BRR = ratio of brain inside cranium; BRV = brain volume; CLV = center of low variability; GMV = grey matter volume; ICV = intracranial volume; ISIM = pair-wise correlation of voxel-wise intensities in image space; SA = pair-wise overlap of vertices in sulcal areas; SSIM = pair-wise vertex-vertex distance WMV = white matter volume. (13 on the right) hemisphere ( Fig. 2 ). For further details, refer to Kruggel and Solodkin (2019) . The set of individual basins that comprise a community is found by graduated assignment from the cohort-level community map.

Structural Features
A comprehensive set of 11 methods ( Table 1 ) was compiled to extract structural features at different structural resolution levels: (a) the brain and (b) its hemispheres; (c) cortical communities; (d) vertices of the cortical surface mesh. Methods in column 3 quantify individual features, while methods in column 4 rate the within-pair similarity. For each method, we describe data and processing, and the reasoning for inclusion in this battery in the paragraphs below. Results compiled in the Tables and Figures of the next section refer to method numbers and abbreviation for derived measures.
Some redundancy is intended in this battery: Firstly, some methods assess the same properties at different resolution levels (e.g., grey matter volume, cortical thickness at the regional and vertex level). Thus, the number of observations that yield an individual measure varied widely. An increase in spatial resolution comes with a higher error variance. In this analytical context, larger random effects E lead to lower heritability estimates A . Thus, a careful evaluation across resolution levels is needed to compare heritability estimates across levels. Secondly, the same properties were used to characterize different structural features (e.g., curvature and geodesic depth to segment basins and sulcal areas), as we wanted to assess the "latent variable " heritability from different observations.
Method 1 -Compartment volumes: The probabilistic segmentation of the extracted brain into classes CSF, GM, and WM ( Section 2.2 , step 2) was used to compute compartment volumes. The probability value in each voxel was interpreted as a volume fraction of the corresponding material, and summed up over the whole domain. The brain volume (BRV) corresponds to the sum of GM and WM volumes, the intracranial volume (ICV) to the sum of the brain and CSF volume, while the brain ratio (BRR) is given by the quotient BRV/ICV. Note that all measures include the cerebellum and brain stem. Except for BRR, values are given in cm 3 . These individual measures are easily obtained and very robust but bear little structural information, as the same volume can be occupied by highly different structures.
Method 2 -Hemispheric similarity in image space: The extracted brain ( Section 2.2 , step 1) of member A of each pair was linearly registered with member B using correlation as similarity metric. The within-pair similarity of the left (resp. right) hemisphere (ISIM LH, RH) is represented by the correlation coefficient ( Bartley et al., 1997 ). This correlation coefficient provides a simple and robust assessment of global structural pair-wise similarity in image space.
Method 3 -Hemispheric surfaces: Surface areas of both hemispheres (AR LH, RH) were computed from the meshes optimally adapted to the GM/WM interface ( Section 2.2 , step 3). Results are given in cm 2 . Note that surface optimization depends on critical parameter choices. While absolute area measures depend on these settings, relative differences as employed here are robust. Of note, structure is not examined here, as highly different structures can have the same area.
Method 4 -Hemispheric similarity in surface space: Triangle meshes representing the GM/WM interface ( Section 2.2 , step 3) were used to compute the shortest Euclidean distance between vertices on the left (resp. right) side (SSIM LH, RH). Results are given as the mean distance in mm. Similar to the image-based correlation coefficient, this measure assesses structural pair-wise similarity of surface shape at the hemispheric level.
Method 5 -Community-level statistics: Previous studies (e.g., Im et al., 2010 ;Meng et al., 2014 ;Regis et al., 2005 ) defined a consistent set of homologue basins across subjects. Depending on definitions, their set size varied between 40 and 70, and left out shallow structures that occupy a considerable portion of the surface. Because we aimed at comparing the patterning of the whole surface, we refrained from establishing basinlevel correspondences, and examined similarity at the community level. The individual set of basin labels was determined for a given community ( Section 2.2 , step 6). For each community, we computed: (a) the number of basins, (b) the size (i.e., the number of vertices), (c) the maximum geodesic depth (in mm), (d) the fraction of vertices in sulcal areas, (e) the mean cortical thickness (in mm), and (f) the mean myelin fraction, based on data obtained in ( Section 2.2 , step 4). These measures quantify individual, regional properties of cortical shape. Method 6 -Community-level similarity of basins: A matrix of the cooccurrence of basin labels in both members of a pair was computed. The NMI of this matrix reflects the similarity of basin patterns between homologue communities of a pair. Basin segmentation depends on the choice of surface filter settings and the approach of using curvature and depth in the region growing process. While the absolute number and extent of basins depend on these choices, relative differences as employed here are robust.
Method 7 -Community-level similarity of sulcal areas: Deep sulcal regions develop early, are less variable than gyral crowns and likely under stronger genetic control ( Clouchoux et al., 2012 ;Dubois et al., 2017 ;Kruggel and Solodkin, 2019 ;Le Guen et al., 2018a ). Regions of deep sulcal vertices (shape index < 0 and geodesic depth > 8 mm) were determined in both members ( Section 2.2 , step 4). The overlap of regions between members was expressed by the Dice similarity index. In comparison with the NMI measure, this definition of sulcal areas is trivial and robust but assesses a fraction of the surface only.
Method 8 -Community-level similarity of sulcal bottom lines: Sulcal bottom lines were often used as structural features for establishing intersubject correspondences. Regions of deep sulcal vertices obtained in method 6 were separated into connected components. Each vertex was assigned a weight corresponding to the product of geodesic depth and shape index. Local minima were selected as deep sulcal vertices. In each connected component, neighboring vertices were connected by shortest paths constrained to deep sulcal areas. Paths were found by Dijkstra's algorithm using the reciprocal of vertex weights. All sets of connected vertices were denoted as sulcal bottom lines. The mean distance between lines in subject pairs was considered as the similarity metric. Because vertices reside on a sphere, we used the angle between unit vectors as (angular) distance here. In comparison with the area-based measure, a sparser representation of sulci was chosen, potentially less robust, but likely less variable.
Method 9 -Centers of low variability: Regions of low variability were segmented from the cohort-wise map obtained in ( Section 2.2 , step 5). The most invariant region per community was selected ( Kruggel and Solodkin, 2019 ). For each subject and community, the deepest vertex within the corresponding search region was detected. We recorded the within-pair absolute difference in geodesic depth (in mm) and used the angular distance between the deepest vertices as similarity metric. Here, a single, invariant location was used to represent a community. Communities are understood as the "common structural core " underlying all subjects. Thus, these locations are expected to be genetically constrained across all subjects, but depend less on subject-or pair-wise differences in genotype.
Method 10 -Community-level similarity of basin graphs: Im et al. (2011) developed a measure for the similarity of basin graphs. These graphs consist of nodes corresponding to basins, and edges to their neighborhood relationships. Node attributes are: (a) the position of the deepest point (in mm), (b) the size (number of vertices), and (c) the number of neighboring basins. Edge weights include the distance between neighboring basins and the ridge height of the gyrus that separates neighboring basins. Graphs are constructed from the basin configuration in a community of member A and compared with the corresponding graph of member B of a subject pair. Spectral matching was used to find the optimal node assignment between graphs ( Leordeanu and Hebert, 2005 ). In contrast to Im et al. (2011) , we computed graphs in normalized spherical space ( Section 2.2 , step 4), with weights optimized for this discrimination problem. This measure provides a complex assessment of the local cortical patterning, and was proposed in the context of addressing anatomical labels to sulci. However, results critically depend on the setting of weights used to integrate the different features.
Method 11 -Vertex-wise features: At the most detailed structural level, individual, vertex-wise estimates of geodesic depth, surface curvature, cortical thickness, and myelin fraction ( Section 2.2 , step 4) were assessed. Vertex-wise estimates of additive genetic variance (heritability A ) were color-coded and mapped onto an "inflated " hemispheric surface.
The key issue here is how subject-wise measures are sampled at homologue brain locations within the cohort. Suppose that structural features in monozygotic twins were identical. Still, between-pair differences are expected to compare with those found in unrelated subjects. Most commonly, filtering in the spatial domain (e.g., by Gaussian kernels of 20 mm width) or nonlinear registration was used to reduce individual variability. Unfortunately, the first approach also lowers spatial resolution, while the second one uses the same (or correlated) features for registration than those assessed later as phenotype trait. Thus, we used two approaches here. Firstly, we computed vertex-wise statistics of heritability estimates that were based on a native vertex-wise definition of homology. Thus, maps are confounded by the inter-individual variability of the cortical patterning. Secondly, we corrected for individual detail based on structural filtering. We used a wavelet filter basis on the sphere ( Yeo et al., 2008 ) and reconstructed vertex-wise measures at six resolution levels, and determined heritability estimates at each level. Especially level 3 is considered to represent the cortical shape at gestation week 26, before individual features develop.
Source code for the methods provided here will be made available on request to the communicating author.

Statistical Models
Two statistical models were employed to analyze structural features obtained with the above methods: (a) individual, paired measures were assessed for additive genetic effects (heritability A ); (b) within-pair similarity measures were used to estimate heritability R according to proportion of genes identical by descent (IBD).
Model A -Statistics for individual measures: To analyze individual, paired phenotype data for genetic influences, a generalized DeFries-Fulker (gDF) regression model was employed ( Lazzeroni and Ray, 2013 ): where y i ,1 , y i ,2 correspond to the quantitative phenotypes of pair i , and x i ,1, k , x i ,2, k to covariates x k for sex (female or male), body height (in cm), age (in years) and race , including an intercept , , 1 = 1 .
The classical ACE twin model ( Horst and Scheike, 2015 ; Neale and Maes, 2004 ) uses a linear mixed effects model to decompose the residual variance into components addressed to additive genetic effects A , common environmental effects C , and random errors E , with A ≥ 0, C ≥ 0 and + + = 1 . In this gDF model, components A and C were estimated by regression. The regressor I A (corresponding to coefficient A ) represented the proportion of genes identical by descent (IBD), which was set to 1 for MZ pairs, 0.5 for DZ and SB pairs, and 0 for UN pairs. Similarly, the regressor I C modeled common environment factors C , which was set to 1 for MZ, DZ and SB pairs, and 0 for UN pairs. The gDF model is simpler to implement and computationally more efficient, which is relevant for the vertex-wise computation of heritability estimates. Refer to the appendix for a formal comparison between the ACE and gDF model.
Errors i were assumed to follow  (0 , 2 ) . Significance of models and regression coefficients was assumed at an error level < 0.05. The number of equations in this model was equal to the number of subjects n , using double entry of pairs (i.e., switching indices 1,2 in Eq. 1 ). Standard errors and significance levels were based on a sandwich variances estimator. If a two-sided test C ≠ 0 was not significant, this factor was dropped, resorting to an AE model. Genetic heritability was assessed by one-sided tests of A > 0 ( Lazzeroni and Ray, 2013 ). Communitylevel results were corrected for multiple comparison using the Bonferroni method, with a factor of 12 for the left, and 13 for the right hemisphere. For vertex-wise measures y , resulting p -values were corrected for a false-discovery rate (FDR) of = 0.05 using the ( Benjamini and Hochberg, 1995 ). Generally, we found that heritability estimates A ≥ 0.1 were significant.
Unless explicitly noted, age-related effects were not significant for this cohort of young subjects with a narrow age range. Similarly, the influence of body height was insignificant for all relative measures y (e.g., Table 2 Influence of sex, body height and heritability A on brain-and hemisphere-level measures. The cohortlevel mean is specified in the third column; sex-related differences and effect size (Cohen's d ) are given as relative to females. For GMV, a significant age-related loss was found (-2.43 cm 3 /yr, p = 1.8e-7). Superscripts correspond to methods in Section 2.3 . Measures: (a) ICV = intracranial volume; (b) BRV = brain volume; (c) GMV = grey matter volume; (d) WMV = white matter volume; (e) BRR = ratio BRV/ICV; (f) AR = surface area of left (right) hemisphere. areas and distances on the unit sphere). Note that measures based on a smaller number of observations (e.g., regional vs hemispheric surface area) are expected to have a larger variance at the cohort-level. This larger variance is reflected as an increase in environmental effects E , with consequentially lower values for A and C . For regional measures, it is difficult to decide whether a relatively lower estimate of genetic effects is due to a higher variability or due to a limited precision of the underlying measure.

Model B -Statistics for similarity measures:
For the assessment of within-pair similarity measures, we used ordinary least squares (OLS) regression: where y i corresponds to the similarity measure of pair i , and x i,k to covariates x k for sex (same or different), body height (absolute difference, in cm), age (absolute difference, in years) and race (same or different), including an intercept x i ,1 = 1. Genetic influences were represented by the IBD regressor I A (as above), which corresponded to the genetic similarity within a pair. Thus, we assessed phenotypic similarity y i in terms of genetic similarity I A . The significance of genetic influences was reported as the p -value for the regression coefficient b a , termed heritability R here. Similarly, common environmental effects were modeled via regressor I C . Note that this coefficient is signed: negative for dis-similarity measures (e.g., pair-wise distances) and positive for similarity measures.
Here, the number of equations was equal to the number of subject pairs n /2. Community-level results were corrected for multiple comparison using the Bonferroni method, with a factor of 12 for the left, and 13 for the right hemisphere. Unless explicitly noted, coefficients for age-, height-, and race-related regressors and environmental effects were not significant.

Brain-and Hemisphere-Level
At the highest structural level, results ( Table 2 ) demonstrate a strong influence of genetic effects on the brain, compartment volumes, and hemisphere surface size. Heritability values A are at this level were published recently ( Le Guen et al., 2018a;Patel et al., 2018;Strike et al., 2019 ), other cohorts Lukies et al., 2017;Panizzon et al., 2009 ), along other anatomical features (e.g., body height in adults: Silventoinen et al., 2003 ). Males had between 8-10% larger compartment volumes and surface sizes than females, taking body height and race into account. Except for GM volume, age was not found to significantly influence these measures, presumably due to the relative narrow age range of this cohort of young adults.
Next, we assessed the within-pair similarity of cerebral hemispheres with several metrics: correlation in image space (ISIM), mean distance between ipsilateral surface vertices (SSIM), similarity of the basin segmentation (BLS), Dice overlap of sulcal areas (SA), mean distance between sulcal bottom lines (BLD), and similarity of whole-hemisphere basin graphs (BGS). Results, compiled in Table 3 , demonstrated a strong positive correlation between the IBD regressor I and similarity metrics (ISIM, BLS, SA, BGS), and a strong negative correlation for distance measures (SSIM, BLD), indicating a strong genetic influence on overall surface features. Simple similarity measures such as BLS and SA allow distinguishing MZ from DZ pairs with more than 90% correctness. Overall, heritability was not different between hemispheres. Sex and body height showed weak influences. We also checked for mirror-symmetric brains in twins (i.e., if the LH in member A is similar to the RH in member B: Steinmetz et al., 1995 ;Teplica and Peekna, 2005 ), but did not find striking examples.
While in Table 2 , volumes and areas were assessed for similarity without taking shape into account, Table 3 includes measures that examine voxel-or vertex-wise similarity based on stricter criteria. The relevance of high heritability estimates (e.g., for brain volume) could hence be considered arguable, because structural similarity is implied but actually not assessed.
We also assessed the correlation between measures of Table 2 and 3 (refer to supplementary Fig. S1 , top row). High correlations between volumetric measures were trivial and expected. Moderate correlations were determined between measures BLS, SA, BLD and BGS which all assess different aspects of basins and their local depth and curvature properties. Correlations between volumetric and surface measures were low, and often not significant.

Community-Level Results
At the regional level of cortical communities, we examined number of basins per community, community size, maximum geodesic depth within a community, fraction of vertices in sulcal areas, mean cortical thickness, and mean myelin fraction ( Table 4 ).
We found the strongest influence of genetic effects for cortical thickness ( ~0.60) and myelin fraction ( ~0.45), similar across hemispheres and communities. Weaker influences were determined for the size of sulcal areas ( ~0.40), maximum depth ( ~0.37), basin count ( ~0.30), and community size ( ~0.20). Here influences were similar across hemispheres, but varied by community. Relatively strong genetic influences were found consistently for all parameters in communities CS and IN, and TL on both sides.
Following, within-pair similarity of communities was assessed via: Dice overlap of sulcal areas (SA), similarity of the basin labels (BLS), mean distance between sulcal bottom lines (BLD), and similarity of community-wise basin graphs (BGS), and depth and position of community centers. By construction, all six features represent different aspects of the local structural configuration. Compiled results in Table 5  Table 3 Influence of sex, body height and heritability R on hemisphere-level similarity measures. Superscripts correspond to methods in Section 2.3 : (a) ISIM = correlation coefficient of left (right) hemisphere in image space; (b) SSIM = mean vertex-wise distance between surface vertices of the left (right) hemisphere; (c) BLS = normalized mutual information of basin labels; (d) SA = Dice overlap coefficient of sulcal areas; (e) BLD = mean distance between sulcal bottom lines; (f) BGS = similarity of whole-hemisphere basin graphs. Means were provided in column 3 so that regression coefficients in columns 4, 7, and 9 can be understood in terms of their quantitative influence. Sex-related differences and effect size (Cohen's d ) are given as relative to females.  revealed strong correlation between IBD rating I and similarity of sulcal areas, basin graphs and basin labels, while the distance between bottom lines, and the difference in depth and position of community centers had less genetic influence. Interestingly, the similarity of sulcal areas (SA) is always greater than whole basins (BLD), confirming that deeper areas are more invariant. Likewise, the center depth is more invariant than the center position. In contrast, the similarity of bottom lines may retain too little information about local structural properties. Again, the amount of correlation differed between communities. On the left side, communities CS, FI, TL, and OM showed the strongest effects, on the right side, OM, CS, MP, and TM were most prominent. Overall, heritability was not different between hemispheres.
We also assessed the correlation between measures of Table 4 and 5 (refer to supplementary Fig. S1 , bottom row). Again, moderate correlations were determined between measures that quantify different aspects of basins and their local depth and curvature properties (i.e., basin count and size, maximum depth and sulcal fractions). Correlations between all other measures were low or not significant.

Results at the Vertex-Level
Finally, vertex-wise estimates of geodesic depth, cortical thickness, surface curvature, and myelin fraction were assessed for additive genetic effects as represented by heritability A in model A. Results are shown for depth in Fig. 3 , and for the other measures in supplementary Figs. S2, S3, and S4 . To ease comparison across figures, the same heritability range 0-0.8 was color-coded from magenta to red. Overall, the strongest heritability was found for geodesic depth (0.3-0.8), while myelin fraction and cortical thickness had moderate influences (0.2-0.6). To provide a spatially continuous representation, we did not threshold for statistical significance. Based on the Benjamini-Hochberg procedure, values above 0.1 were addressed as significant.
For geodesic depth, the strongest heritability (~0.9) was found in the insula and central sulcus on both sides, and along the right superior and middle temporal sulcus. Regions with moderate heritability (0.3-0.5) include the inferior prefrontal sulcus, the intra-parietal and lunate sulci, the calcarine and parieto-occipital sulci, and the rostral segment of the calcarine sulcus, similar in extent and magnitude on both sides. A moderate heritability (~0.4) was determined for the myelin fraction, which was rather uniform on the cortical convexity, but less prominent in midline cortices. Similarly, moderate effects were found for cortical thickness, especially in the insula and central sulcus on both sides.
Note that the analyses above assumed that vertex-wise properties were sampled at homologue cortical positions across all hemispheres. Because linear regression was the sole method to correct for global rotation, we employed a second method to reduce inter-subject variability. Germanaud et al. (2012) demonstrated that observable increasing detail of cortical structures during fetal brain development can be quantified by an increase in high frequency components in the spectral domain. Thus, we removed some of the inter-subject variability induced by developmental processes via low-pass filtering features using wavelet banks ( Yeo et al., 2008 ). Results are shown in Fig. 4 for the left hemisphere, where the left column shows transformed results, the right column results at the original resolution. Rows correspond to features geodesic depth, curvature, thickness, and myelin fraction, resp. The reduction of individual variability from filtering led to an overall increase in heritability estimates. Notably, the spatial pattern for depth and curvature remained similar and regionally non-uniform across scales, while patterns for thickness and myelin fraction became regionally more uniform. We argue that a certain proportion of heritability is "covered " by individual processes modulating cortical detail in the previous analysis.
Summarizing, macro-structural features geodesic depth and curvature showed a similar anisotropic pattern, with a strong heritability in an anterior-posterior direction along the insula -parieto-occipital sulcus, and a superior-inferior direction along the central sulcus and superior portion of the callosal sulcus. In comparison, micro-structural features showed a rather uniform pattern of heritability across the cerebral cortex.

Sex-Related Differences
Sex proved a significant influence on most structural measures studied. Males had about 9.6% larger brain volumes and 7.9% larger surface sizes ( Table 2 ), corresponding to large effects (Cohen's d > 0.8), taking body height, age and race into account. For a deeper insight, we assessed results for regressors sex and body height on the hemisphere-and community-level for surface properties in Table 4 . Results were compiled in supplementary Tables S1-S6 . Due to the larger surfaces, males had more basins than females (left: +7.56, p = 2.16e-8, d = 0.482; right: +8.12, p = 3.57e-10, d = 0.527). Similar differences were determined at the community level. Effect sizes were smaller, due to a considerable variability in the individual community structure. Additional basins in males were small ( "dimples ") and rather allocated in "remaining space " on the slightly larger surface. Males had deeper folds (left: +1.31 mm, p = 2.35e-5, d = 0.373, right: +0.858 mm, p = 4.87e-2, d = 0.175), corresponding to larger cortical areas in folds ( "sulcal fraction " left: +0.671%, p = 1.84e-11, d = 0.530; right: +0.480%, p = 5.17e-7, d = 0.394). The myelin fraction was slightly higher in males (left: +0.0290, p = 1.688e-4, d = 0.327; right: +0.0278, p = 4.42e-4, d = Fig. 3. Vertex-wise mapping of estimated heritability A for geodesic depth, as determined by model A. The left and right hemispheres are viewed from the top (row A), side (row B), bottom (row C), and mid-sagittal plane (row D). Strong heritability of geodesic depth was found along an anterior-posterior axis (insulaparieto-occipital sulcus) and a superior-inferior axis (central sulcus and superior portion of the callosal sulcus). Note that regions with A < 0.1 (in magenta) were included to preserve the visual continuity, but were not statistically significant. 0.305). Effects were approximately similar across communities. In contrast, the cortical layer was about 1% thinner in males (left: -0.0229 mm, p = 4.38e-2, d = -0.162; right: -0.0244 mm, p = 3.84e-2, d = -0.172), with regionally stronger differences in communities PA and MP on both sides. Similar values and regions were found in the larger UK Biobank cohort ( Elliott et al., 2018 ;Ritchie et al., 2018 ). Community size was generally independent of sex. Body height showed a positive correlation with all measures except community size and cortical thickness, where no significant dependence was found. Effects were generally an order of magnitude smaller than sex-related differences, and thus, less often statistically significant.

Asymmetries, Handedness and Heritability
For community-level properties assessed by methods 5-10, we computed a laterality measure = ( − )∕( + ) , where and denote the result on the corresponding side. This measure ranges between -1 (completely right lateralized) and 1 (completely left lateralized). Laterality measures were assessed by models B and C, using sex, handedness, and race as covariates. For results, refer to supplementary  Tables S7 and S8 . While brain-level asymmetries were comparatively small, strong structural asymmetries were found for most measures and communities. Consider communities TL and FI in Table S7 : regions were larger on the left side (basin count and size), shallower on the left (max. depth and sulcal fraction), and had a higher myelin fraction for FI, and a lower myelin fraction for TL. Continuing with Table S8 , basin labels and graphs are more similar on the left side, but sulci are deeper on the right side for TL, and on the left side for FI. The opposite pattern of a right-ward asymmetry is found most prominent in communities OL, OM, and FS.

Discussion
This study aimed at assessing the genetic influences on structural properties in the human cerebral cortex. Using a comprehensive measurement battery, a broad range of cortical features were examined at different spatial scales. Novel aspects of this study include: Fig. 4. Vertex-wise mapping of heritability A for geodesic depth, curvature, cortical thickness, and myelin fraction, as determined by model A. Wavelet-filtered analysis at level 2 (left column) vs level 5 (right column). Note that the spatial pattern for depth and curvature remained stable across levels, while the pattern for thickness and myelin became more uniform and considerably stronger at level 2. Note that regions with A < 0.1 (in magenta) were included to preserve the visual continuity, but were not statistically significant.
(1) Inclusion of all univariate features accessible from anatomical MRI (intensity, cortical thickness, geodesic depth, surface curvature and myelin fraction) in a common analytic framework. Results provide an unprecedented view on heritability of these features across the cortical mantle. Of note, we report the strongest heritability estimates reported so far for geodesic depth. (2) As our battery included multivariate measures that assess structural properties of the cortical surface, it provides a better understanding how genetic control shapes structural patterns and segregates regions of the cortex. That is, we quantitatively confirm the hypothesis that deeper regions (i.e., sulcal areas) are under stronger genetic control than gyral crowns. (3) Similar features were assessed at different scales. For all measures, features of (a) MZ twins were more similar than those of (b) DZ twins and siblings, and (c) those of unrelated subjects. Thus, genetic effects significantly influence structural features assessed at all spatial resolutions. However, while strong effects were found for gross compartment volumes, effects tend to weaken at finer resolutions. (4) Regional hemispheric measures provided a straightforward path for analyzing the strong structural asymmetries for heritability. In comparison with the strong heritability of features above, additive genetic effects on laterality were detectable, but considerably smaller. (5) Finally, we consider the consistency of results across different scales as an internal validation of the robustness of the measures provided in this study.
Based on the above, it could be of great interest to sketch a common framework of genetic influences on structural patterning of the human cortex, integrating findings with the current knowledge about its spatiotemporal development.

Genetic Influences on Cortical Properties
The univariate cortical properties (thickness, geodesic depth, curvature and myelin fraction) that were evaluated at different spatial scales, yielded a consistent pattern of genetic control.
Measures related to cortical thickness and myelin fraction showed a moderate to high heritability across spatial scales that was homogeneous throughout the cortex, confirming earlier reports for thickness Jha et al., 2018 ;Joshi et al., 2011 ;Panizzon et al., 2009 ;Strike et al., 2019 ) and myelin fraction Schmitt et al., 2019 ). The cortex of adults shows regional differences in thickness ( Lyall et al., 2016 ) and myelin content ( Grydeland et al., 2013 ) that differentiate them from cortical development in early childhood. In contrast, geodesic depth and curvature showed a region-dependent pattern with genetic control that varied from weak to strong. Geodesic depth showed the strongest heritability at the vertex-level analysis. We suggest that genetic control acts in two main directions: anterior-posterior along the long axis of the insula (IN) to the parietooccipital sulcus (POS); superior-inferior along the central sulcus (CS) to the superior portion of the callosal sulcus (SCS). Needless to say, primary cortices are allocated along these axes.
As demonstrated earlier ( Kruggel and Solodkin, 2019 ), the depth of sulcal centers can be regarded as a proxy for the time point at which they appear in fetal development: IN and POS emerge typically at gestation week (GW) 18, while CS and SCS do so around GW 20 ( Nishikuni, 2006;Nishikuni and Ribas, 2012 ). Our evidence suggests that, after the hemispheres separate into left and right (GW 5), the next development step segregates the cortical mantle into a superior and an inferior portion, followed by the sub-division of the superior portion into an anterior and a posterior sections. While these basic steps are under strong genetic control, any subsequent changes and allocation of cortical space occurs within the primordial sub-divisions, which could explain the regional inhomogeneity of depth and curvature maps. Comparing the depth and position of invariant centers in communities, genetic control appeared to be stronger for depth than position. Hence, using depth as proxy for time, we state the hypothesis that the spatio-temporal process governing the cortical gyrification is likely controlled by a genetic "clock ", with a spatial pattern as outcome ( Chou et al., 2016 ).
In summary, our results point toward differential mechanisms for the control of cortical thickness and myelin fraction on the one hand and geodesic depth and curvature on the other. While control over the former was spatially uniform, control over the latter differed between regions. This suggests that whereas one process controls the fine-structural properties of the cortical layer, the other controls its regional allocation, presumably as a temporal sequence of sub-divisions into hemispheres, superior and inferior, anterior and posterior quadrants.

Genetic Influences on Structural Features
In accordance to previous studies showing that gyral crowns are structurally more variable than sulcal bottoms ( Kruggel, 2018 ), we found that our structural measures assessing deep areas are under stronger genetic control than measures associated with gyral crowns. Regional differences in heritability of these structural features parallel those described for geodesic depth and curvature above. Using depth as proxy for time, sulcal regions are segregated early in development and remain relatively constant, while gyral regions add detail even after birth ( Lenroot et al., 2009 ;Lyall et al., 2016 ;Im and Grant, 2019 ;Teeuw et al., 2019 ). It is interesting to note that our map of cortical communities is largely similar to a map of 12 regions defined by the clustering of genetic correlations ( Chen et al., 2013 ).
We assessed the heritability of some features at the hemispheric and regional level (BLS, SA, BLD, and BGS). Consistently, stronger effects were found for the hemispheric than for the regional analysis. It is likely that epigenetic and developmental influences lead to an individual variation in regional size and shape, especially in cortical regions allocated later in development, thus, lowering local pair-wise similarity measures. However, the much stronger similarities in MZ vs. DZ twins for wholehemisphere measures indicated that MZ twins are more similar in the combination of all regions, underlining that the common genetic control in MZ twins is present throughout the whole cortex.
In contrast, the pair-wise similarity in position and depth of invariant centers depended much less on the assumed proportion of genes identical by descent. We assume that the allocation of centers (in terms of their position and time point during gestation) is under genetic control common to the whole cohort, thus, similar for groups MZ, DZ, SB, and UN.
Strong sex-related differences were found for all structural measures. Approximating a hemisphere by a half-sphere, the 10% larger brain volume in males would only correspond to a 5.6% increase in surface size.
Thus, the 9% larger surface size must be accommodated by deeper folds in males, which we found in this cohort. The relatively larger interface between GM and CSF is also reflected in the over-proportionally larger CSF volume in males, and the higher brain ratio (BRR) in females. Overall, cortical thickness was 1% smaller in males, likely related to the finding that males had a larger proportion of (thinner) cortex in folds. The larger surface area on the convexity of male brains was apparently addressed to a slightly larger number of small basins ( "dimples "), with their neuro-biological significance undefined at this time. Similar in regional pattern and extent, we found strong structural asymmetries ( Geschwind et al., 2002;Habas et al., 2012;Im et al., 2010;Steinmetz et al., 1995 ) to the left in communities OR, FI, and TL, and to the right in communities FS, OL, and OM. In spite of their strong laterality, statistical models provided comparatively little evidence for genetic influence suggesting that the variance could rather be addressed to environmental factors. Furthermore, we did not encounter clear evidence for mirror-symmetric twin pairs in this cohort.

Mechanisms of Cortical Folding
Although general mechanisms of early cortical development are largely similar in mammals (reviews: Bystron et al., 2008 ;Silbereis et al., 2016 ), distinct differences in genetic patterns of later development stages are determined between lis-encephalic (e.g., mice) and gyr-encephalic species (e.g., ferret, humans). While the former show graded patterns of transcription factors (e.g., Emx2, Pax6, Coup-TFI, and Sp8) along the anterior-posterior and medio-lateral directions ( O'Leary et al., 2007 ), the latter show modulated expression levels (e.g., Cdh1 and Trnp1) concurrent with prospective locations of gyri and sulci ( De Juan et al., 2015;De Juan and Borrell, 2017;Martinez-Martinez et al., 2016 ). Especially the outer sub-ventricular zone (oSVG), which is very prominent in humans, contains outer radial glia (oRG) that extend very long processes to the pial surface. These cells produce neurons in layers II and III, increasing the cortical thickness. It was proposed that patterned differences in neurogenesis lead to a differential expansion of the cortical layer, and thus, to gyri and sulci ( Borrell, 2018 ). In-vitro and in-silico modeling ( Bayly et al., 2013;Ronan et al., 2014;Tallinen et al., 2016;Toro and Burnod, 2005 ) demonstrated that differential expansion (thus, different mechanical properties) of a cortical layer over a (homogeneous) white matter core result in folding reminiscent of the human brain. We found that the cortex was on average 0.6 mm thicker in gyral crowns than sulcal bottoms ( ≈ 20%) which may support the hypothesis of differential expansion of the cortical layer. In contrast, our results indicate that genetic control over cortical features thickness and myelin fraction is spatially rather uniform, which points to a uniform mechanism for cortical development that nonetheless allows for regional variation, i.e., by modulated expression levels of transcription factors.
However, differential expansion of the cortical layer fails to explain the consistent location of primary gyri and sulci across individuals, and the growth of the sub-cortical white matter. A more recent study by Rash et al. (2019) found that neurogenesis in oSVG is followed by gliogenesis before gyri develop, indicating that a patterned growth of cortico-cortical connections induces gyrification. In addition, there are reciprocal structural and functional relationships between primary cortices and thalamic nuclei: the ventro-lateral (VL) nucleus with the motor cortex, the ventro-posterior (VP) nucleus with the somatosensory cortex, the dorsal lateral geniculate nucleus (dLGN) with the visual cortex, and the ventral part of the medial geniculate nucleus (MGv) with the auditory cortex ( O'Leary et al., 2007 ). These thalamo-cortical axons grow into the cortex during GW 20-24 ( Reillo et al., 2011;Silbereis et al., 2016 ), a period during which the corresponding sulci develop. In contrast to the uniform genetic influence on cortical measures, our heritability maps for shape-related features geodesic depth and curvature are spatially non-uniform, especially in the vicinity of primary cortices. This indicates a second mechanism that allocates cortical space via the formation of subcortical connections which is under strong genetic control in primary regions that develop earlier, compared with a weaker control in secondary areas. Notably, regional transcriptional profiles were largely bilaterally symmetrical across the two hemispheres ( Nowakowski et al., 2017;Pletikos et al., 2014 ), in contrast to the observation that macroscopic structural asymmetries were detectable during gestation ( Chi et al., 1977 ). This indicates either that any genetic mechanisms are currently undetected or that non-transcriptional mechanisms play a critical role ( Silbereis et al., 2016 ). Our morphometric study is in line with this conclusion: While we found strong lateralization in several structural measures, statistical modeling yielded only weak genetic influences.
Summarizing, it is likely that two mechanisms govern cortical pattern formation: (1) a gyrification process, driven by differential expansion and the formation of cortico-cortical connections, (2) a spatial allocation process, guided by the formation of sub-cortical connections, that introduces regularity and functional specialization to pattern formation.
While we focused on the study of structural features of cortex, we deliberately disregarded an additional mechanism that shapes the cortex: the development of cortico-cortical connectivity and additional projection fibers that likely contribute to the individual embodiment of gyral crowns. More work is necessary to understand the shape-forming properties of short-and medium-range U-fibers, and their genetic determinants . A deeper understanding of bio-mechanic and genetic constraints governing pattern formation may lead to (simplified) parametric models of sulcal shape, and in turn, to a better definition of homologue locations across cortices. In spite of the fact that large similarities between brains exist, we begin to recognize that individual structural detail is also relevant. Attention to this detail will help to better understand the implementation of functional processes in the human brain, and their variation in development during life.
In contrast to the ACE model, coefficients A, C estimated by gDF regression can be negative. However, A is non-negative and identical to the classical model if the within-pair variance of DZ twins is larger than those of MZ twins ( ar ( , 1 − , 2 ) ≥ ar ( , 1 − , 2 ) ). Likewise, C is non-negative and identical to the classical model if A is non-negative and the within-pair variance of UN pairs is larger than those of DZ twins ( ar ≥ ar ). We compared results between our implementation and the classical ACE model on simulated and real data from our study. Given the constraints above, estimates for A, C never differed by more than 5%.

Supplementary material
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.neuroimage.2020.117169 .