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.
Owing to its origins in the physics of foams [62], 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). 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 [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 E V M = ∑ i = 1 N K i 2 A i - A i 0 2 + ∑ i = 1 N Γ i 2 P i 2 + 2 ∑ μ, ν Λ μ ν l μ ν, (1) where N is the total number of cells, Ai is the area of the cell i, while A i 0 is its reference area. Ki 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. Pi 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 Pi. 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) 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 ∑ i Γ ˜ i / 2 (P i - P i 0) 2, where P i 0 is a reference perimeter, and omit the last term in Eq (1) or completely omit the P2 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 ∑〈μ, ν〉 lμ, ν = ΛiPi. Therefore, if we identify Λ = - Γ ˜ P i 0, 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 / 2 Γ ˜ i (P i 0) 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 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). 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 Ki, Γ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 γ d r μ d t = F μ, (2) 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μ EVM. 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 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 [58].
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 [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 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), 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).
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. We now explore different modes of collective behaviour, i. e. , phases, of the tissue based on the values of parameters of the original VM (Ki, Γi and Λμ, ν), and AVM-specific parameters such as the activity fa, 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 Ki = 1 and set A i 0 = π 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 E V M = ∑ i = 1 N 1 2 A i - A 0 2 + ∑ i = 1 N Γ 2 P i - P 0 2, (20) where P0 = −Λ/Γ and A0 = π. 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 A0 to be generically compatible with the preferred perimeter P0. This sets up a competition between the two terms in Eq (20), giving a natural scale that is determined by the relative ratio of Γ P 0 2 to K A 0 2. In other words, if K A 0 2 > Γ 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 A 0 2 < Γ P 0 2. Bi, et al. [59] introduced the dimensionless shape factor p 0 = P 0 A 0, which controls the ratio of the cell’s perimeter to its area through the target area A0 and target perimeter P0. The value of p0 then determines whether the area or the perimeter term in Eq (20) wins and effectively sets the preferred shape of each cell: cells of different shapes have different values of p0. For example, regular hexagons, pentagons, squares and triangles correspond to p0 = 3. 722, p0 = 3. 812, p0 = 4. 0 and p0 = 4. 559, respectively. 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 p0 = 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 p0 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〉 = A0. 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 liquid-like phase). The unit of time is set by γ/Ka2, where a ≡ 1 is the unit of length. We note that Bi, et al. use a = A 0 as the unit of length. This is possible as long as cells are not allowed to grow, i. e. , when A0 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 p0, we find a system that prefers to be in a state with mostly hexagonal cells, unless the active driving fa is very high. Open systems will shrink at this point so that all cells are close to their target P0, 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 p0 < 3. 81 and low values of driving fa, 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 = exp i q · r (t) - r (0). (21) 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 p0 and fa, 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 p0 ≈ 3. 81 at small fa to a maximum activity fa beyond which the system is fluid at all p0. This is qualitatively, but not quantitatively consistent with the results of Bi, et al. , who find a transition line at roughly twice our fa 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 τ r - 1 = 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 fa for the transition, fa = 0. 03, a factor of 10 compared to the fixed case. We note that the effect also persists at τ 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 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 p0 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 p0 = 3. 95, near the relative minimum τα for a given fa. The distribution of cell centres appears random. In contrast, as can be seen in Fig 6g, at very high p0 = 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 p0, i. e. for p0 ≳ 4. 6. In parts of the phase diagram, we observe a fingering instability [90–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, τ 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 fa = 0. 03. At and below fa = 0. 1, the boundary of the system is stable, and above this threshold, the instability becomes more pronounced with increasing fa and smaller τα. The physical mechanism responsible for the instability involves a subtle interplay of fa, boundary line tension (stronger line tension suppresses the instability), the noise level τ r - 1 (lower τ r - 1 enhances the instability), and p0. 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. 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, p0 = 3. 10, corresponding to Γ = 1 and Λ = −5. 5 and no active driving, i. e. , we set fa = 0. This puts the system into the solid-like phase where T1 events are absent. Our simulation runs for 106 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 p0 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 p0. 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] σ ^ w H = 1 2 ∑ α, β α ≠ β ∫ s = 0 1 - f α β ω (1 - s) r α + s r β - r ⊗ r α - r β d s, (22) 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, ω (r) = 1 / A D if r ∈ D 0 otherwise (23) 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 α β = - ∂ E V M ∂ r α β r ^ α β, where EVM is defined in Eq (1) and r ^ α β 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 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. 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 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. 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 [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 p0 in the range 3. 58–3. 93. We, however, note that p0 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, fa, 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. 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 ni of cell i aligns with its velocity vi. The torque term in Eq (14) is then given by τi = −Jv vi × ni, where Jv 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 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. 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. 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. 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–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–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.
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 model is implemented into the SAMoS simulation software package and is publically available under a non-restrictive open source licence.
cell physiology
cell motility
medicine and health sciences
classical mechanics
cell cycle and cell division
cell processes
cell polarity
simulation and modeling
physiological processes
developmental biology
embryos
wound healing
research and analysis methods
mechanical tension
embryology
tissue repair
biophysics
physics
cell biology
physiology
cell migration
biology and life sciences
physical sciences
computational biology
biophysical simulations
Active Vertex Model for cell-resolution description of epithelial tissue mechanics