Impact of environmental asymmetry on epithelial morphogenesis

Epithelial folding is a universal biological phenomenon in morphogenesis, typical examples being brain gyri, villi of the intestinal tract, and imaginal discs in invertebrates. During epithelial morphogenesis, the physical constraints imposed by the surrounding microenvironment on epithelial tissue play critical roles in folding morphology. In this study, we focused on the asymmetry of the environmental constraints sandwiching the epithelial sheet and introduced the degree of asymmetry, which indicates whether the basal or apical side of the epithelium is closer to the constraint wall. Then, we investigated the relationship between the degree of asymmetry and epithelial folding morphology using three-dimensional vertex simulations. The results show that the folding patterns of the epithelial sheets change from spot patterns to labyrinth patterns and then to hole patterns as the degree of asymmetry changes. Furthermore, we examined the pattern formation in terms of the equation of out-of-plane displacement of the sheet derived from the mechanical energy functional.


Methods
Mathematical models for simulating multicellular dynamics. In the 3D vertex model, a cell is represented as a polyhedron. Because epithelial tissue comprises a group of similar cells closely attached to their neighbors, it can be modeled as a network of connected polyhedrons whose vertices and edges are shared by adjacent ones. The kinematics of tissue sheet deformation can be described based on the locations and movements of the vertices of individual polyhedral units constituting the sheet. The movement of the i-th vertex can be described by where r i represents its position vector, v loc i represents the local velocity vector determined from its current location and those of its neighboring vertices, η is the friction coefficient, U is the potential energy function, and V represents the velocity of the system's center of gravity (CoG). Based on previous studies [23][24][25]28 , the local velocity vector v loc i given in Eq. (1) can be defined as the mean velocity vector of the surrounding vertices: where χ V i j is the indicator function for a subset V i that is the set of all vertices directly connected to the i-th vertex by edges. Here, the number of vertices connected to vertex i is expressed as vertex j χ V i j . In our simulations, this sum for vertex i is always equal to four because the tissues are monolayers and there are three cells surrounding a vertex. This local velocity vector is introduced to satisfy Galilean invariance 28 . There is no noise term in this vertex motion, but there is randomness in cell growth, described below, which is reflected in the equation of motion through the energy function. During morphogenesis, tissue deformation occurs over a longer timeframe than during cell displacements. Therefore, Eq. (1) neglects the effects of inertia on cellular dynamics and predominantly accounts for the effects of viscosity.
Vertex behavior defined by Eq. (1) depends on the potential energy function. This study assumes that individual cells have the following types of potential energy: volume elasticity energy U VE , surface elasticity energy U SE , height elasticity energy U HE , and environmental constraint energy U EN . To account for environmental effects, the tissue kinematics model includes two elastic walls lying parallel to the tissue sheet, one at a distance of l a from the apical surface and the other at a distance of l b from the basal surface (Fig. 1). The total energy function U is defined as follows: Each of these types of energy is defined as follows:

Apical
Basal Wall Wall l a l b Figure 1. Schematic diagram of epithelial tissue and the elastic wall that provides physical constraint to simulate cell proliferation dynamics in our model. The elastic walls are located at the apical and basal sides of the epithelium, and the distance between them is denoted by l a and l b , respectively. where Ʃ i cell represents summation across all cells, and l api i and l bsl i represent the out-of-plane displacement of the center of gravity of the i-th cell from the apical and basal surfaces, respectively. Moreover, k api EN i and k bsl EN i indicate the elasticity coefficients of the wall on the apical and basal sides, respectively. These variables are defined as follows: where k EN represents the characteristic elasticity coefficient of the wall, and S api i and S bsl i represent the apical and basal surface areas of the i-th cell, respectively. In this study, tissue growth is represented by the cell proliferation model 27 . Cell division is represented by the division of polyhedral, and cell growth is represented by changes in V c, eq and corresponding changes in S c,eq ( H c,eq is set constant). Each cell is assumed to divide along the long axis of the cross sectional cell shape normal to apicobasal axis. Each cell is constrained to adopt the shape of prism and to always have a basal and apical side. The cell cycle is represented by the mean cycle period τ cycle and its standard deviation σ cycle . The percentages of the duration of G1, S, G2, and M phases within the cell cycle are set as I , II , III and IV , respectively. This model does not take into account cell removal in this simulation. Table 1 lists all model constants used in this study.
Focusing on the folding structures induced by cell proliferation, we used a flat, homogeneous epithelial monolayer sheet for the initial condition. In the initial state, the tissue consists of 40 × 40 hexagonal prism shaped cells aligned in a regular hexagonal lattice. The model constants are set so that the initial state is a steady state without cell proliferation. In order to analyze the folding structure, we performed cell proliferation simulation under the periodic boundary conditions at the edges of the tissue.

Results
Proximity of the elastic wall and peak-to-peak folding distance. This section presents the results of the 3D sheet morphology simulations obtained by altering the degree of constraint on the out-of-surface deformation. In this run, the same degrees of constraint were assumed on both sides of the epithelial sheet. Specifically, each of the following values was entered simultaneously into both l a and l b in Eq. (7): 0.20, 0.30, 0.40, 0.50, 1.0, 1.5, and 2.0. www.nature.com/scientificreports/ The initial shape of the epithelial sheet and the results of the simulations are presented graphically in Fig. 2. The initial shape is planar as shown in Fig. 2a. As cells proliferate under periodic boundary conditions, buckling of the tissue is induced. This causes the out-of-plane deformed epithelial tissue to collide with the wall. This process is shown in Supplementary Video 1. The sequence in Fig. 2b corresponds to the order of proximity between the elastic wall and the epithelial sheet (i.e., 0.20-1.5). Closer proximity of the elastic wall resulted in a smaller peak-to-peak folding distance.
The degree of asymmetry relative to the XY plane and peak-to-peak folding distance. The simulation results shown in the previous section represent the cases where the same degrees of environmental constraint were applied to the apical and basal sides of the sheet (i.e., l a = l b ). In this study, we evaluated the impact of asymmetric environmental constraint on the folding patterns by changing the degree of asymmetry relative to the xy plane. For this purpose, two variables were introduced: the total distance of the elastic walls from the epithelial sheet l sum (= l a + l b ) and the degree of asymmetry relative to the xy plane Λ (= l b /l sum ). Specifically, l sum had one of the following values: 0.60, 1.0, 1.4, 2.0, or 3.0. For each l sum value, the asymmetry indicator Λ ranged from 0 to 1.
The representative 3D morphology of the apical surface simulated using l sum = 1.0 is illustrated graphically using the emboss effect in Fig. 3a; the sequence represents increasing values of Λ (i.e., 0-1.0). Similarly, the 3D structures of the basal surface simulated using l sum = 1.0 are presented graphically in Fig. 3b, where the sequence represents increasing values of Λ (i.e., 0-1.0). These results demonstrate that increases in Λ from 0 to 0.5 (i.e., from maximum asymmetry to symmetry) resulted in longer (and narrower) folds. Moreover, the folding patterns were similar between the apical side at Λ = p (p representing a given value between 0 and 1) and the basal side at Λ = 1 − p, supporting the mathematical validity of the model.

Environmental constraints and wavenumbers of folds.
The previous section showed that the presence of elastic walls in the sheet kinematics model induced folding of the growing tissue (i.e., emergence of grooves and ridges), and that the folding patterns were dependent on the distances between the walls and the apical and basal surfaces. These findings allow further mathematical considerations of epithelial morphogenesis. In this section, we introduce an indicator for the peak-to-peak folding distance and provide a detailed quantitative analysis of this relationship.
The folding distance indicator I(u x ,u y ) is derived as follows. First, the data regarding the vertices and CoGs of the cells on the apical surface were extracted from the simulation results described above. The area of computation was divided into unit grid cells of the same size, the total number of which was N x × N y , with N x and N y denoting the number of unit grid cells on the x and y axes, respectively. Here, N x and N y were set to 35. The z-axis coordinates of the vertices and CoGs included in a given grid cell were determined, and their average value was defined to represent the out-of-surface displacement. A discrete second-order displacement field z(x,y) was derived by linking the displacements and coordinates (x,y) of the center of the individual unit grid cells. A discrete Fourier transform of the displacement field z(x,y) and its power spectrum I can be expressed as follows: The closer the distance between the epithelial sheet and the elastic wall, the smaller the peak-to-peak distance of folding. where u x and u y represent the wavenumbers over the system's x -and y-axis lengths, X and Y , respectively. The wavenumbers u x and u y were normalized to the characteristic length of the system L = √ X 2 + Y 2 , and u x and u y were defined as follows: The distribution of the power spectrum I(u x ,u y ) thus obtained was used as an indicator of the peak-to-peak folding distance. A study revealed that the distribution of the power spectrum reflects the peak-to-peak folding distance and the orientation of the folds, demonstrating its utility for analyzing folding distance 25 .
Moreover, on the basis of the distribution of the power spectrum derived in this study, characteristic wavenumbers were determined to investigate the relationship between the wall-to-sheet distance and folding distance. Specifically, changes in the characteristic wavenumber resulting from changes in l a and l b were analyzed. Equation (13) provides the average wavenumber u avg based on the power spectrum I(x,y), and u avg was defined as the characteristic wavenumber: Figure 4a graphically presents the apical-side results of the average wavenumber derived using the data described in the previous sections. Each point in Fig. 4 corresponds to a single simulation result. These results supported the finding that a smaller wall-to-sheet distance resulted in a smaller peak-to-peak folding distance. The wavenumbers on the basal side were determined similarly to those on the apical side, and the results are presented in Fig. 4b. The average wavenumbers on the apical and basal sides were generally comparable, with minor differences.

Creating black-and-white images by binarizing Z-axis coordinate values. This section describes
an indicator for the average longitudinal length of folds and provides a detailed quantitative analysis of its relationship with degree of asymmetry Λ. Here, the number of folds formed was chosen as a surrogate parameter inversely proportional to the total sum of the longitudinal lengths of the folds. We adopted this approach because determining the longitudinal length of the fold directly was not technically feasible. One major technical challenge was related to defining the fold's ridge line, which is a key determinant of its longitudinal length. Another major reason for choosing this approach was that we used periodic boundary conditions in this study, and the longitudinal length of the fold was strongly dependent on the size of the simulation system. Consequently, the number of folds was determined as described below. www.nature.com/scientificreports/ First, the area of computation was divided into unit grid cells of the same size, the total number of which was N x × N y , where N x and N y denote the numbers of unit grid cells on the x-and y-axes, respectively (N x = N y = 35). A discrete second-order displacement field z(x,y) was derived by linking the displacements and coordinates (x, y) of the individual grid centers. For each unit grid cell, the Z value was compared with the threshold value max (z(x,y))+min (z(x,y)) 2 , and the unit grid cells whose z values were greater and smaller than the threshold were presented as white and black, respectively. The image processing library openCV was used to delineate and count the white and black areas. Figure 5a presents the binary black-and-white images obtained from the simulations described above, and Fig. 5b presents the number of folds identified based on the data presented in Fig. 3. The graph shows that the number of folds was the smallest at Λ = 0.5 for all simulation runs. In other words, increases in Λ from 0 to 0.5 (i.e., maximum asymmetry to symmetry) led to longer (and narrower) folds.

Mechanical model derived from energy functional.
In order to discuss the pattern formation depending on the degree of asymmetry Λ obtained from the 3D vertex simulation in terms of energy, we defined an energy functional of the cell sheet. Because buckling was caused by compression due to cell proliferation under periodic boundary conditions in the 3D vertex simulations, we considered the cell sheet compressed by force N from x, and y directions. Based on the variational principle, the overdamped relaxation dynamics of the dimensionless displacement field w(t, x, y ) is derived as follows.
where ˜ denotes the dimensionless Laplacian and α is a dimensionless constant. The term with the ramp function R(x) comes from the fact that the stiffness of the elastic wall is expressed in terms of the step function in Eq. (8). The definition of the energy functional and details of the derivation of Eq. (14) are provided in the Supplementary Information.
We numerically integrated this differential equation with a second order Runge-Kutta scheme under periodic boundary conditions with �t = 1 × 10 −3 , �x = �y = 0.6, N x = N y = 128, α = 1000 . The initial state was set with randomly small displacements as shown in Fig. 6a. The displacement field w at t = 2, 000, 000 is shown in Fig. 6b for varying the asymmetry degree Λ from 0.0 to 1.0. The results confirm a transition from a dot pattern to a hole pattern via a labyrinth pattern with increasing Λ as well as that observed in the 3D vertex simulation. Furthermore, to quantitatively evaluate this pattern transition, the displacement field was binarized (Fig. 6c) and the number of outlines was counted, as in the 3D vertex simulation. The results are shown in Fig. 6d. As it shows, the number of contours was large when Λ = 0.0, decreased closer to Λ = 0.5, and increased as Λ was further increased. This trend is comparable to the results of the 3D vertex simulation.

Discussion
We examined the relationship between the epithelial folding patterns and the asymmetry of environmental constraints using the 3D vertex model, and in order to discuss the pattern formation depending on the degree of asymmetry Λ obtained from the 3D vertex simulation, we derived the mechanical model from the energy functional and confirmed that our derived equation could reproduce the results of the 3D vertex model. (14) ∂w ∂t = −� 2w −�w − N(w) ,  According to a previous study 14   www.nature.com/scientificreports/ Examining the stability of the trivial solution u = 0 , we obtain k = −k 4 + k 2 , which means that there is the range of wave number in which order parameter u is not damped, while the linear and cubic terms (third and fourth terms in r.h.s. of Eq. (16)) suppress the increase of order parameter u . In our mechanical model, the terms with the offset ramp functions in Eq. (15) are considered to play the role in suppression of the amplitude of w . Due to its strong nonlinearity term including offset ramp functions in our model, analytical solution of pattern selection is a future challenging work.
In this study, degree of asymmetry was the parameter for pattern selection. From a mathematical point of view, there are other phenomena of pattern formation that can be explained by the SH equation, and the parameters of pattern selection do not necessarily correspond to asymmetry. For example, in the system of spherically shaped www.nature.com/scientificreports/ elastic bilayer materials 14 , the parameter of pattern selection is the effective radius R/h , which is the ratio of the thickness R of the substrate layer to the thickness h of the film layer.
In summary, we investigated the relationship between epithelial folding pattern and environmental constraint on cell displacement using a 3D vertex model that describes morphological changes resulting from cell growth and division. The results revealed that the wall-to-sheet distance was a major determinant of the peak-to-peak folding distance. Furthermore, using a 3D vertex model and a mechanical model derived from the energy functional, numerical simulations showed that the degree of asymmetry with respect to the location of the upper and lower walls relative to the epithelial sheet resulted in different morphological patterns, such as dot patterns, labyrinth patterns, and hole patterns.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.