All relevant data are within the paper and its Supporting Information files. In addition, the source code of SAMoS package can be found at <https://github.com/sknepneklab/SAMoS>.

Introduction {#sec001}
============

Collective cell migration \[[@pcbi.1005569.ref001], [@pcbi.1005569.ref002]\] in epithelial tissues is one of the key mechanisms behind many biological processes, such as the development of an embryo \[[@pcbi.1005569.ref003]\], wound healing \[[@pcbi.1005569.ref004], [@pcbi.1005569.ref005]\], tumour metastasis and invasion \[[@pcbi.1005569.ref006]\]. Due to their layered, tightly connected structure \[[@pcbi.1005569.ref007]\], epithelial tissues also serve as an excellent model system to study cell migration processes. Over several decades \[[@pcbi.1005569.ref008]\] extensive research efforts have been devoted to understanding molecular processes that lead to cell migration \[[@pcbi.1005569.ref009]\] and, at larger scales, on how cell migration drives complex processes at the level of the entire tissue, such as morphogenesis. With recent advances in various microscopy techniques combined with the development of sophisticated automatic cell tracking methods, it is now possible to study collective migration patterns of a large number of cells over extended periods of time with cell-level resolution, both *in vitro* and *in vivo*. Traction force microscopy \[[@pcbi.1005569.ref010]\] experiments revealed that collective cell motion is far more complex than expected \[[@pcbi.1005569.ref011]--[@pcbi.1005569.ref013]\]. It is often useful to draw parallels between the collective behaviour of tissues and systems studied in the physics of colloids, granular materials and foams as these can provide powerful tools for understanding complex cell interactions in biological systems. For example, in a homogenous cell monolayer, one observes large spatial and temporal fluctuations of inter-cellular forces that cannot be pinpointed to a specific location, but cover regions that extend over several cells \[[@pcbi.1005569.ref014]\]. These are reminiscent of the fluctuations observed in supercooled colloidal and molecular liquids approaching the glass transition \[[@pcbi.1005569.ref011]\] and include characteristic features of the dynamical and mechanical response, such as dynamical heterogeneities and heterogeneous stress patterns, that were first observed in glasses, colloids and granular materials and that have extensively been studied in soft condensed matter physics \[[@pcbi.1005569.ref015]\]. It has also been argued that the migration patterns are sensitive to the expression of different adhesion proteins \[[@pcbi.1005569.ref016]\] as well as to the properties of the extracellular environment \[[@pcbi.1005569.ref017]\], such as the stiffness of the substrate \[[@pcbi.1005569.ref018], [@pcbi.1005569.ref019]\]. These observations lead to the development of the notion of *plithotaxis* \[[@pcbi.1005569.ref012]\], a universal mechanical principle akin to the more familiar chemotaxis, which states that each cell tends to move in a way that maintains minimal local intercellular shear stress. While plausible, it is yet to be determined whether plithotaxis is indeed a generic feature in all epithelial tissues.

Equally fascinating are the experiments on model systems that study cell migration in settings designed to mimic wound healing \[[@pcbi.1005569.ref005], [@pcbi.1005569.ref020]--[@pcbi.1005569.ref023]\]. For example, the existence of mechanical waves that span the entire tissue and generate long-range cell-guidance have been established in Madin-Darby Canine Kidney (MDCK) epithelial cell monolayers \[[@pcbi.1005569.ref023]\]. Subtle correlations between purse-string contractility and large-scale remodelling of the tissue while closing circular gaps have also been identified \[[@pcbi.1005569.ref022]\]. Finally, a mechanism dubbed *kenotaxis* has been proposed \[[@pcbi.1005569.ref020]\], which suggests that there is a robust tendency of a collection of migrating cells to generate local tractions that systematically and cooperatively pull towards the empty regions of the substrate.

On the developmental side, in pioneering work, Keller *et al*. \[[@pcbi.1005569.ref024]\] constructed a light-sheet microscope that enabled them to track *in vivo* positions of each individual cell in a zebrafish embryo over a period of 24h. A quantitive analysis \[[@pcbi.1005569.ref025]\] of the zebrafish embryo was also able to relate mechanical energy and geometry to the shapes of the aggregate surface cells. Another extensively studied system that allows detailed tracking of individual cells is the *Drosophila* embryo \[[@pcbi.1005569.ref026]--[@pcbi.1005569.ref030]\]. In recent studies that combined experiments with advanced data analysis, it was possible to quantitatively account for shape change of the wing blade by decomposing it into cell divisions, cell rearrangements and cell shape changes \[[@pcbi.1005569.ref031], [@pcbi.1005569.ref032]\]. Finally, it has recently become possible to track more than 100,000 individual cells in a chick embryo over a period exceeding 24 hours \[[@pcbi.1005569.ref033]\]. This was achieved by developing an advanced light-sheet microscope and state-of-the-art data analysis techniques designed to automatically track individual cells in a transgenic chick embryo line with the cell membranes of all cells in the embryonic and extra embryonic tissues labelled with a green fluorescent protein tag. All these experiments and advanced data analysis techniques provide unprecedented insights into the early stages of embryonic development, making it possible to connect processes at the level of individual cells with embryo-scale collective cell motion patterns.

While there have been great advances in our understanding of how cells regulate force generation and transmission between each other and with the extracellular matrix in order to control their shape and cell-cell contacts \[[@pcbi.1005569.ref009]\], it is still not clear how these processes are coordinated at the tissue-level to drive tissue morphogenesis or allow the tissue to maintain its function once it reaches homeostasis. Computer models of various levels of complexity have played an essential role in helping to answer many of those questions \[[@pcbi.1005569.ref034]\]. One of the early yet successful approaches has been based on the cellular Potts model (CPM) \[[@pcbi.1005569.ref035], [@pcbi.1005569.ref036]\]. In the CPM, cells are represented as groups of "pixels" of a given type. Pixels are updated one at a time following a set of probabilistic rules. Pixel updates in the CPM are computationally inexpensive \[[@pcbi.1005569.ref037]\], which allows for simulations of large systems. In addition, the CPM extends naturally to three dimensions \[[@pcbi.1005569.ref038]\]. While very successful in describing cell sorting as well as certain aspects of tumour growth \[[@pcbi.1005569.ref039]\], CPM has several limitations, the main one being that the dynamics of pixel updates is somewhat artificial and very hard to relate to the dynamics of actual cells. This problem has been to some extent alleviated by the introduction of the subcellular element model (ScEM) \[[@pcbi.1005569.ref040], [@pcbi.1005569.ref041]\]. ScEM is an off-lattice model with each cell being represented as a collection of 100--200 elements---spherical particles interacting with their immediate neighbours via short-range soft-core potentials. Therefore, ScEM is able to model cells of arbitrary shapes that grow, divide and move in 3D. The main disadvantage of ScEM is that it is computationally expensive (off-lattice methods in general require more computations per time step compared to their lattice counterparts), and without a highly optimised parallel implementation, applications of the ScEM are limited to a few hundred cells at most, which is not sufficient to study effects that span long length- and time-scales.

Particle-based models have also been very successful at capturing many aspects of cell migration in tissues \[[@pcbi.1005569.ref042]--[@pcbi.1005569.ref046]\]. However, when it comes to modelling confluent epithelial layers with the resolution of individual cells, one of the most widely and successfully used models is the Vertex Model (VM). The VM originated in studies of the physics of foams in the 1970s and was first applied to model monolayer cell sheets in the early 1980s \[[@pcbi.1005569.ref047]\]. Over the past 35 years it has been implemented and extended numerous times and used to study a wide variety of different systems \[[@pcbi.1005569.ref048]--[@pcbi.1005569.ref053]\]. The VM is at the core of the cell-based \[[@pcbi.1005569.ref054]\] CHASTE \[[@pcbi.1005569.ref055]\], a versatile and widely used software package. The VM assumes that all cells in the epithelium are roughly the same height and thus that the entire system can be well approximated as a two-dimensional sheet. The conformation of the tissue in the VM is computed as a configuration that simultaneously optimises area and perimeter of all cells. While the model is quasi-static in nature, it captures remarkably well many properties of actual epithelial tissues. There have been numerous attempts to introduce dynamics into the vertex model \[[@pcbi.1005569.ref047], [@pcbi.1005569.ref051], [@pcbi.1005569.ref053], [@pcbi.1005569.ref056], [@pcbi.1005569.ref057]\]. However, there are limitations associated with each approach. To the best of our knowledge, most dynamical versions of the vertex model seem to neglect fluctuations with a notable recent exception \[[@pcbi.1005569.ref053]\]. In contrast, recent traction microscopy experiments \[[@pcbi.1005569.ref014]\] suggest that these fluctuations might be a crucial ingredient in understanding collective cell migration. Finally, we point out a technical point that makes the implementation of VM somewhat challenging. Namely, in order to capture topology changing moves, such as cell neighbour exchanges, i.e., T1 transitions, one has to perform rather complex mesh restructuring operations \[[@pcbi.1005569.ref051], [@pcbi.1005569.ref058]\] that require complex data structures and algorithms and that inevitably add to the computational complexity of the model.

Building upon the recently introduced Self-Propelled Voronoi (SPV) model \[[@pcbi.1005569.ref059]\], in this paper we apply the ideas introduced in the physics of active matter systems \[[@pcbi.1005569.ref060]\] to the VM. This allows us to construct a hybrid Active Vertex Model (AVM) that is able to accurately describe the collective migration dynamics of a large number of cells. The AVM is implemented within *SAMoS* \[[@pcbi.1005569.ref061]\], an off-lattice particle-based simulation software developed specifically to study active matter systems. One of the key advantages of the hybrid approach presented in this study is that it not only enables studies of very large systems, but also introduces a very natural way to handle the T1 transitions, thus removing the need for complex mesh manipulations that are of uncertain physical and biological meaning.

Overview of the Vertex Model {#sec002}
----------------------------

Owing to its origins in the physics of foams \[[@pcbi.1005569.ref062]\], in the VM cells are modelled as two-dimensional convex polygons that cover the plane with no holes and overlaps, i.e., the epithelial tissue is represented as a convex polygonal partitioning of the plane ([Fig 1](#pcbi.1005569.g001){ref-type="fig"}). The main simplification compared to models of foams is that most implementations of the VM assume that contacts between neighbouring cells are straight lines. However, we note that there have been several recent studies where this assumption has been removed and cell-cell junctions were allowed to be curved \[[@pcbi.1005569.ref063], [@pcbi.1005569.ref064]\]. In addition, neighbouring cells are also assumed to share a single edge, which is a simplification compared to real tissues, where junctions between two neighbouring cells consist of two separate cell membranes that can be independently regulated. Typically, three junction lines meet at a vertex, although vertices with a higher number of contacts are also possible \[[@pcbi.1005569.ref058]\]. The model tissue is therefore a mesh consisting of polygons (i.e., cells), edges (i.e., cell junctions), and vertices (i.e., meeting points of three or more cells).

![In the Vertex Model (VM), a confluent epithelial sheet is represented as a polygonal tiling of the plane with no holes or overlaps.\
Each cell is represented by an *n*−sided polygon. Neighbouring cells share an edge, which models the cell junction as a straight line. Three edges meet at a vertex (dark blue dots). The behaviour of cell *i* is described by three parameters: 1) reference area $A_{i}^{0}$, 2) area modulus *K*~*i*~, and 3) perimeter modulus Γ~*i*~. In addition, a junction connecting vertices *μ* and *ν* has tension Λ~*μν*~.](pcbi.1005569.g001){#pcbi.1005569.g001}

An energy is associated to each configuration of the mesh. It can be written as $$E_{VM} = \sum\limits_{i = 1}^{N}\mspace{180mu}\frac{K_{i}}{2}\left( {A_{i} - A_{i}^{0}} \right)^{2} + \sum\limits_{i = 1}^{N}\mspace{180mu}\frac{\Gamma_{i}}{2}P_{i}^{2} + 2\sum\limits_{\langle{\mu,\nu}\rangle}\mspace{180mu}\Lambda_{\mu\nu}l_{\mu\nu},$$ where *N* is the total number of cells, *A*~*i*~ is the area of the cell *i*, while $A_{i}^{0}$ is its reference area. *K*~*i*~ is the area modulus, i.e. a constant with units of energy per area squared measuring how hard it is to change the cell's area. *P*~*i*~ is the cell perimeter and Γ~*i*~ (with units of energy per length squared) is the perimeter modulus that determines how hard it is to change perimeter *P*~*i*~. *l*~*μν*~ is the length of the junction between vertices *μ* and *ν* and Λ~*μν*~ is the tension of that junction (with units of energy per length). 〈*μ*, *ν*〉 in the last term denotes the sum is over all pairs of vertices that share a junction. Note that the model allows for different cells to have different area and perimeter moduli as well as reference areas, allowing for modelling of tissues containing different cell types. Finally, for convenience we introduced a prefactor 2 in the last term in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} to compensate for double counting of cell junctions when switching from a sum over junctions to a sum over cells in the force calculation discussed below.

Some authors write the perimeter term as $\sum_{i}{\widetilde{\Gamma}}_{i}/2\mspace{180mu}\left( P_{i} - P_{i}^{0} \right)^{2}$, where $P_{i}^{0}$ is a reference perimeter, and omit the last term in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} or completely omit the *P*^2^ term \[[@pcbi.1005569.ref047], [@pcbi.1005569.ref058]\]. Under the assumption that the values of Λ~*μν*~ for all junctions of the cell *i* are the same, i.e., Λ~*μ*,*ν*~ ≡ Λ~*i*~, the last term in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} becomes Λ~*i*~ ∑~〈*μ*,*ν*〉~ *l*~*μ*,*ν*~ = Λ~*i*~*P*~*i*~. Therefore, if we identify $\Lambda = - \widetilde{\Gamma}P_{i}^{0}$, it immediately becomes clear that the descriptions in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} and the model with the preferred perimeter are identical to each other. Note that this is true up to the constant term $1/2{\widetilde{\Gamma}}_{i}\mspace{180mu}\left( P_{i}^{0} \right)^{2}$, which is unimportant as it only shifts the overall energy and does not contribute to the force on the cell (see below). The description in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} is slightly more general as it allows for the junctions to have different properties depending on the types of cells that are in contact. Retaining the $P_{i}^{2}$ term is also advisable in order to prevent the model from becoming unstable if the area modulus is too small.

It is straightforward to express cell area and cell perimeter in terms of vertex coordinates. Therefore, vertex positions together with their connectivities uniquely determine the energy of the epithelial sheet, hence the name Vertex Model. The main assumption of the VM is that the tissue will always be in a configuration which minimises the total energy in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"}. Determining the minimum energy configuration is a non-trivial multidimensional optimisation problem and, with the exception of a few very simple cases, it can only be solved numerically. A basic implementation of the VM therefore needs to use advanced multidimensional numerical minimisation algorithms to determine the positions of vertices that minimise [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} for a given set of parameters *K*~*i*~, Γ~*i*~ and Λ~*μν*~. Most implementations \[[@pcbi.1005569.ref051], [@pcbi.1005569.ref053], [@pcbi.1005569.ref058]\], also introduce topology (connectivity) changing moves to model events such as cell rearrangements.

There have been several attempts to introduce dynamics into the VM \[[@pcbi.1005569.ref047], [@pcbi.1005569.ref051], [@pcbi.1005569.ref056]\], including a recent study that introduced stochasticity into the junction tension \[[@pcbi.1005569.ref053]\]. The idea behind such approaches it to write equations of motion for each vertex as $$\gamma\frac{d\mathbf{r}_{\mu}}{dt} = \mathbf{F}_{\mu},$$ where *γ* is a friction coefficient and **r**~*μ*~ is the position vector of vertex *μ*. **F**~*μ*~ is the total force on vertex *μ* computed as the negative gradient of [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} with respect to **r**~*μ*~, i.e., **F**~*μ*~ = −∇~**r**~*μ*~~ *E*~*VM*~. We note that the exact meaning of friction in confluent epithelial tissues is the subject of an ongoing debate that is beyond the scope of this study. Here, as in the case of most models to date, we assume that all effects of friction (i.e., between neighbouring cells as well as between cells and the substrate and the extracellular matrix) can be modelled by a single constant. While this may appear to be a major simplification, as we will show below, the model is capable of capturing many key features of real epithelial tissues. [Eq (2)](#pcbi.1005569.e009){ref-type="disp-formula"} is a first order equation since the mass terms have been omitted. This so-called overdamped limit is commonly applied to biological systems, since the inertial effects are typically several orders of magnitude smaller than the effects arising from the cell-cell interactions or random fluctuations produced by the environment. Note that the force on vertex *μ* depends on the position of its neighbouring vertices, resulting in a set of coupled non-linear ordinary differential equations. In most cases those equations can only be solved numerically.

While the introduction of dynamics alleviates some of the problems related to the quasi-static nature of the VM, one still has to implement topology changing moves if the model is to be applicable to describing cell intercalation events. This can lead to unphysical back and forth flips of the same junction and has only recently been analysed in full detail \[[@pcbi.1005569.ref058]\].

Model {#sec003}
=====

Active Vertex Model {#sec004}
-------------------

It is important to note that the VM in its original form is a *quasi-static* model. In other words, it assumes that at every instant in time, the tissue is in a state of mechanical equilibrium. This is a strong assumption, which is in line with many biological systems, especially in the case of embryos where cells actively grow, divide and rearrange. As a matter of fact, biological systems are among the most common examples of systems out of equilibrium. Therefore, while it is able to capture many of the mechanical properties of the tissue, the VM is unable to fully describe the effects that are inherently related to being out of equilibrium. Many such effects are believed to be behind the collective migration patterns observed in recent experiments. In addition, in many dynamical implementations of the VM the effects of both thermal and non-thermal random fluctuations originating in complex intercellular processes and interactions with the environment are either completely omitted or not very clear. While for a system out of equilibrium the fluctuation-dissipation theorem \[[@pcbi.1005569.ref065]\] does not hold, and the relation between random fluctuations and friction is not simple, it is even more important to note that fluctuations can have non-trivial effects on the collective motion patterns \[[@pcbi.1005569.ref053]\].

Here we take an alternative approach and build a description similar to the recently introduced SPV model \[[@pcbi.1005569.ref059]\]. The idea behind the SPV is that instead of treating vertices as the degrees of freedom, one tracks positions of cell centres. Forces on cell centres are, however, computed from the energy of the VM, [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"}. The core assumption of the model is that the tissue conformations correspond to the Voronoi tessellations of the plane with cell centres acting as Voronoi seeds. We recall that a Voronoi tessellation is a polygonal tiling of the plane based on distances to a set of points, called *seeds*. For each seed point there is a corresponding polygon consisting of all points closer to that seed point than to any other. This imposes some restrictions onto possible tissue conformations, i.e., not all convex polygonal tessellations of the plane are necessarily Voronoi, but it has recently been argued that Voronoi tessellations can predict the diverse cell shape distributions of many tissues \[[@pcbi.1005569.ref066]\]. Furthermore, the exact details of the tessellation are not expected to play a significant role in the large scale behaviour of the tissue, which this model aims to describe. We, however, note that recently an interesting model based on a generalised Voronoi description has been proposed \[[@pcbi.1005569.ref067]\].

In the original implementation of Bi, *et al*. \[[@pcbi.1005569.ref059]\] the Voronoi tessellation is computed at every time step. The vertices of the tessellation are then used to evaluate forces at all cell centres, that are, in turn, moved in accordance to those forces and the entire process is repeated. While conceptually clear, this procedure is numerically expensive as it requires computation of the entire Voronoi diagram at each time step. This limits the accessible system size to several hundred cells.

Here, we instead propose an alternative approach based on the Delaunay triangulation. The Delaunay triangulation for a set of points *P* in the plane is a triangular tiling, *DT* (*P*), of the plane with the property that there are no points of *P* inside the circumcircle of any of the triangles in *DT* (*P*) \[[@pcbi.1005569.ref068]\]. A property of a Delaunay triangulation that is key for this work is that it is possible to construct a so-called dual Voronoi tessellation by connecting circumcenters of its triangles. This establishes a mathematical duality between Delaunay and Voronoi descriptions. This duality is exact and quantities, such as the force, expressed on the Voronoi tiling have an exact map onto quantities expressed on its dual Delaunay triangulation. Although being non-linear (see Sec. "Force on the cell centre"), this map is relatively simple, and therefore fast to compute. An important property of the Voronoi-Delaunay duality is that continuous deformations of one map into continuous deformations of the other. In other words, smooth motion of a cell's centre will correspond to a smooth change in that cell's shape. This is crucial to ensuring that during the dynamical evolution the cell connectivity changes continuously, a feature that is essential for accurately modelling T1 transitions. The main advantage of working with the Delaunay description is that while the Voronoi tessellation has to be recomputed each time cell centres move, it is possible to retain the Delaunay character of a triangulation via local edge flip moves ([Fig 2c](#pcbi.1005569.g002){ref-type="fig"}), which drastically increases the efficiency of the Delaunay based approach and enables us to simulate systems containing tens of thousands of cells.

![(a) Two coordinate representations dual to each other. In the AVM we track particles that correspond to cell centres (red dots). Particles form a Delaunay triangulation (red lines) and their positions are labeled with vectors carrying Latin indices. The dual of the Delaunay triangulation is a Voronoi tessellation, with each Voronoi cell (a shaded polygon) representing an actual cell. Cell edges are marked by blue lines. Positions of the vertices (blue dots) of the dual mesh are denoted by Greek indices. (b) Position of the circumcenter **r**~*μ*~ of the triangle with corners **r**~*i*~, **r**~*j*~ and **r**~*k*~. (c) The edge flip move is at the core of the equiangulation procedure. If the sum of the angles opposite to the red edge exceeds *α* + *β* \> 180° the edge is "flipped". As the result, the sum of the angles opposite to the new edge is less than 180°, i.e., *γ* + *δ* \< 180°. Note that this edge flip is local and only affects one Voronoi edge (i.e., cell junction). Therefore, edge flips can only affect the local connectivity of four cells (see also, [Fig 5](#pcbi.1005569.g005){ref-type="fig"}).](pcbi.1005569.g002){#pcbi.1005569.g002}

Before we introduce the Active Vertex Model (AVM), we pause to make a comment about the notation. In the following, we will always use Latin letters to denote cells, i.e. positions of their centres, and Greek letters to denote vertices of the dual Voronoi tessellation, i.e., meeting points of three or more cells. Therefore, vertices of the VM will always carry Greek indices ([Fig 2a](#pcbi.1005569.g002){ref-type="fig"}).

### Force on the cell centre {#sec005}

In the AVM, the area *A*~*i*~ in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} corresponding to the cell (particle) *i* is the area of the polygon (face) of the Voronoi polygon, Ω~*i*~, associated with *i*. It is given as $$A_{i} = \frac{1}{2}\mspace{180mu}\sum\limits_{\mu \in \Omega_{i}}\mspace{180mu}\left( {\mathbf{r}_{\mu} \times \mathbf{r}_{\mu + 1}} \right) \cdot \mathbf{N}_{i},$$ where **r**~*μ*~ is the position of vertex *μ* and **N**~*i*~ ≡ **e**~*z*~ is a unit-length vector perpendicular to the surface of the polygon---which does not depend on the position of the vertices in the plane. The sum is over all vertices of the Voronoi cell and we close the loop with *μ* + 1 ≡ 1 for *μ* = *N*~Ω~*i*~~, where *N*~Ω~*i*~~ is the total number of vertices in the cell Ω~*i*~. This expression is just a discrete version of Green's formula. Similarly, the cell perimeter is $$P_{i} = \sum\limits_{\mu \in \Omega_{i}}\mspace{180mu}\left| {\mathbf{r}_{\mu + 1} - \mathbf{r}_{\mu}} \right|,$$ with the same rules for the next neighbour of the "last" vertex. \|⋅\| represents the Euclidean length of the vector.

In order to make a connection between cell centres and indices of the Voronoi tessellation we recall that for a given triangle in the Delaunay triangulation, the position of its dual Voronoi vertex coincides with the centre of the circumscribed circle. The position of the circumcenter is given by (see [Fig 2b](#pcbi.1005569.g002){ref-type="fig"}) $$\mathbf{r}_{\mu} = \frac{\lambda_{1}}{\Lambda}\mathbf{r}_{i} + \frac{\lambda_{2}}{\Lambda}\mathbf{r}_{j} + \frac{\lambda_{3}}{\Lambda}\mathbf{r}_{k},$$ where **r**~*i*~, **r**~*j*~ and **r**~*k*~ are position vectors of the corners of the triangle and *λ*~1~, *λ*~2~ and *λ*~3~ are the barycentric coordinates (see [S1 Appendix](#pcbi.1005569.s001){ref-type="supplementary-material"}), with Λ = *λ*~1~ + *λ*~2~ + *λ*~3~.

Armed with the mapping in [Eq (5)](#pcbi.1005569.e012){ref-type="disp-formula"}, we can proceed to compute the forces acting on each cell centre based on the VM energy functional in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"}. This is done directly by computing the negative gradient of the energy. Therefore, the force on the centre of cell *i* is computed as $$\mathbf{F}_{i} = - \nabla_{\mathbf{r}_{i}}E_{VM}.$$ The expression for the gradient is made somewhat complicated by the fact that the derivative is taken with respect to the position of the cell centre, while the energy of the VM is naturally written in terms of the positions of the dual vertices. After a lengthy but straightforward calculation (see [S1 Appendix](#pcbi.1005569.s001){ref-type="supplementary-material"}) we obtain: $$\begin{array}{cl}
\mathbf{F}_{i} & {= - \sum\limits_{k = 1}^{N}\mspace{180mu}\frac{K_{k}}{2}\left( {A_{k} - A_{k}^{0}} \right)\sum\limits_{\nu \in \Omega_{k}}\mspace{180mu}\left\lbrack {\mathbf{r}_{\nu + 1,\nu - 1} \times \mathbf{N}_{k}} \right\rbrack^{T}\mspace{180mu}\left\lbrack \frac{\partial\mathbf{r}_{\nu}}{\partial\mathbf{r}_{i}} \right\rbrack} \\
 & {- \sum\limits_{k = 1}^{N}\mspace{180mu}\Gamma_{k}P_{k}\sum\limits_{\nu \in \Omega_{k}}\mspace{180mu}\left( {{\hat{\mathbf{r}}}_{\nu,\nu - 1} - {\hat{\mathbf{r}}}_{\nu + 1,\nu}} \right)^{T}\mspace{180mu}\left\lbrack \frac{\partial\mathbf{r}_{\nu}}{\partial\mathbf{r}_{i}} \right\rbrack} \\
 & {- \sum\limits_{k = 1}^{N}\mspace{180mu}\sum\limits_{\nu \in \Omega_{k}}\mspace{180mu}\left\lbrack {\Lambda_{\nu - 1,\nu}{\hat{\mathbf{r}}}_{\nu,\nu - 1} - \Lambda_{\nu,\nu + 1}{\hat{\mathbf{r}}}_{\nu + 1,\nu}} \right\rbrack^{T}\mspace{180mu}\left\lbrack \frac{\partial\mathbf{r}_{\nu}}{\partial\mathbf{r}_{i}} \right\rbrack.} \\
\end{array}$$ In the last expression, $\left\lbrack \frac{\partial\mathbf{r}_{\nu}}{\partial\mathbf{r}_{i}} \right\rbrack$ is the 3 × 3 Jacobian matrix (see [S1 Appendix](#pcbi.1005569.s001){ref-type="supplementary-material"}) connecting coordinates of cell centres with coordinates of the dual Voronoi tessellation and \[⋅\]^*T*^ \[⋅\] denotes a row-matrix product producing a 3 × 1 column vector. Note that this product does not commute, i.e., the order in which terms appear in the expression above is important. *N* is the total number of cells in the system.

We note that most terms in the *k* sum in [Eq (7)](#pcbi.1005569.e014){ref-type="disp-formula"} are equal to zero, since each vertex coordinate **r**~*ν*~ depends only the cell centres **r**~*k*~ associated to its Delaunay triangle (see [Fig 2b](#pcbi.1005569.g002){ref-type="fig"}). In other words, we only need to consider cell *i* and its immediate neighbours. For clarity, we outline the algorithm for computing the area term and note that perimeter and junction terms can be treated in a similar fashion.

1.  For particle *i* compute $\frac{K_{i}}{2}\mspace{180mu}\left( A_{i} - A_{i}^{0} \right)$ and multiply it by the sum $\sum_{\nu \in \Omega_{i}}\mspace{180mu}\left\lbrack \left( \mathbf{r}_{\nu + 1} - \mathbf{r}_{\nu - 1} \right) \times \mathbf{N}_{i} \right\rbrack^{T}\mspace{180mu}\left\lbrack \frac{\partial\mathbf{r}_{\nu}}{\partial\mathbf{r}_{i}} \right\rbrack$. This sum is over all vertices (corners) *ν* of cell *i*.

2.  For all immediate neighbours *j* of cell *i* compute $\frac{K_{j}}{2}\mspace{180mu}\left( A_{j} - A_{j}^{0} \right)$ and multiply it by the sum $$\sum\limits_{\nu \in \Omega_{i} \cap \Omega_{j}}\mspace{180mu}\left\lbrack {\left( {\mathbf{r}_{\nu + 1} - \mathbf{r}_{\nu - 1}} \right) \times \mathbf{N}_{j}} \right\rbrack^{T}\mspace{180mu}\left\lbrack \frac{\partial\mathbf{r}_{\nu}}{\partial\mathbf{r}_{i}} \right\rbrack.$$ Note that *ν* ∈ Ω~*i*~ ∩ Ω~*j*~ ensures that vertices *ν* surrounding *j* are taken into account only if they are affected by (i.e., also belong to) cell *i*.

From the previous discussion it is clear that the expression for the force is local, i.e., computing the force does not require including cell centre positions beyond the immediate neighbourhood of a given cell. This is extremely beneficial from a computational point of view as one can readily utilise standard force cutoff techniques, such as cell and neighbour lists \[[@pcbi.1005569.ref069]\] in order to speed up force computations. It is also evident that the force in [Eq (7)](#pcbi.1005569.e014){ref-type="disp-formula"} is not pairwise, i.e., it cannot be written as a sum of forces acting on a pair of cells, or in mathematical terms, **F**~*i*~ ≠ ∑~*j*~ **F**~*ij*~. This is not surprising since the position of the Voronoi vertices depends on the position of three cell centres.

### Cell alignment {#sec006}

The AVM also includes the *cell polarity* \[[@pcbi.1005569.ref070]\] modelled as a unit length vector, **n**~*i*~, laying in the *xy* plane. This vector determines the direction of the cell's motion and division and should not be confused with apical-basal polarity, which is not included in this model. In the simplest case, the vector **n**~*i*~ does not take any preferred direction, but instead randomly fluctuates under the influence of uncorrelated random noise originating from the intercellular processes and the environment. This simple model can be augmented to include several different models for cell-cell alignment. For example, a cell can be given a tendency to align its polarity to its immediate neighbours by minimising the alignment energy, $$E_{polar\, align} = - J_{p}\sum\limits_{j\, n.n\, i}\mspace{180mu}\mathbf{n}_{i} \cdot \mathbf{n}_{j},$$ where *J*~*p*~ \> 0 is the alignment strength and the sum is over all nearest neighbours of *i*. This model is similar to the Vicsek model, which is used in studies of flocking, e.g. of swarms of birds \[[@pcbi.1005569.ref071]\]. The molecular details of polar alignment are not yet understood in detail. It as been argued that in some tissue systems this involves components of the planar cell polarity signalling cascade \[[@pcbi.1005569.ref072], [@pcbi.1005569.ref073]\], but even here their precise coordination and involvement are as yet not resolved.

Alternative alignment mechanisms consider only internal processes in the cell and do not directly depend on the polarity of the neighbouring cells. One such mechanism introduced in physics of active matter \[[@pcbi.1005569.ref074], [@pcbi.1005569.ref075]\] assumes that cell polarisation vector aligns with the migration direction of the cell, i.e., $$E_{vel\, align} = - J_{v}\ \mathbf{n}_{i} \cdot {\hat{\mathbf{v}}}_{i},$$ where ${\hat{\mathbf{v}}}_{i}$ is the normalised cell velocity vector and *J*~*v*~ \> 0 is the alignment strength.

Finally, **n**~*i*~ can also align to the direction of the eigenvector of the cell shape tensor, i.e. the cell seeks align direction of **n**~*i*~ along its long axis. This is achieved by minimising the value of the energy $$E_{shape\, align} = - J_{s}\mathbf{n}_{i} \cdot \mathbf{p}_{i},$$ where *J*~*s*~ \> 0 is the alignment strength and **p**~*i*~ is the eigenvector of the matrix \[[@pcbi.1005569.ref076]\] $${\hat{P}}_{i} = \frac{1}{N_{\Omega_{i}}}\mspace{180mu}\sum\limits_{\mu \in \Omega_{i}}\mspace{180mu}\mathbf{r}_{\mu} \otimes \mathbf{r}_{\mu + 1} = \frac{1}{N_{\Omega_{i}}}\mspace{180mu}\sum\limits_{\mu \in \Omega_{i}}\mspace{180mu}\begin{pmatrix}
\left( {x_{\mu} - x_{\mu + 1}} \right)^{2} & {\left( {x_{\mu} - x_{\mu + 1}} \right)\left( {y_{\mu} - y_{\mu + 1}} \right)} \\
{\left( {x_{\mu} - x_{\mu + 1}} \right)\left( {y_{\mu} - y_{\mu + 1}} \right)} & \left( {y_{\mu} - y_{\mu + 1}} \right)^{2} \\
\end{pmatrix},$$ corresponding to its largest eigenvalue. This mechanism is one possible pathway to obtain *plithotaxis*.

Each of these alignment mechanisms causes a torque on **n**~*i*~, given by $$\tau_{i} = - \mathbf{n}_{i} \times \nabla_{\mathbf{n}_{i}}E_{align}.$$ In the *SAMoS* implementation of the AVM, one can select which of these alignment mechanisms are to be included in an actual simulation.

### Equations of motion {#sec007}

We proceed to write equations of motion for the position and orientation of cell *i* in the AVM. As discussed in Sec. "Overview of the Vertex Model", at the cellular scale inertial effects are negligible. In this overdamped limit, the equation of motion for the position of the cell centre is just force balance between friction and driving forces. If we assume that the friction of the cell with its surrounding and the substrate is isotropic and can be modelled by the single friction coefficient *γ*, the equation of motion for the position of cell centre *i* becomes $$\gamma\frac{d\mathbf{r}_{i}}{dt} = f_{\text{a}}\mathbf{n}_{i} + \mathbf{F}_{i} + \mathbf{\nu}_{i}\left( t \right),$$ where **F**~*i*~ is the sum of all forces acting on the cell, i.e., a sum of the VM model forces given in [Eq (7)](#pcbi.1005569.e014){ref-type="disp-formula"} and the soft-core repulsion defined in Eq. (35) in [S1 Appendix](#pcbi.1005569.s001){ref-type="supplementary-material"}. $\mathbf{\nu}_{i}\left( t \right)$ is an uncorrelated stochastic force that models the effects of random fluctuations. These fluctuations originate from the intracellular processes and interactions with the environment. In its simplest possible form, one can assume that $\mathbf{\nu}_{i}\left( t \right)$ has no time correlations, i.e., each cell experiences a random "kick" at a given instant of time whose magnitude and direction do not depend on any past "kicks" on that cell. Formally, we write 〈*ν*~*i*,*α*~ (*t*)〉 = 0 and 〈*ν*~*i*,*α*~ (*t*) *ν*~*j*,*β*~ (*t*′)〉 = 2*Dδ*~*ij*~ *δ*~*αβ*~ *δ* (*t* − *t*′), where *D* measures the average magnitude of the stochastic force and is interpreted as an effective translational diffusion coefficient, *α*, *β* ∈ {*x*, *y*} and 〈...〉 denotes the statistical average. We note that in some cases, such as in the early *Drosophila* embryo, or the cultured monolayer system \[[@pcbi.1005569.ref077]\] there is little or no extracellular matrix or substrate and this friction model would have to be modified. While the inclusion of cell-cell friction is beyond the scope of this study, in the future, we plan to extend the AVM to include such an effect. One important application of the cell-cell friction will be to model the chick embryo \[[@pcbi.1005569.ref033]\], in which a near monolayer of epithelial tissue floats on liquid in the subgerminal space that provides much less friction than a solid substrate.

Finally, the term *f*~a~**n**~*i*~ has its origins in the physics of active matter systems \[[@pcbi.1005569.ref060]\]. It models *activity*, i.e., internal cellular processes that drive it to move in the direction of its polarity. The active force strength *f*~a~ controls the magnitude of this activity. While the biological meaning of *f*~a~ may appear unclear, it quantifies a cell's ability to move on its own due to the complex molecular machinery within it. In many models of individual cells crawling on a substrate with a prominent lamellipodium, the resultant active velocity *v*~0~ = *f*~a~/*γ* is due to either actin tread-milling \[[@pcbi.1005569.ref078]\], differential friction or differential contractility. In an epithelium, on the other hand, *f*~a~ should be understood as an effective parameter that models the non-balanced active force due to junction contractions and internal cell contractility \[[@pcbi.1005569.ref079]\]. Although the molecular origin of *f*~*a*~ is at present not fully understood, even this simple description of the activity combined with the cell polarity alignment models discussed in the previous subsection leads to collective behaviour that resembles behaviour observed in experiments.

We then write equations for motion for the polarity vector. If we define the angle of vector **n**~*i*~ with the *x*−axis as *θ*~*i*~, then **n**~*i*~ = (cos *θ*~*i*~, sin *θ*~*i*~) and we have $$\gamma_{r}\frac{d\theta_{i}}{dt} = \mathbf{\tau}_{i} \cdot \mathbf{N}_{i} + \nu_{i}^{r}\mspace{180mu}\left( t \right),$$ where ***τ***~*i*~ is the torque acting on **n**~*i*~, given in [Eq (12)](#pcbi.1005569.e025){ref-type="disp-formula"}, **N**~*i*~ is the local normal to the cell surface (i.e., simply the unit length vector in the *z*-direction), *γ*~*r*~ is the orientational friction, and $\nu_{i}^{r}\mspace{180mu}\left( t \right)$ is responsible for introducing an orientational randomness. Akin to the stochastic force $\mathbf{\nu}_{i}\left( t \right)$ in [Eq (13)](#pcbi.1005569.e026){ref-type="disp-formula"}, $\nu_{i}^{r}\mspace{180mu}\left( t \right)$ has properties $\left\langle \nu_{i}^{r}\mspace{180mu}\left( t \right) \right\rangle = 0$ and $\left\langle \nu_{i}^{r}\mspace{180mu}\left( t \right)\nu_{j}^{r}\mspace{180mu}\left( t^{\prime} \right) \right\rangle = 2D_{r}\delta_{ij}\delta\mspace{180mu}\left( t - t^{\prime} \right)$, where *D*~*r*~ is an rotational diffusion constant. The related time scale *τ*~*r*~ = *γ*~*r*~/2*D*~*r*~ measures the time necessary for the cell polarisation direction to decorrelate due to fluctuations produced by cellular processes and the environment. The unit-length vector **n**~*i*~ then evolves according to $\frac{d\mathbf{n}_{i}}{dt} = \frac{d\theta_{i}}{dt}\mathbf{N}_{i} \times \mathbf{n}_{i}$ \[[@pcbi.1005569.ref080]\].

Eqs ([13](#pcbi.1005569.e026){ref-type="disp-formula"}) and ([14](#pcbi.1005569.e029){ref-type="disp-formula"}) describe the cell dynamics in the AVM. These equations are solved using the standard Euler-Maruyama method for the numerical solution of stochastic differential equations.

### Cell growth, division and death {#sec008}

Another contribution to the activity comes from cell growth, division and death or extrusion from the cell sheet, present in many epithelial tissues. We will refer to these three processes as *population control*. Cell growth is modelled as a constant-rate increase of the native area, i.e., $$A_{0}\mspace{180mu}\left( {t + \delta t} \right) = \left( {1 + \eta\delta t} \right)\mspace{180mu} A_{0}\left( t \right),$$ where *η* is the growth rate per unit time and *δt* is the simulation time step. We note that this is a very simple, generic model of cell growth that is clearly not applicable to all tissues. For example, in the case of the extensively studied wing imaginal disc in *Drosophila* an alternative cell-growth model has been proposed \[[@pcbi.1005569.ref081]\]. Due to the highly modular, object oriented design used in our implementation of the AVM, extending the model to include more sophisticated growth models would be straightforward (see [S3 Appendix](#pcbi.1005569.s003){ref-type="supplementary-material"}).

Cell death is modelled by tracking the cell's "age". Age is increased at every time step by *δt*. Once the age reaches a critical value, the cell is removed from the system. The cell is removed instantaneously. Upon removal of a cell, the Delaunay triangulation and its dual Voronoi diagram have to be recomputed, causing instantaneous changes of shapes of all neighbouring cells. While we did not explore in detail the effects of this instantaneous removal of cells, a set of simulations performed in different regions of parameter space suggests that there are no appreciable artefacts in the long-time and large-scale behaviour of the system, i.e., the "disturbance" caused by the instantaneous removal of a cell quickly relaxes. This is valid in the regime when the cell death rate is sufficiently low that the system locally relaxes between cell removal events. We note that onset of the cell's death (or alternatively its extrusion from the sheet) may also depend on the stresses or forces exerted on the cell. In the current implementation of the AVM we do not include such effects. Adding such effects to the model would be straightforward.

Finally, cell division is modelled based on the ideas of Bell and Anderson \[[@pcbi.1005569.ref082]\]. The cell can only divide once its area is larger then some critical area *A*~*c*~, in which case it divides with probability proportional to *A* − *A*~*c*~, i.e., $$p_{div} = \begin{cases}
{\chi\mspace{180mu}\left( {A - A_{c}} \right)} & {\,\,\text{for}\,\, A > A_{c}} \\
0 & {\,\,\text{otherwise}} \\
\end{cases},$$ where *χ* is a constant with units of inverse area times time that does not depend on a particular choice of the simulation time step. Upon division the two daughter cells are placed along the direction of the vector **n**~*i*~ and their ages are reset to zero. The Delaunay triangulation is rebuilt after every division or death event.

The cell cycle is clearly a far more complex phenomenon than what is captured by this simple model. The simplest possible extension would be to build upon the ideas of Smith and Martin \[[@pcbi.1005569.ref083]\], or use many other more sophisticated models available in the literature. Adding such extensions to the AVM would be straightforward and will be included in later versions of the model.

### Maintaining a Delaunay triangulation/Voronoi tessellation {#sec009}

In an actual simulation, we start either from carefully constructed initial positions of cell centres or choose cell positions from a particular experimental system. Those positions are then used to build the initial Delaunay triangulation and its corresponding dual Voronoi tessellation. However, as the cell centres move, there is no guarantee that the Delaunay character of the triangulation is preserved. For the model to be able to properly capture cell dynamics we, however, need to ensure that the triangulation is indeed Delaunay at each time step. Building a new triangulation in every time step would be computationally costly. Instead, we apply the so-called equiangulation procedure \[[@pcbi.1005569.ref084]\]: for every edge in the triangulation we compute the sum of the angles opposite to it. If the sum is larger than 180° we "flip" the edge (see [Fig 2c](#pcbi.1005569.g002){ref-type="fig"}). The procedure is repeated until there are no more edges left to flip. One can show \[[@pcbi.1005569.ref084]\] that this procedure always converges and leads to a Delaunay triangulation. While equiangulation is not the most efficient way to build a Delaunay triangulation "from scratch", if one starts with a triangulation that is nearly Delaunay, in practice only a handful of flip moves are required to recover the Delaunay triangulation. Given that cell centres move continuously, in the AVM the equiangulation procedure significantly increases the efficiency of maintaining the Delaunay triangulation throughout the simulation.

It is important to note that edge flips are local and flipping one edge in the Delaunay triangulation only affects one edge of the dual Voronoi tessellation. This means that a single edge flip can only affect junctions between four cells. The centres of those four cells coincide with the locations of the four corners of the polygon formed by two triangles sharing the flipped edge (see, [Fig 2c](#pcbi.1005569.g002){ref-type="fig"}). This is precisely the mechanism behind the T1 transition discussed in Sec. "T1 transitions" below. No other cell contacts are affected by the flip.

### Handling boundaries {#sec010}

The implementation of the SPV by Bi, *et al*. \[[@pcbi.1005569.ref059]\] assumes periodic boundary conditions. This assumption is reasonable if one studies a relatively small region of a much larger epithelial tissue. Many experiments, especially those of cell migration on elastic substrates, however, are performed with a relatively small number of cells where the effects of boundaries cannot be neglected or are even the main focus of the study. Therefore, in the AVM we include an open, flexible boundary.

Changes of the topology (connectivity) of the boundary would require complex, computationally costly geometric manipulations. At present, it is not even clear what appropriate set of update rules would be required in order to cover all possible cases that would involve changes of the boundary topology. Therefore, we assume that the connectivity of the boundary is maintained throughout the entire simulation. This does not mean that the boundary is fixed. It can grow, shrink and change its shape, but it cannot change its topology, e.g. it is not possible to transition from a disk to an annulus. Examples of allowed and disallowed changes of the boundary are shown in [Fig 3](#pcbi.1005569.g003){ref-type="fig"}. While fixing the boundary topology may appear restrictive, we will see below that in practice a model with fixed topology of the boundary allows for detailed studies of a broad range of problems directly applicable to many current experiments.

![Schematic representation of the allowed and disallowed changes of the boundary in the AVM.\
An initial configuration with the topology of a disk (a) is allowed to develop pronounced fingers (b). However, it is not possible to split into two domains (c) or develop a hole (d), which would both lead to the introduction of a new boundary lines, and therefore lead to changes in the topology.](pcbi.1005569.g003){#pcbi.1005569.g003}

In the AVM, the degrees of freedom are cell centres, represented as particles. These particles serve as sites of the Delaunay triangulation. In order to handle boundaries we introduce a special type of *boundary* particle and in addition, we also specify the connectivity between the boundary particles. Boundary particles together with their connectivity information form a *boundary line*. This line sets the topology of the boundary, which is preserved throughout the simulation, and delineates between the tissue and its surrounding. In addition, the boundary line can have an energy associated to it. Biologically, this boundary energy corresponds to the complex molecular machinery, such as actin cables, that is known to affect the behaviour of the free edge of an epithelial sheet. Specifically, we introduce the boundary line tension, $$E_{lt} = \frac{1}{2}\mspace{180mu}\sum\limits_{\langle{i,j}\rangle}\mspace{180mu}\lambda_{i,j}\mspace{180mu}\left( {l_{ij} - l_{0}} \right)^{2}$$ where *λ*~*i*,*j*~ is the line tension of the edge connecting boundary vertices *i* and *j*, *l*~*ij*~ = \|**r**~*i*~ − **r**~*j*~\| is the length of that edge and *l*~0~ is its native (preferred) length. In addition, we also introduce the boundary bending stiffness $$E_{bbend} = \frac{1}{2}\mspace{180mu}\sum\limits_{i}\mspace{180mu}\zeta_{i}\left( {\theta_{i} - \pi} \right)^{2},$$ where *ζ*~*i*~ is the stiffness of angle *θ*~*i*~ at the boundary particle *i*, determined as $$\theta_{i} = \arccos\frac{\mathbf{r}_{ji} \cdot \mathbf{r}_{ki}}{\left| \mathbf{r}_{ji} \right|\left| \mathbf{r}_{ki} \right|},$$ where **r**~*j*~ and **r**~*k*~ are the positions of boundary particles to the left and to the right of particle **r**~*i*~. For simplicity, we assume that each boundary particle has exactly two boundary neighbours. This prevents somewhat pathological, "cross-like" configurations where two otherwise disjoint domains would hinge on a single boundary site.

It is important to note that boundary particles do not represent centres of actual cells. The Voronoi polygon dual to a boundary particle in a Delaunay triangulation extends out to infinity. Consequently, it is not possible to unambiguously assign quantities such as the associated area or perimeter to boundary particles. Therefore, only internal particles correspond to cell centres, while boundary particles should be thought of as "ghosts" that serve to mark the edge of the cell sheet. These particles, however, experience forces from the interior of the tissue as well as from the interactions with their neighbours on the boundary. The boundary line is able to dynamically adjust its shape and length in this manner. We also allow for boundary particles to be added to or removed from the boundary line. Details of the algorithms used to dynamically update the boundary are given in [S2 Appendix](#pcbi.1005569.s002){ref-type="supplementary-material"}.

This concludes the description of the AVM. In [Fig 4](#pcbi.1005569.g004){ref-type="fig"} we present the flow chart of the main steps in computing the time evolution of cell centres and polarity vectors. Technical details of the implementation are given in [S3 Appendix](#pcbi.1005569.s003){ref-type="supplementary-material"} and in the on-line documentation provided with the *SAMoS* code.

![Flow chart of the force and torque calculations in each time step based on Eqs ([7](#pcbi.1005569.e014){ref-type="disp-formula"}) and ([12](#pcbi.1005569.e025){ref-type="disp-formula"}).\
The calculation speed is significantly improved by using the fast equiangulation moves shown in [Fig 2c](#pcbi.1005569.g002){ref-type="fig"}. The Delaunay triangulation is recomputed only if there are cell division and/or death events, which do not occur at every time step.](pcbi.1005569.g004){#pcbi.1005569.g004}

Results and discussion {#sec011}
======================

In order to validate the model and compare it with the results of similar models proposed in the literature, as well as to show its use in modelling actual biological tissues, we now apply the AVM to several problems relevant to the mechanics of epithelial tissue layers.

T1 transitions {#sec012}
--------------

We start by illustrating one of the key processes observed in epithelial tissues, the T1 transition. As detailed in [Fig 2](#pcbi.1005569.g002){ref-type="fig"}, in the AVM T1 transitions are handled through an edge flip in the Delaunay triangulation. An edge flip only happens when, in the notation of [Fig 2c](#pcbi.1005569.g002){ref-type="fig"}, we have *α* + *β* = *δ* + *γ* = 180°. Then both triangles are circumscribed by the same circle passing through its combined four vertices. The location of the T1 transition coincides with the centre of this circle. Due to the continuous connection between the position of sites of the Delaunay triangulation and its dual Voronoi tessellation, we always approach this point smoothly, i.e. a junction between two cells will smoothly shrink to a point, the T1 transition will occur, and then it will expand in a new direction. This process arises naturally in the AVM model, in stark contrast with many currently available implementations \[[@pcbi.1005569.ref085]\] that require a cut-off criterion on the edge length of a cell before a T1 transition can occur. It also avoids discontinuous jumps at finite edge length, bypassing the T1 point altogether, also a feature of a number of models, notably those based on sequential energy minimisation \[[@pcbi.1005569.ref049], [@pcbi.1005569.ref081], [@pcbi.1005569.ref086]\]. Next to its high computational efficiency, the ability to smoothly go through a T1 transition without the need for any additional manipulations of either the Delaunay triangulation or the Voronoi tessellation is one of the key advantages of the AVM approach.

In [Fig 5](#pcbi.1005569.g005){ref-type="fig"}, we illustrate a T1 transition in the bulk, in a region of phase space where the system exhibits liquid-like behaviour, but with very slow dynamics (see next section). The edge linking cells 2 and 4 (in red) slowly shrinks to a point, and then rapidly expands in the opposite direction. This feature points to dynamics akin to certain models of sheared materials \[[@pcbi.1005569.ref087]\], where the active driving pulls the material over an energy barrier from one minimum to the next. It is somewhat different from the activated dynamics which has been proposed for the SPV \[[@pcbi.1005569.ref086]\], which would predict a series of fluctuations through which the barrier between minima is ultimately crossed.

![Time lapse of a T1 transition, each label is the time (in units of *γ*/*Ka*^2^, with *a* = 1) with respect to the T1 event.\
Initially, cells 2 and 4 are in contact (red line). Approaching the transition, the connecting line slowly contracts, until it becomes a point at the transition. Cells 1 and 3 make a new contact which then rapidly expands. The arrows represent the force on the cell centre resulting from the Vertex potential (green), which partially compensates the active self-propulsion force (yellow) to give the resultant total force (red). The Voronoi tessellation is outlined in black, and the Delaunay triangulation is in white.](pcbi.1005569.g005){#pcbi.1005569.g005}

Activity driven fluidisation phase diagram {#sec013}
------------------------------------------

We now explore different modes of collective behaviour, i.e., *phases*, of the tissue based on the values of parameters of the original VM (*K*~*i*~, Γ~*i*~ and Λ~*μ*,*ν*~), and AVM-specific parameters such as the activity *f*~*a*~, the orientational correlation time *τ*~*r*~, and the boundary line tension *λ*. In order to keep the number of independent parameters to a minimum, it is again convenient to rewrite the energy of the VM, [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"}, in a scaled form \[[@pcbi.1005569.ref049], [@pcbi.1005569.ref059]\]. We first choose *K*~*i*~ = 1 and set $A_{i}^{0} = \pi$ as an area scale. For simplicity, we assume that all perimeter and junction tensions are the same, i.e., we set Γ~*i*~ ≡ Γ and Λ~*μ*,*ν*~ ≡ Λ for all *i*, *μ* and *ν*. Then, as discussed below [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"}, we can complete the square on the second and third terms in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} and obtain the scaled VM potential $$E_{VM} = \sum\limits_{i = 1}^{N}\mspace{180mu}\frac{1}{2}\mspace{180mu}\left( {A_{i} - A_{0}} \right)^{2} + \sum\limits_{i = 1}^{N}\mspace{180mu}\frac{\Gamma}{2}\mspace{180mu}\left( {P_{i} - P_{0}} \right)^{2},$$ where *P*~0~ = −Λ/Γ and *A*~0~ = *π*. The first term in [Eq (20)](#pcbi.1005569.e042){ref-type="disp-formula"} penalises changes in the cell area, while the second term penalises changes of the perimeter. There is no reason for the preferred area *A*~0~ to be generically compatible with the preferred perimeter *P*~0~. This sets up a competition between the two terms in [Eq (20)](#pcbi.1005569.e042){ref-type="disp-formula"}, giving a natural scale that is determined by the relative ratio of $\Gamma\mspace{180mu} P_{0}^{2}$ to $K\mspace{180mu} A_{0}^{2}$. In other words, if $K\mspace{180mu} A_{0}^{2} > \Gamma\mspace{180mu} P_{0}^{2}$, the cell will try to optimise its area at the expense of paying a penalty for not having the most optimal perimeter, and the opposite if $K\mspace{180mu} A_{0}^{2} < \Gamma\mspace{180mu} P_{0}^{2}$.

Bi, *et al*. \[[@pcbi.1005569.ref059]\] introduced the dimensionless *shape factor* $p_{0} = \frac{P_{0}}{\sqrt{A_{0}}}$, which controls the ratio of the cell's perimeter to its area through the target area *A*~0~ and target perimeter *P*~0~. The value of *p*~0~ then determines whether the area or the perimeter term in [Eq (20)](#pcbi.1005569.e042){ref-type="disp-formula"} wins and effectively sets the preferred shape of each cell: cells of different shapes have different values of *p*~0~. For example, regular hexagons, pentagons, squares and triangles correspond to *p*~0~ = 3.722, *p*~0~ = 3.812, *p*~0~ = 4.0 and *p*~0~ = 4.559, respectively.

Remarkably, one observes \[[@pcbi.1005569.ref049], [@pcbi.1005569.ref059], [@pcbi.1005569.ref086]\] a transition between a solid-like behaviour of the tissue, where cells do not exchange neighbours, and liquid-like behaviour, where neighbour exchanges do occur, at *p*~0~ = 3.812, a value that corresponds to a regular pentagon. At present, the biological significance of this observation is not clear, but it appears to be a robust feature of many experimental systems \[[@pcbi.1005569.ref088]\]. In order to make the comparison between the AVM and the SPV model, we also adopt *p*~0~ as a main parameter that controls the preferred cell shape.

In order to initialise the simulation, in each run we start by placing soft spheres with slightly polydisperse radii in a circular region. We then use *SAMoS* to minimise the energy of a soft sphere packing in the presence of a fixed boundary. This ensures that initially, cells are evenly spaced without being on a grid. We also fix the packing fraction to *ϕ* = 1, ensuring that the average cell area of the initial configuration is 〈*A*〉 = *A*~0~. The boundary is either fixed (referred to as "fixed system"), or allowed to fluctuate freely ("open system"). [Fig 6](#pcbi.1005569.g006){ref-type="fig"} shows a representative set of the states that we observe. We run the simulation for either 100,000 time steps with step size *δt* = 0.01 in the unstable region (e.g., [Fig 6e](#pcbi.1005569.g006){ref-type="fig"}), or 250,000 time steps with *δt* = 0.025 in the solid-like region (e.g., [Fig 6c](#pcbi.1005569.g006){ref-type="fig"}). For these systems with *N* = 1000 cells in the interior, this takes between 10--40 minutes on a single core of a modern Intel Xeon processor depending on the number of rebuilds of the Delaunay triangulation that are necessary (more in the liquid-like phase).

![Representative snapshots of states seen with fixed and flexible boundaries in the AVM.\
Contact lines between cells are outlined in white, and cells are coloured according to their area. The line connecting tissue boundary points is blue for flexible boundaries (panels a, d, e, f and g), and red for fixed ones (panels b and c). In panels e, f and g, cell centres are denoted by white spheres, and in panels e and f, the Delaunay triangulation is also shown as grey lines. (a) Shrinking cells at *p*~0~ = 2.48, Γ = 1.0, $\tau_{r}^{- 1} = 0.01$ and *f*~a~ = 0.1, no boundary line tension. (b) Same as (a), but for a fixed boundary. (c) Solid-like (glassy) state at *p*~0~ = 3.39, Γ = 1.0, $\tau_{r}^{- 1} = 0.1$ and *f*~a~ = 0.03. (d) Liquid-like state with a fluctuating boundary at *p*~0~ = 3.39, Γ = 1.0, $\tau_{r}^{- 1} = 0.1$ and *f*~a~ = 0.3, no boundary line tension (e) Fingering instability at *p*~0~ = 3.72, Γ = 0.1, $\tau_{r}^{- 1} = 0.01$ and *f*~a~ = 0.3, boundary line tension *λ* = 0.1. (f) Fluid state at *p*~0~ = 3.95, Γ = 0.1, $\tau_{r}^{- 1} = 0.1$ and *f*~a~ = 0.1, boundary line tension *λ* = 0.3. (g) Rosette formation at *p*~0~ = 4.85, Γ = 0.1, $\tau_{r}^{- 1} = 0.1$ and *f*~a~ = 0.03, boundary line tension *λ* = 0.3.](pcbi.1005569.g006){#pcbi.1005569.g006}

The unit of time is set by *γ*/*Ka*^2^, where *a* ≡ 1 is the unit of length. We note that Bi, *et al*. use $a = \sqrt{A_{0}}$ as the unit of length. This is possible as long as cells are not allowed to grow, i.e., when *A*~0~ does not change in time. The AVM allows for the cells to change their size and therefore we need to choose a different unit of length. In our case, *a* is the range of the of soft-core repulsion between cell centres (see [S1 Appendix](#pcbi.1005569.s001){ref-type="supplementary-material"}, Eq. (35)).

At low values of *p*~0~, we find a system that prefers to be in a state with mostly hexagonal cells, unless the active driving *f*~a~ is very high. Open systems will shrink at this point so that all cells are close to their target *P*~0~, as shown in [Fig 6a](#pcbi.1005569.g006){ref-type="fig"}. Consistent with this, larger values of the perimeter modulus Γ lead to stronger shrinking. For fixed systems, this route is blocked, and instead there is a strong inward tension on the boundaries and a gradient in local density, as shown in [Fig 6b](#pcbi.1005569.g006){ref-type="fig"}.

In agreement with the results of Bi, *et al*. \[[@pcbi.1005569.ref059]\], we find that at low *p*~0~ \< 3.81 and low values of driving *f*~a~, cells do not take an organised pattern and do not exchange neighbours. Recast in the language of solid state physics, the tissue is in an *amorphous solid* or *glassy* state. In [Fig 6c](#pcbi.1005569.g006){ref-type="fig"} we show such a state for a fixed boundary. In order to characterise the physical properties of this state, we measure the dynamical time scale of cell rearrangements through a standard tool of the physics of glassy systems, the self-intermediate scattering function \[[@pcbi.1005569.ref015]\] $$F\left( {q,t} \right) = \left\langle {\exp\left\lbrack {i\mathbf{q} \cdot \left( {\mathbf{r}\left( t \right) - \mathbf{r}\left( 0 \right)} \right)} \right\rbrack} \right\rangle.$$ *F* (*q*, *t*) measures the decay of the autocorrelation of cell-centre positions **r**(*t*) at a particular wave vector, **q**, taken usually to be the inverse cell size *q* ≡ \|**q**\| = 2*π*/*a*. The long-time decay of *F* (*q*, *t*) is characterised by the so-called *alpha-relaxation time* *τ*~*α*~ at which *F* (*q*, *t*) has decayed by half. When the system solidifies, i.e. when neighbour exchanges stop, **r**(*t*) remains constant and hence *τ*~*α*~ diverges \[[@pcbi.1005569.ref015]\], and stays infinite within the solid phase. In [Fig 7a--7c](#pcbi.1005569.g007){ref-type="fig"}, we show the phase diagram of *τ*~*α*~ as a function of *p*~0~ and *f*~a~, for several systems with different boundary conditions.

![Phase diagrams of the AVM.\
Panels (a), (b) and (c): *α*−relaxation time *τ*~*α*~ determined from the self-intermediate scattering function, [Eq (21)](#pcbi.1005569.e055){ref-type="disp-formula"}. These plots indicate that it is possible to initiate cell intercalation events by changing values of *p*~0~ or active driving. High values of *τ*~*α*~ (dark regions) correspond to the solid-like (glassy) phase where T1 events are suppressed. (a) Fixed system with Γ = 1 and $\tau_{r}^{- 1} = 0.01$. (b) Open system with Γ = 1, $\tau_{r}^{- 1} = 0.01$ and boundary line tension *λ* = 0. (c) Open system with Γ = 0.1, $\tau_{r}^{- 1} = 0.1$ and boundary line tension *λ* = 0.1. All systems are solid-like at low *p*~0~ and low driving *f*~a~. The critical *p*~0~ = 3.81 is the same in (a), (b) and (c), but the open system becomes liquid-like at much lower values of active driving. Lowering Γ also lowers this transition point. The dashed white line represents a rough boundary between solid-like and fluid-like behaviour. (d) Characteristic time scale *τ*~inst~ needed to reach the boundary instability, determined from reaching a threshold boundary length (see text), for the same parameters as (c). Sufficiently high driving always leads to an unstable system.](pcbi.1005569.g007){#pcbi.1005569.g007}

In [Fig 7a](#pcbi.1005569.g007){ref-type="fig"} we show regions of solid-like and liquid-like phases in a system with fixed boundaries, at Γ = 1, and a low noise value of *τ*~*r*~ = 0.01. We find a boundary of the solid-like phase that stretches from *p*~0~ ≈ 3.81 at small *f*~a~ to a maximum activity *f*~a~ beyond which the system is fluid at all *p*~0~. This is qualitatively, but not quantitatively consistent with the results of Bi, *et al*., who find a transition line at roughly twice our *f*~a~ values. Several factors are likely implicated in this discrepancy. Our systems, at *N* = 1000 cells, are more than twice as large as the *N* = 400 systems considered by Bi, *et al*., and finite system size effects seem to play an important role, as shown below. We measure *τ*~*α*~ at a value of 1/*q* corresponding to displacements of one cell size. However, even though displacements are large, we have evidence that this may not be sufficient to induce T1 transitions and therefore fluidise the system. Finally, fixed boundaries were used here and the periodic boundaries of Bi, *et al*. are likely not strictly equivalent.

The influence of the type of boundary conditions is very significant. In [Fig 7b](#pcbi.1005569.g007){ref-type="fig"}, we show the phase diagram for the same Γ = 1 and $\tau_{r}^{- 1} = 0.01$ as in [Fig 7a](#pcbi.1005569.g007){ref-type="fig"}, except with open boundary conditions and boundary line tension *λ* = 0. Separately, for Γ = 0.1, we have also confirmed that the value of the boundary line tension does not significantly affect the onset of the solid-like regime. We find a significantly lower maximum *f*~*a*~ for the transition, *f*~a~ = 0.03, a factor of 10 compared to the fixed case. We note that the effect also persists at $\tau_{r}^{- 1} = 0.1$, but is less pronounced. While we do not have a full explanation for this result, we do note that fluctuations of the boundary allow for rearrangements that are otherwise strongly suppressed by the fixed boundary. For example, the system in [Fig 6d](#pcbi.1005569.g006){ref-type="fig"} shows significant boundary fluctuations. It is liquid-like with *τ*~*α*~ ≈ 10, whereas the equivalent fixed system has *τ*~*α*~ ≈ 100. In view of the significant role of the boundary, we expect a strong system-size dependence.

At very high *p*~0~ and low active driving, we observe a systematic increase of *τ*~*α*~ (especially visible in [Fig 7c](#pcbi.1005569.g007){ref-type="fig"}). This unexpected result is accompanied by structural changes in the cell patterns that we observe. [Fig 6f](#pcbi.1005569.g006){ref-type="fig"} shows a liquid system at *p*~0~ = 3.95, near the relative minimum *τ*~*α*~ for a given *f*~a~. The distribution of cell centres appears random. In contrast, as can be seen in [Fig 6g](#pcbi.1005569.g006){ref-type="fig"}, at very high *p*~0~ = 4.85, cells arrange themselves into rosette shapes, where many vertices meet in a point. Rosettes are a feature of many developmental systems \[[@pcbi.1005569.ref089]\], so it is interesting to see that they do appear in the AVM context. Cell centres also arrange themselves in equidistant chains, hinting at a connection to one of the various pattern-formation instabilities studied in nonlinear dynamics. We note that this regime is numerically delicate, and the addition of the soft repulsive core between cell centres (see [S1 Appendix](#pcbi.1005569.s001){ref-type="supplementary-material"}, Eq. (35)) is necessary to make simulations stable. This simulation is performed deep inside the liquid regime where cell junctions are not necessarily accurately represented by straight lines and the applicability of the model is questionable. At present it is, therefore, not clear whether the configuration shown in [Fig 6g](#pcbi.1005569.g006){ref-type="fig"} is an artefact of the model or it indeed has biological significance. One should proceed with caution when assigning biological interpretation to any configuration with very large values of *p*~0~, i.e. for p~0~ ≳ 4.6.

In parts of the phase diagram, we observe a fingering instability \[[@pcbi.1005569.ref090]--[@pcbi.1005569.ref093]\] where regions a few cells wide migrate outward from the centre, as shown in [Fig 6e](#pcbi.1005569.g006){ref-type="fig"}. When *τ*~*α*~ drops below approximately 10, we observe that the fluctuations of the boundary already present in [Fig 6d](#pcbi.1005569.g006){ref-type="fig"} become unconstrained. This is a mechanically unstable regime: Eventually, these cells will detach, a process we are not yet able to model due to the topological change that it would imply (see, [Fig 3](#pcbi.1005569.g003){ref-type="fig"}). We have observed that fluctuations need to reach a threshold of approximately \>5% of a length increase in the boundary to break through to an unconstrained growth, otherwise the system remains stable, see e.g. [Fig 6d--6g](#pcbi.1005569.g006){ref-type="fig"}. We then associate a time scale *τ*~inst~ with reaching this threshold and use it to measure the degree of instability: a small time scale denotes a rapid growth rate of fluctuations. [Fig 7c](#pcbi.1005569.g007){ref-type="fig"} shows *τ*~*α*~, and [Fig 7d](#pcbi.1005569.g007){ref-type="fig"} shows *τ*~inst~ for the same open system with Γ = 0.1, $\tau_{r}^{- 1} = 0.1$ and boundary line tension *λ* = 0.1. We note that the transition line between solid-like and fluid-like states is low, at *f*~a~ = 0.03. At and below *f*~a~ = 0.1, the boundary of the system is stable, and above this threshold, the instability becomes more pronounced with increasing *f*~a~ and smaller *τ*~*α*~. The physical mechanism responsible for the instability involves a subtle interplay of *f*~a~, boundary line tension (stronger line tension suppresses the instability), the noise level $\tau_{r}^{- 1}$ (lower $\tau_{r}^{- 1}$ enhances the instability), and *p*~0~. The instability resembles observations of finger formation in MDCK monolayers \[[@pcbi.1005569.ref090], [@pcbi.1005569.ref091]\]. Existing models link it to either leader cells \[[@pcbi.1005569.ref093], [@pcbi.1005569.ref094]\], a bending instability \[[@pcbi.1005569.ref092]\], or an active growth feedback loop \[[@pcbi.1005569.ref095]\], while here it emerges naturally. A detailed account of this phenomenon will be published elsewhere.

Growth and division {#sec014}
-------------------

Division and death processes are important in any living tissue, for example, cell division and ingression processes play essential roles during development. Therefore, as noted in Sec. "Cell Growth, Division and Death", the AVM is equipped to handle such processes. It is important to note, however, that the removal of one cell during apoptosis or ingression and the addition of two new cells during division in the AVM causes a discrete change in the Voronoi tessellation which implies a discontinuous change of the local forces derived from the VM. We have simulated the growth of a small cluster of cells to assess whether this discrete change in geometry can lead to any instabilities in the model. These test runs did not reveal any artefacts due to discontinuities in the force caused by the division events.

In order to illustrate the growth process, we choose a shape factor, *p*~0~ = 3.10, corresponding to Γ = 1 and Λ = −5.5 and no active driving, i.e., we set *f*~a~ = 0. This puts the system into the solid-like phase where T1 events are absent. Our simulation runs for 10^6^ time steps at *δt* = 0.005, corresponding to 5000 time units, starting from 37 cells and stopping at about 24,000 cells. To balance computational efficiency with a smooth rate of division, cells are checked for division every 25 time steps. We show snapshots of different stages of the tissue growth in [Fig 8a--8e](#pcbi.1005569.g008){ref-type="fig"}.

![Snapshots of a growing epithelial tissue.\
Frames (a), (b), (c), (d) and (e) have 37, 63, 124, 4633 and 23787 cells and are at times 50, 500, 1000, 3500 and 5000, respectively. Cells have a chance to divide if their area is greater than a critical area *A*~*c*~ = 2.8*a*^2^ after which the probability of a cell to divide increases linearly with its area (see [Eq (16)](#pcbi.1005569.e037){ref-type="disp-formula"}). In this simulation, the shape factor was set to *p*~0~ = 3.10. (f) Log-linear plot of the total number of cells as a function of simulation time. The growth rate of the patch is initially exponential but starts to slow at around 3000 cells. This is due to cells in the centre of the cluster being prevented from expanding by the surrounding tissue. (g) Tissue after 5000 time units (10^6^ time steps) with each cell coloured by pressure. Pressure has built up in the centre of the tissue while close to the edge the average pressure is low. (h) Average pressure (averaged over the polar angle) as a function of the radial distance from the centre of the tissue. (i) Distribution of the number of neighbours for cells in the system shown in (g).](pcbi.1005569.g008){#pcbi.1005569.g008}

We note that the numerical stability of the simulation that involves growth is quite sensitive to the values of the parameters used in the AVM. For example, divisions of highly irregularly shaped cells, as commonly observed in the high *p*~0~ regime, can put a significant strain on the simulation and even cause a crash. Helpfully, some of these problems are alleviated by the soft repulsive potential between cell centres (see [S1 Appendix](#pcbi.1005569.s001){ref-type="supplementary-material"}, Eq. (35)). This in addition to a smaller time step is used to mediate the impact of cell divisions for growing systems with high *p*~0~.

In [Fig 8f](#pcbi.1005569.g008){ref-type="fig"} we show the tissue size as a function of the simulation time. In this simulation there are no apoptosis or cell ingression events and, as expected, the tissue size grows exponentially. However, at long times, the growth slows down and deviates from exponential growth. This is easy to understand, as the centre of the tissue is prevented from expanding by the surrounding cells. The effect can be seen in [Fig 8e](#pcbi.1005569.g008){ref-type="fig"}, where cells located towards the centre have, on average, smaller areas and in [Fig 8g](#pcbi.1005569.g008){ref-type="fig"}, which shows a clear pressure buildup in the centre. This suggests that in the later stages, the simulated tissue is not in mechanical equilibrium any more. The pressure is computed as the trace of the Hardy stress tensor, defined as \[[@pcbi.1005569.ref096]\] $${\hat{\sigma}}_{w}^{H} = \frac{1}{2}\mspace{180mu}\sum\limits_{\binom{\alpha,\beta}{\alpha \neq \beta}}\mspace{180mu}\int_{s = 0}^{1}\mspace{180mu}\left\lbrack {- \mathbf{f}_{\alpha\beta}\omega\mspace{180mu}\left( {\left( 1 - s \right)\mathbf{r}_{\alpha} + s\mathbf{r}_{\beta} - \mathbf{r}} \right) \otimes \left( {\mathbf{r}_{\alpha} - \mathbf{r}_{\beta}} \right)} \right\rbrack\mspace{180mu} ds,$$ where *ω* is a smoothing function and the sum is over all junctions and triangulation edges. The reason is that we can determine the shape and area of each cell using just the scalar lengths, {*r*~*αβ*~}, of these two sets of edges. We eliminate the integral by choosing, $$\omega\left( \mathbf{r} \right) = \begin{cases}
{1/A_{\mathcal{D}}} & {\text{if}\mathbf{r} \in \mathcal{D}} \\
0 & \text{otherwise} \\
\end{cases}$$ where $\mathcal{D}$ is the region for which the stress is calculated and $A_{\mathcal{D}}$ is the area of $\mathcal{D}$. For convenience we choose $\mathcal{D}$ to coincide with individual cells. Finally, $\mathbf{f}_{\alpha\beta} = - \frac{\partial E_{VM}}{\partial r_{\alpha\beta}}{\hat{\mathbf{r}}}_{\alpha\beta}$, where *E*~*VM*~ is defined in [Eq (1)](#pcbi.1005569.e002){ref-type="disp-formula"} and ${\hat{\mathbf{r}}}_{\alpha\beta}$ is the unit-length vector along the junction or triangulation edge denoted by (*α*, *β*). This particular choice of *ω* and $\mathcal{D}$ reduces our calculation to that of the more common virial stress \[[@pcbi.1005569.ref096]\], but we note that the Hardy stress is a useful tool for coarse-graining. A full discussion regarding coarse-graining of stress calculations in the AVM will be published elsewhere. We also see clear heterogeneities in the local pressure shown in [Fig 8g](#pcbi.1005569.g008){ref-type="fig"}. In [Fig 8h](#pcbi.1005569.g008){ref-type="fig"}, we show the radial pressure profile in the tissue at the end of the simulation. From the figure it is also evident that the there is a substantial pressure buildup close to the centre of the tissue as well as that angular averaging substantially reduces local pressure fluctuation notable in [Fig 8g](#pcbi.1005569.g008){ref-type="fig"}. The origin of these effects warrants a detailed investigation and will be addressed in a later publication, we note however that stress inhomogeneities are a persistent feature of the epithelial cell monolayers that have been investigated by traction force microscopy \[[@pcbi.1005569.ref012], [@pcbi.1005569.ref097]\]. A similar pressure buildup has also been investigated in a growing cell aggregate \[[@pcbi.1005569.ref098]\] and in a mathematical model of nonuniform growth in a layer of tissue \[[@pcbi.1005569.ref099]\].

Finally, in [Fig 8i](#pcbi.1005569.g008){ref-type="fig"} we show the distribution of the number of neighbours for this model system. The observed distribution is in qualitative agreement with observations on several tissues both in animals and plants (e.g., *Drosophila* wing disc, *Xenopus* tail epidermis, *Hydra* outer epidermis, *Anagallis arvensis* meristem, cucumber epidermis) reported in the literature \[[@pcbi.1005569.ref100], [@pcbi.1005569.ref101]\]. We note that while the relative values of the percentages of nearest neighbours are dependent on the parameters used in the model, we observed that the overall shape of the distribution is maintained over a several sets of parameters chosen in different regions of the parameter space. A detailed study of this distribution and its precise dependence on the model parameters is, however, beyond the scope of this study.

Modelling mechanically heterogeneous tissues {#sec015}
--------------------------------------------

The AVM is equipped to allow for cell-specific parameters, which enables us to investigate tissues with locally varying mechanical properties. A commonly studied example of the effects such heterogeneities is cell sorting. As an example we show simulations that display sorting of two distinct cell types. We achieve this by setting the junction tension Λ for each pair of cell-cell and cell-boundary contacts. All our simulations consist of 1000 cells with half chosen randomly to be of the "red" type and the others being of the "blue" type. In these simulations, boundaries have been kept fixed. We observe sorting behaviour akin to that found in other commonly used tissue models \[[@pcbi.1005569.ref035], [@pcbi.1005569.ref102]\]. Using *r*, *b* and *M* to denote red, blue and the boundary respectively, we start by fixing *K* = 1 and Γ = 1 and set −6.8 = Λ~*rr*~ \< Λ~*rb*~ \< Λ~*bb*~ = −6.2, corresponding to *p*~0~ in the range 3.58--3.93. We, however, note that *p*~0~ is a quantity defined "per cell" and one should understand it in this context only as a rough estimate whether a given cell type is in the solid or fluid phase. All cells are subject to small random fluctuations of their position which allows for T1 transitions that can bring initially distant cells into contact. We set the active driving, *f*~a~, to zero. We have chosen different values for Λ~*rr*~ and Λ~*bb*~ to reflect the idea that the surfaces of these cells have different adhesive properties \[[@pcbi.1005569.ref103]\]. Note that the Λ parameter for a particular contact is proportional to its energy per unit length. Sorting of cells into groups of the same type occurs when the energy of two red-blue contacts is greater than the energy of one red-red contact and one blue-blue contact, corresponding to Λ~*rb*~ \> (Λ~*rr*~ + Λ~*bb*~)/2, see [Fig 9a--9c](#pcbi.1005569.g009){ref-type="fig"}. In this regime, for cells of the same type it is energetically favourable for the new contact to elongate while local red-blue contacts are shortened. Conversely, if Λ~*rb*~ \< (Λ~*rr*~ + Λ~*bb*~)/2, then cells maximise their red-blue contacts forming a "checkerboard" pattern ([Fig 9d](#pcbi.1005569.g009){ref-type="fig"}). The final pattern is not without defects, the number and location of which depend on the initial conditions. The tissue boundary consists of contacts between cell centres and boundary particles so Λ~*rM*~ and Λ~*bM*~ need also to be specified to reflect the way in which the cell types interact with the extracellular matrix or surrounding medium. Initially we set Λ~*rM*~ = Λ~*bM*~ = −6.2 and observe that blue cells cover the boundary enveloping red cells because this facilitates lower energy red-blue and red-red contacts being formed. If we incentivise red-boundary contacts by setting Λ~*rM*~ \< Λ~*rr*~ + Λ~*bM*~ − Λ~*rb*~ = −6.6 we make red-boundary contacts preferable \[[@pcbi.1005569.ref102]\]. This case is shown in [Fig 9e](#pcbi.1005569.g009){ref-type="fig"} for Λ~*rM*~ = −6.8.

![(a-c) Snapshots of a system of two cell types at times 10, 500 and 5000 with Λ~*rb*~ = −6.4. (d) A "checkerboard" pattern formes immediately (at time 10) when red-blue cell-cell contacts are energetically favourable compared with pairs of red-red and blue-blue contacts, Λ~*rb*~ = −6.7. (e) Same initial system as (a-c) but with red-boundary contacts slightly favoured over blue-boundary contacts. The system gradually separates into compartments of each cell type. The uncorrelated random fluctuations are sufficient to drive neighbour exchanges within the bulk of both the red and blue cell compartments. Cells on the compartment boundary can sometimes move parallel to it but meet strong resistance when trying to move across it.](pcbi.1005569.g009){#pcbi.1005569.g009}

Finally, we note that all simulations for mechanically heterogeneous systems were performed without any active driving. While a detailed study of the effects of active driving on cell sorting is beyond the scope of this study, we note that a limited set of test runs suggest, as expected, that including activity can reduce the time over which cells sort but that very high activity leads to mixing.

Effects of cellular alignment {#sec016}
-----------------------------

We now briefly turn our attention to the effects of several models of cell polarity alignment. So far, we have assumed torque ***τ***~*i*~ = 0 in [Eq (14)](#pcbi.1005569.e029){ref-type="disp-formula"}, i.e., we are in the situation where the the polarisation vector of each cell is independent of the surrounding cells and its direction diffuses randomly over time. In biological systems, it is known that a cell's polarity responds to the surrounding and many forms of polarity alignment have been proposed. Here, we highlight two alignment mechanisms that are compatible with the current understanding of cell mechanics. In [Fig 10a--10c](#pcbi.1005569.g010){ref-type="fig"}, we have used the alignment model defined in [Eq (9)](#pcbi.1005569.e021){ref-type="disp-formula"} that assumes that the polarity vector **n**~*i*~ of cell *i* aligns with its velocity **v**~*i*~. The torque term in [Eq (14)](#pcbi.1005569.e029){ref-type="disp-formula"} is then given by ***τ***~*i*~ = −*J*~*v*~ **v**~*i*~ × **n**~*i*~, where *J*~*v*~ is the alignment strength. This model was first developed for collectively migrating cells (modelled as particles) \[[@pcbi.1005569.ref074]\], and it exhibits global polar migration, i.e. a state in which all particles align their velocities and travel as a flock. In dense systems of active particles confined to a finite region, velocity alignment has been shown to be intimately linked to collective elastic oscillations \[[@pcbi.1005569.ref075]\]. It is remarkable that the main hallmarks of this active matter dynamics are also observed in the model tissue. In [Fig 10a](#pcbi.1005569.g010){ref-type="fig"}, we show velocity alignment dynamics for a confined system in the solid-like phase; here the collective oscillations are very apparent. They are strikingly reminiscent of the collective displacement modes observed in confined MDCK cell layers \[[@pcbi.1005569.ref104], [@pcbi.1005569.ref105]\].

![Effects of alignment.\
(a-c) are self-alignment of the cell polarisation with the cell velocity. (d) is self-alignment of the cell polarisation with the long axis of the cell. (a) A confined system in the solid-like region (*f*~a~ = 0.03, *p*~0~ = 3.385, $\tau_{r}^{- 1} = 0.01$), at alignment strength *J*~*v*~ = 1.0 shows oscillating collective modes. (b) In the liquid region (*f*~a~ = 0.1, *p*~0~ = 4.4, $\tau_{r}^{- 1} = 0.1$), the system enters a vortex rotation state instead (*J*~*v*~ = 1.0). (c) An open system at *f*~a~ = 0.1, *p*~0~ = 4.4, $\tau_{r}^{- 1} = 0.1$, boundary line tension *λ* = 0.1 and *J*~*v*~ = 0.1 migrates collectively. (d) An open system at *f*~a~ = 0.1, *p*~0~ = 3.95, $\tau_{r}^{- 1} = 0.1$, boundary line tension *λ* = 0.1 and alignment with cell shape with *J*~*s*~ = 1 migrates collectively with complex fluctuations. Arrows are cell velocity vectors **v**~*i*~; they are coloured by the magnitude of *v*~*x*~ for (a) and *v*~*y*~ for (b-d). Orange is positive (pointing right/upwards) and blue is negative (pointing left/downward).](pcbi.1005569.g010){#pcbi.1005569.g010}

In [Fig 10b](#pcbi.1005569.g010){ref-type="fig"}, we apply the same dynamics, but now to a system that is in the liquid-like phase, with fixed boundaries. Here, the collective migration wins, but the confinement to a disk with fixed boundaries forces the cells into a vortex-shaped migration pattern. Finally, in [Fig 10c](#pcbi.1005569.g010){ref-type="fig"}, in the absence of confinement, we recover the collective polar migration of the cell patch.

In [Fig 10d](#pcbi.1005569.g010){ref-type="fig"}, we show the effects of aligning the cell's polarity to the largest principal axis of the cell shape tensor, defined in [Eq (10)](#pcbi.1005569.e023){ref-type="disp-formula"}. This type of alignment also leads to collective motion in an unconfined system, however there are significant fluctuations as the allowed cell patterns are highly frustrated by the constraint to remain in a Voronoi tesselation.

These preliminary results serve as a showcase of the non-trivial effects of cell-cell alignment on the collective behaviour of the entire tissue. A more detailed account of the effects of different alignment models will be published elsewhere.

Modelling non-circular tissue shapes {#sec017}
------------------------------------

In the previous discussions, all examples assumed a patch of cells with the topology of a disk. However, the AVM is not restricted to the circular geometry and can be applied to systems of arbitrary shapes, including domains with complex connectivity. Such situations often arise when modelling experimental systems where cells surround an obstacle, or in studies of wound healing. In [Fig 11](#pcbi.1005569.g011){ref-type="fig"} we present a gallery of non-circular shapes that can be readily studied using the AVM as it is implemented in *SAMoS*. The annular geometry shown in [Fig 11a](#pcbi.1005569.g011){ref-type="fig"} would be suited for modelling wound healing problems as well as situations where cells migrate in order to fill a void. A common experimental setup where cell colonies are prepared as rectangular strips \[[@pcbi.1005569.ref090]\] is shown in [Fig 11b](#pcbi.1005569.g011){ref-type="fig"}, where three separate patches grow towards each other. Finally, in [Fig 11c](#pcbi.1005569.g011){ref-type="fig"} we show an example of yet another very interesting situation \[[@pcbi.1005569.ref106]\], where cells are grown in a confined region of space.

![(a) Snapshots of the simulation of a cell sheet with an annular geometry used to illustrate how cells divide and fill the circular void in the centre. The system configuration was recorded at times 0, 1000 and 1700. The initial configuration has *p*~0~ = 3.46, i.e., it was in the solid-like phase. While it is not simple to define *p*~0~ for a growing system due to a constant change in the cell target area *A*~0~, we note that throughout the simulation, cell shapes remain regular. (b) Illustration of a common experimental setup where cells are grown in rectangular "moulds" for a system with initial *p*~0~ = 3.35. Once the entire region is filled with cells, the mould is removed and the colony is allowed to freely grow. Images on the right show two of the strips about to merge. In (c), we model tissue growth in confinement using a system with initial *p*~0~ = 3.46. Grey beads form the boundary of the confinement region that constrain the cell growth. Initially, cells do not touch the wall and freely grow. As the colony reaches the wall, one starts to notice pressure buildup. Finally, the entire cavity is filled with cells and any subsequent division leads to increase of the pressure in the system. Snapshots were recorded at times 1200, 1480 and 1600.](pcbi.1005569.g011){#pcbi.1005569.g011}

Summary and conclusions {#sec018}
-----------------------

In this paper we have introduced the Active Vertex Model. It is a hybrid model that combines ideas from the physics of active matter with the Vertex Model, a widely used description for modelling confluent epithelial tissues. Active matter physics is a rapidly growing field of research in soft condensed matter physics, and it is emerging as a natural framework for describing many biophysical processes, in particular those that occur at mesoscales, i.e., at the scales that span multiple cells to the entire tissue. Our approach is complementary to the recently introduced Self-Propelled Voronoi model \[[@pcbi.1005569.ref059]\], for it allows modelling of systems with fixed and open, i.e. dynamically changing boundaries as well as cell-cell alignment, cell growth, division and death. The AVM has been implemented into the *SAMoS* software package and is publicly available under a non-restrictive open source license.

The AVM utilises a mathematical duality between Delaunay and Voronoi tessellation in order to relate forces on cell centres to the positions of the vertices of the dual lattice, i.e. meeting points of three of more cells---a natural description of a confluent epithelial tissue. This not only allows for a straightforward and efficient implementation using standard algorithms for particle-based simulations, but provides a natural framework for modelling topological changes in the tissue, such as intercalation and ingression events. In other words, in the AVM T1 transition events arise spontaneously and it is not necessary to perform any additional steps in order to ensure that cells are able to exchange neighbours.

Furthermore, our implementation of the AVM is very efficient, allowing for simulations of systems containing tens of thousands of cells on a single CPU core, thus enabling one to probe collective features, such as global cell flow patterns that span length-scales of several millimetres. In addition, the AVM is also able to handle multiple cell types and type specific cell contacts, which allows simulations of mechanically heterogeneous systems.

All these features make the AVM a strong candidate model to address many interesting biological and biophysical problems related to the mechanical response of epithelial tissues, especially those that occur at large length and time scales that are typically only accessible to continuum models. Unlike in the case of continuum models where relating parameters of the model to the experimental systems is often difficult and unclear, the AVM retains the cell-level resolution, making it simpler to connect it to the processes that occur at scales of single cells. The AVM is, however, not designed to replace continuum models, but to serve as the important layer that connects the complex molecular processes that occur at the cell level with the global collective behaviour observed at the level of the entire tissue.

A natural question to ask is: What is the advantage of the AVM compared to other particle based models, e.g. \[[@pcbi.1005569.ref043], [@pcbi.1005569.ref044]\], that are computationally far less demanding? There is clearly a tradeoff between computational costs and the required level of details necessary to address a specific biological question. It is, however, not a priori clear where the right balance between the two lies. The key advantage of the AVM compared to other cell-centre based models is that keeping track of both cell centres and cells themselves allows for straightforward extensions of the model to include effects such as active forces on vertices, active non-linear response of cell junctions, etc. that would be impossible to implement and parametrise in purely cell-centre based models.

With that in mind, there are, of course, still many ways the model can be improved. For examples, it would be very interesting to augment the AVM to include the effects of biochemical signalling. This would require solving a set of differential equations for signals in each time step, and then supplying those solutions to the mechanics part of the model. Adding such functionality would substantially increase the computational cost of simulations, however at the same time it would allow for detailed studies of the coupling between chemical and mechanical signalling. These are believed to be essential for developing a full understanding of the mechanics of epithelial tissues. Given the layout of the AVM and its implementation, implementing such functionality would be straightforward.

Furthermore, in the current version of the AVM activity is introduced in a very rudimentary manner, via assuming that cells self-propel in the direction of their polarity vector. This is a very strong assumption that would need a much stronger experimental support than currently available. It is also possible that a far more sophisticated model would be required to fully capture cell motility. However, one also needs to keep in mind that there is a tradeoff between being as biologically accurate as possible and retaining a sufficient level of simplicity to be able to efficiently perform simulations of large systems. With all this in mind, we argue that the our simulations clearly show that even this simplistic model of self-propulsion is capable of capturing many features of real systems, and that it can serve as a good starting point for building biologically more accurate descriptions.

Another potentially very interesting feature that is currently not supported would be splitting and merging of the boundaries, that is, allowing for topological changes of the entire sheet such as those depicted in [Fig 3c](#pcbi.1005569.g003){ref-type="fig"}. This would allow us to study detachment of a part of the tissue or opening apertures as well as the opposite problem of closing holes and gaps. The latter is of great importance for studying problems related to wound healing. Unfortunately, setting up a set of general rules on how to automatically split a boundary line or merge two boundaries into a single one is not a simple task from a point of view of computational geometry. The problem is further complicated if those rules are also to be made biologically plausible, which is essential for the model to be relevant to actual experiments.

We conclude by noting that recently there have been several interesting attempts at extending the Vertex Model description to three dimensions \[[@pcbi.1005569.ref085], [@pcbi.1005569.ref107]--[@pcbi.1005569.ref112]\]. While in principle possible, extending the AVM to a curved surface or making it fully three-dimensional would be quite involved. Being able to study curved epithelial sheets would be of great interest to systems where curvature clearly cannot be ignored, e.g., as is the case in modelling intestinal crypts \[[@pcbi.1005569.ref113]--[@pcbi.1005569.ref116]\]. While there is nothing in the description of the AVM that is unique to the planar geometry, there are several technical challenges associated with directly porting it onto a curved surface. Most notably, building a Delaunay triangulation on an arbitrary curved surface is not a simple task. In addition, quantities such as bending rigidity that are naturally defined on triangles would have to be properly mapped onto contacts lines between neighbouring cells. This is not straightforward to do. Developing a fully three-dimensional version of the AVM would be an even a greater challenge since the duality between Delaunay and Voronoi descriptions central to this model has no analogue in three dimensions.

We hope that the AVM will provide a useful and complementary tool for probing the aspects of the epithelial tissue mechanics that are not available to other methods, as well as serve as an independent validation for the results obtained by other methods.

Supporting information {#sec019}
======================

###### Detailed derivation of [Eq (7)](#pcbi.1005569.e014){ref-type="disp-formula"}, the AVM force on the cell centres.

(PDF)

###### 

Click here for additional data file.

###### Details of the algorithm for handling boundaries of tissues and cell monolayers.

(PDF)

###### 

Click here for additional data file.

###### Overview of *SAMoS* organisation and AVM implementation.

(PDF)

###### 

Click here for additional data file.

###### Shrinking cells with shape index *p*~0~ = 2.48.

Parameters: Γ = 1.0, $\tau_{r}^{- 1} = 0.01$ and *f*~a~ = 0.1, no boundary line tension. Video associated with [Fig 6a](#pcbi.1005569.g006){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### Same as S1 at *p*~0~ = 2.48 but with a fixed boundary, preventing cells from shrinking.

The red arrows represent the tensile force on the boundary. Parameters: Γ = 1.0, $\tau_{r}^{- 1} = 0.01$ and *f*~a~ = 0.1. Video associated with [Fig 6b](#pcbi.1005569.g006){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### Solid-like (glassy) state at *p*~0~ = 3.39, Γ = 1.0, $\tau_{r}^{- 1} = 0.1$, *f*~a~ = 0.03.

Boundary is fixed. Video associated with [Fig 6c](#pcbi.1005569.g006){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### Liquid-like state with a fluctuating boundary at *p*~0~ = 3.39, Γ = 1.0, $\tau_{r}^{- 1} = 0.1$ and *f*~a~ = 0.3.

The boundary line tension *λ* is zero here. Video associated with [Fig 6d](#pcbi.1005569.g006){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### Fingering instability at *p*~0~ = 3.72.

This instability occurs when correlated fluctuations in the bulk at high driving overcome small line tension. Parameters: Γ = 0.1, $\tau_{r}^{- 1} = 0.01$ and *f*~a~ = 0.3, boundary line tension *λ* = 0.1. Video associated with [Fig 6e](#pcbi.1005569.g006){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### Fluid state at *p*~0~ = 3.95, Γ = 0.1, $\tau_{r}^{- 1} = 0.1$ and *f*~a~ = 0.1, boundary line tension *λ* = 0.3.

The distribution of cell Voronoi centres, marked by white spheres, appears to be random. Video associated with [Fig 6f](#pcbi.1005569.g006){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### Rosette formation at (very high) *p*~0~ = 4.85.

Cell centres form chain and ring structures, and multiple vertices converge in a single point, leading to rosette-like structures. Parameters: Γ = 0.1, $\tau_{r}^{- 1} = 0.1$ and *f*~a~ = 0.03, boundary line tension *λ* = 0.3. Cell Voronoi centres are marked by white spheres. Video associated with [Fig 6g](#pcbi.1005569.g006){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### Simulation of cells which are growing and dividing (starting with about 1000 cells here).

The likelihood that a cell will divide increases in proportion to its size. Cells have an initial *p*~0~ = 3.1 and the few intercalation events between cells are primarily caused by cell divisions. Parameters: Γ = 1.0 and *f*~a~ = 0. Video associated with [Fig 8](#pcbi.1005569.g008){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### Two distinct cell types, intially randomly dispersed, attempt to minimise their contact with opposing-type cells by sorting into compartments.

The sorting behaviour is assisted by small uncorrelated random forces applied to all cells which acts to prevent the system from settling into a local minimum. In other words it can bring separate patches of same-type cells into contact. Video associated with [Fig 9a--9c](#pcbi.1005569.g009){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### A confined system in the solid-like region showing oscillating collective modes.

Cells align their active force *f*~a~ with their velocity at a rate defined by the alignment strength *J*~*v*~ = 1.0. Parameters: *f*~a~ = 0.03, *p*~0~ = 3.385, $\tau_{r}^{- 1} = 0.01$. Video associated with [Fig 10a](#pcbi.1005569.g010){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### For high *p*~0~ (fluid phase) the system enters a vortex rotation state when placed in confinement.

Parameters: *f*~a~ = 0.1, *p*~0~ = 4.4, $\tau_{r}^{- 1} = 0.1$, alignment strength *J*~*v*~ = 1.0. Video associated with [Fig 10b](#pcbi.1005569.g010){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### System with a free boundary that migrates collectively.

Parameters: *f*~a~ = 0.1, *p*~0~ = 4.4, $\tau_{r}^{- 1} = 0.1$, boundary line tension *λ* = 0.1, alignment strength *J*~*v*~ = 1.0. Video associated with [Fig 10c](#pcbi.1005569.g010){ref-type="fig"}.

(MP4)

###### 

Click here for additional data file.

###### System with shape alignment in free boundary that migrates collectively but with complex fluctuations in the bulk.

Here, each cell aligns its active force *f*~a~ with the principal axis of its geometric shape. Parameters: *f*~a~ = 0.1, *p*~0~ = 3.95, $\tau_{r}^{- 1} = 0.1$, boundary line tension *λ* = 0.1 and the rate of the alignment with cell shape is *J*~*s*~ = 1. Video associated with [Fig 10d](#pcbi.1005569.g010){ref-type="fig"}.

(AVI)

###### 

Click here for additional data file.

###### A growing tissue patch with a hole cut from the centre to form an annular geometry which mimics those used in wound healing studies.

Intial *p*~0~ = 3.46.

(MP4)

###### 

Click here for additional data file.

###### Example of a common experimental setup where cells are grown in separated rectangular regions and then allowed to grow freely towards eachother.

Inital *p*~0~ = 3.35.

(MP4)

###### 

Click here for additional data file.

###### Cells growing and dividing in a confined space.

In this case the wall is fixed in place and pressure in the system builds up quickly after the cell colony fills the available space. Initial *p*~0~ = 3.46.

(MP4)

###### 

Click here for additional data file.

###### Compressed folder containing input data and configuration files for reproducing all figures and supplementary movies in this paper.

(ZIP)

###### 

Click here for additional data file.

We would like to acknowledge many valuable discussions with Geoff Barton, Dapeng (Max) Bi, Luke Coburn, Amit Das, Tamal Das, Alexander Fletcher, François Graner, M. Lisa Manning, M. Cristina Marchetti, Kirsten Martens, Inke Näthke, and K. R. Prathyusha. SH acknowledges hospitality of the Kavli Institute of Theoretical Physics during the workshop "From Genes to Growth and Form".

[^1]: The authors have declared that no competing interests exist.
