Understanding the mechanisms causing buckling of epithelial monolayers

Epithelial monolayers form the building blocks of many tissues and organs in the body. Deformation and buckling of these monolayers is an important process in early development and in tissue renewal. It also plays an important role in the progression of cancer within epithelial tissues. In this study we apply the recently developed Rigid Body Multi–Cellular Framework (RBMCF) to investigate the mechanisms of buckling for an epithelial monolayer attached to a basement membrane and the surrounding stromal tissue. The epithelial monolayer and supporting basement membrane and stromal tissue are modelled using two separate vertex dynamics models and interactions between the two are considered using the RBMCF to ensure biologically realistic interactions. Model simulations are used to investigate the effects of cell–stromal attachment and membrane rigidity on buckling behaviour. We demonstrate that there are two competing modes of buckling, stromal deformation and stromal separation.


Introduction
Epithelial layers can be described as a sheet of tightly packed cells that line the surfaces and cavities of organs, forming an interface with the environment [27].They can be found in numerous parts of the body, including ⇤ jmosborne@unimelb.edu.au,ORCID:0000-0002-5622-0104 April 1, 2024 the skin, the bladder, the kidneys, the digestive tract, and reproductive system [27].A functioning epithelial layer is made up of three distinct parts [19] (Figure 1):

Epithelial Cells Basement Membrane
• a stromal layer, made up of connective tissue, vessels and nerves which forms the structure for the cells to sit upon; • a thin basement membrane layer made up of a web of collagen and glycoprotein molecules which separates the stroma from the cell; and, • the epithelial cells themselves.
Blood vessels found in the stromal layer transport nutrients to (and in many cases, from) the epithelial cells via di↵usion through the connective tissue and basement membrane [27].The basement membrane provides a substructure that the cells attach to, which persists even when cells are lost.
As long as the membrane remains intact, epithelial layers can regenerate [27].
The epithelial cells themselves come in many di↵erent types, depending on the particular function the epithelium needs to perform.
Epithelial tissues, particularly monolayers, are often the objects of study of mathematical modelling.Due to their simplicity; they often lend themselves very favourably to implementation using multicellular modelling techniques, where individual cells, and their interactions are considered.
Epithelia can be categorised into two broad classifications: stratified and simple [27].A stratified epithelium is made up of multiple layers of cells, where the bottom-most layer is attached to a basement membrane, and the higher layers sit on top of other cells.Cells proliferate in the bottom layer and transit up through the layers until they die and are shed into the environment [21].
A simple epithelium on the other hand, is characterised by a monolayer of cells where each cell maintains connection to the basement membrane [27].
The single layer functions as a means to absorb and filter material from the environment into the organ [32], as well as secrete material (mucus, enzymes) into the environment [7].Unlike the stratified epithelium, cells that lose connection to the basement membrane will normally die through a process called anoikis [20].Loss of this behaviour is understood to be a contributing factor to cancer [20,22].By definition, cells in a simple epithelium maintain a single layer, so dividing cells naturally create migratory pressures parallel to the basement membrane [15].
The lining of the the small and large intestines, is a key example of a simple columnar epithelium (so called because of the tall narrow shape of the monolayer of cells).Within the surface a vast number test-tube-shaped indentations know as the crypts of Lieberkühn can be found.The intestinal epithelium attracts a lot of research attention.This is due to the association with bowel cancer, where a loss of proper function of the epithelial self renewal process leads to pre-cancerous lesions known as polyps, which left untreated can become malignant [6].It is not currently well understood how new crypts are formed, but it is thought that they can emerge via a surface invagination process (in a developing embryo) or via a process where an existing crypt splits in two, known as fission [23].Modelling has been used to examine a range of questions about the crypt.One example is clonal conversion [35].
Which has important applications for carcinogenesis, as it helps explain how a single genetic mutation can spread.Other modelling studies focus on the developmental perspective, for example, aiming to understand how crypt-like structures are formed in colonic organoids grown in vitro [4].
Another example that attracts significany attention is the imagnial wing disc of the Drosophila melanogaster (fruit-fly) lava.This is a popular biological system to study for wide range of reasons, including the short life cycle of the animal, and the ease with which it can be cultivated [33].In the early stages of growth, a feature appears from the ectoderm made up of cuboidal epithelial cells, which flattens into a disc shape, and through a sequence of changes, eventually becomes a fully formed wing [10].There are many questions surrounding the processes that regulate the changes, and some of these have been investigated using vertex models [18].Among the questions of interest are: how organ size is regulated during growth [2], how cell growth is regulated by mechanical forces [1], and how the final cell topology is formed [16].
A common modelling approach is to treat the surrounding tissue as static and treat the monolayer as a 2D planar sheet [35] or 2D surface [13,14] This kind of modelling is useful to examine things such as clonal conversion, and cell proliferation, however it has limitations due to the fixed structure.An alternative view is to consider a deformable membrane.The common way to achieve this take a cross-section of the layer, representing it as a single row of cells.This approach is best suited to examining movements that happen perpendicular to the layer, including buckling, and cells being ejected from the layer [4,11,12].There are examples where deformable substrates have been implemented in three dimensions [9], however due to their node centric implementation, monolayer deformation is limited [8].
Continuum modelling has also been used to examine the buckling of epithelial layers.In [15], Edwards and Chapman investigated a flat beam of cells viscoelastically attached to a supporting stromal layer.They see that each of: the applied force; level of cell attatchment; and cell proliferation can independently trigger buckling.Jones and Chapman [24] employ shell theory to explore the role of apical constriction of epithelial cells in gastrulation during embryonic development.Additionally, Nelson and colleagues published two studies that examine the buckling of a one-dimensional beam into crypt-like structures [30], as well as the patterns created by the buckling of a two-dimensional epithelial sheet [31].Using such models relies on being able to approximate the tissue as a continuous medium.This can be an unrealistic assumption at certain scales and in certain situations depending on the behaviours that are of interest, but when the conditions are right, it can provide valuable insights.
However, none of these models have been able to model the epithelial layer, and the supporting sub-structure as separate objects, allowing an investigation of their interactions.The Rigid Body Multi-Cellular Framework (RBMCF), from [8], makes this possible.The RBMCF provides us with a range of modelling capabilities, including the ability to represent interacting deformable surfaces.Broadly speaking the RBMCF represents cells, and membranes, as collections of edges, which exert forces on each other.The RBMCF has been used to simulate: a tumour spheroid with decognal cells; a non-intersecting epithelial layer consisting of connected quadrilateral cells; rod based bacterial cells; and spherical cells proliferating in a closed polygonal membrane [8].One of the key advances of the RBMCF is its ability to represent interactions between surfaces in a mechanically consistent way.
In this paper we exploit the RBMCF to develop a model of an epithelial monolayer, attached to a basement membrane and stromal tissue, in order to understand the mechanisms behind buckling of the epithelial monolayers.
The remainder of this paper is structured as follows.First in Section 2 we present our model of an attatched epithelial monolayer.In Section 3 we present simulations of our model and use it to demonstrate multiple mechanisms for buckling.Finally in Section 4 we discuss the impact of our results and present some avenues for future work.

Methods
Our model couples the edge based epithelial model (with connected quadrilateral cells) from [8] with an edge based model for the basement membrane and tissue stroma.In the RBMCF vertices in the model, x i (for i = 1, . . ., N, where N is the total number of vertices), move based on the balance of forces ( The model equations presented in this study are all non-dimensionalised by a typical cell diameter (taken to be 10µm [27]) therefore all positions and lengths are given in terms of Cell Diameters (CD).
This section is organised as follows.We present the model for the epithlial layer in Section 2.1 and the model for the basement membrane and tissue stroma in Section 2.2.The coupling between these models is described in Section 2.3 and, finally, details of model implementation are provided in Section 2.4.

Modelling epithelial monolayer
Following the epithelial layer example in [8] the epithelial layer is defined by a set of connected quadrilaterals (shown in Figure 2 (a)).Let the vertices be denoted by x i for i = 1, . . ., N C , where N C is the number of vertices in the epithelial layer.For a given vertex x i , the force due to the embodied energy of the cells is calculated by Figure 2: (a) schematic of epithelial monolayer (b) the division process for epitheial cells.When a cell is ready to divide, a new edge is introduced that splits the parent cell into two adjacent daughter cells.The top and bottom edges are split in two by adding a node to their centre.The new edge then joins the new top node to the new bottom node.Modified from [8].(c) the stromal tissue is represented by a single large object defined by the nodes and edges that make up its boundary.It has a target area and target perimeter and forces are applied to the nodes via energy methods to push it into its preferred shape.The bottom and sides are fixed in place (red and green vertices), while the top surface is free to move (green and black vertices).
(d) interaction between the epithelial layer and the stromal tissue.The node from one tissue (in black) interacts with the edges of the other (at the red points). .
where U Cell j is the energy embodied in cell j is, r i is the gradient with respect to the coordinates of vertex i and M i is the set of all cells that contain vertex i.
Following [8] we use energy methods to drive the area and perimeter, of our polygonal cells, in a similar manner to the vertex dynamics model [29].
where a j represents the current area, and A Cell j (t) the target area (which changes based on the cells age), and where p j , and likewise P Cell The cell cycle (which controls cell division) is commonly broken up into four phases [3], and occasionally as many as 11 [26].However, we are mainly concerned with the growth and division processes therefore, following [8], we define two phases: a "growth" phase G where the cell actively increases in size, and a "pause" phase P where the cell maintains its current size.A new cell, immediately after division enters the pause phase where it sits for an amount of time t P which is drawn from a U ( tP 2, tP + 2) distribution where tP is the mean pause duration (minimum 5 1 ).At the end of the P phase, the cell then enters the G phase, where it grows in size over a period of time, t G drawn from a U ( tG 2, tG + 2) distribution where tG is the mean growth duration (minimum 5).After growth has completed, the cell is able to immediately divide into two child cells by splitting the top and bottom edges in half, and adding a new edge between the new top and bottom nodes see Figure 2 (b).
The actual growth of cell j is controlled by setting the target area, A Cell j (t), and perimeter, P Cell j (t), and letting the internal forces cause it to expand.At the end of the growth period the cell will divide into two child cells Here, we decouple the target perimeter from the target area, and specify 1 minimum specified to ensure computational stability.
their values by the equations where t j is the age of cell j at time t.In this model we set the initial perimeter to P 0 = 3.4, the grown perimeter to P G = 4, the initial area to A 0 = 0.55, and the grown area A G = 1.The "grown" values were chosen to have square cells in isolation, and the "initial" values were empirically chosen to minimise pinching 2 , hence preventing artificial layer detachment.

Modelling the basement membrane and tissue stroma
We extend the epithelial ring exemplar from [8], to include the basement membrane and stromal tissues that provides the supporting structure for the epithelial layer.To that end, we model the basement membrane and tissue stroma (referred to here as the stroma) as a large polygonal object, illustrated in Figure 2 (c).The stroma is modelled as a large rectangular body, with the bottom and vertical sides represented by a single edge each, and the top surface made up of multiple edges.The bottom and side edges are fixed in space and the top is allowed to move freely to create a dynamic surface.
Let x i for i = N C + 1, . . ., N C + N S where N S is the number of vertices in the stroma.
Internal forces on the stroma are calculated using the same energy methods as used for cells but with one large unconnected polygon.For a given vertex x i , the force due to the embodied energy of the stroma is calculated by where again r i is the gradient with respect to the coordinates of vertex i and U Stroma is the embodied energy of the stroma.A target area and target 2 Pinching occurs because of the sudden change of internal area and outer perimeter when a cell divides.This results in unbalanced forces being introduced to the new edge between the two daughter cells, causing it to contract sharply.Pinching can artificially cause layer detachment, but fortunately it can be mitigated by carefully controlling the growth of the cell target area and perimeter.perimeter are specified, and the di↵erence between those and their measured values is used to calculate the embodied energy: Here the parameters ↵ S and S are the energy factors specifically for the stroma, and they normally will be di↵erent from the respective energy factors for the epithelial cells.The target area (A Stroma = 39) and perimeter (B Stroma = 28.8)will remain constant throughout a simulation, set to the values for the initial rectangular stroma.In this study we will investigate how the stromal mechanical parameters influence the behaviour of the epithelial monolayer.
Note that here, we are assuming that the basement membrane and stromal tissue can be appropriately approximated by an elastic membrane filled with homogeneous fluid.A normal stromal tissue is made up of a wide range of components (blood vessels, ducts, other cells etc.), held together in a connective tissue and future work will use a more complicated model for the stromal tissue.
We fix the bottom vertices of the stroma in place and fix the x-coordinate of the upper edge vertices, see the red and green vertices in Figure 2 (c) respectively.Specifically, the edges are held at x = 0 and x = 10CD.The edges being fixed in place implies a rigid boundary between the tissue of interest and the rest of the biological system.In various tissues (for example the intestinal crypt), we can see features that may play a similar role to a fixed boundary (e.g. the lamina propria [12]), although these likely cannot be considered completely rigid.At some point we can expect that the transfer of force through a supporting tissue has dispersed over a wide enough area/volume that any deformation is orders of magnitude smaller than the scale we are observing.Hence, as a model feature when approximating stromal tissue, the bottom rigid boundary is a fair option when used appropriately.

Coupling the epithelium and stroma
To couple the epithelium with the stroma we set the horizontal width of the stroma to 10CD, accommodating an initial population of epithelial 20 cells placed on top of the stroma.Epithelial cells divide and push adjacent cells along the stroma and are removed once they extrude past the horizontal boundaries of the stroma, see Figure 2 (a).
The epithelial layer will sit on the stroma and interactions between the two will be determined using node-edge interactions (Figure 2 (d)).This allows the layer to remain separate from the stroma, but still a↵ect its shape.
Following [8], the interaction force, F Interaction , between vertices and edges will be calculated according to In the analysis to follow, we will only be examining cases up until the point of buckling, meaning that interactions between epithelial cells or between the stroma and itself are unlikely to occur however we include repulsion here in order to avoid overlapping of edges.
When this interaction force acts on the vertices directly (the black circles in Figure 2 (d)) we can calculate the force applied to the vertex directly as where x is the separation between vertex i and neighbouring edge k (the length of the dashed line in Figure 2 (d)), and xki is the unit vector in the direction of the vertex i from the point of intersection on the edge k (the direction of the dashed line shown in Figure 2 (d)).
When this interaction force acts on an edge k (the red circles in Fig- , following [8] there are two components, one from the direct interaction with the centre of mass on the edge and one from the perpendicular force which results in a torque on the edge.For a rigid body with a single externally applied force F that does not pass through through the centre of drag, it can be shown that F can be replaced by a force acting at the centre of drag and a moment [28].Therefore, following [8], the motion of the body can be completely determined by evaluating the two equations Where x D i is the centre of drag of the edge, ✓ is the angle the edge makes with the x-axis, ⌘ D is the total drag on the edge, I D is the moment of drag of the edge, and d is a vector from the centre of drag to the line of action of F such that the two vectors are perpendicular, see [8] for more detail and formal definitions.We can use Equation ( 10) and a Forward Euler approximation to calculate the motion of the edge due to the applied force F in a small timestep t.This can be used to calculate two virtual forces, which when applied to each end of the edge, k 1 and k 2 , would give the equivalent motion, , and These virtual moves can be combined for each edge that vertex i is in to give.
where E i is the set of edges containing vertex i and i k is the index (1 or 2) that vertex i is in edge k.
Finally, from Equations ( 2) and ( 6), the force applied to the vertex i in the model, due to internal interactions, is given by.
Taking Equations ( 9), ( 12) and ( 13) and balancing the applied force with the drag experienced by the vertex (Equation ( 1)) allows us to derive the following equations of motion for vertex i.
We use a Forward Euler approximation for the derivative to simulate the evolution of the vertex locations [8].

Implementation
This

Epithelial layer exhibits buckling
Running the model, we observe di↵erent types of behaviour depending on the input parameters.Specifically we set the mechanical parameters ↵ S = 10, S = 6, and  Attract = 10 and vary the proliferation rate.In Figure 3 (a) (and Supplementary Movie 1) we present a simulation where the input parameters intuitively would be conducive to a stable layer (slow proliferation, tG = tP = 11), we observe that the epithelium and stroma maintain contact.
The layer may exhibit some small deviations, but over long simulation periods, there are no signs of buckling or detachment.In Figure 3 (b) (and Supplementary Movie 2) we present a simulation in the reverse situation (faster proliferation, tG = tP = 7).We see buckling and eventually detachment in simulations.
We can examine the causes of buckling systematically by performing a parameter sweep across a range of suitable input values.But first we need a way to determine when buckling has occurred.When an epithelial layer in our model buckles, it goes from having a predictable shape with a bounded number of cells, to some unspecified shape with an exponentially growing hump.A simply way to quantify the whether a layer has buckled is to compare the path length of the layer to the value expected for a stable layer.
To that end, we define the buckling ratio, given by where s is the measured path length of the underside of the epithelial cells, and w is the horizontal width the layer covers (here 10 CD).By definition, the ratio will be positive and greater than 1, with values very close to 1 representing a near flat layer, and numbers greater than 1 representing increased undulations in the layer.We can calculate this value for both the epithelial layer r E and the top surface of the stroma r S .
For the examples in Figure 3 we see that the stable layer (a) has maximum buckling ratios of r Max Buckling ratio  1. Simulations are stopped at t = 200 hours or when r E = 1.2.Simulations are shown in Supplementary Movies 1 and 2 respectively.

Buckling is driven by increased proliferation
A common observation from simulations of epithelial layers is that increased proliferation can drive buckling [15,5].Therefore, in our first experiment, we test buckling over a range of cell cycle lengths by varying the pause phase and growth phase duration's.All parameters are fixed to the same value in simulations (See Table 1), except for the values of tP and tG which we vary here.Simulations are run until buckling occurs (defined by r E 1.05), or until the simulation reaches 500 hours.Each parameter set is run 50 times, each starting with a di↵erent random number generator seed.
Figure 4 plots the proportion of buckled simulations for a range of mean pause and growth phase durations between 5 and 12 hours.This range captures interesting features in the parameter space, and is roughly in line with the reported cell cycle lengths for crypt cells in mice [34].Each point in the plot represents a fixed parameter set where 50 simulations were run.
The colour of the point indicates the proportion of those simulations where buckling occurred, with dark red corresponding to all buckling, and dark blue corresponding to all simulations remaining stable for the full 500 hours.
Here we can see that faster cell cycles do indeed increase the likelihood of buckling, confirming that the model can reproduce established results.The colours between red and dark blue (i.e.orange, yellow, green, teal) represent parameter sets where some simulations buckled, and some remained unbuckled for the full 500 hours.
In Figure 4 we can observe distinct bands of each colour, forming a "transition front" that can reasonably be fitted with a straight line.If we fit a level curve for parameter sets with specific proportions of buckled simulations, we produce a family of lines with gradient approximately 15  9 .Specifically the level curve for parameter sets where half the simulations buckled forms roughly a straight line with equation tP = 15 9 tG + 20.By a simple geometric argument, it is straight forward to see that passing through the transition region from all buckled to none buckled is much faster when tG is varied compared to when tP is varied.This implies that changing the length of the growth phase can more rapidly alter the buckling behaviour of the layer.This logically makes sense, since the growth phase is when additional internal stresses are added to the layer, and if these stresses are introduced more rapidly (i.e. with a shorter tG ), one would expect buckling to be more likely.Reducing the pause phase while keeping the growth phase constant will not change how quickly internal stresses are introduced, but it will alter the total number of cells that are introducing internal stress at the same time.This again, can cause buckling, but as demonstrated, is slower to take e↵ect.

Buckling is independent of stromal compression
A novel feature of this model is its ability to separate the e↵ects of the stromal tissue surface behaviour, the stromal tissue volume behaviour and the strength of adhesion to the surface.We now perform some experiments to investigate how these e↵ects interact with regards to buckling.
As a first experiment, we vary the parameters for area energy, ↵ S , and perimeter energy, S , for the stromal tissue.The design of the energy method (Equation ( 7)) comes from the vertex model [17] where the properties of a cell boundary are likened to the surface tension of soap films, and the internal properties are modelled as a compressible elastic body.In that regard, ↵ S can be thought of as a measure of the tissue's resistance to compression, while S can be thought of as the spring sti↵ness of the upper surface which is covered by the basement membrane.By sweeping across the two energy parameters, we can investigate the e↵ect of the stromal tissue compared against the e↵ect of the basement membrane (with all other parameters held fixed, as in Table 1, with tG = 5, tP = 10, and  Attract = 10).As with the previous experiment, we run 50 simulations for each distinct parameter set and count the proportion of simulations where buckling occurs (r E 1.05) before 500 simulation hours have elapsed.In the context of the model this is a sensible result.The bottom and side boundaries of the stroma are fixed in place, so the only way the area or perimeter can change is by the top free surface deforming.In its equilibrium state, the surface forms a straight line from corner to corner.Any deviation from this state will lengthen the perimeter, since the shortest distance between two points in Euclidean space is a straight line.However, for a given small deviation, it is possible to introduce an additional deviation that restores the area.In other words, there is only one spatial configuration that results in zero perimeter energy, but there are many that result in zero area energy.The assumption of a homogeneous fluid filling the stromal volume Figure 5: E↵ect of stromal parameters on buckling.Proportion of simulations where buckling was observed when varying the area energy ↵ S and the perimeter energy S .Each point represents 50 simulations, and the colour represents the proportion of those that buckled before the end of 500 simulation hours.tP = 10, tG = 5,  Attract = 10, and all other parameters as in Table 1.
implies that all ways of generating the same measured area are equal.In reality, this is unlikely since stromal tissue is not a perfectly homogeneous fluid, rather it is a solid tissue made up of numerous components that exhibits visco-elastic properties as discussed in Section 1.That said, there will be some degree of comparable behaviour in the physical situation since much of a tissue is made up of fluid.To the extent that a real tissue behaves like a fluid, this result suggests the basement membrane influences the chance of buckling for a epithelial monolayer to a much greater degree.

Buckling can occur through multiple mechanisms
We now turn our attention to the interaction force, between the stroma and epithelial layer, which is given by Equation ( 8).The parameter  Attract captures the strength of attraction between the water-bed stroma and the rectangular cell model of the epithelial layer, hence by increasing or decreasing  Attract we can simulate a stronger or weaker attachment between basement membrane and epithelial cell.Intuitively, a weaker interaction should lead to increased surface detachment and buckling.However, in Figure 5 we also demonstrated that buckling is driven by the surface tension in the membrane, S .We now perform an experiment to test the interaction S and  Attract .Once again, we perform a parameter sweep across the two variables keeping all others constant (as in Table 1, with tP = 10, tG = 5, and ↵ S = 10).We run 50 simulations for each parameter set (each starting from di↵erent random seeds), and track the proportion of simulations that end in buckling (r E 1.05) before completing 500 simulation hours.1.
either by the layer breaking away from a flat stroma, or losing connectivity at the peak of a local undulation.
We can investigate where the di↵erent modes occur by examining the buckling ratio for the stromal surface r S .We have defined buckling to occur when r E 1.05, so we can use the value of r S at this time to determine which buckling mode has occurred.Looking at Figures 6 (d) and (e) we can see that if r S ⇠ r E = 1.05, then we know the layer has stayed attached to the stroma up to the point of buckling (Figure 6 (e)), but if r S << r E = 1.05, then detachment has occurred to trigger buckling (Figure 6 (d)).
Figure 7 illustrates the average value of r S when buckling occurred for each of the parameter sets.Not all repetitions (i.e.random seeds) exhibit buckling for a given parameter set, particularly in the transition region, so the average r S is calculated only from those simulations that did end by buckling.Red indicates the average stroma buckling ratio is equal to the layer buckling ratio (i.e.r E ⇠ r S ), and blue indicates the stroma is close to flat (r S ⇠ 1).We can clearly see that there are two distinct regions in parameter space, roughly divided by the line  Attract = S .To the bottom right of the plot we see that the stroma has a buckling ratio close to 1, indicating that it has remained flat while the layer has buckled, implying separation.Above the line, we see that the stroma is more buckled as S (the stromal perimeter energy) is decreased, or  Attract (the layer adhesion) is increased.These observations fall in line with intuition: a surface that in some sense is less sti↵ will more easily deform under load, and a stronger adhesion force is more likely to pull a sti↵ surface with it as it buckles.

Discussion
In this study we have applied the RBMCF [8] to a problem that, from a multi-cellular modelling perspective, has seen little attention [5].This is perhaps due to the di culty in modelling surface interactions with pre-existing tools.Using our model of a stroma-epithelial interaction, we were able to reproduce the established result, that increased proliferation is associated with increased buckling [15].More importantly, the novel construction of the model allowed us to investigate separately the impact of basement membrane sti↵ness and epithelial layer adhesion on the likelihood of buckling.
We were able to show that buckling can occur via two modes: one where the stroma itself is weakened, allowing the layer to buckle while still attached, and the other where internal stresses cause the layer to break away from the basement membrane.This result would have not been possible without a modelling approach that can handle two distinct interacting bodies.
As detailed during the description of the model there are a number of implicit assumptions in the construction of our model impacting how closely it will approximate stromal tissue.We assumed that the basement membrane and stromal tissue can be approximated by an elastic membrane filled with homogeneous fluid.We also considered the system to be represented by a cross sectional geometry.Despite these limitations, the model does allow us a means of investigating how a proliferating monolayer of cells would interact with a deformable supporting structure, and it serves as a good first approximation.
Future work will look to work on relaxing these assumptions.First, we will implement an explicit model for the membrane and use a multicellular model (like the spheroid exemplar from [8]) for the surrounding stromal tissue.This will enable us to look at realistic crypt geometries as in [12].
Further work will extend the model to three dimensions, so we have a two dimensional epithelial layer.This will enable us to look at the e↵ect of cell sorting of the heterogeneous cells in the layer which can be seen to induce buckling in, for example, crypt development [25].
This work represents a step change in the mechanical detail available for simulations of epithelial monolayers which has the potential to make an impact in multiple areas of developmental biology.

Figure 1 :
Figure 1: Schematic of epithelial monolayer.Example of simple epithelial monolayer.Showing the epithelial cells basement membrane and stromal tissue, containing: stromal cells; connective tissue; and blood vessels.

.
due to: the internal model for the cell (or other structures), F Internal i ; the interactions between cells or between cells and other structures, F Interaction i ; and the viscous drag on the vertices, F Drag i for i = 1, . . ., N.

j
(t), represent the current perimeter and its target value.The parameters ↵ C and C are energy density factors controlling how strongly the cell wants to return to its preferred area and perimeter.Here, following [8], we use the values of ↵ C = 20 and C = 10.
) Where x is the perpendicular distance between edge k and vertex i,  Attract is the spring sti↵ness for attraction (and is zero for cell-cell or stroma-stroma interactions), and  Repel is the spring sti↵ness for repulsion.Here we choose a fixed value of  Repel = 10 and vary  Attract .The typical interaction separation and the maximum interaction distance are given as d Sep = 0.1 and d Lim = 0.2 respectively.

E ⇡ 1 .⇡ 1 .= 1 .2 and r Max S ⇡ 1 .
02 and r MaxS 02 where as the buckled layer has r Max E 025 when the simulations is stopped.Note we set these simulations to stop whenever r E = 1.2.The stable simulations maintain r E < 1.025 for however long we run them (Figure3 (a)) however once a layer moves beyond r E increases to r E = 1.05 and keeps going.We therefore choose r E = 1.05 to define the point where buckling has occurred.This is shown by the dashed orange line in Figures3 (c) and (d).

Figure 3 :
Figure 3: Epithelial monolayer can exhibit buckling.(a) simulation with slow proliferation tP = tG = 11.(b) simulation with fast proliferation tP = tG = 7. Snapshots are given over time for three indicative times t 1 < t 2 < t 3 (di↵erent for (a) and(b)).(c) and (d) show buckling ratios, r E (blue) and r S (red), over time for the slow proliferation rate, and (d) for fast proliferation rate.r Buckled = 1.05 threshold is shown by the orange dashed line.Mechanical parameters are ↵ S = 10, S = 6, and  Attract = 10.All other parameters as inTable 1. Simulations are stopped at t = 200 hours or when r E = 1.2.Simulations are shown in Supplementary Movies 1 and 2 respectively.

Figure 4 :
Figure 4: Proliferation drives buckling.Proportion of simulations where buckling occurred for varying values of the pause phase duration and the growth phase duration.Red indicates all 50 simulations exhibited buckling and dark blue indicates all simulations remained stable for 500 hours.As expected, faster cell proliferation (i.e.smaller values of tP and tG ) results in a higher chance of buckling.Mechanical parameters are ↵ S = 10, S = 6, and  Attract = 10 All other parameters as in Table 1.Simulations from Figure 3 are shown by the white and black +'s for Figure 3 (a) ( tP = tG = 11) and Figure 3 (b) ( tP = tG = 7) respectively.

Figure 5
Figure 5 illustrates the parameter sweep for 1  ↵ S , S  20, with all other parameters held constant.Here we can see a clear vertical band indicating almost zero e↵ect from varying the ↵ S , and a pronounced e↵ect from varying the S over the same range.This leads us to the conclusion that the compression resistance of the stromal tissue has no impact on the buckling of the epithelial layer.

Figure 6 (Figure 6 :Figure 6 :
Figure6(a) illustrates the parameter sweep for varying stromal sti↵ness and membrane adhesion (1  S ,  Attract  20).As anticipated, the plot demonstrates that both weakening the layer interaction with the surface and lowering the perimeter energy increase the likelihood of buckling.What is not clear however, is whether the buckling in these two cases is caused by the same mechanism.Figures6 (b)-(c) illustrate snapshots from two di↵erent modes observed from simulations.Figures6 (d)-(e) illustrate the corresponding buckling ratio over time.For high perimeter energy and low adhesion ( S = 15,  Attract = 5), the stroma remains relatively flat, but when the conditions are right, the epithelial layer buckles and breaks away, Figures6 (b) and (d).In the opposite case (low S = 1 and high 6 = 1), buckling occurs while the layer maintains attachment to the stroma, Figures6 (c) and (e).This creates regions of high local curvature, which eventually lead to the layer once again pulling away.Buckling can occur in two di↵erent ways,