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

Abstract

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.

Fig 1. In the Vertex Model (VM), a confluent epithelial sheet is represented as a polygonal tiling of the plane with no holes or overlaps.

Each cell is represented by an n−sided polygon. Neighbouring cells share an edge, which models the cell junction as a straight line. Three edges meet at a vertex (dark blue dots). The behaviour of cell i is described by three parameters: 1) reference area $A_i^0$, 2) area modulus K_i, and 3) perimeter modulus Γ_i. In addition, a junction connecting vertices μ and ν has tension Λ_μν.

Fig 2

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

Fig 3. Schematic representation of the allowed and disallowed changes of the boundary in the AVM.

An initial configuration with the topology of a disk (a) is allowed to develop pronounced fingers (b). However, it is not possible to split into two domains (c) or develop a hole (d), which would both lead to the introduction of a new boundary lines, and therefore lead to changes in the topology.

Fig 4. Flow chart of the force and torque calculations in each time step based on Eqs (7) and (12).

The calculation speed is significantly improved by using the fast equiangulation moves shown in Fig 2c. The Delaunay triangulation is recomputed only if there are cell division and/or death events, which do not occur at every time step.

Fig 5. Time lapse of a T1 transition, each label is the time (in units of γ/Ka², with a = 1) with respect to the T1 event.

Initially, cells 2 and 4 are in contact (red line). Approaching the transition, the connecting line slowly contracts, until it becomes a point at the transition. Cells 1 and 3 make a new contact which then rapidly expands. The arrows represent the force on the cell centre resulting from the Vertex potential (green), which partially compensates the active self-propulsion force (yellow) to give the resultant total force (red). The Voronoi tessellation is outlined in black, and the Delaunay triangulation is in white.

Fig 6. Representative snapshots of states seen with fixed and flexible boundaries in the AVM.

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

Fig 7. Phase diagrams of the AVM.

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

Fig 8. Snapshots of a growing epithelial tissue.

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

Fig 9

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

Fig 10. Effects of alignment.

(a-c) are self-alignment of the cell polarisation with the cell velocity. (d) is self-alignment of the cell polarisation with the long axis of the cell. (a) A confined system in the solid-like region (f_a = 0.03, p₀ = 3.385, ${\tau }_r^{-1}=0.01$), at alignment strength J_v = 1.0 shows oscillating collective modes. (b) In the liquid region (f_a = 0.1, p₀ = 4.4, ${\tau }_r^{-1}=0.1$), the system enters a vortex rotation state instead (J_v = 1.0). (c) An open system at f_a = 0.1, p₀ = 4.4, ${\tau }_r^{-1}=0.1$, boundary line tension λ = 0.1 and J_v = 0.1 migrates collectively. (d) An open system at f_a = 0.1, p₀ = 3.95, ${\tau }_r^{-1}=0.1$, boundary line tension λ = 0.1 and alignment with cell shape with J_s = 1 migrates collectively with complex fluctuations. Arrows are cell velocity vectors v_i; they are coloured by the magnitude of v_x for (a) and v_y for (b-d). Orange is positive (pointing right/upwards) and blue is negative (pointing left/downward).

Fig 11

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