Active Vertex Model for cell-resolution description of epithelial tissue mechanics

We introduce an Active Vertex Model (AVM) for cell-resolution studies of the mechanics of confluent epithelial tissues consisting of tens of thousands of cells, with a level of detail inaccessible to similar methods. The AVM combines the Vertex Model for confluent epithelial tissues with active matter dynamics. This introduces a natural description of the cell motion and accounts for motion patterns observed on multiple scales. Furthermore, cell contacts are generated dynamically from positions of cell centres. This not only enables efficient numerical implementation, but provides a natural description of the T1 transition events responsible for local tissue rearrangements. The AVM also includes cell alignment, cell-specific mechanical properties, cell growth, division and apoptosis. In addition, the AVM introduces a flexible, dynamically changing boundary of the epithelial sheet allowing for studies of phenomena such as the fingering instability or wound healing. We illustrate these capabilities with a number of case studies.


Author summary
We present a detailed analysis of the Active Vertex Model to study the mechanics of confluent epithelial tissues and cell monolayers. The model combines the commonly used Vertex Model for describing epithelial tissue mechanics with the active matter dynamics extensively studied in soft matter physics. We utlise an exact mathematical mapping that enables a very efficient numerical implementation using standard methods for simulating particle-based models. System sizes accessible to this model allow us to probe the dynamical motion patterns that occur in tissues over a range of length-and time-scales previously inaccessible to available simulation tools. Our model also includes a number of essential features required to properly describe actual biological systems such as cell growth, cell division and aptotsis, as well as the dynamic boundary of the epithelial sheet. This allows us to study phenomena such as the finger-like protrusions in cell monolayers and processes related to wound healing. The

Introduction
Collective cell migration [1,2] in epithelial tissues is one of the key mechanisms behind many biological processes, such as the development of an embryo [3], wound healing [4,5], tumour metastasis and invasion [6]. Due to their layered, tightly connected structure [7], epithelial tissues also serve as an excellent model system to study cell migration processes. Over several decades [8] extensive research efforts have been devoted to understanding molecular processes that lead to cell migration [9] 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 [10] experiments revealed that collective cell motion is far more complex than expected [11][12][13]. 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 [14]. These are reminiscent of the fluctuations observed in supercooled colloidal and molecular liquids approaching the glass transition [11] 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 [15]. It has also been argued that the migration patterns are sensitive to the expression of different adhesion proteins [16] as well as to the properties of the extracellular environment [17], such as the stiffness of the substrate [18,19]. These observations lead to the development of the notion of plithotaxis [12], 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 [5,[20][21][22][23]. 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 [23]. Subtle correlations between purse-string contractility and large-scale remodelling of the tissue while closing circular gaps have also been identified [22]. Finally, a mechanism dubbed kenotaxis has been proposed [20], 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. [24] 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 [25] 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 [26][27][28][29][30]. 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 [31,32]. Finally, it has recently become possible to track more than 100,000 individual cells in a chick embryo over a period exceeding 24 hours [33]. 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 [9], 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 [34]. One of the early yet successful approaches has been based on the cellular Potts model (CPM) [35,36]. 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 [37], which allows for simulations of large systems. In addition, the CPM extends naturally to three dimensions [38]. While very successful in describing cell sorting as well as certain aspects of tumour growth [39], 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) [40,41]. 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 [42][43][44][45][46]. 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 [47]. Over the past 35 years it has been implemented and extended numerous times and used to study a wide variety of different systems [48][49][50][51][52][53]. The VM is at the core of the cell-based [54] CHASTE [55], 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 twodimensional 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 [47,51,53,56,57]. 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 [53]. In contrast, recent traction microscopy experiments [14] 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 [51,58] 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 [59], in this paper we apply the ideas introduced in the physics of active matter systems [60] 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 [61], 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
Owing to its origins in the physics of foams [62], in the VM cells are modelled as twodimensional 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). The main simplification compared to models of foams is that most implementations of the VM assume that 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 0 i , 2) area modulus K i , and 3) perimeter modulus Γ i . In addition, a junction connecting vertices μ and ν has tension Λ μν . 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 [63,64]. 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 [58]. 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).
An energy is associated to each configuration of the mesh. It can be written as where N is the total number of cells, A i is the area of the cell i, while A 0 i 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). hμ, νi 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) 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 where P 0 i is a reference perimeter, and omit the last term in Eq (1) or completely omit the P 2 term [47,58]. 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) becomes Λ i ∑ hμ,νi l μ,ν = Λ i P i . Therefore, if we identify L ¼ À GP 0 i , it immediately becomes clear that the descriptions in Eq (1) and the model with the preferred perimeter are identical to each other. Note that this is true up to the constant term 1=2G i ðP 0 i Þ 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) 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 2 i 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). 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) for a given set of parameters K i , Γ i and Λ μν . Most implementations [51,53,58], also introduce topology (connectivity) changing moves to model events such as cell rearrangements.
There have been several attempts to introduce dynamics into the VM [47,51,56], including a recent study that introduced stochasticity into the junction tension [53]. The idea behind such approaches it to write equations of motion for each vertex as 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) with respect to r μ , i.e., F μ = −r 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) 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 quasistatic 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 [58].

Active Vertex Model
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 nonthermal 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 [65] 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 [53].
Here we take an alternative approach and build a description similar to the recently introduced SPV model [59]. 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). 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 [66]. 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 [67].
In the original implementation of Bi, et al. [59] 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) [68]. 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 nonlinear (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), which drastically increases the efficiency of the Delaunay based approach and enables us to simulate systems containing tens of thousands of cells.
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).
Force on the cell centre. In the AVM, the area A i in Eq (1) corresponding to the cell (particle) i is the area of the polygon (face) of the Voronoi polygon, O i , associated with i. It is given as 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 This expression is just a discrete version of Green's formula. Similarly, the cell perimeter is 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) 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), with Λ = λ 1 + λ 2 + λ 3 . Armed with the mapping in Eq (5), we can proceed to compute the forces acting on each cell centre based on the VM energy functional in Eq (1). This is done directly by computing the negative gradient of the energy. Therefore, the force on the centre of cell i is computed as 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) we obtain: In the last expression, @r n @r i h i is the 3 × 3 Jacobian matrix (see S1 Appendix) 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) are equal to zero, since each vertex coordinate r ν depends only the cell centres r k associated to its Delaunay triangle (see Fig 2b). 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.
This sum is over all vertices (corners) ν of cell i.

For all immediate neighbours
and multiply it by the sum Note that ν 2 O i \ O 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 [69] in order to speed up force computations. It is also evident that the force in Eq (7) is not pairwise, i.e., it cannot be written as a sum of forces acting on a pair of cells, or in mathematical terms, This is not surprising since the position of the Voronoi vertices depends on the position of three cell centres. Cell alignment. The AVM also includes the cell polarity [70] 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, 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 [71]. 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 [72,73], 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 [74,75] assumes that cell polarisation vector aligns with the migration direction of the cell, i.e., wherev 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 where J s > 0 is the alignment strength and p i is the eigenvector of the matrix [76] 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 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. 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 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) and the soft-core repulsion defined in Eq. (35) in S1 Appendix. ν i ðtÞ 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 ν i ðtÞ 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 hν i,α (t)i = 0 and hν i, where D measures the average magnitude of the stochastic force and is interpreted as an effective translational diffusion coefficient, α, β 2 {x, y} and h. . .i denotes the statistical average. We note that in some cases, such as in the early Drosophila embryo, or the cultured monolayer system [77] 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 [33], 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 [60]. 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 [78], 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 [79]. 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 where τ i is the torque acting on n i , given in Eq (12), 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 n r i ðtÞ is responsible for introducing an orientational randomness. Akin to the stochastic force ν i ðtÞ in Eq (13), n r i ðtÞ has properties hn r i ðtÞi ¼ 0 and hn r i ðtÞn r j ðt 0 Þi ¼ 2D r d ij d ðt À t 0 Þ, where D r is an rotational diffusion constant. The related time scale τ r = γ r /2D 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 [80]. Eqs (13) and (14) 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. 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., 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 [81]. 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). 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 [82]. 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., 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 [83], 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. 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 [84]: 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). The procedure is repeated until there are no more edges left to flip. One can show [84] 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). 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. The implementation of the SPV by Bi, et al. [59] 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. 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.
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, 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 where z i is the stiffness of angle θ i at the boundary particle i, determined as 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.
This concludes the description of the AVM. In

Results and discussion
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
We start by illustrating one of the key processes observed in epithelial tissues, the T1 transition. As detailed in Fig 2, 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, 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 [85] 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 [49,81,86]. 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, 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 [87], 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 [86], which would predict a series of fluctuations through which the barrier between minima is ultimately crossed.

Activity driven fluidisation phase diagram
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), in a scaled form [49,59]. We first choose K i = 1 and set A 0 i ¼ p 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), we can complete the square on the second and third terms in Eq (1) and obtain the scaled VM potential where P 0 = −Λ/Γ and A 0 = π. The first term in Eq (20) 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), giving a natural scale that is determined by the relative ratio of GP 2 0 to K A 2 0 . In other words, if K A 2 0 > GP 2 0 , 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 A 2 0 < GP 2 0 . Bi, et al. [59] introduced the dimensionless shape factor p 0 ¼ P 0 ffiffiffi ffi Remarkably, one observes [49,59,86] 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 [88]. 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 hAi = A 0 . The boundary is either fixed (referred to as "fixed system"), or allowed to fluctuate freely ("open system"). Fig 6 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), or 250,000 time steps with δt = 0.025 in the solid-like region (e.g., Fig 6c). 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 liquidlike phase).
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 ¼ ffiffiffiffiffi A 0 p 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, 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. 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. In agreement with the results of Bi, et al. [59], 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  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 [15] 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 [15], and stays infinite within the solid phase. In Fig 7a-7c, we show the phase diagram of τ α as a function of p 0 and f a , for several systems with different boundary conditions. In Fig 7a 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, we show the phase diagram for the same Γ = 1 and t À 1 r ¼ 0:01 as in Fig 7a, 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 t À 1 r ¼ 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 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). This unexpected result is accompanied by structural changes in the cell patterns that we observe. Fig 6f 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, 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 [89], 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, 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 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 [90][91][92][93] where regions a few cells wide migrate outward from the centre, as shown in Fig 6e. When τ α drops below approximately 10, we observe that the fluctuations of the boundary already present in Fig 6d  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). 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. 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 shows τ α , and Fig 7d shows τ inst for the same open system with Γ = 0.1, t À 1 r ¼ 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 t À 1 r (lower t À 1 r enhances the instability), and p 0 . The instability resembles observations of finger formation in MDCK monolayers [90,91]. Existing models link it to either leader cells [93,94], a bending instability [92], or an active growth feedback loop [95], while here it emerges naturally. A detailed account of this phenomenon will be published elsewhere.

Growth and division
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.
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, 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 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, where cells located towards the centre have, on average, smaller areas and in Fig 8g, 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 [96] 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, where D is the region for which the stress is calculated and A D is the area of D. For convenience we choose D to coincide with individual cells. Finally, f ab ¼ À @E VM @r abr ab , where E VM is defined in Eq (1) andr ab is the unit-length vector along the junction or triangulation edge denoted by (α, β). This particular choice of ω and D reduces our calculation to that of the more common virial stress [96], but we note that the Hardy stress is a useful tool for coarsegraining. 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. In  Fig 8h, 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 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.8a 2 after which the probability of a cell to divide increases linearly with its area (see Eq (16)). 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. in Fig 8g. 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 [12,97]. A similar pressure buildup has also been investigated in a growing cell aggregate [98] and in a mathematical model of nonuniform growth in a layer of tissue [99].
Finally, in Fig 8i 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 [100,101]. 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
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 cellcell 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 [35,102]. 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 [103]. 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. 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). 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 [102]. This case is shown in Fig 9e for Λ rM = −6.8.
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
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), 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, we have used the alignment model defined in Eq (9) that assumes that the polarity vector n i of cell i aligns with its velocity v i . The torque term in Eq (14) 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) [74], 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 [75]. It is remarkable that the main hallmarks of this   9. (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.
active matter dynamics are also observed in the model tissue. In Fig 10a, 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 [104,105].
In Fig 10b, 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, in the absence of confinement, we recover the collective polar migration of the cell patch.
In Fig 10d, we show the effects of aligning the cell's polarity to the largest principal axis of the cell shape tensor, defined in Eq (10). 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
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 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 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 [90] is shown in Fig 11b, where three separate patches grow towards each other. Finally, in Fig 11c we show an example of yet another very interesting situation [106], where cells are grown in a confined region of space.

Summary and conclusions
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 [59], 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. Fig 11. (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. 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. [43,44], 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. 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 [85,[107][108][109][110][111][112]. 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 [113][114][115][116]. 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 S1 Appendix. Detailed derivation of Eq (7), the AVM force on the cell centres. (MP4) S13 Video. 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, t À 1 r ¼ 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. (AVI) S14 Video. 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) S15 Video. 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) S16 Video. 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) S1 Dataset. Compressed folder containing input data and configuration files for reproducing all figures and supplementary movies in this paper. (ZIP)