Shape-tension coupling produces nematic order in an epithelium vertex model

We study the vertex model for epithelial tissue mechanics extended to include coupling between the cell shapes and tensions in cell-cell junctions. This coupling represents an active force which drives the system out of equilibrium and leads to the formation of nematic order interspersed with prominent, long-lived $+1$ defects. The defects in the nematic ordering are coupled to the shape of the cell tiling, affecting cell areas and coordinations. This intricate interplay between cell shape, size, and coordination provides a possible mechanism by which tissues could spontaneously develop long-range polarity through local mechanical forces without resorting to long-range chemical patterning.

Introduction.-Epithelialcell monolayers and tissues are prime examples of dense, active, viscoelastic materials [1][2][3].Understanding their behaviour is also of fundamental importance in biology and medicine since epithelia line all organs and cavities in the body, and the majority of cancers develop in epithelial cells [4].The vertex model [5][6][7] has played an important role in modelling epithelial mechanics.It can naturally capture the solid-to-fluid transition observed in experiments [8] by tuning its geometric parameters [9].It is also straightforward to extend it to include active effects, both as self-propulsion [10,11] and as cell junction activity [12][13][14][15][16][17][18][19].The vertex model is, therefore, a powerful tool to study active processes in tissues at the cell scale.
In many epithelial systems, one also observes nematiclike features at scales that span several cells [20][21][22][23][24].These are a readout of elongated cell shapes, and the defects associated with the nematic state have been argued to play important biological roles, e.g. as sites of cell extrusions [21] or stress-organising centres that drive tissue morphogenesis [25].Theories of active nematics [26][27][28][29] have also been successful in explaining large-scale collective behaviours such as active turbulence both in epithelial cell monolayers [30,31], bacterial suspensions [32] and cell membranes [33].
Cell motion in biological systems is, however, rarely turbulent, and instead one often observes coordinated motion over distances much larger than the typical cell size.Such coordinated movements are key, e.g. during embryonic development [34][35][36].One of the central open questions is how these motions are generated, sustained, and regulated.A closely related question is to what extent such large-scale features require guidance by biochemical patterning, e.g.via spatio-temporal coordination of morphogens, or whether they can spontaneously emerge as a result of cellular behaviours.It is, therefore, important to understand how cell-level processes coordinate to form tissue-and organ-scale structures.
It has recently been shown that coupling between ten-sion and a global nematic field leads to active T1 transitions that drive tissue shape changes and can elongate cells [19].In this paper, we explore how nematic order can emerge in a vertex model by introducing coupling between the local cell shape, a proxy for the nematic director, and the tension on cell-cell junctions.We find that this model also leads to prominent +1 defects in the nematic order, and we primarily focus on the role of these defects in determining local cell shapes and global tissue tiling.
Model details.-We begin with the vertex model for planar epithelia [6].Each cell is represented by a polygon, two cells share a junction, and three or more cells meet at a vertex.The dynamics of the vertices is described by overdamped equations of motion, Here, R (i) is the position of the i−th vertex, the overdot denotes the time derivative, ∇ R (i) indicates the gradient with respect to R (i) , E VM is the elastic energy of the model tissue, η is the viscosity, and act is the active force on the vertex due to coupling between cell shape and junctional tension.The explicit form of The elastic energy of the vertex model takes the form where the sum is over all cells.A (c) and P (c) are, respectively, the area and perimeter of the cell c, whereas A are its target area and perimeter, which we assume to be identical for all cells.K A and K P are the area and perimeter elasticity moduli.
To recast the expression for energy into a dimensionless form, we divide it by K A A 2 0 , which gives where e VM is the dimensionless energy, measured in units of K A A 2 0 , and a (c) = A (c) /A 0 is the dimensionless area of the cell c, measured in units of A 0 .This sets the unit of length as a * = √ A 0 .Finally, p (c) = P (c) / √ A 0 , k P = K P /(K A A 0 ), and the target cell-shape index is defined as p 0 = P 0 / √ A 0 .The parameter p 0 plays an important role in determining the mechanical response of the vertex model.Namely, if p 0 exceeds a critical value, the model undergoes a solid-to-fluid transition [37].For regular hexagonal tilings, this occurs at p c 0 = 8 √ 3 ≈ 3.72, while for a random tiling, this occurs at a critical value p c 0 ≈ 3.81 − 3.92 [38][39][40].If time is measured in units of t * = η/ (K A A 0 ), and force in units of f * = K A A 3/2 0 , the equations of motion (1) can be recast into a dimensionless form, where r i is the dimensionless position of the i−th vertex and the overdot and ∇ r (i) now, respectively, indicate the derivative with respect to dimensionless time and the gradient with respect to r (i) .
The dimensionless active force takes the form where the sum is over the set J of all cell-cell junctions, γ (j) (t) is the tension of the j−th junction at time t, and l j is its length.The tension in a junction evolves according to where τ γ sets a characteristic relaxation time scale and γ (j) 0 is the target tension of the junction.Here, it is selected to couple the tension in the junction with the elongation of its neighbouring cells, by choosing ϑ and ϑ are the angles between junction j and the directors of its neighbouring cells c and c (Fig. 1a), and ζ is a coupling constant.The sign of ζ determines whether the active forces act to extend or contract a junction that is aligned with the cell's director.We define a cell's director to point along the eigenvector corresponding to the largest eigenvalue of the cell's gyration tensor, given as Here the sum is over the set V c of n (c) vertices of the cell c, and r , is the position of the cell's geometric centre.Cells for which the gyration tensor has two identical eigenvalues do not contribute to Eq. (7).
It is important to note that the choice of the coupling between the cell geometry and junction tensions is key for making the system active.Namely, for the above model, the active force cannot be written as a gradient of a line tension contribution to the energy, which would lead to the dynamics of the system corresponding to passive energy minimisation.Instead, the tensions γ (j) are coupled to the instantaneous geometry through Eqs. ( 7) and ( 6), but the resulting forces are only along the junctions.The movement of the vertices they produce can and does change the tissue shape in such a way that the coupling in Eq. ( 7) further increases the energy because angles θ and θ have changed.This is qualitatively different from inserting γ 0 from Eq. ( 7) directly into an energy of the form γ (j) 0 l (j) , in which case the gradient would lead to additional terms that would rotate junctions, and the model tissue would relax towards a local energy minimum.Effectively, in our model, the gradient does not "know" that the angles θ and θ are included in the line tension.As a result, the movement of the vertices does not in fact minimise the energy, rendering the system active.Figure 1c shows how the energy of the model tissue changes in time.
We solve the equation of motion [Eq.( 4)] using the first-order Euler scheme with the time step δt = 0.01.In all simulations, we set k p = 0.02 and τ γ = 1.We start simulations with a hexagonal lattice in which each vertex is perturbed by a displacement vector δr (i) = δr (i) (cos(α (i) ), sin(α (i) )), where δr (i) and α (i) are uniformly distributed random numbers respectively drawn from the intervals [0, 0.1] and [0, 2π].This is necessary to generate non-zero active forces that otherwise vanish for regular hexagonal tiling.Unless otherwise specified, we used a 32 × 32 lattice of cells placed in a periodic simulation box of fixed size.The longest simulations were run until t max = 4 • 10 5 .
Finally, to capture plastic events (i.e.local cell rearrangements), we implemented T1 transitions [41].These are neighbour exchanges facilitated by the shrinking of a junction shared by two cells and then creating a new junction between cells that were previously separated.T1 transitions are implemented as follows: If the length of a junction falls below a threshold value l 0 = 0.01 and has decreased since the previous time step, the junction is merged into a four-fold vertex.The four-fold vertex is maintained for ∆t = 0.03.It is then resolved into the new configuration, separating the two newly created vertices by l = 0.001.Immediately following this resolution, the tension of the new junction is set to 0. Note that our implementation does not allow for a vertex with more than four neighbours, so a junction that is connected to a four-fold vertex cannot undergo a T1 transition.
Results.-We focus on the case ζ < 0, i.e.where junctions are under a higher tension if they are aligned with local cell elongation.Above a threshold value ζ c , cell shapes are irregular polygons, but all cells have the sixfold coordination of the initial configuration.For ζ < ζ c , however, a local alignment of elongated cells starts to develop (Fig. 1b).The absolute value of the threshold ζ c at which local nematic order emerges decreases with increasing p 0 .
To quantify the emerging nematic order, we calculated the order parameter −0.1, p 0 = 3.2, shown in blue in Fig. 2a) cells do not move significantly and defects in the nematic ordering do not disappear from the tissue.
To illustrate this further, we calculated the mean square displacement of cells from their initial positions as a function of time, where N C is the total number of cells and r (c) 0 (t) is the position of the centre [defined below Eq. ( 8)] of cell c at time t, whereas r (c) (0) is its initial position at the start of the simulation.The resulting plot is shown in Fig. 2b.While the movement of the tissue is for the most part arrested once the fully nematically ordered state develops, it is important to note that due to active tension coupling this state is not necessarily the result of energy minimisation.
This nematic-like state also features very high values of the cell-shape index, as shown in Fig. 2c.In all cases, we find average cellshape indices well above 3.81, which is generally associated with a fluid-like tissue.This remains the case even after cell movement has been arrested.
A striking feature of the cell configurations is the presence of prominent, vortex-like +1 defects where cells form concentric rings around the defect core (Fig. 1b).These defects either remain for the entire duration of the simulation or vanish over time, for parameter values that lead to the global nematic order.Figure 2d shows the number of +1 defects as a function of time.
While +1 defects are featured most prominently, defects with different topological charges are also present and can therefore be involved in the annihilation of +1 defects.It should, however, be stressed that cells in the vertex model are not hard rods and can deform in such a way that the director changes discontinuously, e.g.going through an isotropic intermediary shape for which the gyration tensor [Eq.( 8)] has two identical eigenvalues and therefore the cell has no director.Therefore, topological charge is not necessarily strictly conserved.
As noted, the negative values of ζ favour high tension in junctions that are aligned with the cell director.This leads to the formation of concentric circles of high-tension junctions around vortex-like +1 defects (Fig. 3a) that constrict the defect core.To quantify this, we define the angle φ between the direction of the junction and the radial vector connecting its centre to the defect core.In Fig. 4a we show how the correlation of the angle φ and the junction tension changes as a function of distance from the core Appendix B for details).result of the concentric circles of high-tension junctions, cells close to the core of a nematic vortex-like +1 defect are compressed compared to the cells further away, and the area elasticity balances the junctional tensions (Figs.3b  and 4b).This is a plausible explanation for the long life of +1 defects, as the circles of high-tension junctions prevent a +1 defect from splitting into a pair of +1/2 defects.
The measured cell-shape index [Eq.(11)] also depends on the position of +1 defects (Figs.3c and 4c).In particular, we find a decrease in q in the defect core, followed by a maximum moving outwards from the defect.While cells in the defect core generally have smaller areas, they are also more isotropic, leading to smaller perimeters and, therefore, lower shape indices.In turn, the cells near, but outside, the defect core has a more pronounced elongation compared to those further away, leading to a higher q .Directly measuring how elongation changes with distance from the defect core gives a similar profile to that seen in Fig. 4c (see Appendix B).
The nematic order couples back to the cell tiling through these +1 defects.In particular, we find that cells with five neighbours are very common in the vicinity of the +1 defect cores (Fig. 3d).To quantify this, we computed the distribution function, where the sum is over all n d +1 defects in all 20 simulations used, and n (r) is the number of n-sided cells at distances between r and r + ∆r from +1 defect d.The normalisation was chosen such that the two-dimensional integral of the function would give the number of n−sided cells in the model tissue if f n (r) were based on a single simulation run.Interestingly, both +1 defects [25,42] and five-neighbour cells [43][44][45] in the tiling have previously been reported as connected to the formation of budded structures in 3D tissues, so their co-location in our model may be relevant, though typically the experimentally reported +1 defects leading to budding are of the aster type.
Summary and Discussion.-In this paper, we analysed a vertex model extended to include coupling between the elongation of cells and junctional tensions, leading to an active force on the vertices.For appropriate values of parameters, the model tissue forms nematic ordering of elongated cells which, surprisingly, features prominent vortex-like +1 defects.Experimentally, +1 defects are less common in biological systems than ±1/2 defects, but they have been reported both in vitro [46] and in vivo [47,48].It has also been recently argued that +1 defects play a role in morphogenesis [25,42].
Moreover, defects in the nematic order are coupled to the tiling of the confluent tissue.Cells around the +1 defects arrange themselves in such a way that the defects are surrounded by nearly concentric circles of hightension junctions.This leads to the compression of cells near the defect cores.Moreover, cells that form the defect cores typically have five neighbours.
An alternative way to introduce active dipolar forces into the vertex model has been proposed in Ref. [19].
The key difference to our approach is that tensions in our model are coupled directly to the local elongation of cells, rather than to an external, uniform nematic field.This allows the active tensions to reorient the nematic field to which they are coupled and does not require global patterning to drive the formation of nematic order.Furthermore, recent studies have proposed mechanisms by which extensile stresses can arise in a purely contractile nematic system due to polar fluctuating forces [49,50] or active inter-cellular interactions [51] Such mechanisms could provide insights into the behaviour of epithelial monolayers on a substrate.
We also note that the emergence of nematic order due to activity arising from coupling with the local cell elongation has been reported in Ref. [52].That study introduces activity via an active stress term proportional to a Q tensor, constructed from the positions of the junctions of each cell.This results in different active forces on vertices.It leads to both tissue fluidisation and the emergence of nematic order with ±1/2 defects, for sufficiently high activity.
Regarding the effects of noise that are inherently present in all biological systems, the movement of cells in our model arises directly from the shape-tension coupling [Eqs.(7) and ( 6)] and does not rely on noise mod-elled as an Ornstein-Uhlenbeck process used, e.g. in Refs.[12,15,19].A noise term was, therefore, omitted from our model.
Finally, the model studied here assumes a specific form of coupling between cell shape and junctional tension.As discussed, this results in an active force that leads to rich collective behaviours.It is, therefore, important to ask how such local coupling could arise in real tissues.There is evidence that cells can sense their shape [53].In, e.g. the fly, large-scale chemical patterning of cytoskeletal molecules is observed [54][55][56] that gives global directionality to the tissue, making the model of Duclut, et al. [19] applicable.On the other hand, in systems such as earlystage avian embryos [35], there is no such global patterning, yet local anisotropy of cell shapes and actomyosin orientation is apparent, albeit with no clear nematic order.While the molecular details of such coupling are likely to be very intricate, this suggests that it is plausible to consider a scenario in which a cell can locally inform its junctions about its current direction.result may slightly depend on the order in which the vertices are considered for this step).We found that this approach is generally more accurate for finding +1 defects compared to the grid-based approaches implemented in, e.g.Refs.[50,52].
This algorithm does not distinguish between vortexand aster-shaped defects.Both are therefore included in calculations for Figs. 4 and 5. Aster-shaped defects are, however, not common, so the statistics pertain chiefly to vortex-like +1 defects.

FIG. 1 .
FIG. 1.(a) Schematic representation of the model.Thick red lines indicate the orientations of cell directors.The tension of an edge j (blue line) shared between the two cells depends on its alignment with the cell directors.(b) Snapshots of a simulation of the vertex model with ζ = −0.18 at time t = 10 3 (left) and t = 10 4 (right).Cell directors are shown in red, and cores of +1 nematic defects are shown as blue disks.In the long time limit, defects can disappear.In panel (b), for clarity, only a quarter of the actual simulation box is shown.(c) Tissue energy as a function of time, including the γ(t)l term.

FIG. 2 .
FIG.2.Time-dependence of (a) nematic order parameter S, (b) mean square displacement of cell centres, (c) average cell-shape index q of the tissue, (d) number of +1 nematic defects.All panels show data for the same three parameter sets, and the values obtained are averaged over 10 simulation runs.

FIG. 3 .FIG. 4 .
FIG.3.Zoom in on a +1 nematic defect from Fig.1b(defect core shown in red).Cells are coloured as follows: (a) with edges coloured according to the tension arising from the shape-tension coupling, (b) according to the cell area, (c) according to the cell-shape index [Eq.(11)],(d) with cells with five neighbours in cyan and others in white.Cell directors are shown in black on all panels.

FIG. 5 .
FIG. 5. (a) Schematic representation of the ingredients used to construct the correlation function in Fig. 4a.Junctions are coloured as in Fig. 3a.(b) Average cell elongation as a function of distance from +1 nematic defects, combining data from 20 simulations with p0 = 3.2 and ζ = −0.18 at time t = 10 3 , for a 64 × 64 lattice.