Epithelial tissue folding pattern in confined geometry.

The primordium of the exoskeleton of an insect is epithelial tissue with characteristic patterns of folds. As the insect develops from larva to pupa, the spreading of these folds produces the three-dimensional shape of the exoskeleton of the insect. It is known that the three-dimensional exoskeleton shape has already been encoded in characteristic patterns of folds in the primordium; however, a description of how the epithelial tissue forms with the characteristic patterns of folds remains elusive. The present paper suggests a possible mechanism for the formation of the folding pattern. During the primordium development, because of the epithelial tissue is surrounded by other tissues, cell proliferation proceeds within a confined geometry. To elucidate the mechanics of the folding of the epithelial tissue in the confined geometry, we employ a three-dimensional vertex model that expresses tissue deformations based on cell mechanical behaviors and apply the model to examine the effects of cell divisions and the confined geometry on epithelial folding. Our simulation results suggest that the orientation of the axis of cell division is sufficient to cause different folding patterns in silico and that the restraint of out-of-plane deformation due to the confined geometry determines the interspacing of the folds.

Introduction

Holometabolous insects undergo a remarkable change in surface shape as they grow from larvae to pupae. The mechanism enabling this change is hidden in the structure of imaginal primordia in the larval body. The imaginal primordium comprises epithelial monolayer tissues that are highly folded in a compact manner. Upon the onset of pupation, these folds spread out to form the three-dimensional (3D) shape of the exoskeletal body including appendages of the insect, such as legs and horns.

Folds of epithelial tissue are not limited to the imaginal primordium but are also found in organs of vertebrates, such as brains (Tallinen and Biggins 2015), oviducts (Koyama et al. 2016), and intestines (Hannezo et al. 2011; Krajnc and Ziherl 2015). Tangential expansion of gray matter constrained by white matter is important to the formation of brain folds (Tallinen and Biggins 2015; Tallinen et al. 2016). Intestinal villi are formed by cell proliferation in a monolayer constrained by the underlying stroma. Such a constrained growth of tissue results in compression of the tissue that leads to a mechanical instability such as buckling. Folded shapes are observed in the typical patterns of engineering materials that have buckled (Chan and Crosby 2006; Guvendiren et al. 2010). Buckling is therefore anticipated to be a mechanism of folding epithelial tissues from the mechanical viewpoint. In fact, the tissue folding mentioned above has been explained by mechanical instability arising from constrained tissue growth (Drasdo 2000; Hannezo et al. 2011; Krajnc and Ziherl 2015; Tallinen and Biggins 2015; Tallinen et al. 2016).

Although there is potentially a mechanism of folding epithelial tissues that is common to insects and vertebrate organs, the main difference between the primordium and vertebrate organs is the function of the folded structure. The epithelial tissue of the primordium is premised to unfold to establish 3D exoskeletal shapes at a subsequent stage while vertebrate organs are designed to remain folded for their functionality. Therefore, the folded structure of the imaginal primordium is not the final form of development. There is thus one more mechanism hidden in the structure of imaginal primordia in the larval body that establishes the final 3D shape from the primordium.

To reveal the mechanism of encoding the final 3D shape in the primordium, the relationship between the 3D shape of the pupa and primordia has been clarified in experiments in which the beetle horn primordium of Trypoxylus dichotomus has been extracted and observed to spread out using several methods (Matsuda et al. 2017). Those experimental results show that there are characteristic folding patterns, such as labyrinth, stripe, and concentric patterns, and revealed what part of the 3D surface shape of the horn is encoded in these characteristic patterns. In addition, because of the imaginal primordia of the Drosophila leg (Morata 2001) and the tip of the beetle horn (Matsuda et al. 2017) have similar (concentric) patterns, there might be a robustness of the relationship between folding patterns and their 3D shapes. However, it is not well understood how these folding patterns of the primordia form in the larval body.

In terms of the development of primordia, because of imaginal primordia are surrounded by other tissues, such as the peripodial membrane (Beira and Paro 2016; Milner et al. 1984), cell proliferation proceeds within a confined geometry. The tissue deformation driven by cell proliferation is thus restrained by the surrounding tissue. Observations of the imaginal primordium of the wings of Drosophila have revealed that cells proliferate with a certain orientation of cell division (Mao et al. 2011, 2013). The regulation of the orientation of cell division is the planar cell polarity signaling pathway. In Drosophila, planar cell polarity molecules such as dachsous molecules function as global direction cues. In response to the dachsous gradient, Dachs, an atypical myosin, is localized in a planar-polarized manner. Dachs controls the long axis of the cell shape on the apical side and thus regulates the orientation of cell division (Mao et al. 2011). Mutation of the dachsous gene results in a wing of Drosophila that is shorter and wider than that of the wild type (Baena-López et al. 2005). However, it has been confirmed that randomization of the axis of cell division does not appreciably affect the adult form of Drosophila (Zhou et al. 2019). In growing beetle horn primordia, the global shape of the horn primordia is mushroom like, with dense local furrows. The global shape is affected in the dachsous-gene knocked-down beetle, in which the direction of cell division has been altered randomly, while dense local furrows are not appreciably affected (Adachi et al. 2018). The exact role of the orientation of cell divisions is therefore not clear in vivo.

During the last decade, mathematical models for the computational simulation of epithelial tissue mechanics have been well established (Hannezo et al. 2011; Heisenberg and Bellaïche 2013; Honda et al. 2004, 2008; Okuda et al. 2013a, b; Takeda et al. 2019) and applied to test hypotheses motivated by biological experiments (Amar and Jia 2013; Eiraku et al. 2011; Fletcher et al. 2014; Inoue et al. 2016, 2017; Okamoto et al. 2013). To suggest a possible role of cell proliferation with a certain orientation of cell division in the confined geometry, we examine whether the above mechanical effects are sufficient to form a proper pattern of folds in silico. We here employ a 3D vertex model (Okuda et al. 2013a) that considers the behavior of an individual cell to express epithelial tissue deformations driven by the activity of each cell. Using the model, we perform computer simulations for several orientations of cell division and magnitudes of the restraint of tissue deformation to examine effects of the orientation of cell division and confined geometry on the formation of folding patterns.

Three-dimensional vertex model expressing multicellular dynamics

The 3D vertex model expresses the shape of a cell as a polyhedron comprising vertices and edges. The tissue is represented as an aggregate of multiple cells (Fig. 1a, b), where the vertices and edges of each cell are shared with neighboring cells. The vertices and edges thus constitute a network that represents the entire shape of the tissue (Fig. 1c, d).

Fig. 1

Shapes of a epithelial tissue and b cells respectively modeled as c a network and d polyhedrons using the 3D vertex model

The equation of motion of vertex i with position vector at time t iswhere is the friction coefficient and is the mean velocity vector of vertex i. To satisfy the Galilean invariance of the equation of motion, the friction force expressed in Eq. (1) is calculated on a local inertial frame. In the 3D vertex model, vertex i is directly connected to four adjacent vertices by edges. Indexing the directly connected vertices as j(i), the mean velocity vector is defined asAssuming that all vertices have the same friction coefficient, Eqs. (1) and (2) are derived from the previous 3D vertex model (Okuda et al. 2015).

The right-hand side of Eq. (1) represents a force acting on vertex i derived from the total energy function U, which represents the mechanical properties and morphogenetic events of the cells. U is expressed aswhere indicates a summation over all cells. The energy function U includes the cell volume elastic energy , cell surface elastic energy , cell height elastic energy , restraint energy with respect to out-of-plane deformation , and self-collision energy . Mathematical expressions for these energy functions areEquations (4–6) are derived from a 3D vertex model of epithelial cells (Inoue et al. 2016, 2017). The volume , surface area , and height of the th cell are represented as variables. The superscript eq used for several variables in Eqs. (4–6) refers to the value in the stress-free state. The constants , , and are, respectively, the volume elasticity, surface elasticity, and height elasticity. The cell height is defined as the distance between the centroid of the apical face and that of the basal face of the cell.

The restraint of out-of-plane deformation is introduced in terms of the harmonic potential that is responsible for the restoring force acting on the th cell with respect to its out-of-plane displacement, , along the z-axis. The mathematical form of is the same as that introduced in (Brau et al. 2013), in which a small out-of-plane displacement of a thin sheet is considered. We assume that growing tissue is in contact with the surrounding tissue on the apical side. We thus use the out-of-plane displacement of the centroid of the apical face of the th cell, , as . Because the restraint energy is defined for the unit area (Brau et al. 2013), the total restraint energy for the cell is obtained by multiplying by the apical surface area. This is why the restraint coefficient is expressed as the product of the restraint constant, , and the surface area of the apical face of the th cell, . The self-collision energy is introduced to prevent simulated tissue from penetrating in the direction of collision using a penalty-based method (Tang et al. 2012), where , , and are the self-collision energy constant, distance between vertices i and n, and threshold distance of the collision, respectively.

All model constants are listed in Table 1. Here, to focus on folded structures created by cell proliferation, we ignore the asymmetry of physical properties dependent on epithelial polarity and assume a flat, homogeneous, epithelial monolayer sheet (Fig. 2) for the initial condition. The size of the sheet is corresponding to cells arranged in a regular hexagonal lattice in the plane under initial conditions. Although we have not performed a pre-simulation for equilibration of the system, the initial conditions have been confirmed to be stable in a planar and hexagonal packing of the columnar cells (animation shown in Supplementary Material 1). The apical surface of the sheet is set to face to the +z direction. Periodic boundary conditions are adopted for x and y directions. To express cell proliferation, we employ a polyhedron-division model (Okuda et al. 2013b), where the timing of cell division is determined by a mean cell cycle time and standard deviation .

Table 1

Model constants

Symbol	Value	Descriptions
k^v	20	Volume elasticity
k^s	0.256	Area elasticity
k^h	0.1	Height elasticity
k^r	0.0–1.6×10^-2	Restraint constant
k^sc	40	Self-collision energy constant
v^eq	1.0	Cell volume at stress free state
s^eq	2v^eq/h^eq+223v^eqh^eq	Cell surface area (hexagonal prism) at stress free state
h^eq	1.0	Cell height at stress free state
σ	1.0	Collision threshold distance
τ_cc	1000	Mean cell cycle time
σ_cc	10	Standard deviation of the cell cycle time
η	0.25	Friction coefficient of vertex
Δt	2.0×10^-4	Time step size for numerical integration of Eq. (1)

Fig. 2

Initial shape of the monolayer cell sheet. The x-axis and y-axis are defined on the plane of the initial sheet, and the z-axis is defined normal to the plane. Periodic boundary conditions are applied for x and y directions. Cells on the boundaries are not visualized

Results

Folding pattern obtained for the orientation of cell division in a confined geometry

We perform 3D vertex simulations of epithelial tissue growth to examine whether tissue growth in a confined geometry produces a folded tissue structure. Simulation results in Fig. 3 show that the restraint of out-of-plane deformation results in a short interspacing of the folds (Fig. 3a, animation shown in Supplementary Material 2) as compared with the interspacing under the condition of no restraint (Fig. 3b, Suppl. Mov. 3).

Fig. 3

Snapshots of epithelial folding simulated under the conditions of a restraint of out-of-plane deformation () and b no restraint of out-of-plane deformation (). All snapshots show simulation results at time

We next investigate how the orientation of the cell division axis affects the folded structure. We consider three orientations of the cell division axis: the longitudinal direction of each cell shape, uniaxial direction in global coordinates, and radial direction toward the centroid of the tissue. The simulation results in Fig. 4 show that (a) random patterns of folds form for the longitudinal orientation (animation shown in Supplementary Material 2), (b) folds align in one direction like stripes for the uniaxial orientation (Suppl. Mov. 4), and (c) folds form a concentric pattern for the radial orientation of cell division (Suppl. Mov. 5). The interspacing of folds is almost the same regardless of the orientation of the axis of cell division (with detailed analysis presented in the next section), suggesting that the orientation of the cell division axis plays a role in determining the pattern of folds.

Fig. 4

Snapshots of epithelial folding simulated for three orientations of the cell division axis: a longest axis of each cell, b uniaxial direction (x-axis), and c radial axis (). All snapshots show simulation results at time

Length scale of folds determined by the restraint of out-of-plane deformation

We examine the effect of the restraint constant, , on the interspacing of folds to determine how the length scale of the folds is affected. The simulation results in Fig. 5 show that an increase in leads to shorter interspacing of the folds. To quantify the interspacing, we analyze the wavenumber, u, of the fold pattern using the Fourier transformation and obtain the mean wavenumber, , as a weighted average of u using the Fourier coefficient.

Fig. 5

Snapshots of epithelial folding simulated under the conditions of a , b, c, d, e, and f. The orientation of cell division is along the longest axis of each cell. All snapshots show simulation results at time

Figure 6 shows the power law relationship between the mean wavenumber and restraint constant. All data points obey a unique power law regardless of the orientation of the axis of cell division, suggesting that interspacing of the folds can be determined by the restraint of the out-of-plane deformation due to the surrounds of the growing tissue.

Fig. 6

Relationship between the mean wavenumber and restraint constant. The legend gives the orientation of cell division. The line shows the least-squares fitting line:

The characteristic wavenumber of folds has been derived theoretically whereby the balance of normal forces under uniaxially compressed sheet registing on a soft substrate is considered (Brau et al. 2013). In the theory (Cerda and Mahadevan 2003; Brau et al. 2013), assuming that the bending stiffness of the tissue and the restraint constant are independent of the size of deformation, the power exponent could be 0.25, while the power exponent is approximately 0.20 in our simulations. Because of the restraint constant in our simulation is independent of the size of deformation, we speculate that a discrepancy between the exponent in simulation and that in theoretical prediction is due to discreetness originating from the cell size. Because of the tissue bending will not be realized on a length scale shorter than the cell width, the wavenumber is upper bounded by the corresponding wavenumber of the cell width. In fact, applying least-squares fitting to data in the small wavenumber range (), we obtained the power exponent of 0.25 (Appendix).

In addition to the discreetness originating from the cell size, because of the bending stiffness of the tissue is not given directly in the simulation but induced by the combination of the energy functions of the cell, there is the possibility that the bending stiffness effectively depends on the wavenumber of the deformation. To explain the power exponent of 0.20 in this simulation using the theory (Brau et al. 2013), it is speculated that the bending stiffness is proportional to the wavenumber if we consider the independence of the restraint constant on the wavenumber. It is future work to use a 3D vertex model to derive the mechanical properties of the tissue from cell energy functions.

Discussion

We showed that the orientation of cell division and the restraint of out-of-plane deformation are respectively sufficient to cause the characteristic pattern and interspacing of folds in silico. The three patterns of folds obtained in this study coincide with the characteristic pattern of folds observed in beetle horns (Matsuda et al. 2017). The concentric pattern of folds is also found in the imaginal disk of Drosophila legs (Beira and Paro 2016; Schubiger et al. 2012).

However, we do not insist that the mechanism proposed in this study is the only mechanism determining the folding pattern. The actual beetle horn is cylindrical and there is thus a possibility that there is a mechanical difference in the ease of forming folds in the direction of the cylinder axis and in the circumferential direction. In fact, curvature-inducing wrinkles have been pointed out in stiff thin film on curved soft substrates (Stoop et al. 2015), and a similar mechanism might be available for actual beetles.

Although the global shape of the beetle horn primordia is affected in the dachsous-gene knocked-down beetle, in which the direction of cell division has been altered randomly, dense local furrows are not appreciably affected (Adachi et al. 2018). In the development of beetle primordia, actin accumulation, indicating apical constrictions of cells, has been observed at the position of a future furrow even before the local furrow forms, suggesting pre-patterning for the fine furrow pattern (Adachi et al. 2018).

Furthermore, it has been confirmed that perturbation of the axis of cell division does not appreciably affect the adult form of Drosophila (Zhou et al. 2019). Although a spatial pattern of actomyosin accumulation at a supracellular scale has not been reported for Drosophila primordia development, the spatial pattern of actomyosin accumulation or a still unknown mechanism may be considered for the primordia development of such small-sized bodies and tips of appendages. We speculate that there are at least two mechanisms that determine global and local folding patterns. The relationship between the spatial pattern of apical constriction and the 3D shape of deformed tissues has been investigated using a 3D vertex model (Inoue et al. 2017) and there are thus future opportunities to combine the present model with the pattern of apical constriction and to apply the model to clarify how these two mechanisms do or do not crosstalk.

In this study, the effect of confinement was introduced by the restraint of out-of-plane deformation on the apical side solely using Eq. (7). In addition to the apical side, it is possible to introduce the effect of confinement on the basal side, such as in the case of an extracellular matrix, using Eq. (7) and replacing the apical surface area term with a basal surface area term. However, the simulation results obtained when restraining both apical and basal sides coincide with those obtained when restraining the apical side solely in terms of there being an effective restraint constant that is the sum of values for the restraint constants on apical and basal sides (Appendix). With respect to introducing the effect of confinement, another possible model is tissue bounded by rigid walls on apical and basal sides. Because the rigid wall cannot be displaced or deformed, the interspace between the tissue and wall under initial conditions is an important parameter with which to determine tissue deformation. There is future opportunity to clarify how the interspace determines the tissue folding.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary material 1 (pdf 151 KB)

Supplementary material 2 (mp4 4100 KB)

Supplementary material 3 (mp4 7638 KB)

Supplementary material 4 (mp4 7644 KB)

Supplementary material 5 (mp4 7747 KB)

Supplementary material 6 (mp4 7679 KB)