A Data-Driven Approach to Morphogenesis under Structural Instability

We construct digital libraries of morphological patterns that are then used to identify abnormal patterns, predict their development, and assist in diagnosis and risk reduction (Fig. 1b,c). To model brain morphologies, a representative core-shell geometry is used for its simplicity [28, 16, 29, 3], as implemented to explore the mechanical instability in cortical folding [15, 30, 13, 31, 32, 14]. The outer spherical shell represents the cerebral cortex (gray matter), and the inner core for the white matter. Following the experimental evidences [33, 34, 35, 36], the cortical thickness ranges from $0.03-1.63$ $\rm{mm}$ according to the abnormal and normal human cerebral cortex measurements and the scale factor [33, 34], and the relative shear modulus ($\mu_{\rm grey}/\mu_{\rm white}$) ranges from $0.65-1$ [35, 36]. The tangential growth (TG) model is used to simulate the cellular mechanisms that create the growth stresses and lead to the pattern evolution (see Methods for details) [28, 15, 9].

In aerospace engineering design, thin plates and cylindrical shells are commonly used in primary structures of aircraft and rockets for structural efficiency [23, 24]. The Johnson-Cook model of plasticity is used for the materials [37, 38]. Axial compression is applied to trigger buckling of the thin-walled structures [39].

To resolve the morphological complexity and for diagnosis purposes, the first step is to recognize the geometrical features during brain development, such as the GI and mean curvatures. We discretize the brain surfaces into a tetrahedral mesh with a density of $10^{6}$ tetrahedra/cm³ [28, 15, 16]. The convergence of calculated descriptors is confirmed. The value of GI measures the depth and width of the cortical folding [40] and is defined as the areal ratio of the cortical surface to the smooth convex hull (Fig. 2a),

$${\rm GI}=\frac{S_{\rm cortex}}{S_{\rm convex}}.$$

(1)

The mean curvatures are defined as the average of the two principal curvatures ($k_{\rm max}$ and $k_{\rm min}$) (Supplementary Fig. 2a),

$$\kappa=\frac{k_{\rm max}+k_{\rm min}}{2},$$

(2)

which quantify the local corrugation of cortical folding. The average value of absolute curvature, $\lvert\kappa\rvert$, is defined as a key geometrical indicator [16]. The variance in $\kappa$ increases as the cortical convolution becomes stronger. In addition to the geometrical measures, topological features are also crucial for neuronal message-passing through chemical and electrical signals across the gyri-to-gyri contact, and vital for brain morphologies and functions [41]. A localization factor ($L_{\rm F}$) is thus defined as the mean distance of contact mapped into the reference configuration to characterize the topology of contact in cortical convolution [42], that is

$$L_{\rm F}=\frac{\sum_{i}\alpha_{i}}{N},$$

(3)

$$\alpha_{i}=\frac{\sum_{j}D_{ij}}{M_{i}},$$

(4)

where $D_{ij}$ is the distance between the nodes $i$ and $j$ in the mesh grid, $M_{i}$ is the number of contacts for the node $i$, $j$ presents all nodes in contact with $i$, $\alpha_{i}$ represents the average contact distance of the node $i$, and $N$ is the number of all nodes that are in contact with others (Fig. 2b).

The combined set of parameters (GI, $\kappa$, $L_{\rm F}$) is used to characterize the brain morphology. Unsupervised clustering of simulated brain morphologies in the space of GI and $\kappa$ is summarized in Fig. 2d. Four distinct types are identified and labeled as ‘convolution-free’, ‘small folds’, ‘thick convolutions’, and ‘normal convolutions’, which aligns well with the experimental evidence of lissencephaly, polymicrogyria, pachygyria, and normal brain, respectively (Fig. 2d) [4]. The labels are then fed to supervised learning using PointNet. The confusion matrix of the well-trained classifier which is summarized in Fig. 2e,f, Supplementary Fig. 3 shows that the accuracy of morphological identification exceeds $95\%$.

Interestingly, we find that the complexity of pattern development can be coarse-grained according to low-dimensional features such as the network of grooves and ridges in the brain morphology [43]. Compared to the original point-cloud representation, the reduced models constructed from the boundaries between grooves and ridges can be used to develop evolutionary dynamics models to understand the underlying biophysical principles. Local patterns such as the nodes are identified in the groove-and-ridge network (Supplementary Fig. 1b-d), which are nucleated, grow, and merge, showing distinct bifurcation behaviors.

To efficiently represent the pattern development, we adopt the PointNet autoencoder by considering rotation and permutation invariance. The learned representations are then fed to the LSTM network to include the historical dependency in the latent space (Fig. 1c). The predicted morphologies at specific stages of development are shown in (Supplementary Fig. 2e,f) for two representative samples with different cortical thicknesses. The model predictions are evaluated by the Chamfer distance, which measures the discrepancy between predicted morphologies at different values of observation time $t_{\rm obs}$ and the ground truth [44]  (Fig. 2h, Supplementary Fig. 2b, c). As $t_{\rm obs}$ increases, more information is extracted and the model corrects earlier predictions, e.g., from cortical convolutions to smooth appearance which forecasts the risk of lissencephaly (Supplementary Fig. 2b,e). The models also modify the local patterns continuously, demonstrating the capability to resolve the fine morphological features during the evolution (Fig. 2g,h, Supplementary Fig. 2f).

Pattern development from smooth to highly curled surfaces involving buckling of the thin shells, which is related to the malformation of the brain (Supplementary Fig. 1b) resulting from migration and growth of neurons as simulated by FEA using the TG model. Identifying the bifurcation behaviors in the morphological changes (e.g., the cortical convolutions) is of vital importance for prognosis purposes (Fig. 2i, Supplementary Fig. 2d). Our results show that bifurcation behaviors are identified from the evolution of the mean value of $\lvert\kappa\rvert$ (Supplementary Fig. 2a). Their features are further analyzed through the descriptors such as GI and local $L_{\rm F}$, as well as the DEVs, which characterize the occurrence and nature of bifurcation (Fig. 2a-c, Supplementary Fig. 4,5, see Methods for details) [45]. Primary and secondary bifurcation are identified in Supplementary Fig. 5. It should be noted that the mapping between the simulation time of brain development and the gestational age is non-linear [16]. The time series of the global indicator, GI, and DEVs analysis indicate a primary fold bifurcation at $t=0.725$ (Fig. 2a, Supplementary Fig. 5a,b, see Methods for details). On the other hand, the time series of the local indicator, $L_{\rm F}$, and DEVs analysis indicate a secondary fold bifurcation at $t=0.937$ (Fig. 2b, Supplementary Fig. 5c,d). Impressively, we find that our framework forecasts the bifurcation behaviors at $t=0.428$ through the Chamfer distance between the prediction and the ground truth (Fig. 2h, Supplementary Fig. 2b,c).

Moreover, similar to the coarse-grained representation of morphologies in the configurational space, the complexity of developmental patterns in the time domain can be addressed by dimensionality reduction according to the occurrence of bifurcation. Determining the nature of branching could advance our understanding of the physics that governs the spatiotemporal evolution of complex systems (e.g., for brain development).

The well-trained models (recognizer and predictor) based on the digital libraries of the brain can be applied to clinical diagnosis of malformations and abnormal development (Fig. 1c, Supplementary Fig. 3a,b, Supplementary Fig. 6a). The point clouds reconstructed from MRI can be fed to the models to assess the current state and predict its future evolution (Fig. 1c, Supplementary Fig. 6a). Once transition into a critical branch in the bifurcation diagram is predicted (such as normal brain to lissencephaly, polymicrogyria, or pachygyria), medical intervention may need to be imposed in conjunction with diagnosis by medical experts, forecasting the event and type of bifurcation is thus highly valuable for the diagnosis purposes (Fig. 4a,c).

Predicting post-buckling behaviors allows to develop instability-tolerant design of aerospace structures (Fig. 3, Supplementary Fig. 6b). Once the load increases beyond the buckling limit, plates and cylindrical shells lose their stability through the primary bifurcation, and weak wrinkles form on their surfaces (Fig. 3a). As the load continues to increase, a secondary bifurcation takes place before the collapse. More significant and irregular crumples form and the structures can no longer withstand the loads (Fig. 3b, Supplementary Fig. 6b). In conservative design, limit and ultimate loads are set to be lower than the first buckling load to avoid structural failure [7]. In practice, one can alleviate the design constraint by permitting post-buckling, which reduces the weight and prolong the lifespan of structural components. Based on accurate and efficient simulations, the European project Improved POstbuckling SImulation for Design of Fibre COmposite Stiffened Fuselage Structures (POSICOSS) lift the limit and ultimate loads above the first buckling load for the fuselage structures, which were validated by comprehensive experimental databases and derived design guidelines [26]. The follow-up project Improved MATerial Exploitation at Safe Design of COmposite Airframe Structures by Accurate Simulation of COllapse (COCOMAT) further accounts for the collapse of structures by considering material degradation [27]. However, these strategies are enforced in the design stage. An instability-tolerant design proposed in this work can allow for bifurcations with weak effects on the load-bearing capacity (Supplementary Fig. 6b), and the tolerance of subcritical wrinkles in service reduces the time and cost of maintenance (Fig. 3a).

In addition to the morphological evolution, internal stress also arises in the structures as driven by mechanical instability. A recent work studying 2D macromolecular membranes shows that including the lattice strain and topological contact in the machine-learning process significantly improves the performance of morphological classification and identification [42]. In contrast, the outstanding performance ($>$ $95$%) reported here using the geometric information only (Fig. 2e and Supplementary Fig. 3c) suggests that there is an intrinsic relation between morphology and internal stress, while the morphology of 2D macromolecules is sensitive to thermal fluctuation, which is externally perturbed [46].

The capability to identify bifurcation behaviors helps to construct phase diagrams of structural stability. Morphological development in complex structures can switch from one branch to another, which is prohibited in traditional path-following methods (Fig. 4c) [47]. Fold bifurcation leads to a catastrophic shift, period-doubling results in fluctuation, and Hopf bifurcation drives the system to an oscillatory state (Fig. 4a, Supplementary Fig. 6c). Identifying the nature of bifurcation based on dimensionality reduction and LSTM can thus guide the diagnosis and structural health management (SHM) to minimize potential losses [48]. The phase diagram of structural stability explored in the current model is limited by the constructed digital libraries and can be extended by using generative models [49, 50].

Pattern development under external forcing often experiences uncertain or rare events. For example, the initial geometry and intrinsic wrinkles may vary from sample to sample. The cortical thickness may change as a result of neuronal cell migration. Genes or genomic mutations also modulate the molecular and cellular processes and brain morphology [51, 52]. For aerospace engineering, bird strikes, conflicts between airplanes, collisions between space debris and satellites, and extreme weather conditions such as lighting strikes or wake vortex circulations also result in uncertain factors that may lead to unexpected morphological pattern evolution and structural failures [53]. These rare events  (Fig. 4b) can be handled in the framework by enriching the datasets, using active learning to improve the efficiency of sampling, or reweighting the samples in calculating the loss function [54]. However, the relatively smooth predictions from data-driven models can underestimate the impact of rare events [55].

Brain development is regulated by genetic, molecular, cellular, and mechanical factors across multiple spatiotemporal scales [8, 9], and the differential tangential growth hypothesis is commonly used [15, 8, 9]. Finite element analysis (FEA) can model morphological evolution during brain growth at the continuum level [15, 28, 13, 31, 16]. In the TG model, the tangential growth of the outer gray matter is faster than the inner white matter [15]. Compression resulting from the mismatch in deformation may then leads to mechanical instability of the brain surface, forming characteristic sulci and gyri structures [28, 15, 30, 13, 16, 31, 32].

In continuum modeling, the reference configuration can be mapped to the current one through the deformation gradient tensor as

$$\textbf{F}=\textbf{F}^{\rm e}\cdot\textbf{G},$$

(5)

where $\textbf{F}^{\rm e}$ is the elastic deformation gradient and G is the growth term. In the TG model, the growth tensor G is

$$\textbf{G}=g\textbf{I}+(1-g)\hat{\textbf{n}}\otimes\hat{\textbf{n}},$$

(6)

where $\hat{\textbf{n}}$ is the surface normal of the reference configuration, I is the unit tensor, and

$$g=1+\frac{\alpha_{t}}{1+e^{10(\frac{y}{T}-1)}},$$

(7)

is the growth coefficient, where $\alpha_{t}$ controls the magnitude of local cortical expansion. There is a smooth transition from the surface of the gray matter layer to the white matter layer with a gradually decreasing growth coefficient. $y$ is the distance to the surface, and $T$ is the thickness of the cortex. The brain is modeled as a nonlinear neo-Hookean hyperelastic material, where the strain energy density is

$$W=\frac{\mu}{2}[\rm{Tr}(\textbf{F}^{\rm e}\textbf{F}^{\rm eT})\emph{J}^{-2/3}-3]+\frac{\emph{K}}{2}(\emph{J}-1)^{2},$$

(8)

where $\mu$ is the shear modulus, $J$ is the determinant of Jacobian matrix, $K$ is the bulk modulus. The interaction between surfaces is modeled with an energy penalty via vertex-triangle contact, which prevents the nodes from penetrating the faces of elements [56]. The explicit solver is used to minimize the total (elastic and contact) energy of the quasi-static system. The time step $\Delta t=0.05a\sqrt{\rho/\emph{K}}$ is set to ensure the convergence, where $a$ is mesh size and $\rho$ is mass density [57].

The morphological data obtained from FEA are stored in point clouds, which possess high-dimensional characteristics and redundant information. Experimentally, point-cloud data can be collected from 3D MRI images [58]. Feature extraction and dimensionality reduction are carried out to improve the performance of downstream tasks (Fig. 1c, Supplementary Fig. 7). PointNet is used to extract key features and train models for morphological phase classification and recognition [59]. The K-means algorithm is used for clustering, and the labels are fed into PointNet (Fig. 1c). The historical dependency of data is addressed by using an autoencoder design based on LSTM (Fig. 1c, Supplementary Fig. 8) [20, 60].

DEVs analysis is performed to characterize the bifurcation behaviors based on Refs. [45, 61, 62], which makes use of principal eigenvalues (PEs) of the Jacobian matrix in the reconstructed state space of a univariate time series [45]. A dynamical system approaches the bifurcation point as the Euclidean norm of PEs approaches $1$ (Fig. 4a, Supplementary Fig. 6c) [45, 63]. According to the characteristics of PEs in the complex plane, the type of bifurcation can be determined (Fig. 4a, Supplementary Fig. 6c) [45, 63]. Specifically, a fold bifurcation has a real part of $1$ and an imaginary part of $0$. The period-doubling bifurcation can be recognized by real and imaginary parts of $-1$ and $0$, respectively. The Hopf bifurcation has a non-zero imaginary part and the Euclidean norm of PEs approaches $1$ [45, 63].