An integrated finite element method and machine learning algorithm for brain morphology prediction

Abstract

The human brain development experiences a complex evolving cortical folding from a smooth surface to a convoluted ensemble of folds. Computational modeling of brain development has played an essential role in better understanding the process of cortical folding, but still leaves many questions to be answered. A major challenge faced by computational models is how to create massive brain developmental simulations with affordable computational sources to complement neuroimaging data and provide reliable predictions for brain folding. In this study, we leveraged the power of machine learning in data augmentation and prediction to develop a machine-learning-based finite element surrogate model to expedite brain computational simulations, predict brain folding morphology, and explore the underlying folding mechanism. To do so, massive finite element method (FEM) mechanical models were run to simulate brain development using the predefined brain patch growth models with adjustable surface curvature. Then, a GAN-based machine learning model was trained and validated with these produced computational data to predict brain folding morphology given a predefined initial configuration. The results indicate that the machine learning models can predict the complex morphology of folding patterns, including 3-hinge gyral folds. The close agreement between the folding patterns observed in FEM results and those predicted by machine learning models validate the feasibility of the proposed approach, offering a promising avenue to predict the brain development with given fetal brain configurations.

brain development, cortical folding, computational modeling, machine learning, surrogate model

Introduction

The development of human brain entails a complex and intricate series of processes that culminates in the transformation of a small and smooth brain into a folded structure between 25 and 40 gestational weeks. During development, the cerebral cortex undergoes a significant expansion in volume and surface area accompanied by cortical folding, which results in the formation of convex hills known as gyri and concave valleys called sulci on the cerebral cortical surface (Vasung et al. 2016). Primary folds are the largest, deepest, and earliest to emerge, whereas secondary and tertiary folds emerge through the subdivision of primary folds (Huang et al. 2009). Despite interindividual variation in cortical morphology, the primary gyro-sulcal cortical pattern of species is remarkably preserved across individuals (Lohmann et al. 2008; Dubois et al. 2008b; Li et al. 2016; Zhang et al. 2016). The size and morphology of these folds exert a profound influence on the brain’s performance (Fernández et al. 2016). Evidence suggest that numerous neurological disorders, including schizophrenia, autism spectrum disorder, lissencephaly, and polymicrogyria, are linked to the deviation from the normal gyrification of the brain (Landrieu et al. 1998; Sallet et al. 2003; Tortori-Donati et al. 2005; Nordahl et al. 2007; Dubois et al. 2008a; Giedd and Rapoport 2010; Li et al. 2014). Thus, a comprehensive understanding of the cortical folding could be a vital key for early detection of cognitive impairments as well as neurodevelopmental and psychiatric disorders.

Folding patterns of the human cerebral cortex show quite different morphologies between subjects. Evidence also has shown that the folding patterns of the human cerebral cortex can predict its cytoarchitecture (Fischl et al. 2008). Therefore, understanding of the underlying mechanism and quantitative description of folding patterns are 2 important research goals. In the past decade, scientists think the underlying mechanism stems from a complex interplay of biological and mechanical processes. Several possible internal and external factors, including radial constraint (Rakic 1995), cranial constraint (Chen et al. 2010; Garcia et al. 2018), differential growth due to cellular factors (Richman et al. 1975; Shatz 1992; Cartwright 2002; Tallinen et al. 2014), and axon maturation (Essen 1997; Sur and Rubenstein 2005; Hilgetag and Barbas 2006; Nie et al. 2012), have been proposed to contribute to the brain gyrification. For example, one of the earliest theories posited that external constraints and forces, such as the cranial constraint and cerebrospinal fluid pressure acting against the expanding brain, were responsible for cortical folding (Jones and Peters 2012). However, experimental studies have indicated that internal forces play a dominant role in this phenomenon. It is believed that differential growth within the cortex is a possible stimulus for secondary and tertiary folding (Richman et al. 1975; Tallinen et al. 2016). According to the tangential differential growth hypothesis, the outer layers of brain show more rapid growth than the inner layers. This growth rate mismatch acts as the driving mechanism for brain structure instability and gyrification. The tangential differential growth hypothesis has been supported by recent experimental and computational analyses (Ronan et al. 2014; Jalil Razavi et al. 2015; Razavi et al. 2015; Tallinen et al. 2016), although it is believed that other factors such as axon fiber position and growth also play a non-negligible role in the regulation of morphology (Holland et al. 2015; Zhang et al. 2016, 2017).

As we also know, primary folding is notably preserved among individuals, whereas secondary and tertiary foldings evolve after primary folding is completed, and they can vary widely across individuals (Lohmann et al. 2008; Gilles et al. 2013). Folding patterns of the human cerebral cortex show quite different morphologies between subjects; this is true even for major gyri and sulci. Many neuroscientists still widely use Brodmann map as a common reference for mapping neuroimaging data to study the brain structure. However, it does not provide a quantitative description of brain structural well. Different researchers could define their own Brodmann’s map for each individual brain. More recently, researchers have employed a straightforward way to account for huge variability among individual brains by “averaging” a population’s brains (Glasser et al. 2016) and defining more fine-grained cortical regional maps on these averaged, much smoothed, and regularized brain. A basic assumption of image registration methodology is that the brain images under consideration are similar and can be matched well. However, this assumption is fundamentally questionable, considering the enormous variability of cortical anatomy (Brett et al. 2002; Derrfuss and Mar 2009), which substantially limits the accuracy and reliability of image registration algorithms. Meanwhile, a variety of folding pattern descriptors (Thompson and Toga 1996; Van Essen et al. 1998; Woods et al. 1998; Fischl et al. 1999, 2002; Mangin et al. 2004; Li et al. 2008; Yang et al. 2008; Awate et al. 2009; Boucher et al. 2009), ranging from vertex curvature to a small patch, to gyrus and sulcus, to a cortical lobe, and to a whole-brain cortical surface, have been proposed to quantify cortical folding. Notably, our recent work has conceptualized a new type of brain folding descriptor termed 3-hinge gyral junctions (Li et al. 2010, 2016), which is a conjunction of 3 gyral crestlines coming from different directions in cortical folding. Exemplar 3-hinges on a human brain are represented by colored bubbles and curves in Fig. 1. Our prior study shows that 4 gyral crests rarely meet to form a 4-hinge gyrus, whereas 3-hinges are abundant; and they are able to act as an identity tag for individuals. Moreover, these gyral hinges behave more like cortical hubs in the cortico-cortical networks and compose a majority portion of the network’s “core.” Quantitative characterization of gyral folding patterns via hinge numbers with cortical surfaces constructed from MRI data has been used to identify 6 common 3-hinge gyral folds that exhibit consistent anatomical locations across humans, chimpanzees, and macaques, as well as 2 unique 3-hinge patterns in macaques, 6 in chimpanzees, and 14 in humans (Li et al. 2016). It is not surprising that the number of 3-hinge patterns identified in the human brain is 2.5 and 7.8 times greater than the number found in chimpanzee and macaque brains, respectively. Therefore, the shape and number of 3-hinge patterns in a growing brain could be a new metric to characterize the folding of a primate brain.

Fig. 1

An example of 3-hinge gyrus. The centers are represented by 3 color bubbles (blue, orange, and purple). The spokes are represented by solid curves.

Meanwhile, computational modeling has emerged as a powerful tool for investigating the complex and nonlinear interactions that occur during the rapid growth and folding of the human brain, given the difficulties of in vivo experiments, longitudinal imaging, and close observation (Fernández et al. 2016; Wang et al. 2019, 2021; Costa Campos et al. 2021; Garcia et al. 2021; Zarzor et al. 2021; Darayi et al. 2022). Such models enable researchers to explore different hypotheses and examine the effect of various parameters including material properties (Razavi et al. 2015; Tallinen and Biggins 2015), cortical thickness (Tallinen et al. 2014), curvature (Toro and Burnod 2005), axonal fibers distribution (Holland et al. 2015; Razavi et al. 2017; Chavoshnejad et al. 2021; Garcia et al. 2021), in predicting realistic folded morphologies. For example, these studies have demonstrated that a higher cortical growth rate increases gyrification, whereas higher cortical stiffness or thickness reduces gyrification (Jalil Razavi et al. 2015). Moreover, even slight variations in geometry, material properties, loading, and boundary conditions can significantly alter the final morphology of the brain, highlighting the intricacy and uniqueness of cortical morphology among individuals (Naidich et al. 2012).

Generally, the computational models can be classified into 2 categories: mechanistic models and statistic models. Mechanistic models, such as the finite element method (FEM), are developed based on prior knowledge of the mechanical behavior of a given tissue or process. These methods utilize the governing equations of mechanics to solve the problem explicitly or approximately (Bakhaty et al. 2017; Madani et al. 2017, 2019). FEM involves dividing the subject into a large number of elements and nodes, and applying mechanical governing equations to each of these elements. Although this approach yields accurate results, the computational cost associated with FEM can be prohibitively high, making it difficult to use in patient-specific treatments. In contrast, statistical modeling, particularly machine learning (ML), treats the system as a black box without any preconceived behavior assumption, allowing ML to solve the problem in real-time. Promising ML models have been developed that can predict the mechanical response of brain tissues with comparable accuracy to FEM, while dramatically reducing the computational cost (Shim et al. 2020; Wu et al. 2020; Menichetti et al. 2021). However, to ensure accuracy and interpretability, ML models requires large amount of data for training, which can often be obtained through simulations rather than experiments. As a result, the combination of FEM and ML has emerged as a hot topic in the field due to its potential to leverage the advantages of both methods (Martin et al. 2016; Taghizadeh et al. 2016; Martínez-Martínez et al. 2017; Tonutti et al. 2017; Karami et al. 2018; Madani et al. 2019; Calka et al. 2021).

Fig. 2

a) Initial model before applying geometrical undulations. b) Initial model after applying geometrical undulations. The patch model consists of gray and white matters with a thickness of 1.5 and 58.5 mm, respectively. Symmetric boundary conditions are imposed on all 4 sides and the bottom of the model is fixed.

In this study, we present an ML-based FEM surrogate model to predict the morphology of brain and investigate the underlying folding mechanism. Initially, we conduct extensive computational simulations to obtain multiple FEM results that represent brain gyrifications. These results are further used as the inputs for training and validating GAN-based ML models. Subsequently, we implement the trained surrogate model to predict the brain morphology based on the testing data. Finally, we compare the predictions of ML models with results from FEM simulation to evaluate the prediction performance of the ML-based finite element surrogate model.

Theoretical methods and computational models

Constitutive relationship for brain tissue

In continuum mechanics, we introduce the deformation map |$\boldsymbol{\varphi}$| that maps a material particle from original position |$\boldsymbol{X}$| in the undeformed (reference) configuration |${\mathcal{B}}_0$| to its new position |$\boldsymbol{x}=\boldsymbol{\varphi} \left(\boldsymbol{X}\right)$| in the deformed (spatial) configuration |${\mathcal{B}}_t$|⁠. In addition, we introduce the deformation gradient |$\boldsymbol{F}={\nabla}_{\mathbf{X}}\boldsymbol{\varphi}$| to measure the mapping of the line element from reference to spatial configuration and Jacobian J to describe the local volume alternation from reference to spatial configuration. Both |$\boldsymbol{F}$| and J can be further decomposed multiplicatively into an elastic part and a growth part (Garikipati 2009),

$$\begin{equation} \boldsymbol{F}=\boldsymbol{A}\cdot \boldsymbol{G}\ \mathrm{and}\ J=\det \left(\boldsymbol{F}\right)={J}^e{J}^g \end{equation}$$

(1)

where A and |$\boldsymbol{G}$| are elastic and growth deformation tensor, respectively. |${J}^e$|and |${J}^g$| are Jacobian corresponding to the elastic and growth tensor, with the relation described as |${J}^e=\det \left(\boldsymbol{A}\right)$| and |${J}^g=\det \left(\boldsymbol{G}\right)$|⁠. Following previous studies (Budday et al. 2014; Garcia et al. 2021), the growth tensor |$\boldsymbol{G}$| only contains the diagonal components as:

$$\begin{equation} {\displaystyle \begin{array}{c}\left[\boldsymbol{G}\right]=\operatorname{diag}\left({G}_1,{G}_2,{G}_3\right)\end{array}} \end{equation}$$

(2)

here |${G}_1$|⁠, |${G}_2$|⁠, |${G}_3$| are 3 orthogonal components of the growth tensor in the local coordinate. For isotropic growth, |${G}_1={G}_2={G}_3=\gamma$| and |${J}^g={\gamma}^3$|⁠, where |$\gamma$| is the growth rate; for the tangential growth in the x–y plane, as illustrated in Fig. 2, |${G}_3=1$|⁠. The elastic part of brain tissue including the gray matter and white matter is modeled as an isotropic and perfectly incompressible neo-Hookean material with a strain energy function |$W$|⁠, parameterized exclusively in terms of the elastic tensor A and its Jacobian |${J}^e$|⁠,

$$\begin{equation} \begin{array}{c}W(\boldsymbol{A})=\frac{1}{2}\lambda{{ln}}^2{J}^e+\frac{1}{2}\mu [\boldsymbol{A}:\boldsymbol{A}-3-2\ln ({J}^e)]\end{array} \end{equation}$$

(3)

where λ and μ are Lame constants. Following conventional arguments of thermodynamics, the Piola stress |$\boldsymbol{P}$| is energetically conjugate to the deformation gradient (Budday et al. 2014):

$$\begin{equation} {\displaystyle \begin{array}{c}\boldsymbol{P}=\frac{\partial W\left(\boldsymbol{A}\right)}{\partial \boldsymbol{F}}=\frac{\partial W\left(\boldsymbol{A}\right)}{\partial \boldsymbol{A}}:\frac{\partial \boldsymbol{A}}{\partial \boldsymbol{F}}={\boldsymbol{P}}^e\cdot{\boldsymbol{G}}^{-T}\end{array}} \end{equation}$$

(4)

where |${\boldsymbol{P}}^e$| denotes the elastic Piola stress, which can be expressed more in detail for neo-Hookean material,

$$\begin{equation} {\displaystyle \begin{array}{c}{\boldsymbol{P}}^e=\frac{\partial W\left(\boldsymbol{A}\right)}{\partial \boldsymbol{A}}=\mu \boldsymbol{A}+\left[\lambda \ln \left({J}^e\right)-\mu \right]{\boldsymbol{A}}^{-T}\end{array}} \end{equation}$$

(5)

In the absence of body forces, the balance of linear momentum is simplified by expressing it as the vanishing divergence of the Piola stress:

$$\begin{equation} {\displaystyle \begin{array}{c}\mathrm{Div}\left(\boldsymbol{P}\right)=0\end{array}} \end{equation}$$

(6)

Computational models for brain tissue

A 3D cuboid model measuring 60 × 60 × 60 mm is created to represent a small portion of the brain to study the gyrification of the brain in the presence of geometrical undulations. This model has previously been used to study consistent gyrus formation and also fiber density effect on 3-hinge formation (Ge et al. 2018). The model consists of 2 distinct layers: a thin layer of gray matter situated above a white matter substrate. The bottom layer is the core, which is supposed to be a simple representation of the subplate, intermediate zone, and ventricular zone. In human brains, the cortex (gray matter) is a thin (⁠|$2\hbox{--} 4\ \mathrm{mm}$|⁠) layer (Bayly et al. 2014), in contrast to the core that has a much greater thickness of around |$50\ \mathrm{mm}$| (Tallinen et al. 2014). Therefore, the initial thickness of gray matter prior to folding is set to 1.5 mm to be in the range of mature thickness after folding. The dimension of the model was selected based on experimental data gathered from small pieces of a brain. The dimensions and boundary conditions of the model are illustrated in Fig. 2, where symmetric boundary conditions are imposed on all 4 sides of the model. A fixed boundary condition is applied to the bottom surface of the white matter. The dimension of the model was large enough in comparison with the wavelength of folded patterns observed in experiments to prevent boundary effects. In a bilayer model with an isotropic growth for both cortex and subcortex without stress-dependent growth and no boundary confinement, only their relative growth ratio is a determinant factor for instability and folding. Also, in a human premature brain, the volume of the cortical plate increases by 4-fold in the first 2 months (22–30 gestational weeks), whereas the subplate plus intermediate zone (SP + IZ) increases approximately by 3-fold. In our recent analytical-computational study, we have already shown that the growth ratio of 4/3 between 2 distinct layers is sufficient enough to trigger the structural instability. This ratio is independent of the growth ratios’ absolute values. Without loss of generality, we set the cortex grow gradually with no growth in the core (Tallinen et al. 2014; Zhang et al. 2016). The growth of the cortex was assumed to be linear in time.

Moreover, a self-contact constraint is applied on the top surface of gray matter to avoid self-penetration of the outer surface. Brain growth is simulated using relative tangential growth. Research has demonstrated that gyrification in the brain primarily stems from the disparity in growth rates between gray matter and white matter (Budday et al. 2014; Ronan et al. 2014; Tallinen et al. 2016). To capture this, we assume that growth occurs solely in the gray matter layer as the important parameter in the initiation and progression of cortical folding is the growth rate difference between layers not their individual rates. The tangential expansion along the x and y-axes is regulated to approximately achieve a gray matter thickness of 3 mm after folding, with growth being monitored by adjusting the thermal expansion coefficient of the cortex and temperature increase in a dynamic explicit simulation. The analogy between volumetric growth and thermal expansion has been discussed in the literature (Cao et al. 2012b). The growth scenario is the same for all the models and only their initial surface morphologies are different.

Both gray matter and white matter are meshed using an 8-node linear brick (C3D8R) in the ABAQUS FE software. A smaller mesh size is employed in proximity to the top layer to reduce the computational cost. Additionally, a mesh convergence study is performed to eliminate the effect of the mesh size. A python code is developed to implement the geometrical undulations by relocating nodes in the meshed model to create different initial surface morphologies. The geometrical undulations are achieved using the following Gaussian distribution function:

$$\begin{equation} {\displaystyle \begin{array}{c}z\left(x,y\right)=K\exp \left(-\left(\frac{{\left(x-{x}_0\right)}^2}{2{\sigma}_x^2}+\frac{{\left(y-{y}_0\right)}^2}{2{\sigma}_y^2}\right)\right)\end{array}} \end{equation}$$

(7)

where |$K$| is the amplitude, |${x}_0$| and |${y}_0$| are the coordinates of the center, |${\mathrm{\sigma}}_x$| and |${\mathrm{\sigma}}_y$| are the spread of blob along the x and y-axes, respectively. To create an irregular surface of the cortical model, we employed a randomized approach. Specifically, we employed a varying number of Gaussian distributions, ranging from 5 to 15, with amplitudes ranging between −2.5 and 2.5 mm. These undulations mimic the surface irregularities observed in the fetal smooth brain, such as the central sulcus, distinguish the initial surface morphologies of various models, and ultimately facilitate folding (Cao et al. 2012a). For each model, we randomly select the number, amplitude, and position of Gaussian distributions. We apply the undulations to the surface of the cortex prior to the growth process, making their specifications independent of the cortex’s growth. It should be noted that the thickness of the cortex is uniform at the undulation locations and consistent with that at other locations. This results in the maximum undulation being exhibited on the top surface of the gray matter, whereas the bottom surface of the white matter remains a perfect plane. Figure 2 depicts a model before and after the application of geometrical undulations.

Detection of 3-hinges

Our previously developed gyral net pipeline (Chen et al. 2017) is used to extract the number and shape of the 3-hinge folds in the models after folding. Briefly, the pipeline consists of 2 major steps: (i) Gyral crest segmentation: we define the “mid-surface” as a line that separates gyri and sulci. This “mid-surface” is chosen so that the mean of the displacements of all surface vertices from their original locations is zero. Gyral altitude for a vertex is defined as the movement from its original location to the “mid-surface” in the surface normal direction. Based on these definitions, the watershed algorithm is adopted to segment the gyral crest (regions over an altitude level). More details and effects of the watershed algorithm can be found in our previous works. (ii) Gyral skeleton extraction: this step is to skeletonize gyral crests. Gyral skeletons are defined as the crest curves located in the central parts of gyral crest. First, a distance transformation algorithm is conducted on the segmented gyral crests to highlight their central regions. Then, a tree marching algorithm is adopted to successively connect the vertices to form multiple tree-shape graphs. After the redundant branches are pruned, the major branches are left and taken as the skeleton of the gyral crests. Detailed information regarding the extraction of 3-hinge folds can be found in our previous publications (Chen et al. 2017). Regarding the construction and classification of 3-hinges from the cortical surface, we provide a compact description, as shown in Fig. S1 in the Supplemental Materials, based on our previous work (Zhang et al. 2018). For more details, such as the rationales and parameter settings, please refer to our previous publication (Zhang et al. 2018).

ML algorithm and data preparation

ML algorithm

Although the FEM could successfully predict brain development, it is still very time-consuming and computationally expensive. To efficiently generate a smooth and accurate model to demonstrate brain development, we take the results of FEM as “ground-truth” developed models and aim to propose an end-to-end GAN-based ML framework to utilize the initial model to estimate the developed models to be similar to the ground-truth. First, we formulate the problem as |$\boldsymbol{t}=\boldsymbol{s}\odot \boldsymbol{f}$|⁠, where |$\boldsymbol{t}$| is the developed model obtained by FEM and obeys the distribution |${p}_{\boldsymbol{t}}$|⁠. |$\boldsymbol{s}$| is the initial model. |$\boldsymbol{f}$| is the deformation field to be estimated, which should be smooth. |$\odot$| denotes the pixel-wise multiplication. Therefore, during training, our GAN-based framework is aimed to learn the mappings between the input initial models |$\boldsymbol{s}$| and the “ground-truth” developed models |$\boldsymbol{t}$|⁠. During testing, we use the trained framework to predict the developed model |$\boldsymbol{t}$| for any unseen initial model |$\boldsymbol{s}$|⁠. The pipeline of our framework is illustrated in Fig. 3.

Fig. 3

Architecture of the proposed ML network. (i) During the training phase, the initial model and FEM generated models are introduced to guarantee the reliable prediction of the deformation field and (ii) during the testing phase, the testing initial model is input into the trained generator to estimate the deformation field, which further multiplies the testing initial model to achieve the predict developed model to exhibit the brain development. Notably, the discriminator is removed in the testing phase.

As a GAN-based method, the proposed end-to-end framework includes a generator and a discriminator, adopts the initial models as the input, and predicts the corresponding developed models. The detailed ML network architecture can be found in Figs S3 and S4 for the generator and discriminator, respectively, in the Supplemental Materials.

Generator

To reliably estimate the developed model |$\boldsymbol{t}$|⁠, the generator (G) is built based on a UNet-like architecture with an encoder–decoder structure. The 3D coordinates of the initial model |$\boldsymbol{s}$| are formed into 3 channels, concatenated as input to the encoder. The encoder includes 5 encoding blocks, each of which is constructed in a form of Conv-IN-LeakReLU, where Conv denotes convolution layer, IN is the instance normalization, and LeakReLU is the activate function. Similarly, the decoder also contains 5 decoding blocks with a form of Upconv-IN-ReLU, where Upconv denotes the transposed convolution layer and ReLU is the activate function. To alleviate checkerboard artifacts (Isola et al. 2017),the number of training epochs is increased for sufficient training. Moreover, the skip connections are performed between each corresponding pair of encoding and decoding blocks to directly provide the feature representations at different levels, helping preserve essential details that may be lost in the down-sampling operations of the encoder.

Discriminator

The discriminator (D) is designed to force the similarity of the distributions between the estimated models and the ground-truth models. The discriminator is built with 5 convolutional layers. The first layer is the input layer with the form of Conv-LeakReLU. The following 4 layers are used for feature extracting, which keep the consistent form with encoding blocks in generator. The last output layer is a convolutional block (Conv) to output an identity matrix. The average of the identity matrix is the real/fake estimation of the input model.

Training and testing

During training, the generator is trained to learn the longitudinal brain development based on the initial model and predict the deformation field to obtain the developed model. Meanwhile, the discriminator is trained to identify similarity between the FEM generated “ground-truth” and the models estimated by the generator. In practical applications, for each unseen model, we only utilize the trained generator to predict the development.

Loss functions

The proposed generator and discriminator are, respectively, trained by alternately minimizing a generative loss and a discriminative loss.

Generative loss function

To estimate the deformation field from the initial models under a premise of smoothness, we design a generative loss function (⁠|${\mathbb{L}}_G$|⁠), which includes 3 complementary loss terms, i.e. adversarial loss term, similarity loss term, and smoothness loss term.

$$\begin{equation} {\displaystyle \begin{array}{c}{\mathbb{L}}_G={L}_{Adversarial}+{\lambda}_1\ast{L}_{similar}+{\lambda}_2\ast{L}_{\mathrm{Smooth}}\end{array}} \end{equation}$$

(8)

where |${\lambda}_1$| and |${\lambda}_2$| are the tunable parameters to help balance the contribution of each loss term.

Adversarial loss term

The adversarial loss term is introduced to force the predicted models to follow the same distribution as the “ground-truth” models generated by FEM. The original GANs commonly used sigmoid cross-entropy as the adversarial loss, which suffers from the vanishing gradient problem as it can only penalize the objects at the wrong side of the decision boundary. When the fake samples are at the right side of the decision boundary, the gradient is vanishing, whereas the fake samples could still be far away from the real data, making a limited accuracy and an unstable training procedure. To handle these challenges, many loss functions have been proposed, i.e. least-square adversarial loss (Mao et al. 2017) and Wasserstein distance-based loss (Arjovsky et al. 2017), which can help continue pushing the fake samples to the real data even they are already at the right side of decision boundary to provide enhanced performance and stability. To mitigate the vanishing gradient problem, we use the least-square adversarial loss, which illustrated an improved training stabilization compared with the cross-entropy loss and less computational costs compared with Wasserstein distance-based loss (Lucic et al. 2018) is defined as |${L}_{Adversarial}=\frac{1}{2}\mathbb{E}\left[\left(D\left(\boldsymbol{s}\circ G\left(\boldsymbol{s}\right)\right)-{1}^2\right.\right]$|⁠, where |$\mathbb{E}$| is the expectation.

Similarity loss term

To stabilize the training of the generator and improve the convergence of the adversarial network, we perform a similarity loss based on the typically used mean square error, which minimizes the difference between the estimated model and the “ground-truth” model generated by FEM as |${L}_{similar}=\parallel \boldsymbol{s}\circ G\left(\boldsymbol{s}\right)-\boldsymbol{t}{\parallel}_2^2$|⁠.

Smoothness loss term

To guarantee the smoothness of the predicted models, the estimated deformation field should be smooth as well. Therefore, we add a Laplacian operator-based smoothness loss term to supervise the smoothness of the estimated deformation field. The Laplacian operator measures the average value of this deformation field over a small ball centered at each point, and thus is commonly used in estimating the continuous displacement fields in registration tasks to provide explicit guidance of the smoothness constraint. During training, the framework learned to minimize the following loss to ensure the smoothness of the estimated deformation field G(s): |${L}_{Smooth}=\parallel \varDelta G\left(\boldsymbol{s}\right){\parallel}_2^2$|⁠, where |$\varDelta$| is the Laplacian operator.

By minimizing this joint loss function, the proposed framework can be stably trained to estimate the developed models, under the premise of smoothness.

Discriminative loss function

The discriminative loss is implemented to train the discriminator to distinguish the estimated developed models from the “ground-truth” models, which is based on the least square loss: |${\mathbb{L}}_{\mathrm{D}}=\frac{1}{2}\mathbb{E}{\left(D\left(\boldsymbol{f}\right)-1\right)}^2+\frac{1}{2}\mathbb{E}{\left(D\left(G\left(\boldsymbol{v}\right)\right)\right)}^2$|⁠.

Implementation details and data set

We trained the proposed framework on a Dell Alienware PC with a 6 core Intel Core I7-8700K 3.7 GHz CPU, 16 GB NVIDIA RTX 2080Ti GPU, and 64 GB of memory. The training parameters were empirically set as (based on the validation data set): learning rate of Adam optimizer = 0.0002, 1 = 10, and 2 = 1. The 4 × 4 × 4 convolutional kernels with double stride were employed for the generator and discriminator. The input models have a size of 151 × 151 × 151. We padded them to 160 × 160 × 160 to fit for the framework input dimensions as preprocessing. To investigate the performance of the proposed framework, we randomly selected 1,000 models from the simulated database into 5 parts to perform stratified, 5-fold, cross-validation in each experiment. Additionally, we choose the remaining 70 cases as the validation data set. For each patch growth model, the FE computational cost approximately was 5 SUs (1 SU = running a serial job on 1 CPU core for 1 h). The computational cost for the ML part is almost 13 h for each training based on the above computational resource. The trained network can be used to perform the inference, and each inference approximately costs 9 s. To facilitate the reproducibility and further exploration of our proposed ML framework, we have made the source code available on GitHub via the link: https://github.com/cljun27/3-Hinge.

Results and discussion

Our developed computational models are capable of accurately predicting the complex morphology of the folded human brain. Our previous computational-imaging study demonstrated that 3-hinge patterns are present in the convoluted human brain. The 3-hinge patterns predicted by the mechanical models closely resemble those observed in the human brain. Furthermore, our previous work has shown that our models can even capture the dominant shape of 3-hinge patterns observed in MR images. A visual and quantitative comparison between the predicted 3-hinge patterns from our mechanical models and MRI images can be found in Razavi et al. (2015) and Fig. S2 in the Supplemental Materials. In this section, we first analyze the finite element results and study the impact of geometrical undulation on the morphology of the brain model after folding, without considering the ML part. Then, we compare and evaluate the outcomes of both ML and finite element modeling approaches.

Effect of initial surface curvature on brain folding pattern and 3-hinge gyri

Figure 4 depicts the convolution process of a developing brain model at various stages of growth. In the model, the gray matter grows tangentially on the white matter, resulting in compression stress in the gray matter that triggers instability and amplifies any initial undulations. In the state of instability, the system releases its strain energy by transitioning into a post-instability morphology. As the model continues to grow, the peaks and valleys begin to wrinkle and fold, forming gyri and sulci. As documented in the literature, the emergence of instability in the rapidly growing multilayer biological tissues/organs is attributed to the mismatch in the growth rates of different layers (Richman et al. 1975; Tallinen et al. 2016). This phenomenon is valid for the human brain, where the gray matter volume increases 4-fold, whereas the white matter volume increases 3-fold during the 22–30 gestational weeks.

Fig. 4

Morphological evolution steps for a growing bilayer cubic brain model from the top view. The contour is rendered with displacement along the z-axis.

In Fig. 4, the contour shows displacement along the z direction. It is evident that, from the beginning to 60% of the growth, the amplitude of the undulations increases, but there is no gyrification in the gray matter. However, beyond 60% growth, the model loses stability and enters an unstable state that eventually leads to developed folds. With the continued growth, the emerged folds become more convoluted, resembling those of a real mature brain. By comparing the initial and final states of the models, a clear relationship can be observed between the initial geometrical undulations and the location of gyri and sulci after gyrification. Figure 5 further demonstrates this by comparing the initial and final state of 5 cases with different initial undulations. To facilitate comparison, the initial and final state of each model is plotted into a 2D image with topographical contour along the z-axis. It is noteworthy that each pit of the initial undulation transforms into a sulcus, whereas every peak on the initial surface becomes a gyrus. As previously discussed, the amplitude of the undulations increases during early stages of growth. Therefore, every small negative or positive undulation turns into a valley or peak with high amplitude at the beginning. These valleys and peaks subsequently become more convoluted after certain amount of growth, leading to deep sulci and elevated gyri. To demonstrate the consistency of this approach, 4 additional cases with different initial undulations are included in Fig. 5b and e. By tracking the position of valleys and peaks in the initial undulation of each model, the same trend can be observed.

Fig. 5

a) Comparison between the initial and final state of a model regarding the effect of initial undulations on the gyrification of gray matter. Six initial positive (outward) and negative (inward) undulations are highlighted with dashed circles, the dashed arrow indicates the approximate correspondence between these 2 states. b–e) Growth and folding of the brain models for 4 different initial undulations. The upper row displays the models at initial state applied to different undulations; the bottom row shows models in the final state. The contour is rendered with displacement along the z-axis.

By analyzing the location of the positive initial undulations after folding and comparing it with the location of 3-hinges, our investigation indicated that positive undulations not only turn into gyri, but also exhibit a clear association with the location of 3-hinges. In Fig. 6a, the initial state of a model is depicted from a top-down view, with the contour indicating the altitude of gray matter surface. The red regions represent the positive undulations, whereas blue regions correspond to the negative undulations. The red and blue dots indicate the centers of each peak and pit, respectively. Figure 6b shows the final stage of the model, with the contour showing the altitude of the gray matter surface. The red and green dots indicate the location of peaks after folding and the location of the nearest 3-hinges. Our analysis revealed that 3-hinges develop in the vicinity of 95% of positive undulations. In this paper, we considered a 3-hinge is paired with a positive undulation if the distance between them is <6 mm. Using this criterion, we found that 5 out of 6 positive undulations in the initial state had corresponding 3-hinge pairings, whereas 1 was classified as a single. Figure 6c displays the results of pair and single positive undulations in 61 cases. It is clearly demonstrated that most positive undulations have a 3-hinge pair, and some even exhibit a 3-hinge pairing for all of their undulations. The right column of the figure shows the average single and pair positive undulations across all cases.

Fig. 6

a) The initial state of a model with 6 positive undulations (red dots) and 3 negative undulations (blue dots). b) Final state of the model after gyrification. Location of positive undulations after folding (red dots) and location of 3-hinges (green dots). Dashed green circles represent the 3-hinges pairings c) statistical data indicating the correlation between positive undulations and formation of 3-hinges. The contour is rendered with displacement along the z-axis.

ML prediction of surface morphology and 3-hinges

In this subsection, we present a comparison of the ML model predictions with FEM results for 4 different cases. Figure 7 illustrates the altitude contour of the gray matter surface for each case from the FEM and ML models. The first row of the figure shows the initial state of the gray matter surface, which was used for both FEM analysis and ML prediction. The second and third rows display the FEM and ML results, respectively. The comparison between the ML predictions and FEM results demonstrates its ability to accurately predict the position of gyri and sulci based on the initial state, though the trained algorithm may not produce identical results. Comparing the results on a case-by-case basis reveals that the location of high altitudes depicted with red color (gyri) and low altitudes depicted with blue color (sulci) are notably similar in FEM and ML results.

Fig. 7

Comparison of FEM results and ML surface morphology prediction based on the same initial undulations. The contour is rendered with displacement along the z-axis.

To study the correlation between surface initial undulations and gyrification, we clustered all initial undulations from each case into 4 groups based on their amplitude: positive low, positive high, negative low, and negative high. Subsequently, we tracked their location after folding to determine if they were located in gyri or sulci. Figure 8a shows the statistical results for both FEM and ML. From the figure, 99.5 and 100% of the initial positive undulations turn into gyri after folding for the FEM cases with low and high amplitudes, respectively. Similarly, 99.5 and 100% of negative undulations become sulci for the FEM cases with the low and high amplitudes, respectively. Furthermore, in cases with the high amplitudes, almost all positive undulations become gyri and almost all the negative undulations become sulci after folding. The predictions of the ML models show a similar trend to the finite element modeling, with 97.1 and 99.5% of positive undulations with the low and high amplitude, respectively, turn into gyri and 98 and 99% of the negative undulations with the low and high amplitude, respectively, turn into sulci.

Fig. 8

a) Statistical analysis indicating the correlation between the initial undulations and formation of gyri and sulci. b, c) Linear regression of the coordinates of peaks/pits after deformation between FEM and ML results.

Figure 8b and c depicts a comparison between the final position of initial peaks and pits obtained through FEM and ML. To enable comprehensive comparison of all 3 coordinates (“X,” “Y,” “Z”), all the data are normalized based on their minimum and maximum values. The presented results for both peaks and pits show excellent agreement between FEM and ML in “X” and “Y” directions. The slope of the trend line is very close to 1 for both “X” and “Y” data, with r-squared values ranging from 0.96 to 0.99, indicating that the prediction results are highly consistent with the FEM results. Although the accuracy of the prediction data along the “Z” direction is not as precise as that in the “X” and “Y” directions, the trend line and r-squared values still demonstrate good agreement between FEM and ML results.

To evaluate the validity of ML predictions, we examined the local peaks of gyri without considering the initial undulations. Comparison of the local peaks between the FEM and ML results reveals that the ML is capable of predicting the local peaks of gyri. Figure 9a and b shows results of FEM and ML for one case with the local peaks indicated by red dots in both images. Figure 9c compares the local peaks in FEM and ML results, where the blue dots represent local peaks in ML that have a paired peak in FEM, shown as green dots. The red dots show local peaks in ML that do not have a corresponding pair in FEM. We considered a pair when the distance between the closest local peak in FEM and the corresponding local peak in ML is <6 mm. The comparison of the local peak locations indicates that the majority of local peaks in ML have a corresponding peak in FEM. In the presented case, out of the 33 local peaks in ML, 27 have a corresponding peak in FEM, whereas only 6 peaks are single. We also performed a statistical analysis, as presented in Fig. 9d. The blue section of each column indicates the proportion of local peaks that have a corresponding pair in FEM for each case, and the last column shows the average data of all cases. The results show that more than 70% of the local peaks in ML predictions have a corresponding pair in FEM, indicating the ability of ML models to predict the final state of folding.

Fig. 9

a) Local peaks (red dots) of an FEM result after gyrification. b) Local peaks (red dots) of the corresponding ML result after gyrification. The contour is rendered with displacement along the z-axis. c) Pairs of the local peaks in FEM and ML results. Green dots show the local peaks in the FEM result that has a pair in the ML result. Blue dots are the corresponding pairs in ML. Red dots are the local peaks in the FEM result that does not have a pair in the ML result. d) The statistical analysis of pair in 70 cases.

We also performed an ablation study to examine the impact of both adversarial loss and smoothness loss on the prediction of local peak in brain morphology. Figure S5a in the Supplemental Materials displays the results considering both adversarial and smoothness loss function terms, revealing a statistical confidence of 74% in predicting pairs of the local peaks in a brain surface when comparing FEM simulations with ML predictions. In Fig. S5b, we observe the results without the adversarial loss function term, which yield a statistical confidence of 76% in predicting the pairs of local peaks. Similarly, Fig. S5c presents the outcomes without the smoothness loss function term, demonstrating a statistical confidence of 73% in predicting the pairs of local peaks. Our analysis of Fig. S5 underscores that ML models trained without the adversarial or smoothness loss terms are incapable of accurately predicting the morphology of folding patterns. Specifically, the absence of the adversarial loss leads to the generation of results that deviate from the FEM model, displaying inadequate fold complexity. On the other hand, ML models trained without the smoothness loss solely generate scattered and non-smooth surface patterns. Moreover, it is worth noting that the model without the adversarial loss exhibits slightly improved performance in predicting gyral peaks. However, when compared with the FEM model, it falls short in accurately capturing the size and degree of convolution of folds with insufficient complexity. Conversely, the results obtained with the adversarial loss term are more realistic and capture a greater level of detail in the folding patterns. This aligns with previous study (Isola et al. 2017) that the inclusion of adversarial loss was found to enforce similar distributions between the outputs and the ground-truth, thus enhancing the representation of specific details within the ground-truth, rather than producing ambiguous results with globally low loss values.

Fig. 10

a) 3-hinges (red dots) of an FEM result after gyrification. b) 3-hinges (red dots) of the corresponding ML result after gyrification. c) Pairs of the 3-hinges in FEM and ML results. Green dots show the 3-hinges in the FEM result, and blue dots are the corresponding pairs in ML. Red dots are the 3-hinges in the FEM result without a pair in the ML result. d) Statistical analysis in pairs of the 3-hinges in 70 cases, in which blue represents paired 3-hinges and red represents singled 3-hinges.

We further examined the accuracy of ML predictions in determining the location of 3-hinges in comparison to FEM results. The location of the hinges was analyzed without considering the initial undulations. Figure 10a and b illustrates the location of 3-hinges in FEM and ML for a specific case, whereas Fig. 10c shows the comparison between the 2 methods. The green dots indicate the location of 3-hinges in ML that have corresponding 3-hinges pair in FEM, represented by the blue dots. Red dots represent the location of 3-hinges in ML that do not have any corresponding pair in FEM. Similar to the previous analysis, we considered a pair to exist when the distance between the 3-hinge locations is <6 mm. The statistical data in Fig. 10d showed the position of 3-hinges that have a corresponding pair in blue and single 3-hinges in red. The last column shows the average data of all cases. The results show that 76% of the 3-hinges in ML have a corresponding 3-hinge in FEM, indicating the satisfactory prediction of the location of the 3-hinges by the ML model.

Limitation of the model

Our study serves as a pilot investigation, innovatively linking the morphogenesis of the brain with surface undulations. However, due to the simplified nature of the models utilized, they are unable to fully capture the complex interplay between cortical folds and surface undulations that is observed in actual brains. Therefore, there are limitations to our study that should be considered. These limitations provide ample opportunities for further research, some of which are discussed below. One of the limitations of our study is the use of a simple cubic model for constructing 3D models. In reality, the developing brain possesses an irregular curved shape, particularly in the early stages of growth. Therefore, our patch models may not be suitable for investigating the relationship between surface curvature and the shapes of 3-hinge patterns. However, based on our current and previous findings, we have observed the formation of 3-hinge patterns on both flat and curved surfaces (Razavi et al. 2015, 2021; Chavoshnejad et al. 2021), indicating that the mechanism behind their formation is independent of surface curvature. To address this limitation, future research could focus on constructing initial geometries of the model by converting MRI data from developing brains into 3D geometries. This would allow for more realistic and accurate representations of brain development. Another limitation of our study pertains to the use of a fixed thickness for the cortex and specific material properties for both the cortex and the white matter. In reality, the thickness of the cortex varies across different regions and stages of development, and the material properties of these brain tissues are complex and heterogeneous. Therefore, a more realistic approach would be to incorporate nonuniform thickness and varying material properties into the models used for studying brain development. Despite the abovementioned limitations, the current study represents a significant step toward the application of artificial intelligence in predicting the growth and folding of the human brain. By using computational models to simulate the development of the brain, we were able to gain insights into the complex interplay between different factors that contribute to cortical folding. Our findings provide a foundation for future research aimed at refining these models and incorporating additional factors to further improve our understanding of brain development. Ultimately, the use of artificial intelligence in this context holds great promise for enhancing our ability to predict and diagnose brain disorders, and for developing more effective treatments for these conditions.

Conclusion

Data-driven computational simulations have emerged as a valuable tool for studying the neurodevelopment of the human brain by complementing neuroimaging data and providing a flexible and versatile approach for exploring the mechanical factors involved in cortical folding. However, the efficiency of this approach is limited by computational inefficiencies and resource constraints, particularly when dealing with massive simulation data. In this paper, we presented a robust and high-throughput ML-based finite element method (FEM) surrogate model that can predict brain development and its relevant 3-hinge gyral characteristics for any given fetal brain surface model. With massive data from FEM simulations, we trained a GAN-based ML algorithm to predict the brain surface morphology. Our statistical analysis indicates that the location of high and low altitudes in brain patch models are comparable between FEM results and ML predictions, and the ML algorithm is capable of predicting the location of 3-hinges with more than 75% accuracy compared with FEM. This study represents a significant step toward developing a surrogate computational model for better understanding brain development and other tissue modeling, and highlights the potential of ML in addressing the computational challenges faced in this field.

Author contributions

Poorya Chavoshnejad (Data curation, Formal analysis, Methodology, Writing—original draft), Liangjun Chen (Data curation, Investigation, Methodology, Writing—original draft), Xiaowei Yu (Data curation, Formal analysis, Software), Jixin Hou (Methodology, Writing—original draft), Nicholas Filla (Data curation, Visualization), Dajiang Zhu (Funding acquisition, Investigation, Supervision, Writing—review & editing), Tianming Liu (Conceptualization, Funding acquisition, Writing—review & editing), Gang Li (Conceptualization, Funding acquisition, Writing—review & editing), Mir J. Razavi (Conceptualization, Data curation, Funding acquisition, Writing—review & editing), and Xianqiao Wang (Conceptualization, Funding acquisition, Investigation, Supervision, Writing—original draft, Writing—review & editing)

Funding

This work is partially supported by National Science Foundation (IIS-2011369 and CMMI-2123061) and National Institutes of Health (R01AG075582, RF1NS128534, RF1MH123202, and R01MH116225).

Conflict of interest statement

The authors declare no competing financial interests.

Author notes

Poorya Chavoshnejad, Liangjun Chen with equal contribution.