Biophysical modeling of anatomically realistic prenatal cortical folding development

TL;DR

This paper introduces a new biophysical model that realistically simulates how the human fetal brain develops its complex folded structure during gestation.

Contribution

The novel framework integrates data-driven growth laws with anatomically accurate geometry to model cortical folding with unprecedented realism.

Findings

Region-specific growth fields derived from prenatal MRI data produce realistic cortical folding patterns.

The model replicates clinically atypical brain phenotypes like lissencephaly and polymicrogyria through controlled perturbations.

The framework enables high-fidelity synthetic brain datasets for AI-driven developmental neuroscience.

Abstract

Cortical folds encode the architecture of human cognition, yet the mechanisms that transform the smooth fetal cortex into its convoluted geometry remain elusive. Biophysical modeling enables mechanistic insight into cortical morphogenesis, but existing models often lack anatomical realism and fail to capture key hallmarks and morphometrics of dynamic cortical folding in the developing human brain. Here, we introduce a novel whole-brain developmental framework that integrates region-specific, data-driven growth laws with anatomically accurate cortical geometry to enable realistic and biologically interpretable modeling of cortical morphogenesis during gestation. Growth fields derived from large-scale prenatal magnetic resonance imaging data capture spatiotemporal variations in cortical expansion and thickness across parcellated regions. Incorporating this heterogeneous growth yields anatomically faithful folding patterns that closely match qualitative landmarks and quantitative morphometrics from human imaging. Systematic perturbations of geometry and growth attributes delineate control parameters that produce realistic morphological variability and replicate clinically atypical brain phenotypes consistent with lissencephaly, pachygyria, and polymicrogyria. This framework provides a quantitative foundation for elucidating the mechanisms of typical and atypical fetal brain development and can serve as a promising generative engine for high-fidelity, longitudinal synthetic brain datasets to advance AI-driven developmental neuroscience and clinical translation.

Figures

d](#F2)). To further characterize these differences, we summarized the regional growth ratios at mid-gestation (31 GWs) and at term (40 GWs) based on the identified growth models ([Fig. 2e](#F2)). Tangential growth ratios are consistently higher than their radial counterparts, and their near-linear trajectories preserve the overall spatial distribution from mid-gestation to term. Notably, the temporal, parietal, and peri-Sylvian regions exhibit relatively higher tangential growth ratios. This observation aligns with neuroimaging evidence showing accelerated prenatal development in areas implic

c](#F3)). In the modeled brains, the relative positions of major folds stabilize by approximately 32 GWs; beyond this stage, the sulci primarily elongate and widen, while secondary and tertiary branches gradually emerge to further enrich folding complexity. By contrast, sulcal configurations in real brains appear less consistent across gestational ages, reflecting inter-individual variability inherent in cross-sectional datasets. Importantly, the modeling framework allows continuous tracking of folding evolution, generating brain morphologies at any developmental stage (Movie 1). This capabili

d](#F4)). In both modeled and developing brains, vertex-wise shape index distributions consistently show a characteristic double-peak profile centered around ± 0.5. This reflects the increasing prominence of ridge and rut features, consistent with previous reports^[38](#R38),[39](#R39)^ ([Fig. 4e](#F4)). Marked shifts in the distribution are observed after the initial stage at 22 GWs, coinciding with the onset of surface buckling that initiates fold formation. As development proceeds, the distributions gradually converge in parallel with the maturation of major folds during late gestation. To

d](#F5), difference in the sulcus counts among five growth conditions (Real, OGR, TGR, TGU, and IGU; n = 50, 10 per group) were first evaluated using the Kruskal–Wallis test to evaluate the null hypothesis that group medians were equal. When the overall test was significant (*P* < 0.05), post hoc pairwise comparisons were carried out using Dunn’s test with Bonferroni correction for multiple testing. In Fig. S5e, paired comparisons of sulcal counts between left and right hemispheres (n = 20, 10 per hemisphere) were assessed using the Wilcoxon signed-rank test. The null hypothesis was that the m

b](#F6)). Two levels of offset were introduced: modest adjustments (ICT1, increased thickness; RCT1, reduced thickness) and larger adjustments (ICT2 and RCT2), corresponding to 40% and 80% of the mean initial cortical thickness. Increasing the relative thickness produces fewer folds, characterized by bulkier gyri and deeper sulci, whereas decreasing it generates more densely distributed folds with smaller gyri and shallower sulci. This observation aligns with core–shell buckling theory, which reveals that greater shell thickness leads to less complex folding patterns^[41](#R41)^. To quantify t

h](#F7)). Relative deviations again reveal that regional alterations also propagate to neighboring areas ([Fig. 7i](#F7)). In summary, these analyses indicate that growth ratio magnitude predominantly governs the formation of major folding patterns, while growth trajectories modulate localized folding features, as reflected by local quantitative deviations. On this basis, we investigated hemispherical symmetry in folding patterns, as the two hemispheres can differ in both geometry and growth profiles even for the same subject^[42](#R42)^(Fig. S5a). The modeled brain showed asymmetric folding p

d](#F8)). To replicate this phenotype, we exclusively reduced the tangential growth ratios by half in the affected regions (R1: dorsal precentral & posterior cingulate, R4: lateral precentral & postcentral, and R11: caudal middle frontal) (Movie 5). Conversely, the regional polymicrogyria case, which exhibits abnormally dense but shallow folds in the right frontal lobe ([Fig. 8e](#F8)), was recreated by reducing initial cortical thickness by 40% in selected frontal regions while keeping all other parameters unchanged (Movie 6). Another polymicrogyria case, accompanied by volume atrophy in the

Full text

Introduction

The human brain undergoes remarkable morphological changes during gestation, giving rise to the highly convoluted folds of the cerebral cortex^1^. This folding process, which intensifies during the third trimester, coincides with rapid volumetric expansion and surface growth that foster neuronal proliferation and promote efficient information processing. Alongside the emergence of intricate folding morphology, cortical development exhibits pronounced spatiotemporal heterogeneity^2^. Cross-sectional neuroimaging analyses reveal distinct regional trajectories of cortical thickness and surface expansion across gestational ages^3,4^. This regionalization is further supported by transcriptomic and cellular evidence showing that spatially patterned gene expression regulates neural progenitors proliferate, differentiate, and form cortical layers, leading to divergent growth dynamics across cortical regions^5–7^ (Fig. S1). However, despite increasing evidence of heterogeneous growth signatures, how they physically translate into the complex folding patterns of the cortex remains poorly understood.

Cortical folding displays both consistent and variable features across individuals. Primary folds arise at stereotypical locations around mid-gestation, forming early landmarks such as the central sulcus and Sylvian fissure that are conserved across individuals. In contrast, secondary and tertiary folds emerge later, enriching local surface complexity and contributing to inter-individual variability^8,9^. This coexistence of regularity and variability points to a potential interaction between regional growth dynamics and geometric constraints^10^, yet their specific roles in shaping cortical folding patterns remain unclear. When normal developmental processes are disrupted, the resulting cortical malformations, characterized as atypical folding geometries, are clinically associated with neurodevelopmental disorders such as epilepsy, schizophrenia and autism spectrum disorder^11–13^. Clarifying the developmental mechanisms underlying these malformations is therefore crucial for uncovering their pathological origins and for informing early diagnosis and risk assessment.

With advances in computational and biophysical techniques, in silico modeling has emerged as powerful approach to simulating the complex spatiotemporal development of cortical folding within a physically grounded framework. To this end, various hypotheses and theories have been proposed to incorporate distinct biological mechanisms^10,14–18^. Among them, the prevailing explanation is the differential tangential growth hypothesis, which posits that deformation mismatch between the expanding cortical plate and the underlying white matter drives cortical gyrification^19–23^. While this paradigm has yielded important insight into folding mechanics, current models often yield simplified or generic patterns that fail to capture the anatomical complexity observed in human fetal magnetic resonance imaging (MRI) data. This discrepancy largely arises from the use of conceptual geometrical models and oversimplified growth representations. To approximate brain geometry, many models adopt generalized shapes such as spheres or ellipsoids^24–27^, which capture the brain gross morphology but lack regional anatomical details. Although some studies have employed reconstructed brain surfaces as simulation domains, they generally apply globally uniform growth assumptions derived from bulk volume approximations, which fail to account for region-specific developmental dynamics^10,22,23^. This raises a key question of whether uniting anatomically realistic brain geometries with biologically grounded growth representations can more faithfully capture cortical morphogenesis.

In this study, we present an advanced whole-brain computational modeling framework that integrates heterogeneous developmental growth laws to simulate the intricate prenatal cortical development with anatomically realistic folding patterns. Built upon prenatal brain geometries at 22 gestational weeks (GWs), the framework incorporates, for the first time, atlas-based regional growth laws derived from large-scale imaging-based measurements of cortical expansion and thickness spanning 22–40 GWs. We show that heterogeneous growth generates folding patterns with high anatomical fidelity, closely matching both cortical landmark features and quantitative morphological measures. Systematic variation of geometric factors, including the initial cortical surface configuration and cortical thickness, together with growth factors such as trajectory and magnitude, further demonstrates how their interplay shapes distinct folding morphologies. Beyond accurately modeling cortical development, our framework provides a generative platform for generating numerous anatomically realistic synthetic fetal brain datasets, overcoming the scarcity of longitudinal prenatal imaging. The framework also sheds insight into the mechanistic, evolutionary, and clinical aspects of brain morphogenesis while providing a scalable foundation for developing high-fidelity synthetic brain datasets that advance AI-based neuroscience research.

Results

Growth trajectories are heterogeneous across brain regions

We first quantified regional cortical thickness and surface area in 120 fetal brains (22–40 GWs) curated from the developing Human Connectome Project (dHCP), using a developmental atlas specifically designed to parcellate the prenatal and neonatal cortex (Figs. 1a–d, 2a). Across 22–40 gestational weeks, these regions followed distinct morphogenetic trajectories, reflected by divergent patterns of cortical thickness and surface area (Figs. 1c–d). To facilitate cross-regional comparison, all measurements were normalized to their earliest-stage mean values (22 GWs), yielding relative growth ratios. Using a data-driven approach based on symbolic regression (Fig. S2), we derived mathematic models that capture regional growth dynamics, distinguishing radial growth (changes in cortical thickness; Fig. 2b) from tangential growth (changes in surface area; Fig. 2c).

Both radial and tangential growth exhibit marked heterogeneity in their developmental trajectories. Cortical thickness either increases steadily (e.g., region 1: dorsal precentral & posterior cingulate), reaches an early plateau (region 10: dorsal superior frontal), or follows a biphasic course. In the latter pattern, thickness peaks around mid-gestation (~ 31 GWs) and subsequently declines modestly (region 18: lateral occipital) or stabilizes at a plateau (region 17: middle and inferior temporal). These characteristic trends are consistent with previous fetal MRI observations^28–30^. The non-monotonic decline in thickness likely reflects localized microstructural remodeling processes, including synaptic pruning and dendritic compaction, that reduce apparent cortical thickness^31^. In contrast, tangential growth consistently follows linear trajectory, with rates substantially exceeding those of radial growth. This rapid tangential expansion during the prenatal period drives pronounced cortical enlargement, which in turn induces the emergence of folds. The resulting increase in cortical surface area accommodates proliferating neurons and establishes the structural foundation for nascent brain function^5^.

These regional distinctions demonstrate that cortical development varies across both space and time, resulting in spatiotemporal heterogeneity across brain regions (Fig. 2d). To further characterize these differences, we summarized the regional growth ratios at mid-gestation (31 GWs) and at term (40 GWs) based on the identified growth models (Fig. 2e). Tangential growth ratios are consistently higher than their radial counterparts, and their near-linear trajectories preserve the overall spatial distribution from mid-gestation to term. Notably, the temporal, parietal, and peri-Sylvian regions exhibit relatively higher tangential growth ratios. This observation aligns with neuroimaging evidence showing accelerated prenatal development in areas implicated in auditory and language functions^6,32,33^. Radial growth displays a distinct pattern: cortical thickness increases more rapidly around the central sulcus regions during late gestation, in line with the earlier maturation of primary sensorimotor cortices that support motor control and somatosensory processing^34^. Growth patterns are also not entirely symmetrical between hemispheres (Fig. S3). Although the exact trajectories differ, both hemispheres follow broadly similar trends, especially in tangential growth. Collectively, these findings reveal the pronounced spatiotemporal heterogeneity in cortical development, reflecting the inherent complexity underlying the emergence of human brain folding.

Modeled folding patterns align with those observed in the developing brain

To examine how heterogeneous growth shapes cortical folding, we employed an in-silico modeling approach. A whole-brain computational framework was constructed with anatomically realistic geometry (Fig. S4), into which the identified regional growth models were incorporated as deformation drivers to recapitulate cortical folding development across mid-to-late gestation (22–40 GWs). Starting from an initial smooth geometry, the model progressively generated convoluted folds that conform to realistic brain morphology (Fig. 3).

The resulting brains closely resemble developing human brains not only in overall volumetric growth but, more importantly, in the emergence and evolution of major folding patterns across corresponding developmental stages (Fig. 3a). We note that reference brains shown here represent cross-sectional samples from different individuals and serve as developmental benchmarks rather than direct subject-level comparisons. The modeled trajectories indicate that primary folding patterns are largely established by ~ 32GWs, after which existing sulci mainly deepen and widen in concert with overall brain expansion. This developmental timeline aligns with in utero imaging observations showing that primary folds emerge by 32–34 GWs, while subsequent complexity arises through secondary and tertiary fold elaboration^35^. Notably, the framework reproduced key anatomical folding features, accurately capturing six major sulcal landmarks in both location and relative dimension (Fig. 3b). These sulci are well-recognized markers of cortical organization due to their structural and functional relevance. For example, the central sulcus demarcates the boundary between primary motor and somatosensory cortices, which are essential regions for motor control and sensory perception^36^, whereas the superior temporal sulcus represents a prominent landmark associated with higher-order auditory and language processing, as well as social cognition^37^.

To visualize and compare cortical folding patterns, we projected the brain morphologies onto a unit sphere with overlaid major sulcal contours (Fig. 3c). In the modeled brains, the relative positions of major folds stabilize by approximately 32 GWs; beyond this stage, the sulci primarily elongate and widen, while secondary and tertiary branches gradually emerge to further enrich folding complexity. By contrast, sulcal configurations in real brains appear less consistent across gestational ages, reflecting inter-individual variability inherent in cross-sectional datasets. Importantly, the modeling framework allows continuous tracking of folding evolution, generating brain morphologies at any developmental stage (Movie 1). This capability underscores a key strength of computational modeling in reconstructing the full longitudinal development of cortical folding from a single MRI scan.

To evaluate the accuracy of folding reproduction, we quantified multiple morphometrical parameters and compared results across global, regional, and vertex scales (Fig. 4). At the whole-brain level, folding complexity was assessed using the global gyrification index (GI), which measures the cortical surface convolution relative to the enclosing convex hull. During development, the global GI exhibits an accelerated near-linear increase after ~ 27 GWs, marking the onset of rapid fold formation (Fig. 3). Moreover, the global GI trajectories of the modeled brains closely mirror those of developing brains, reflecting consistent global folding dynamics across development (Fig. 4a). At the regional level, we evaluated the mean local GI and sulcal depth. The local GI refines the global metric by quantifying folding complexity within specific cortical territories, whereas sulcal depth directly measures the depth of individual folds (Fig. 4b). At term, local GI peaked prominently in the insular (region 2) and isthmus cingulate (region 3) areas, consistently with their persistent deep concavities throughout development. Overall, regional local GI values in modeled brains closely match those in developing brains, with relatively larger deviations in the peri-Sylvian and isthmus cingulate regions. Such discrepancies likely stem from geometric complexity that complicates both modeling and morphometric quantification. We further tracked regional sulcal depth trajectories across gestation for both modeled and real brains (Fig. 4c). Sulcal depth exhibited a biphasic growth pattern, characterized by an initial slow increase followed by accelerated deepening, with a transition occurring near 29 GWs (t ≈ 0.4). In most cortical regions, the modeled trajectories faithfully captured both the trend and magnitude of sulcal deepening observed in developing brains.

Subsequently, we evaluated folding accuracy at the vertex level using a curvature-based measure, the shape index. The shape index encodes local surface geometry, with positive values corresponding to gyri and negative values to sulci (Fig. 4d). In both modeled and developing brains, vertex-wise shape index distributions consistently show a characteristic double-peak profile centered around ± 0.5. This reflects the increasing prominence of ridge and rut features, consistent with previous reports^38,39^ (Fig. 4e). Marked shifts in the distribution are observed after the initial stage at 22 GWs, coinciding with the onset of surface buckling that initiates fold formation. As development proceeds, the distributions gradually converge in parallel with the maturation of major folds during late gestation. To quantify vertex-wise similarity, we computed a similarity index between each modeled brain and the average of ten real brains at the corresponding gestational stage (Fig. 4f). Across all stages, similarity indices exceed 0.8 for all brain subjects, indicating high vertex-level fidelity in the modeled brains. Together, these qualitative and quantitative comparisons demonstrate that the framework reliably generates anatomically realistic cortical folding patterns and serves as a powerful complement to longitudinal imaging datasets.

Heterogeneous growth produces more realistic folding outcomes

While the imaging-informed framework reproduces folding patterns with anatomical fidelity, the specific contribution of regional growth heterogeneity to this process remains to be systematically elucidated. To isolate this effect, we examined alternative growth modes by introducing isotropic (same in all directions) and tangential (restricted to tangential planes) formulation in addition to the orthotropic growth assumed in the cortical gray matter (Fig. 5a). Two growth profiles, uniform and regionally varying, were further applied across the cortical gray matter, yielding four modeling cases: isotropic growth with a uniform profile (IGU), tangential growth with uniform (TGU) or regional (TGR) profiles, and orthotropic growth with a regional profile (OGR). Among these, IGU and TGU represent the simplified growth assumptions commonly adopted in previous studies.

These growth conditions produce distinct folding outcomes (Fig. 5b). Under isotropic growth, the enhanced radial expansion yields less complex morphologies, characterized by bulky gyri and widely spaced sulci. Tangential and orthotropic growth generate broadly similar folding patterns, with TGR and OGR producing nearly indistinguishable morphologies. In contrast, TGR and TGU differ more noticeably, especially around the central sulcus, which was less developed in TGU. This observation suggests that regionalization of tangential growth contributes substantially to shaping cortical folding patterns. To quantitatively compare these effects, we analyzed regional averages of local GI and cortical thickness (Fig. 5c). Deviations under different growth conditions were assessed against real-brain metrics using a similarity score, where values closer to one indicate stronger correspondence. Among the tested conditions, OGR and TGR show the closest agreement with real-brain measures, consistent with their folding observations in Fig. 5b. Their close correspondence, despite differing only in the radial growth component, implies that thickness variation plays a relatively minor role in driving major fold formation^40^. By contrast, uniform growth conditions, particularly isotropic growth, exhibit poorer alignment across most regions, with pronounced discrepancies in cortical thickness. To further evaluate growth effects, we compared the number of major sulci in the left hemisphere, identified using empirical curvature and patch-area thresholds (Fig. 5d). As expected, IGU produces substantially fewer sulci than both real brain and other growth models, consistent with its sparse folding patterns. The remaining three conditions generated sulcus counts comparable to real brains, with TGU showing slightly elevated median values, though these differences were not statistically significant. Overall, these findings demonstrate that heterogeneous tangential growth is indispensable for reproducing realistic cortical folding development^19^.

Anatomical geometry and regional growth shape cortical folding patterns

In this section, we investigate how anatomical configuration and growth models influence cortical folding outcomes. By systematically examining these factors, we demonstrate the framework’s capacity to generate diverse folding morphologies under varying control parameters and to reveal the developmental mechanisms underlying cortical folding.

We first examined the role of anatomical configuration, focusing on initial brain geometry and cortical thickness (Fig. 6). Owing to individual developmental variation, brain morphologies exhibit noticeable geometrical differences while remaining largely smooth at 22GWs (Fig. 6a). From a mechanical perspective, these variations translate into distinct initial geometries that influence mesh generation and introduce varying undulations capable of triggering different folding patterns. This effect is evident in the modeled folding outcomes across subjects (Fig. 6a). While localized folding patterns differ, major features such as primary sulci and gyri are consistently reproduced. This correspondence mirrors imaging observations and underscores the robustness of our framework in generating anatomically realistic brain morphologies.

We further assessed the influence of initial cortical thickness on folding patterns. Although this factor has been widely explored in theoretical analyses using simplified brain geometries such as spheres or ellipsoids^20,24,25^, it has not been systematically investigated in realistic brain models. To address this, we varied cortical layer thickness by offsetting the cortical surface inward or outward, thereby altering the initial cortical thickness (Fig. 6b). Two levels of offset were introduced: modest adjustments (ICT1, increased thickness; RCT1, reduced thickness) and larger adjustments (ICT2 and RCT2), corresponding to 40% and 80% of the mean initial cortical thickness. Increasing the relative thickness produces fewer folds, characterized by bulkier gyri and deeper sulci, whereas decreasing it generates more densely distributed folds with smaller gyri and shallower sulci. This observation aligns with core–shell buckling theory, which reveals that greater shell thickness leads to less complex folding patterns^41^. To quantify these effects, we measured regional mean sulcal depth and counted the number of major sulci in the left hemisphere. Sulcal depth increases systematically with initial cortical thickness across all regions, with larger relative thickness producing proportionally greater changes (Fig. 6c). Conversely, major sulci enrich as the relative thickness decreased (Fig. 6d), mirroring the observed morphological patterns. Together, these findings demonstrate that variations in either initial geometry or cortical thickness can generate distinct folding morphologies while preserving consistent major folding features observed in early developing brains.

Next, we investigated the role of growth models in directing regional deformation during development (Fig. 7). This analysis is pertinent to developmental delays observed in neurodevelopmental disorders such as autism spectrum disorder^11^ and schizophrenia^12^, both linked to atypical cortical growth trajectories during prenatal stages. Unlike the growth conditions examined in the previous section, here we focused specifically on growth functions, considering two attributes: growth trajectories and growth values. Because thickness variation plays only a minor role in shaping major folding patterns, all modifications were applied exclusively to tangential growth models.

To examine the influence of growth trajectories, we introduced three representative temporal profiles, namely the hyperbolic tangent (HYT), Gompertz distribution (GOM), and polynomials (POL) distributions, while maintaining identical final growth ratios (Fig. 7a). These trajectory alterations were applied globally across all cortical regions. Altering the growth trajectories produces noticeable changes in the timing of fold emergence (Fig. 7b): the onset of folding varies among cases, following the temporal characteristics of each trajectory, yet the final folding morphologies remain largely similar (Movie 2). To quantify local differences, we calculated relative percentage errors in local GI with respect to the baseline trajectory (SYR) (Figs. 7c–d). Substantial deviations were observed for HYT and POL in both vertex-wise and regional comparisons, with HYT generally overestimating and POL underestimating local GI values, consistent with their growth divergence from the baseline growth trajectory. In contrast, the number of major sulci was comparable across all cases, though POL tends to yield slightly fewer sulci (Fig. 7e). These findings suggest that growth trajectories primarily modulate localized folding characteristics while exerting relatively minor influence on global folding organization^25^.

Because the trajectory alterations above were applied globally, we next examined the effects of regional perturbations on folding evolution. Two cases were considered: alteration within a single region (SR) and alteration across multiple adjacent regions (MR), with all other regions following the baseline trajectory. Folding morphologies were projected onto a spherical surface to visualize sulcal regions, defined by negative mean curvature (Fig. 7f). Regional alterations produce distinct sulcal patterns compared to the global alteration case (AR), particularly under single-region modification. Quantitative analyses reveal that trajectory alterations influenced not only the perturbed regions, but also their immediate neighbors, while distant regions (e.g., R14: medial occipital, R17: middle and inferior temporal, and R18: lateral occipital) show minimal changes in sulcal depth (Fig. 7g).

We then explored the effect of growth ratio magnitude on folding outcomes. Reducing the growth ratios by half, either regionally or globally, results in less developed folds within the affected regions, whereas unaffected areas remain largely unchanged (Fig. 7h). Relative deviations again reveal that regional alterations also propagate to neighboring areas (Fig. 7i). In summary, these analyses indicate that growth ratio magnitude predominantly governs the formation of major folding patterns, while growth trajectories modulate localized folding features, as reflected by local quantitative deviations. On this basis, we investigated hemispherical symmetry in folding patterns, as the two hemispheres can differ in both geometry and growth profiles even for the same subject^42^(Fig. S5a). The modeled brain showed asymmetric folding patterns, while preserving consistent major features across hemispheres (Fig. S5b). Although measurable differences were detected in both folding morphology (Fig. S5c) and regional measures (Fig. S5d), these did not reach statistical significance (Fig. S5e).

Our modeling framework well captures abnormal cortical folding in malformations

Through the preceding analyses, we demonstrated that the framework can generate diverse brain morphologies with realistic folding features by varying geometrical factors and regional growth models. Building on this, we further show that by adjusting these control parameters according to expert case descriptions derived from MRI scans (Table S1), the framework can capture abnormal folding phenotypes of typical malformations (Fig. 8). Representative examples include lissencephaly^43^ and polymicrogyria^44^, which are characterized by atypically reduced or excessive cortical folds compared to normal brain morphology shown in Fig. 8a. For example, Fig. 8b illustrates a severe lissencephaly case presenting a smooth cortical surface and markedly increased initial cortical thickness. This phenotype was reproduced in our framework by doubling initial cortical thickness and reducing regional growth ratios to 30% of their normal values (Movie 3). Similarly, pachygyria, a subtype of lissencephaly in which major sulcal landmarks are preserved despite an overall smoother surface, was recapitulated by applying the same two-fold increase in initial cortical thickness while maintaining normal regional growth models (Fig. 8c) (Movie 4).

Beyond global malformations, our framework can also capture region-specific developmental deficiencies. In one case, pachygyria is localized to the frontal and parietal lobes of the right hemisphere and presented as a smooth cortical surface with preserved thickness (Fig. 8d). To replicate this phenotype, we exclusively reduced the tangential growth ratios by half in the affected regions (R1: dorsal precentral & posterior cingulate, R4: lateral precentral & postcentral, and R11: caudal middle frontal) (Movie 5). Conversely, the regional polymicrogyria case, which exhibits abnormally dense but shallow folds in the right frontal lobe (Fig. 8e), was recreated by reducing initial cortical thickness by 40% in selected frontal regions while keeping all other parameters unchanged (Movie 6). Another polymicrogyria case, accompanied by volume atrophy in the right prefrontal cortex (yellow arrow in Fig. 8f), was reproduced by reducing initial cortical thickness by 40% in regions R6 (inferior frontal) and R16 (dorsal prefrontal) and halving both cortical and subcortical growth ratios (Movie 7). Although these reproductions remain qualitative owing to the limited availability of MRI data, particularly longitudinal scans, the results underscore the framework’s capability to replicate malformation-specific folding phenotypes. By linking growth dynamics and anatomical configuration to pathological morphology, this approach provides a mechanistic basis for interpreting cortical malformations within a unified developmental framework.

Discussion

Cortical folds encode critical information about brain structural organization and functional specialization, yet their formation mechanisms remain poorly understood. While multiple factors such as neuronal migration, axonal tension, and cranial constraints have been implicated in shaping folding, mechanical forces arising from differential growth between gray and white matters are increasingly recognized as key drivers of cortical morphogenesis^50^. Previous computational efforts have provided valuable insight into this process but often lacked anatomical fidelity and fell short of reproducing the spectrum of folding patterns observed in the developing human brain^21,27^. In this study, we introduced a whole-brain computational modeling framework that integrates heterogeneous, region-specific growth models with anatomically realistic prenatal brain geometries, thereby unprecedentedly enabling faithful simulation of developmental folding trajectories documented in MRI datasets. The framework captures the spatiotemporal complexity of cortical morphogenesis, achieving anatomical fidelity in both the qualitative folding features and the quantitative morphometric measures.

Our findings demonstrate that tangential growth is the dominant driver of cortical folding, in agreement with both neuroimaging evidence and theoretical predictions^51^. Importantly, the inclusion of regional heterogeneity in tangential growth is essential for reproducing realistic folding morphologies. In the developing human brain, such heterogeneity is well documented and primarily arises from regionalized gene expression^52^ and molecular specialization^53^ across cortical areas. For example, spatial transcriptomic profiling has revealed that the prenatal cortex exhibits area- and layer-specific gene expression programs as early as mid-gestation, with frontal enrichment of genes such as CBLN2 and CYP26A1 and sharp transcriptional boundaries emerging between adjacent visual areas^7^. These spatially patterned molecular programs give rise to developmental heterogeneity across cortical regions, reflected in tissue-level organization such as axonal connectivity and vascular architecture, which together shape diverse folding morphologies^54^. Despite its biological relevance, growth heterogeneity has been underexplored in computational modeling of cortical folding. This lack of focus largely stems from the scarcity of high-quality MRI data and the inherent challenges of incorporating growth heterogeneity into anatomically realistic simulations^22^. Previous studies often relied on highly simplified geometries or imposed periodic growth fields on planar or patch-based models, limiting their ability to capture the complexity of cortical morphogenesis^25,55^. Such oversimplifications, while mathematically tractable, yield insights that are difficult to generalize to realistic developmental scenarios. Our framework directly addresses these limitations. By constructing anatomically faithful whole-brain models and MRI-based regional parcellations, we are able to assign growth laws derived directly from actual developmental trajectories into biophysical simulations. This approach provides, for the first time, a means to mechanistically examine how heterogeneous growth shapes cortical folding within a realistic developing brain. Notably, the framework offers a flexible platform that can be extended to incorporate other region-specific properties, such as spatial variations in material or cellular composition, provided reliable anatomical measurements are available. Beyond such extensions, it also enables mechanistic exploration of how intrinsic tissue properties, local growth dynamics, and emergent geometric constraints interact during cortical morphogenesis.

Both anatomical configuration and growth models critically shape cortical folding patterns. By systematically varying these factors, we showed that the framework can generate diverse yet realistic folding morphologies. This capability not only advances our understanding of the fundamental determinants of folding development but also, more importantly, provides a reliable platform for generating large-scale, high-fidelity synthetic brains with continuous developmental trajectories. Such synthetic datasets could serve supportive educational purposes, for instance by providing intuitive visualizations of cortical morphogenesis and helping learners or practitioners better understand developmental concepts. Moreover, they can effectively supplement the limited availability of longitudinal brain imaging datasets by supplying anatomically faithful inputs for ML models that predict cortical folding development, a particularly valuable contribution given the computational cost of biophysical simulations. Recent advances in ML models, including GAN-based generative frameworks^56,57^, physics-informed neural networks (PINNs)^58^, transfer-learning approaches^59^, and sequence-learning architectures such as PointNet–LSTM^60^, have shown strong promise for capturing and forecasting cortical morphogenesis across gestation. Our framework provides biologically grounded synthetic data to train and validate these approaches. When extended to simulate patient-specific or atypical developmental patterns, it could generate abnormal developmental datasets to train predictive models for early detection and intervention in neurodevelopmental disorders.

Beyond the forward process of predicting folding patterns from initial brain morphologies, our modeling framework also suggests the potential for reverse inference, thus estimating early brain states from mature cortical configurations. Because the simulations are constrained by biologically informed growth dynamics and biophysical principles, they preserve surface topology and generate smooth geometric transitions, thereby maintaining a coherent developmental trajectory from fetal to adult stages. This continuity mirrors realistic longitudinal development and opens the possibility of reconstructing early anatomy from later morphologies. This strategy resonates with the broader concept of “reverse engineering the brain”^61^, here applied to cortical morphogenesis, where mature folding patterns are used to infer plausible early anatomical states. Although still speculative, this direction could yield extensive physically grounded representations of fetal brain morphology. Such synthetic datasets can alleviate the scarcity of fetal imaging and provide valuable context for studies of neonatal development and adult brain disorders rooted in early developmental deficits. Furthermore, this approach also holds promise for advancing evolutionary neuroscience, for example by informing the interpretation of fossil endocasts^62^ and the reconstruction of extinct or ancestral brain states^63^.

Our framework can qualitatively replicate abnormal folding phenotypes, including lissencephaly, pachygyria, and polymicrogyria, by modifying growth laws or geometric parameters according to clinical case descriptions. These reproductions imply how malformation-specific perturbations give rise to distinct phenotypes. Clinically, such models could support early risk stratification and diagnosis by evaluating malformation progression from prenatal scans, offering a complementary tool for obstetric imaging and counseling. Beyond gross shape changes, certain conditions such as autism are also associated with localized variations in folding patterns^64^. While our framework does not directly capture such subtle cases, it provides a potential mechanical interpretation: alterations in growth trajectories, though not significantly affecting major folding architecture, may modulate local folding characteristics within affected regions and their neighbors. Furthermore, by disentangling the relative contributions of geometry and growth models, we found that primary folding patterns are more strongly governed by geometrical factors, including the initial cortical surface configuration and cortical thickness. In contrast, growth functions act more as dynamic modulators, influencing the magnitude and distribution of tissue deformation through their growth ratios.

While the present framework successfully captures key aspects of cortical folding, several limitations also point toward opportunities for refinement. First, axonal fiber growth was not explicitly represented; its potential contribution was approximated through the imposed spatial growth heterogeneity. Yet, developmental processes such as axonal elongation^65^, pathfinding^66^, and fiber organization^67^ have been implicated in the emergence of cortical folds^50^. Second, white matter was modeled as a homogeneous material with isotropic growth, without accounting for its internal anatomical structures. Incorporating these structural details in future models may improve the accuracy of folding predictions^5^. Additionally, spatial variations in mechanical properties were not considered, despite their recognized influence on brain biomechanics^68^. Nonetheless, the framework is inherently extensible and can readily integrate such regional heterogeneities in growth and material composition, paving the way for increasingly comprehensive and biologically realistic models of cortical development.

Collectively, this work establishes an imaging-informed, whole-brain modeling framework that mechanistically links regional growth heterogeneity to cortical folding. By uniting population MRI with biophysical modeling, the framework offers a quantitative bridge between developmental imaging and physical theory to capture the spatiotemporal complexity of cortical morphogenesis. Beyond reproducing anatomically faithful folding patterns, it provides a quantitative and extensible platform for probing the developmental origins of neurodevelopmental disorders and for generating high-fidelity synthetic datasets that support data-driven neuroscience.

Methods

Brain imaging data

To characterize regional growth models, we assembled longitudinal measurements of brain surface area and cortical thickness from the developing Human Connectome Project (dHCP)^69^, which is a publicly available resource focused on early brain development. To ensure reconstruction and registration fidelity, we selected 120 scans (60 males, 60 females) aged 22–40 GWs to form the finalized study dataset. Sex differences were not considered given the broad developmental focus of this study. All images were processed with the infant dedicated computational pipeline, iBEAT V2^70^, to reconstruct high-quality cortical surfaces, followed by cross-subject surface registration in spherical space to ensure vertex-wise anatomical correspondence. For surface partitioning, we applied our previously developed developmental brain atlas^4^, which divides the cortex into 36 regions (18 per hemisphere). The parcellation was constructed using a data-driven method, non-negative matrix factorization (NMF), that groups cortical vertices into regions based on shared developmental patterns of surface area expansion^4^. Regional measurements of cortical surface area and thickness across gestational age were then extracted. These longitudinal profiles served as the basis for further characterization of regional growth models.

Symbolic regression for growth models characterization

We used symbolic regression to characterize regional growth models. Inspired by Darwinian principles of natural selection, symbolic regression automatedly discovers mathematical relationships directly from data without assuming a predefined functional form, thereby offering high flexibility in model discovery. Symbolic regression is implemented through genetic programming, where candidate functions are efficiently represented as binary tree structures composed of variables, mathematical operators (unary or binary), and constants (Fig. S2a). During training, the algorithm emulates evolutionary processes using tree-based operations, such as evaluation, selection, mutation, and crossover. Mutation introduces diversity by randomly modifying nodes within an expression tree, while crossover generates new candidate functions by recombining components from different parent trees. These processes are repeated iteratively until either an optimal function is identified, or the maximum number of generations is reached.

To depict surface area expansion and cortical thickness variation, we decomposed cortical growth into two independent modes: tangential (in-plane) growth and radial (out-of-plane) growth. For computational implementation, the data were converted into unitless growth ratios and a virtual time variable using the following formulas:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ denote tangential and radial growth ratio, respectively; $[eqn]$ and $[eqn]$ are the initial surface area and cortical thickness measured at 22 GWs ( $[eqn]$ ); and $[eqn]$ corresponds to the maximum gestational age in the dataset (40 GWs). After conversion, $[eqn]$ serves as the single input variable for the symbolic regression algorithm to identify the optimal growth function $[eqn]$ for each region (Fig. S2b).

All symbolic regression analyses were performed using PYSR^71^, an open-source package built on Julia. Training was restricted to a limited operator set: binary operators included addition (+), multiplication (*), and polynomial terms ( $[eqn]$ ), while unary operators were constrained to exponential ( $[eqn]$ ), hyperbolic tangent ( $[eqn]$ ), natural logarithmic ( $[eqn]$ ) functions. Each training run was capped at 30 minutes or 1,000 iterations, whichever occurred first, as most qualified models converged within this range during preliminary testing. The maximum expression tree depth was set to 10, and model complexity was constrained to 100. We used the default convergence criterion (“best”) to guide the model selection process, which efficiently balances fitting accuracy and model simplicity. To satisfy growth ratio normalization, the predicted growth function was constrained to $[eqn]$ , implemented through a penalty term in the customized objective function. For each brain region, training was repeated at least three times, and the best-performing candidate was selected as the final model. All analyses were conducted on a Legion PC equipped with a 12-core Intel Core i9–14900HX 2.2GHz CPU and 32GB of RAM.

Biophysics in modeling brain folding

Cortical folding arises from complex biological processes, but from a mechanical perspective it can be approximated as a buckling and post-buckling problem. To capture growth-driven folding in the developing brain, we adopt the well-established differential growth theory, which attributes folding to deformation mismatch between the rapidly expanding cortical layer and the more slowly growing white matter^72^. The deformation kinematics can be described within the framework of continuum mechanics. Firstly, we consider a one-to-one mapping between the undeformed and deformed continuum body, denoted as $[eqn]$ , which carries a material particle $[eqn]$ in the reference configuration $[eqn]$ to its position $[eqn]$ in the current configuration $[eqn]$ . This transformation is characterized by the deformation gradient, $[eqn]$ . Following the principle of multiplicative decomposition, the deformation process can be interpreted as a combination of stress-free growth and elastic accommodation. Accordingly, both the deformation gradient $[eqn]$ and Jacobian determinant $[eqn]$ , which quantifies volumetric change, are decomposed into elastic and growth components:

$$[eqn]$$

For the growth component $[eqn]$ , we modeled gray matter growth as orthotropic, comprising both the in-plane surface expansion and out-of-plane thickness variation:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the tangential (in-plane) and radial (out-of-plane) growth ratios, respectively, with their expressions $[eqn]$ and $[eqn]$ derived from symbolic regression. For white matter, we assume isotropic growth $[eqn]$ , implying uniform expansion in all directions, a simplification widely supported in the literature^10,23,25^. The vector $[eqn]$ denotes the local surface normal, oriented outward in the radial direction.

For the elastic component $[eqn]$ , we followed previous brain-folding models and adopted a standard neo-Hookean hyperelastic model to describe the constitutive behavior of both gray and white matters^19^. The strain energy density function is expressed in terms of the elastic deformation gradient $[eqn]$ and its Jacobian $[eqn]$ ,

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are bulk and shear moduli, respectively. From thermodynamic principles, the first Piola-Kirchhoff stress $[eqn]$ is work-conjugate to the deformation gradient $[eqn]$ , and specially related to $[eqn]$ :

$$[eqn]$$

The Piola stress enters the standard balance of linear momentum, which governs mechanical equilibrium. In the absence of body forces, this reduces to the vanishing divergence of the first Piola-Kirchhoff stress:

$$[eqn]$$

In summary, gray and white matter growth defines a stress-free configuration, which must be balanced by the elastic deformation arising from material confinement. The neo-Hookean model provides this elastic response, ensuring mechanical compatibility and preventing unphysical deformations.

Computational model of the developing brain

To realistically model cortical folding in the prenatal brains, we used reconstructed cortical surfaces at 22 GWs as the initial simulation geometry. As shown in Fig. S4a, the gray and white matter surfaces were registered and parcellated according to the developmental atlas. Solid models were reconstructed, merged, and meshed to produce a high-quality hexahedral volume. Hexahedral meshes were chosen over tetrahedral meshes for their superior robustness and efficiency in large-scale brain folding simulations. A mesh size of 0.4 mm was used, resulting in ~ 1,600,000 hexahedral elements for a whole-brain model. As seen in Fig. S4b, interior elements formed uniform voxels, while boundary elements were split, reprojected, and smoothed to preserve surface conformity. Elements on both sides of the gray–white interface shared faces and nodes, ensuring geometric and mechanical continuity. This meshing strategy ensured at least five elements across cortical thickness, providing sufficient resolution to capture folding development in finite element simulation^23^. Mesh quality was evaluated using scaled Jacobian and edge ratio metrics, with most elements close to 1, confirming that the mesh satisfied the requirements for reliable simulation (Fig. S4c). After meshing, gray matter elements were subdivided into 19 regions according to cortical surface parcellation for implementing regional growth in the modeling (Fig. S4d). To maintain consistency with anatomical parcellation, the regional convoluted surface was constructed from gray and white matter surfaces sharing the same regional index, and enclosed elements were assigned to the corresponding region. Elements near regional interfaces were allocated to the adjacent region with the higher index. Furthermore, orthotropic growth was enforced by assigning each gray matter element a material orientation, calculated as the normal vector from the element centroid to the nearest vertex on the cortical surface (Fig. S4e).

Brain folding was simulated as thermal expansion in ABAQUS (Dassault Systems, Paris, France) with the dynamic-explicit solver, considering the analogy between the volumetric growth and thermal expansion^73^. Both gray and white matter were modeled as nearly incompressible neo-Hookean materials, with comparable shear stiffness values of 0.65 kPa and 0.68 kPa, respectively^74^. Regional growth models $[eqn]$ and $[eqn]$ , derived from symbolic regression, were applied to the gray matter through a user-defined subroutine VUEXPAN. Following previous studies^19^, we applied a sigmoid shape growth profile across cortical depth to reflect the layered growth patterns observed in cortical development^6,7^ (Fig. S4f). In contrast, the white matter was modeled using an isotropic growth model that captures its volumetric changes during normal brain development (Fig. S6). Free boundary conditions were imposed on the cortical surface, together with a self-contact constraint to prevent self-penetration. A uniform pressure was applied to represent the mechanical confinement of the meninges and developing skull^23^. The total simulation time was set to the maximal rescale gestational time ( $[eqn]$ ). To avoid computational instabilities, the maximum time step was set as $[eqn]$ , where $[eqn]$ is the average mesh size, $[eqn]$ is the mass density, $[eqn]$ is bulk modulus^25^. Detailed mechanical parameters are provided in Table. S2. All simulations were performed on a Dell Alienware desktop equipped with an Intel^®^ Core^™^ Ultra 9 285K processor (24 cores) @ 3.7 GHz, and 64 GB of RAM. For a single subject, the simulation of brain development from 22 to 40 GWs required approximately 12 hours of computation.

Postprocessing analyses and quantitative metrics

In total, we simulated folding development for 10 brains (2 hemispheres per case) at 22 GWs. Simulation output files were exported and analyzed with in-house postprocessing codes. At selected simulation stages, node coordinates were extracted to reconstruct the cortical surfaces. For regional analysis, alignment between simulated and real brain parcellations was achieved by matching vertex positions between the undeformed simulated surface and the real brain geometry at 22 GWs. The assigned regional indices were directly propagated to subsequent deformed states, as the simulation preserved the original cortical surface connectivity and maintained geometry without distortion. Notably, the framework allows cortical surfaces to be reconstructed at any simulation moment, allowing continuous tracking of folding progression that cannot be fully captured by longitudinal imaging. From the reconstructed cortical surfaces, multiple quantitative metrics were calculated to evaluate the accuracy of folding simulations. These include global and local gyrification index, major sulcus count, sulcal depth, curvature measures, enabling systematic comparisons at whole-brain, regional, and vertex levels.

Global gyrification index:

To quantify the folding complexity of the brain surface, we employed the global gyrification index (GI), a widely used whole-brain metric. The global GI is defined as the ratio between the total cortical surface area and the area of its enclosing convex hull^19^. In our implementation, the convex hull was approximated using an alpha-shape surface, which fully encloses the cortex while minimizing the enclosing area to provide a more accurate estimate.

Local gyrification index

In addition to the global GI, we quantified regional folding complexity using the local gyrification index. The local GI complements the global measure by capturing spatial variations in cortical folding and providing a region- or vertex-wise characterization. To enhance accuracy, we applied the anisotropic kernel method^75^, which uses wavefront propagation to construct spatially adaptive kernels that conform to cortical geometry, elongating along gyral crowns and sulcal fundi while remaining uniform over flatter regions. Compared with conventional kernels, this approach provides more reproducible local gyrification estimates that better reflect biologically relevant variations. Herein, the local GI values were computed at each vertex of the cortical surface.

Sulcal depth

Sulcal depth is another widely used quantitative descriptor of cortical folding. Several approaches have been proposed for computing sulcal depth, differing in their choice of reference surface and distance metric. In this study, we defined sulcal depth as the Euclidean distance between the cortical surface and its convex hull, consistent with the definition used for calculating the global GI^27^. For each vertex on the cortical mesh, the nearest projection onto the triangular facets of the convex hull was identified, and the minimal Euclidean distance was taken as the sulcal depth at that vertex. Given the convexity of the hull, the shortest distance is uniquely defined, ensuring stable vertex-wise measurements across the surface.

Curvatures

We consider the Weingarten curvature matrix to quantify curvatures in all directions. Its eigenvalues define the principal curvatures, $[eqn]$ and $[eqn]$ , corresponding to the maximum and minimum normal curvatures. From these, the mean curvature is given by $[eqn]$ , and the Gaussian curvature by $[eqn]$ . Among these descriptors, mean curvature has been widely used in cortical folding studies. However, its magnitude depends on surface size and geometry, making it sensitive to brain scale. To obtain a scale-independent measure, we adopted a dimensionless mean curvature^27^, defined by multiplying $[eqn]$ with a characteristic length $[eqn]$

$$[eqn]$$

where $[eqn]$ is the surface area. Dimensionless mean curvatures were computed at each vertex of the triangular mesh and, for simplicity, are hereafter denoted as “mean curvature”. In addition, we defined the shape index (SI) as another scale-independent descriptor:

$$[eqn]$$

with values confined to [−1, 1]. Unlike raw curvature values, this normalized metric provides a normalized characterization of local surface types. Specifically, negative values correspond to concave features such as sulci, positive values to convex features such as gyri, and values near zero to saddle-like transitional regions^29^.

Sulcus number estimation

To quantify major sulci, we employed a curvature- and size-based filtering approach. Sulcal regions were identified as contiguous patches of the cortical surface with mean curvature value less than − 0.05. To ensure robustness and to exclude spurious indentations, we further required a minimal surface area of 10 mm^2^ and at least 60 vertices per patch. We note that this threshold-based method may underestimate the number of minor or tertiary sulci, which can be functionally relevant but are easily smoothed out by curvature filtering or mesh resolution. While curvature and size thresholds are empirically defined, this method offers a straightforward means of quantifying gross sulcus counts and supporting consistent comparison between real and simulated brains. Similar curvature-based criteria have been applied in surface-based morphometry studies of cortical folding^76^.

Spherical mapping, sulcal landmarks, and similarity assessment:

To facilitate comparison of cortical folding patterns, we projected the brain morphology and associated quantitative measures onto a unit sphere. This was achieved using the FLASH algorithm for spherical mapping^77^. Surface registration is a necessary step to ensure meaningful cross-subject comparison. Since the real brain morphology had already been registered during cortical surface mapping, it was sufficient to build a mapping between the reconstructed brain and the undeformed simulated brain. This mapping, propagated across the computational steps, was obtained through one-to-one distance-based vertex matching once the computational mesh was constructed. For consistency, we defined two element intersections as the north and south poles, respectively: the intersection of regions 6, 15, and 16 as the north pole, and the intersection of regions 9, 17, and 18 as the south pole (Fig. 3c).

The sulcal landmarks (Fig. 3) were manually extracted and visualized using ParaView. We focused on six major sulci, including central sulcus, postcentral sulcus, intraparietal sulcus, superior temporal sulcus, lateral sulcus, and inferior frontal sulcus, which were selected for their anatomical prominence and accessibility of extraction. These sulci also carry important functional relevance, as they delineate boundaries of primary sensorimotor cortices (central and postcentral sulcus), visuospatial and attentional networks (intraparietal sulcus), auditory and language-related regions (superior temporal and lateral sulcus), and higher-order executive areas (inferior frontal sulcus). Although the alignment is not restricted to these six sulci, a more detailed landmark alignment analysis lies beyond the scope of this study. Here, we highlight these six representative sulci to demonstrate the capability of our whole-brain framework to model brain morphology with anatomical fidelity.

After surface registration, the similarity $[eqn]$ between distributions of a quantitative metric $[eqn]$ in real brains ( $[eqn]$ ) and simulated brains ( $[eqn]$ ) was evaluated using the following formula^39^:

$$[eqn]$$

where $[eqn]$ is the total number of vertices; $[eqn]$ is the mapping from the real brain to the simulated brain, accommodating the higher vertex density of the real mesh; $[eqn]$ denotes Euclidean norm, i.e., $[eqn]$ . Noted, to ensure that the similarity measure always falls within the range of [0,1], the quantitative metric $[eqn]$ must be normalized to [−1,1], with a higher value of indicating greater similarity^39^. In this study, we calculate the similarity based on the distribution of the shape index, which inherently lies within the required range. By contrast, the similarity score used in Fig. 5c was defined as $[eqn]$ , where $[eqn]$ and $[eqn]$ denotes the metric measure of simulated and real brain cases, respectively. This similarity score has values between 0 and 1, with 1 indicating perfect agreement. The relative percentage errors in Fig. 7 was defined as $[eqn]$ , where $[eqn]$ and $[eqn]$ denotes the metric measure of studied and reference brain case, respectively.

Growth profile definition

To investigate the influence of growth trajectories on cortical folding simulation, we introduced three additional growth models based on a high-order polynomial (POL), hyperbolic tangent (HYT), and Gompertz distribution (GOM), respectively (Fig. 7a). These functional forms have been used in modeling tissue growth processes such as tumor progression^78^ and the development of myelinated white matter in the human brain^79^. The growth functions are defined as

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ are constants. As an illustrative case, we considered the tangential growth model of region 1 (Fig. 7a), for which the symbolic regression (SYR) predicted model is given by $[eqn]$ . Accordingly, the parameter values were set as: POL, $[eqn]$ , HYT, $[eqn]$ , and GOM, $[eqn]$ . The parameter values were empirically assigned to ensure that all models converge to a consistent growth ratio at the final stage, thereby eliminating differences in absolute growth magnitude and allowing the analysis to focus exclusively on the effect of growth trajectories.

For the growth conditions shown in Fig. 5, the isotropic growth with a uniform profile (IGU) case was defined by linear growth functions $[eqn]$ , whereas the tangential growth with a uniform profile (TGU) case used $[eqn]$ The chosen growth ratio values follow those commonly adopted in previous studies for modeling simplification^10,25^.

Statistical analysis

Statistical analyses were performed using non-parametric tests given the small sample size (n = 10 per group), although data satisfied normality as assessed by the Shapiro–Wilk test. In Fig. 5d, difference in the sulcus counts among five growth conditions (Real, OGR, TGR, TGU, and IGU; n = 50, 10 per group) were first evaluated using the Kruskal–Wallis test to evaluate the null hypothesis that group medians were equal. When the overall test was significant (P < 0.05), post hoc pairwise comparisons were carried out using Dunn’s test with Bonferroni correction for multiple testing. In Fig. S5e, paired comparisons of sulcal counts between left and right hemispheres (n = 20, 10 per hemisphere) were assessed using the Wilcoxon signed-rank test. The null hypothesis was that the median differences between pairs was zero. All tests were two-tailed, with P < 0.05 considered statistically significant. Statistical analyses were performed using OriginPro 2025.

Supplementary Material

Supplementary Files

This is a list of supplementary files associated with this preprint. Click to download.

• Movie3.mp4

• Movie6.mp4

• Movie7.mp4

• Movie5.mp4

• Movie4.mp4

• Movie1.mp4

• Movie2.mp4

• Supplementaryinformation.docx