Active Vertex Model for cell-resolution description of epithelial tissue mechanics
[['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 a
['Introduction', 'Model', 'Results and discussion']
['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 .']
['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']