## Abstract

The formation of three-dimensional structures from patterned epithelial sheets plays a key role in tissue morphogenesis. An important class of morphogenetic mechanisms relies on the spatio-temporal control of apical cell contractility, which can result in the localized bending of cell sheets and in-plane cell rearrangements. We have recently proposed a modified vertex model that can be used to systematically explore the connection between the two-dimensional patterns of cell properties and the emerging three-dimensional structures. Here we review the proposed modelling framework and illustrate it through the computational analysis of the vertex model that captures the salient features of the formation of the dorsal appendages during *Drosophila* oogenesis.

This article is part of the themed issue ‘Systems morphodynamics: understanding the development of tissue hardware’.

## 1. Introduction

It is believed that the organization of the first multicellular animals resembled an envelope, a structure formed by two epithelial layers, with the lower layer adhering to the external surface and top layer exposed to the external environment [1]. Such two-layered organization is observed in *Trichoplax adhaerens*, one of the simplest extant metazoans, an organism which can be studied in the laboratory [2]. This organism moves along the surface in an amoeboid fashion in search of food and can bend in the direction normal to the surface to engulf particles [3]. A similar class of shape transformations, involving in-plane reshaping and out-of-plane bending, is observed in the early embryos of all animals, from ascidians to humans. In particular, at an early stage of development the embryo forms a blastula, which comprises an epithelial sheet enclosing a fluid-filled cavity. As embryogenesis proceeds, more complex structures are formed by progressive bending and in-plane rearrangements of epithelial subdomains [4,5]. Imaging studies across species revealed that these events are driven by controlled shape changes of individual epithelial cells. By regulating the activity of actomyosin networks, these cells can control their apical areas and lengths of individual edges [6,7].

The molecular machinery regulating these processes is incredibly complex, but the general principles of morphogenesis driven by spatio-temporal control of actomyosin and cell adhesion systems have become progressively understood, to the point that we can now recognize the same functional subroutines in diverse experimental models. For instance, localized apical constriction can induce localized bending of epithelial sheets [7]. Among the numerous examples of this well-studied process are mesoderm invagination in *Drosophila* [8], primary invagination in *Ciona intestinalis* [9] and early stages optic cup formation in *Xenopus* [10]. In all of these cases, cells within the invaginating region constrict their apical surfaces (figure 1*a*). A finer mode of control is needed to induce in-plane reshaping of epithelial sheets during the convergent extension movements observed in animal gastrulation [11]. As was shown by studies of the germband extension in *Drosophila* [12] and *Xenopus* gastrulation [13], this type of morphogenesis relies on the regulation of only a subset of cell edges within the epithelium. In these cases, increased contractility of a subset of cell edges across the sheet triggers cell intercalations, leading to in-plane reshaping of an epithelial subdomain (figure 1*b*).

A combination of in-plane reshaping of epithelial subdomains and their out-of-plane bending can generate a wide range of three-dimensional multicellular shapes. Computational models of abstracted mechanisms can be used both to explore the functional capabilities of a class of mechanisms and test the feasibility of mechanisms proposed to explain the dynamics observed in specific developmental systems. Starting from the seminal papers by Odell *et al.* [14], cell-based computational models have been used successfully for both of these purposes. However, most existing models focus either on out-of-plane bending or on in-plane cell shape changes and rearrangements. The first class of models used to mimic epithelial invaginations and evaginations relies on cells with distinct apical and basal surfaces [9,15–19]. Most of these models are essentially two-dimensional (describing the cross-section of a deforming epithelial shell), designed to capture the emergence of the characteristic omega-like shapes during early gastrulation. In the second class of models, designed to explore in-plane deformations, cells are approximated by convex polygons that can deform and exchange neighbours but are always confined to the flat surface [20–22].

In our recent work we have developed a model that can be used to describe both in-plane cell constriction and out-of-plane bending of epithelial sheets [23]. This model can be used to systematically explore the connections between two-dimensional patterns of cell contractility and the resulting three-dimensional tissue deformations. Here we show how the same model can be extended to include spatially controlled cell rearrangements. The paper is organized as follows. In §2, we review the key aspects of the mathematical formulation and numerical implementation of the model. In §3, we demonstrate how the model can be used to explore the multicellular dynamics involved in the three-dimensional morphogenesis of the respiratory appendages on the *Drosophila* eggshell, an experimental system in which morphogenesis can be studied using genetic, imaging and modelling approaches [24,25]. We conclude in §4 by discussing the limitations of the proposed model and outline directions for future work.

## 2. Physical description of an epithelial sheet

The starting point for our computational analysis of three-dimensional epithelial morphogenesis is a model in which the epithelium is modelled as a flat sheet constructed from polygonal cells [20,21,26,27]. The energy of the system depends on the geometrical configuration of the vertices of the polygons through the areas of cells and lengths of cell–cell edges:
2.1The first sum in equation (2.1) runs over all cells *s* and penalizes the deviation of each cell's area *A _{s}* from a target value , with elastic modulus

*μ*. The second sum runs over all edges

_{s}*j*, representing an effective line tension, where

*l*is the length of each edge and

_{j}*σ*is a line tension coefficient. This term models the joint effect of cell–cell adhesion and actomyosin contractility (figure 2

_{j}*a*).

In a uniform tissue at equilibrium, each cell in the model epithelium is a hexagon whose side length is determined by the balance of the area and edge energies (figure 2*b*). However, when the tissue is patterned, cell shapes can change as a function of position (figure 2*c*). When the vertices are allowed to move in three dimensions, some of the two-dimensional patterns of cell properties can cause out-of-plane deformations (figure 2*d*). For instance, contractile contours, which are implemented by increasing the line tension along the edges forming a ring around the patch of cells, can lead to a pitchfork bifurcation, in which a flat configuration loses stability with respect to out-of-plane displacements [28]. Specifically, when all cells within the patch increase their edge tension, cell shapes within the patch and its surroundings change, but the sheet remains flat. This is in contrast to what is observed in a wide range of experimental and computational systems, in which apical constriction of a patch of cells causes localized tissue bending [7].

The bending of an epithelium due to apical constriction is commonly described by models in which cells have distinct apical and basal surfaces and the volume of each cell is conserved [9,14–17,29–31]. In these models, constriction of the apical surface leads to expansion of the basal surface due to the incompressibility of the cell. Each cell thus adopts a wedge-like shape that provides localized curvature to the epithelium, which bends towards the basal side (figure 2*e*,*f*).

In our recent work [23], we found that these effects can be captured by adding a new term to the energy function in the model with polygonal cells (figure 3*a*,*b*). The additional term captures the distinction of apical and basal surface while maintaining a two-dimensional representation for the epithelium. The term mimics the effect of the localized contractile segments that are present on the apical surface of an epithelium. The vertices of these segments are effectively ‘offset’ from the cell centres along the cell normals. The energy for such offset segments can be defined as
2.2where the index *i* represents the edge shared by two adjacent cells with endpoints *s*1(*i*) and *s*2(*i*), *C** _{s}* denotes a cell centre, defined as the mean position of the vertices belonging to cell

*s*,

*N**denotes the unit normal to that cell*

_{s}*s*,

*Γ*characterizes the line tension along the segments joining the cell centres, and

*h*represents the half-thickness of the epithelial monolayer.

In short, is the line tension energy term for the segments joining the points . We only keep the linear truncated form of (equation (2.2)) and add it to the energy defined in equation (2.1) to define the total energy of the epithelial sheet. The truncated form is given by
2.3where *u** _{i}* is the unit vector joining the cell centres:
2.4

The total energy of the system is given by *E* = *E*_{2D} + *E _{Γ}* (equations. (2.1), (2.3)).

As *Γ* is increased, the sheet favours the configuration where the cell centres are closer together. Additionally, when *h* ≠ 0, an increase in *Γ* brings the endpoints of the normal vectors close together (*h* > 0) or further apart (*h* < 0), causing the sheet to bend. This can be appreciated from analysing the following simplified setting, which assumes that *N*_{s}_{1}, *N*_{s}_{2} and *u** _{i}* are coplanar (figure 3

*b*). Consider two cells

*s*1 and

*s*2 adjacent to a particular edge

*i*. Let

*θ*be the angle between the normals

*N*

_{s}_{1}and

*N*

_{s}_{2}and assume that the unit vector

*u**makes an angle (*

_{i}*π*−

*θ*)/2 and (

*π*+

*θ*)/2 with

*N*

_{s}_{1}and

*N*

_{s}_{2}, respectively. In the additional term (equation (2.3)), the contribution from this segment for the case is proportional to

*Γ*(

*l*+

_{s}*hθ*) where

*l*is the distance between the cell centres,

_{s}*C*

_{s}_{1}and

*C*

_{s}_{2}, adjacent to the edge

*i*. From the simplified expression, it is more apparent that to minimize energy in the presence of the additional term, the sheet can now constrict (decrease

*l*) as well as bend (increase or decrease

_{s}*θ*, depending on the sign of

*h*).

The modified two-dimensional model successfully captures the essential behaviour of the three-dimensional model [23]. For example, the epithelium bends only in the apical-to-basal direction in response to the two apical contractility patterns that we considered – a contractile ring or a uniformly constricting patch of cells. The epithelial sheet smoothly deforms into either invaginated or evaginated state in the presence of the apical contractility patterns, depending on the direction of apicobasal polarity. The same effects are preserved when additional features, including natural curvature of the sheet, or the presence of an enclosed fluid, are added to the model.

Finally, the same modelling approach can be used to describe dynamics of cell and tissue deformations. This is done using an overdamped setting, in which the forces acting on each vertex are calculated by taking the partial derivatives of the energy function with respect to the coordinates of the vertex [25,27]. When two vertices become too close and are about to intersect, we change the local connectivity of the tissue. To this end, a discrete event, known as the T1 transition, is defined to resolve an edge that is about to vanish: a new edge is introduced, perpendicular to the vanishing edge and passing through its mid-point. This event changes the local connectivity within the cell sheet (figure 3*b*). If the vertices are free to move in three-dimensions, cells adjacent to the shrinking edge can be non-coplanar. In such cases, an average plane is defined in which the T1 transition takes place (as described in greater detail in the electronic supplementary material, section S1).

## 3. Modelling of dorsal appendage morphogenesis

In this section, we show how the proposed modelling framework can be used to describe morphogenesis of the respiratory appendages during *Drosophila* oogenesis, an established experimental model for studying the mechanisms by which patterned cell sheets give rise to complex three-dimensional structures.

The final product of *Drosophila* oogenesis is a single cell, the oocyte, surrounded by an elaborately patterned eggshell, a proteinaceous structure that houses the oocyte and the future embryo and mediates their interaction with the environment [32]. The eggshell plays a critical role in controlling the respiration of the embryo, ensuring an adequate supply of oxygen and preventing dehydration. In a number of *Drosophila* species, including *Drosophila melanogaster*, the eggshell is adorned with respiratory appendages, which can vary in shape and number, and aid respiration when the egg is buried in a soft oviposition substrate, like a rotting fruit [33,34]. The entire eggshell, including the respiratory appendages, is derived from the epithelial layer that surrounds the developing oocyte [24]. The apical surfaces of epithelial cells in this layer face the oocyte, while their basal surfaces adhere to a common extracellular matrix that surrounds the developing egg follicle. For several decades, the formation of dorsal appendages has served as an important model for genetic studies of tissue patterning and morphogenesis [24,35–38].

Most of our current understanding of this process has been derived from studies in *D. melanogaster*, which has an eggshell with two respiratory appendages. Midway through oogenesis, the follicular epithelium is patterned by an inductive signal that establishes two groups of appendage-producing cells. Each group consists of two cell types—a patch of ‘roof’ cells that eventually form the top of the tube and a single-cell-wide row of ‘floor’ cells that form the lower side of the tube (figure 4*a*). Each appendage primordium is made up of approximately 70 cells and is characterized by a specific pattern of increased apical contractility [25]. This pattern initiates a robust morphogenetic transformation, in which the appendage primordium first bends out from the epithelium and is then transformed, through an ordered sequence of cell rearrangements, into a conical three-dimensional structure (figure 4*b*1–*b*5), [25,39]).

Both steps of this process have been recently simulated using a two-dimensional vertex model [25]. In response to prepatterning, modelled through an increase in the energetic cost of areas and edges in a subset of cells, the primordium first buckles out of the sheet and then undergoes a stereotyped sequence of intercalations. The computational models were based on the vertex description of epithelial sheets that were more complex than equation (2.1), but cannot capture the three-dimensional nature of the epithelium. Here we show how the model, with the added term describing the apicobasal polarity, can be used to describe the dynamics of appendage formation with minimal patterning.

This transformation relies on spatial and temporal patterning of mechanical properties of cells in the model epithelium. The first type of patterning increases the contractility of the roof cells, causing their uniform constriction and localized bending of the primordium in the direction away from the oocyte. In the model, this is implemented with the newly added term (equation (2.3)). This type of patterning is motivated by the observed increased levels of myosin in the roof cells [25]. The second type of patterning increases the line tension along the floor–midline boundary. This is realized through spatially varying the value of *σ* in equation (2.1). Specifically, the value of sigma is peaked at the centre of this boundary, motivated by the observed peaked localization of myosin in this region of the follicular epithelium (figure 5*a*,*b*,*e*). This type of patterning is also supported by experimental observations of myosin localization.

To identify minimal prepatterning needed for ordered cell rearrangements, we start with a small primordium case with a hexagonal array of cells, where the number of appendage-producing cells in the model is less than in the real system. A representative example is shown in figure 5, which demonstrates how the peaked line tension term at the floor/midline border induces repeated cell neighbour exchanges, in which the line of the floor cells is progressively bent, as these cells lose their contacts with the midline cells. As a consequence of these cell rearrangements, the length of the interface between the floor and midline domains is decreased.

Our initial expectation was that these dynamics could be used to generate appendages from primordia independently of their size, as long as one dynamically re-centres the line tension at the midline/floor border. However, we found that the process in larger primordia fails and neighbour exchanges at the midline/floor border stall after the first few cell intercalations. We found that this problem can be repaired if the tissue outside the appendage primordium is made more compliant, which can be implemented by decreasing the values of *μ _{s}* and

*σ*in equation (2.1) in all cells excluding the roof and floor domains. With this additional assumption, the system could form the conical structure with cell–cell connectivity that closely resembles the connectivity revealed by reconstruction of the experimental images of dorsal appendages (figure 6).

In summary, we used our computational model of multicellular dynamics to determine the minimal patterning strategy needed for the transformation of the flat primordium into a three-dimensional conical structure. This strategy combines uniform increase of cell contractility within the primordium, dynamically maintained and peaked profile of cell contractility along the border of the primordium, as well as uniform change of cell properties outside the primordium. Our parametric studies revealed that the resulting morphogenesis is robust with respect to reasonable variations in cell properties and is only sensitive to the parameter G (width of the peaked cable); similar observations were also made in a previous study [25]. Using the model, one can readily explore the differential contributions of the three different prepatterning mechanisms. In the future, some of the model predictions, such as the uniform change in the compliance of the tissue outside the primordium, could be tested experimentally by direct measurement of cell–cell tensile forces. In parallel, one could computationally explore alternative mechanisms, such as directed migration of the roof cells, an effect that plays a key role during the later stages of dorsal appendage morphogenesis in *D. melanogaster* [39,40].

## 4. Discussion

As compared to models based on concepts of classical mechanics [41,42], cell-based models, like the vertex models and the cellular Potts models, are better suited to capture epithelial morphogenesis driven by cell shape changes and rearrangements [26,43]. At the same time, current models of dynamic cell shapes in epithelial sheets are still far behind the models of single cells, which can treat cell shapes using a free boundary value problem formulation [44]. Furthermore, describing some of the key aspects of cell dynamics in epithelial sheets, such as the emergence of wedge-like cell shapes in response to apical cell constriction, requires models that treat cells as three-dimensional objects [23,29]. We believe that the presented model, which is based on a vertex modelling framework, provides the simplest possible mathematical description that can simultaneously handle three-dimensional epithelial deformations and dynamic cell rearrangements. As was shown in our earlier study [23], the two-dimensional model that ignores apical localization of contractile segments is not consistent with the complete three-dimensional description of an epithelial monolayer. Addition of new terms that can effectively offset the contractile segments from the two-dimensional sheet can capture the essential effects of the three-dimensional model, while maintaining the simpler two-dimensional framework [23]. Furthermore, the same model can be used to describe morphogenetic processes that involve ordered cell rearrangements.

Already at this point, this class of vertex models can be used to explore the feasibility of morphogenetic mechanisms proposed on the basis of live imaging studies of epithelial dynamics. For example, in the case of the respiratory appendage formation in *D. melanogaster*, experiments suggested that the process is driven by the localized pattern of cell contractility that induces sequential neighbour exchanges leading to the formation of a conical structure [25]. Our computational exploration of this mechanism suggested that the localized pattern of contractility is not sufficient and must be accompanied by softening of the surrounding tissue. Without this added effect, the process stalls after just a few neighbour exchanges. Another source of the force that has been suggested to play a role in the appendage formation and can get rid of stalling in the simulations is crawling of follicle cells along the overlying basement membrane [39].

Epithelial morphogenesis accompanied by cell rearrangements can be viewed from two perspectives, kinematic and dynamic. First, the transformation between the initial and final structures, such as the flat primordium and fully formed appendage in our example, can be described by specifying an ordered sequence of neighbour exchanges. Second, one can ask about the forces that realize this sequence [45,46]. In our work, the dynamic mechanism is suggested by the spatial patterns of myosin localization and computational modelling. Importantly, the question about the kinematics of transitions between the two states with different cell adjacencies can be viewed separately. When the number of cells in the system is constant, the full lists of cell adjacencies in the initial and final configurations can be represented as graphs and one can ask whether two different graphs can be connected by a path in which each step is realized by a single T1 neighbour exchange. Formalizing this reachability problem and finding a way to solve it efficiently can provide insights into the mechanisms by which flat cell sheets can give rise to a wide range of three-dimensional structures. This purely kinematic approach may be especially useful in the analysis of experimental systems that are not readily amenable to live imaging. Also, when this kinematic path has been identified, the next step could be to use our model and identify a pattern of our model parameters that can drive the system along this path.

Although the present study emphasized mechanisms that rely on ordered cell rearrangements, one should keep in mind that real epithelial sheets can rely on strategies other than T1 exchanges. As an example, recent analysis of eggshell patterning in *Scaptodrosophila pattersoni* revealed that epithelial morphogenesis in this system relies not on ordered cell rearrangements, but on dramatic cell deformations, whereby the lengths of some cell edges become greatly elongated [39]. A modelling framework for describing these changes is yet to be developed.

## Authors' contributions

M.M. and S.Y.S. designed the research; M.M. performed the research. M.M., B.A. and S.Y.S. wrote the article.

## Competing interests

We declare we have no competing interests.

## Funding

S.Y.S. acknowledges support from grant no. 1R01GM107103 from the National Institute of General Medical Sciences. S.Y.S. and M.M. were partially supported by the National Science Foundation (NSF) Science and Technology Center Emergent Behavior of Integrated Cellular Systems (EBICS) grant no. CBET-0939511.

## Acknowledgements

We thank Yannis Kevrekidis, Alexander Fletcher, Chengcheng Gui, Miriam Osterfield, Satoru Okuda and Mona Singh for helpful discussions and comments on the manuscript. We also thank Miriam Osterfield for providing figure 4.

## Footnotes

One contribution of 13 to a theme issue ‘Systems morphodynamics: understanding the development of tissue hardware’.

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3689209.

- Accepted August 26, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.