## Abstract

Over the past decades, the buckling instability of layered materials has been the subject of analytical, experimental, and numerical research. These systems have traditionally been considered with stress-free surfaces, and the influence of surface pressure is understudied. In this study, we developed a finite element model of a bilayer experiencing compression, and found that it behaves differently under surface pressure. We investigated the onset of buckling, the initial wavelength, and the post-buckling behavior of a bilayer system under two modes of compression (externally applied and internally generated by growth). Across a wide range of stiffness ratios, 1 < μ_{f}/μ_{s} < 100, we observed decreased stability in the presence of surface pressure, especially in the low-stiffness-contrast regime, μ_{f}/μ_{s} < 10. Our results suggest the importance of pressure boundary conditions for the stability analysis of bilayered systems, especially in soft and living matter physics, such as folding of the cerebral cortex under cerebrospinal fluid pressure, where pressure may affect morphogenesis and buckling patterns.

## 1 Introduction

Layered materials can become unstable under compressive strain, forming different buckling patterns such as wrinkles [1,2], folds [3,4], and creases [5–8]. The formation of wrinkles and folds is sometimes considered a mode of failure, but in some cases, it is actually a desirable behavior that enhances functionality [9–12]. For example, many biological systems experience instability throughout their development or during some functions, including cell membranes, leaves, skin, intestines, and the lungs [10,13–16].

The brain is another example of a layered material experiencing instability, resulting in prominent peaks, called gyri, and valleys, called sulci. The brain is composed primarily of two tissues, gray matter and white matter. The outer gray matter layer of the brain, also known as the cortex, begins its transition from smooth to highly folded during the third trimester of gestation. Neuron cell bodies, responsible for the brain’s information processing abilities, are primarily housed in the gray matter, while axons that connect neurons spread throughout the white matter. The cortical folding process, also called gyrification, appears to stem from instabilities induced by compressive stresses arising in the gray or white matter. Several theories have been developed to describe the mechanism of cortical folding, including different growth rates in cortical layers [17] and tension in white matter axons [18]. Understanding the process of cortical folding is crucial because folding patterns are associated with the functionality of the brain. Abnormalities in gyrification have been associated with epilepsy, lissencephaly, polymicrogyria, and autism spectrum disorder [19–22].

Instabilities in bilayered structures have been investigated using various approaches including theoretical, experimental, and numerical studies [23–25], including both finite element and incremental methods [26–28]. The influence of compression modes—e.g., due to external constraints, growth, prestretch, or shrinking—on the critical strain has also been studied both analytically [29,30] and numerically [31]. Computational models of the developing brain have led to a comprehensive understanding of internal factors governing the formation of gyri and sulci, and have highlighted the importance of growth rate, cortical thickness, and stiffness ratio between gray and white matter in brain folding [25,32–37].

The main focus of this research has been the instability behavior under a zero-stress boundary condition. However, many layered structures are subjected to a pressure boundary condition, either constant or variable, from huge geological formations (e.g., water pressure on the Earth’s crust) to barely macroscopic biological systems (e.g., air pressure on the lungs and fluid pressure on the digestive tract). Inside the skull, the brain is subject to pressure from the cerebrospinal fluid (CSF). CSF is a body liquid that circulates in the space between the brain surface and meninges, protecting the brain from external damage. Normal CSF pressure is from 0.2 kPa to 0.8 kPa in an infant, and 1.3 kPa to 2 kPa in an adult [38,39], although it can increase up to 6 kPa [40,41] in severe cases. Several neurodevelopmental disorders are associated with abnormalities in CSF circulation, such as hydrocephalus, which is an abnormal accumulation of CSF that results in larger ventricles. Excessive pressure over the brain surface in hydrocephalus causes changes in cortical pattern, including flattened gyri and narrow sulci regions. Other neurodevelopmental disorders associated with abnormal CSF also show variations in cortical thickness and significant changes in folding patterns [42–45].

Most models of cortical folding have primarily concentrated on intrinsic factors, assuming that gray and white matter act as an isolated bilayer system [35,46,47]. Indeed, while an early theory hypothesized that cortical folding happens because of the constraint of the skull [48], experiments demonstrated that the folding of the cortex does not rely on a disproportionate increase in the size of the brain with respect to the capacity of the cranial cavity [49]. Despite this, more recent research has suggested that factors outside the brain itself can also affect cortical folding. Several studies have shown that the skull constraint influences the shape of the cortical layers, flattening them as they make contact with the skull [25,50]. Beyond the skull, the mechanical properties of the pia mater have also been found to influence cortical folding, increasing the wavelength [51]. Finally, including CSF pressure improves the predictions of absolute Gaussian curvature compared to a brain bilayer model in isolation [52].

Recently, some studies on cortical folding have moved beyond a narrow focus on instability characteristics, such as critical strain and wavelength, and have investigated more subtle features of the brain. For instance, during brain development, a consistent thickness difference between the gyri and sulci emerges. This phenomenon has been the object of study for nearly a century [53], has continued to attract interest. Our recent work has shown that this is the result of the mechanical forces acting on the cortex during growth and folding [36,54], and is further affected by increased growth in gyri [55]. Cortical thickness is an important biomarker of neurological health, and small regional variations in cortical thickness have been associated with disorders such as autism spectrum disorder [21].

We have previously studied the influence of surface pressure on instabilities via a linear instability analysis of an inhomogeneous bilayer with pressure applied to the top surface [56]. Our results showed that a bilayer with surface pressure is always more unstable than an equivalent system in isolation. While the effect of external surface pressure on instability has been theoretically investigated in cylindrical [57] and spherical geometries [58], a comprehensive numerical study in flat systems has not been conducted. In this study, we consider the cerebrospinal fluid as an external constraint over the brain surface, and investigate the stability criteria and post-buckling behavior of the brain bilayer using a numerical method. Finally, we aim to understand how the cerebrospinal fluid pressure influences thickness variation throughout the brain.

## 2 Method

### 2.1 Mathematical Model.

**, which undergo a motion**

*X**φ*to their corresponding location in the deformed configuration,

**=**

*x**φ*(

**,**

*X**t*). The deformation gradient is given by

**=**

*F*

*F*^{e}

*F*^{g}. Here, $Fg$ irreversibly maps the initial state to the stress-free grown intermediate configuration representing how the body changes through the addition or removal of material, while

*F*^{e}elastically deforms the intermediate state to the final stressed state to ensure compatibility.

*L*and

*μ*are the Lamé constants, and $Ce=FeTFe$

*C*^{e}=

*F*^{e}

^{T}

*F*^{e}is the elastic right Cauchy–Green deformation tensor. The Cauchy stress is thus given by

*B*^{e}=

*F*^{e}

*F*^{e}

^{T}is the elastic left Cauchy–Green deformation tensor.

### 2.2 Computational Model.

We implemented our computational model in commercial software abaqus/explicit [60] with an explicit dynamic solver for large strains. To ensure the quasi-static equilibrium of the system, the kinetic energy needs to be less than $5%$ of the internal energy. To ensure this, we define a small viscous pressure, 1 × 10^{−5} kPa, compared to the values for the external pressure, at the top surface of the bilayer. Additionally, both film and substrate are given a density of 1 kg/m^{3}. In order to approximate incompressibility, we considered the bulk modulus of the materials to be much larger than the shear modulus (*K* = 1000*μ*).

We created a three-dimensional finite element model of a bilayer rectangle with ℓ = 100 mm, *H* = 20 mm, *W* = 0.25 mm. The initial thickness of the film is *T* = 1 mm, (Fig. 1). For both cases, the film’s top surface is allowed to make sliding, frictionless, contact with itself without penetration. The non-zero stress boundary is added as a follower surface traction in the *X*_{2} direction, and the values are normalized by the film stiffness as $P/\mu f$.

*X*

_{3}direction. To trigger the buckling over the surface, a small imperfection, either in the geometry or the material properties, is needed [55,61]. For our simulations, we introduced a geometric imperfection in a small region (2% of the length) in the middle of the domain by slightly changing the

*X*

_{2}coordinates of nodes on the interface of the two layers according to

*y*coordinate,

*X*

_{2}is the original

*y*coordinate, and

*u*is the magnitude of the imperfection as a fraction of the film thickness. Our results were found to be independent of the magnitude, with

*u*ranging from 0.02 to 0.1, and the location along the length of the model. The following results use a perturbation of 2% in the center of the domain.

We modeled two modes of compression, whole-domain compression and film growth (Fig. 1). For the case of whole-domain compression, there is no growth (i.e., the growth tensor is $Fg=I$ throughout the entire domain). Compression is applied on one side, reducing the length to *λ*ℓ, while the other is fixed in the *X*_{1} direction due to the symmetry of the problem. The front and back faces are fixed in the *X*_{3} direction, and the bottom is fixed in the *X*_{2} direction. The axial strain in the domain is then calculated as $\epsilon =1\u2212\lambda $.

**represent the growth parameter and the unit vector normal to the plane of growth in the reference configuration, respectively. The growth parameter quantifies the area change in the growing film, $detFg=\u03d1g$. To establish a simple correspondence between simulation time and volume growth, we assume the evolution of $\u03d1g$ is constant in time, such that $\u03d1g=1+\u03d1\u02d9gt$. The axial strain in the film is then calculated as $\epsilon =1\u22121/\u03d1g$. The substrate is purely elastic, i.e., $Fg=I$. Finally, roller conditions are assigned to the side, front, and back faces, and the bottom is fixed in all directions.**

*n*_{0}We explored a wide range of stiffness ratios, $\beta =\mu f/\mu s=$$[1,1000]$, along with different values of normalized pressure, $P/\mu f=[0.5,4]$, for both compression modes. We found the critical strain for each stiffness ratio by looking at the total strain energy of the film; the strain energy reaches a peak value when the bilayer becomes unstable, and decreases after the buckling point. The initial wavelength was calculated by counting the number of apparent waves, $kc$, shortly after the buckling point, and then normalized by the film thickness, $Lc/T=L/kc/T$, throughout our simulations. Finally, thickness was calculated at each node on the top surface as the minimum distance to a node on the film-substrate interface. To identify gyral and sulcal points, a horizontal threshold was identified at each time point that divided the nodes into roughly half above and half below. From this, we calculated the average gyral and sulcal thicknesses and the thickness ratio, $tg/ts$.

## 3 Results and Discussion

### 3.1 Critical Strain.

To verify our results, we first compared them against the analytical solution of the zero-pressure case [29], finding good agreement (Fig. 2). To investigate the effect of surface pressure on the stability of our bilayer system, we then further compared the critical strains in the zero-pressure cases with various normalized pressure values (Fig. 2). As expected, adding pressure accelerates the initiation of buckling. This is especially true for softer films; when adding a high pressure ($P/\mu f=4$) to a bilayer with low stiffness ratio (*β* = 3), the critical strain decreases by $28%$ and $13%$ for the whole-domain compression and film growth, respectively. However, increasing the stiffness ratio reduces the pressure effect gradually.

Increasing the pressure makes the system more unstable in both compression modes, but it clearly has a more significant impact on the whole-domain compression (Fig. 2(a)) models while the effect on growth-induced compression models is much smaller (Fig. 2(b)).

Also, we aimed to compare the analytical data from Ref. [56] with our finite element results for the whole-domain compression case. The general trend is similar in both sets of results, with stiffness ratio playing a dominant role, correlating inversely with buckling point, and increasing pressure further reducing the stability. However, that work found a discontinuity in critical strain when the pressure is more than two times the film’s stiffness and suggested further exploration of this phenomenon in numerical simulations. Our simulations did not reveal the discontinuities predicted by Ref. [56]. As those discontinuities occurred when the system was very nearly, but not quite, unstable at strain values that caused instability for nearby parameter combinations, it is likely that, in both numerical simulations and physical systems, natural or prescribed perturbations push the system over into instability.

### 3.2 Critical Wavelength.

To investigate the effect of surface pressure on the resulting buckling characteristics of our bilayer system, we compared the initial buckling wavelengths in the zero-pressure cases with various normalized pressure values (Fig. 3). Our results agree with the analytical results from Ref. [29], showing that a higher stiffness ratio leads to a larger normalized wavelength for both compression cases. Adding the pressure to our bilayer model results in slightly smaller normalized wavelengths compared to the zero-pressure bilayer, in agreement with analytical predictions from Ref. [56]. However, these effects are more modest than those seen in the critical strain; for example, for relatively small stiffness ratios (*β* = 3), high pressure ($P/\mu f=4$) results in a $15%$ and $5%$ decrease in the critical wavelength for whole-domain compression (Fig. 3(a)) and film growth (Fig. 3(b)), respectively.

### 3.3 Thickness Ratio.

Cortical thickness is a biomarker of neurological health [20,22,43]. Generally, gyral regions are thicker than sulcal regions in the brain [36]; this can be measured by the thickness ratio, $tg/ts$, which is generally greater than one. Here, we measured the thickness ratio by averaging the thickness of the apparent gyral peaks and sulcal valleys, based on the vertical position of the nodes. We changed the threshold in each time point and each simulation in order to divide the film into roughly half gyral and half sulcal nodes.

As previous results have shown that thickness ratio decreases with increasing thickness ratio, here we consider a soft film, with *β* = 3, in the case of growth only. As growth takes place, the average thickness increases until buckling occurs and the thickness abruptly decreases (Fig. 4(a)). When applying surface pressure, instability happens earlier; then as pressure increases, the average thickness across the domain consistently decreases, as pressure on the surface will resist the film’s expansion in that direction after buckling. When considering gyral and sulcal thicknesses separately (Fig. 4(b)), we see that both gyri and sulci thin after buckling, but that sulci thin much more. As pressure is added, this trend is exacerbated, with higher pressure causing larger drops in sulcal thicknesses. However, after a secondary instability, gyri thin further and the two thicknesses begin to converge.

As seen previously [54,55], the thickness ratio starts from one (representing a film with uniform thickness) and gradually increases to reach a maximum value before decreasing (Fig. 4(c)). The dramatic thinning of sulci under pressure causes the thickness ratio to increase as pressure increases; throughout the course of growth, higher pressures correspond to higher thickness ratios.

Finally, we investigated the changes in thickness ratio across a large range of stiffness ratios (Fig. 4(d)). Previous theoretical and numerical predictions [36,55] suggested that this ratio decreases with increasing stiffness ratios; in essence, this emergence of thickness inhomogeneities is solely a phenomenon of soft films on soft substrates (*β* < 10). Our results agree with this, showing that the gyral-sulcal thickness ratio decreases as the stiffness ratio increases, even in the presence of surface pressure.

### 3.4 Wrinkling Morphologies.

Several studies have investigated the formation and evolution of various instability patterns, and the effects of stiffness ratio, thickness ratio, growth rate, and other parameters [55,62–64]. Here, we sought to investigate the role of surface pressure in wrinkling patterns for both whole-domain compression (Fig. 5(a)), and film growth (Fig. 5(b)) shortly after the initial bifurcation (critical) point.

In both cases of compression, as the normalized pressure increases, the pattern shifts from sharper creases, with deep sulcal fundi and wide gyral crests, to more symmetric wrinkles. Additionally, the wave number increases with increasing pressure, and the wavelength becomes smaller. Our results suggest that emergent morphologies are sensitive to changes in boundary conditions for low stiffness ratios, and increasing the pressure ratio modifies the instability pattern.

## 4 Conclusion

The instability behavior of bilayer systems of different film-substrate stiffness ratios has been widely investigated; however, few studies have addressed the influence of surface pressure on stability. Our study comprehensively examines the buckling and post-buckling behavior of a film-substrate bilayer subjected to two different compression modes under varying amounts of applied surface pressure. Although we have included a wide range of stiffness contrasts between the film and substrate, we mainly focused on bilayer materials with smaller stiffness regimes to recapitulate the effect of the cerebrospinal fluid pressure on brain tissue instability. We showed that increasing pressure always decreases the stability of the system, especially for bilayer systems with similar stiffnesses in the film and substrate, and also slightly decreases the wavelength. These effects are most pronounced when the contrast between the film and substrate stiffnesses is less than ten, losing their influence when the stiffness ratio increases. In addition, surface pressure decreases the average film thickness, and increases the gyral-sulcal thickness ratio by causing sulci to thin more. Finally, the addition of surface pressure leads to the formation of shorter waves and more symmetric patterns with wider sulcal fundi, sharper gyral peaks.

## Footnote

## Acknowledgment

This work was supported by the National Science Foundation, USA Grant Nos. IIS 1850102 and CMMI 2144412. The first author would like to thank Professor Lisa Oglesbee for giving feedback on writing the first draft of this paper.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The data and information that support the findings of this article are freely available.^{2}