SPV Model.

The SPV model describes an epithelium as a network of polygons. Each cell i is endowed with a position vector 𝒓i, and cell shape is defined by the Voronoi tessellation of all cell positions (Fig. 1), which has been shown to provide a good representation of some real epithelia, such as the blastoderm of the red flour beetle Tribolium castaneum and the fruit fly Drosophila melanogaster (32). Like for vertex models, tissue forces are obtained from an effective energy functional E({𝒓i}) for N cells, given by (7, 15, 26, 33, 34)E=∑i=1NEi,Ei=KA(Ai−A0)2+KP(Pi−P0)2,[1]with Ai and Pi the cross-sectional area and perimeter of the ith cell, respectively. The first term in Eq. 1 arises from incompressibility of the layer in three dimensions and its resistance to height fluctuations, with A0 a preferred cross-sectional area. The second term represents the competition between cortical tensions from the actomyosin network at the apical surface and cell–cell adhesions from adhesive complexes at intercellular junctions (26), with P0 a preferred perimeter resulting from this competition. We simulate N cells in a square box of area AT, with A¯=AT/N the average cell area and with periodic boundary conditions. The system is initialized with a set of N random cell positions, independently drawn from a uniform distribution. Throughout the simulations, we set A¯=A0=1 unless otherwise noted, and KA=KP=1.

Each Voronoi cell is additionally endowed with a constant self-propulsion speed v0 along the direction of polarization 𝒏^i=(cos⁡θi,sin⁡θi) describing cell motility. The dynamics of each Voronoi cell are governed by∂t𝒓i=μ𝑭i+v0𝒏^i,∂tθi=2Drηi(t),[2]where 𝑭i=−∇iE is the force on cell i and μ is the mobility. The direction of cell polarization is randomized by orientational noise of rate Dr, with ⟨ηi(t)⟩=0 and <ηi(t)ηj(t′)>=δijδ(t−t′). The timescale τr=1/Dr controls the persistence of single-cell dynamics. As in self-propelled particle (SPP) models, an isolated cell performs a persistent random walk with a long-time translational diffusivity D0=v02/(2Dr) (29, 35, 36). After each time step, a new Voronoi tessellation is generated based on the updated cell positions. The cell shapes are determined in the process and the exchange of cell neighbors occurs naturally through topological transitions (13).

We showed in ref. 15 that the SPV model exhibits a transition from a solid-like state to a fluid-like state upon increasing the single-cell motility v0, the persistence time τr, or the cell shape parameter P0/A0 that characterizes the competition between cell–cell adhesion and cortical tension. The phase diagram in the (P0,v0) plane is reproduced in Fig. 1B. The transition is identified by setting the shape index q=⟨Pi/Ai⟩ to the value q= 3.813, where ⟨…⟩ denotes the average over all cells. It was shown in ref. 15 that the transition line located by q= 3.813 coincides with the one based on the vanishing of the effective diffusivity obtained from the cellular mean-square displacement. Note that for fixed system size AT and cell number N, the preferred cell area A0 does not affect the interaction forces or cellular shapes. Hence the solid–fluid transition is insensitive to A0, as shown analytically in SI Text. The preferred area A0 only shifts the total pressure of the tissue by a constant.

Developing and Validating Traction-Based Mechanical Inference.

It is well established that in a model tissue described by the energy Eq. 1 , the mechanical state of cell i is characterized by a local stress tensor σαβ(i)int given by (21, 23, 37)σαβ(i)int=−Πiδαβ+12Ai∑ab∈iTabαlabβ,[3a]Πi=−∂E∂Ai, 𝑻ab=∂E∂𝒍ab,[3b]where Πi is the hydrostatic cellular pressure and 𝑻ab=Tab𝒍^ab is the cell-edge tension, with 𝒍^ab=𝒍ab/|𝒍ab|. Here we use Roman indexes i,j,k,… to label cells and a,b,c,… to label vertices, while Greek indexes denote Cartesian components. The summation in Eq. 3a runs over all edges of cell i, and 𝒍ab is the edge vector joining vertices a and b when the perimeter of cell i is traversed clockwise, as shown in Fig. S1A. The factor of 1/2 in the second term on the right-hand side of Eq. 3a is due to the fact that each edge is shared by two cells. We have used the convention that the cellular stress is positive when the cell is contractile and negative when the cell is extensile. Contractile stress means that if the cell is cut off from its neighbors, it tends to contract, consistent with the contractility of the actomyosin network within the cell.

Note that in a vertex model the interaction stress as defined in Eq. 3a is indeed the stress acting on the tissue boundary. This is, however, not the case for the Voronoi model because the Voronoi construction introduces constraints not present in the vertex model. We have verified that the differences are small (SI Text) and in the following use Eq. 3a as a good approximation for the Voronoi model.

Our goal is to obtain the distribution of cellular stresses in a layer of motile cells, where cellular configurations do not minimize the tissue energy, but are governed by the dynamics described by Eq. 2. In this case, as discussed in the Introduction, the local cellular stress can be written as the sum of contributions from interactions and propulsive forces asσαβi=σαβ(i)int+σαβ(i)swim,[4]where, following recent work on active colloids (29, 30),σαβ(i)swim=−v0μAinαirβi[5]describes the flux of propulsive force v0μ𝒏i across a boundary, calculated from a virial expression. The negative sign in Eq. 5 ensures that the swim stress follows the same convention as the interaction stress, i.e., is positive for contractile stress and negative for extensile stress. The swim stress is proportional to the cell motility v0 and vanishes for nonmotile cells. The contribution σαβ(i)int is still given by Eq. 3a, but depends implicitly on cell motility through the instantaneous values of Πi and Tab that are determined not by energy minimization, but by the system dynamics governed by Eq. 2. From the local stresses one can then obtain the total mean stress in the tissue as (SI Text) σαβ=σαβint+σαβswim, withσαβint=1AT∑iAiσαβ(i)int,σαβswim=1AT∑iAiσαβ(i)swim.[6]In simulations where the energy functional is known, it is simple to directly extract the instantaneous pressures and tensions from cell shapes to calculate the interaction contribution σαβ(i)int. The definitions, Eq. 3b, giveΠi=−2KA(Ai−A0),Tab=2KP[(Pj−P0)+(Pk−P0)],[7]where ab is the cell–cell interface that separates cells j and k. Both Ai and Pi are obtained at every time step of the simulation and implicitly depend on cell motility, parameterized by speed v0 and persistence τr. This method directly infers the interaction stress from instantaneous cell shape fluctuations through Eq. 3a, yielding what we call shape-based stresses. While easily implemented numerically, this method is of limited use in experiments where the energy functional is not known and likely more complicated.

For this reason we develop a new mechanical inference method for motile monolayers that attempts to approximate the interaction stresses using only information that is accessible in experiments. Specifically, the proposed traction-based mechanical inference infers tensions and pressures from segmented images of cell boundaries and traction forces obtained by TFM. In the SPV model, we define the traction force at each vertex as the gradient of the tissue energy with respect to the vertex position 𝒕a=−∇aE, which balances the interaction force 𝑭a. Equilibrium mechanical inference methods express the interaction force 𝑭a=−∇aE at each vertex in terms of cellular pressures and edge tensions and the measured geometry of the cellular network (Eqs. S6–S8). Pressures and edge tensions are then obtained by inverting the equations 𝑭a({Πi},{Tab})=0. For a nonequilibrium epithelial layer of motile cells we invert the force balance equations𝑭a({Πi},{Ti})=𝒕a,[8]where the edge tensions Tab have been written as the sum of cortical tensions Ti of adjacent cells (Eq. S7), reducing the number of independent unknowns. The interaction contribution to the local cellular stresses is then calculated using Eq. 3a. A constraint counting yields 4N force balance equations for 2N variables, rendering the system overdetermined, which requires the implementation of a least-squares minimization for the mechanical inference (SI Text).

The equations developed thus far require knowledge of the tractions at each vertex, which is again not realistic in experiments. Therefore, we have developed and implemented a coarse-grained version of this approach that uses experimentally accessible traction forces averaged over a square grid, with a grid spacing of the order of a cell diameter. Pressures and tensions (Πi,Ti) are then calculated by inverting the force balance equation at the center of each grid element,𝑭grid({Πi},{Ti})=𝒕grid.[9]An outline of the coarse-graining procedure is given in SI Text. We refer to stresses inferred from Eq. 9 as traction-based stresses. We emphasize that the traction-based stress from mechanical inference is generally different from the intercellular stress obtained with monolayer stress microscopy (MSM), which rests upon the assumption that the tissue is an isotropic, homogeneous, and linearly elastic material (18). The mechanical inference does not make such assumptions and is compatible with any epithelia whose cell–cell interactions can be decomposed into tensions at cell junctions and pressures within cell bodies. Examples include Drosophila ectoderm (9) and mesoderm (25), Drosophila wing disk (26), and Madin–Darby canine kidney (MDCK) epithelium (12).

Within the framework of the SPV model, we have validated the coarse-grained method by showing that the resulting traction-based stresses agree with the shape-based stresses computed exactly from the simulations (Fig. 2 C and F).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Comparison of shape-based and coarse-grained traction-based stress. (A–C) Solid state at v0=0.5, P0=3.3. (D–F) Liquid state at v0=0.5, P0=3.8. (A and D) Interaction normal stress σn(i)int calculated from the instantaneous cell shapes obtained from Eqs. 3a and 7. Red denotes positive (contractile) stress and blue negative (extensile) stress. (B and E) Interaction normal stress σn(i)int calculated using the coarse-grained traction-based mechanical inference by inverting Eq. 9 and using Eq. 3a. The arrows denote the traction forces. (C and F) The coarse-grained traction-based mechanical inference is validated by plotting the traction-based stress against the shape-based stress in the solid (C) and in the liquid (F) state. The data are for 400 cells in a square box of side L=20 with Dr=0.1 and with periodic boundary condition.

Stress Characterizes Rheological Properties of the Tissue.

To study the mechanical properties of motile confluent tissues, we simulate a confluent cell layer with periodic boundary conditions, using the SPV model. By examining the temporal correlations of the mean stress in the tissue, as defined in Eq. 6, we show that the tissue displays distinct mechanical properties in the liquid and the solid states. Thus, mechanical measurements such as those provided by TFM can be used to characterize the rheological properties of the tissue.

The stress tensor σαβ is symmetric and has three independent components in two dimensions. Both the mean and local stress components are most usefully expressed in terms of normal stress σn, shear stress σs, and normal stress difference σd, withσn,d=12(σxx±σyy),σs=12(σxy+σyx).[10]Each of the interaction and swim contributions can similarly be split in normal, shear, and normal difference components. Below we focus on normal and shear stresses. TFM probes the forces exchanged between tissue and substrate, which by force balance are determined entirely by intercellular forces and hence by interaction stresses. In contrast, the swim components of stress and pressure cannot be probed in TFM, but contribute to the pressure Π=−σn=Πint+Πswim that the tissue would exert laterally on a confining piston. As we show below, the swim contribution dominates the pressure in the liquid state.

Using the expression for the local stress obtained from cell shapes, the mean interaction normal stress of the tissue can be expressed entirely in terms of area and perimeter fluctuations in a virial-like form (SI Text)σnint=1AT∑i[2KAAi(Ai−A0)+KPPi(Pi−P0)].[11]The first term represents the interaction contribution from the pressures within the cells. The second term is the contribution from the competition between actomyosin contractility and cell–cell adhesion that controls the cortical tensions. In our simulation, the cellular pressure is suppressed by setting A¯=A0 and the normal stress comes mainly from the cortical tensions. Eq. 11 then provides a way for extracting mechanical information directly from cell shape based on snapshots of segmented cell images.

Normal stresses are contractile in the solid phase and may become extensile deep in the liquid phase.

We show in Fig. 2 snapshots of the local interaction normal stress in the solid state (A–C) and in the liquid state (D–F). In both the solid and the liquid the interaction normal stress is on average contractile (red), with relatively weaker spatial fluctuations, but much larger mean value in the solid, where contractile cortical tension exceeds cell–cell adhesion.

Fig. 3 displays the total mean normal stress (the separate contributions from interaction and swim stress are shown in Fig. S2) across the solid–liquid transition. The color map shows that the total normal stress is contractile in the solid and across the transition line (Fig. 3A, black crosses), but changes sign and becomes extensile deep in the liquid. While the interaction stress is always positive due to cell contractility and consistent with experimental observations (1⇓–3, 38, the change in sign of the total stress is due to the swim stress that is zero in the solid and always negative in the liquid (Fig. S2), indicating that motility induces extensile stresses, tending to stretch the tissue. The total normal stress is analogous to the stress on a wall confining an active Brownian colloidal fluid (29, 30). We speculate that its change in sign could lead to an expansion of the tissue if released from confinement due to substrate patterning or to surrounding tissue and may contribute to epithelia expansion in wound-healing assays. In our model confinement is provided by the periodic boundary conditions.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Mean total normal stress. (A) Heat map of the mean total normal stress of the tissue in the (v0,P0) plane. The black crosses outline the solid–liquid phase boundary determined by q=3.813. Red indicates contractile stress and blue extensile stress. (B) Mean total normal stress as a function of P0 at various v0, showing a change in sign deep in the liquid state (400 cells for T=1,000 and Dr=0.1 with periodic boundary condition).

The tissue effective shear viscosity diverges at the liquid–solid transition.

While the local shear stress averages to zero in both the liquid and the solid states, its temporal correlations provide a distinctive rheological metric for distinguishing the liquid from the solid and identifying the transition. The time autocorrelation function of the interaction shear stress,Css(τ)=⟨σsint(t0)σsint(t0+τ)⟩t0,[12]where ⟨…⟩t0 denotes the average over the length t0 of the simulation, is shown in Fig. 4A for various v0 across the liquid–solid transition. Shear stress correlations decay in the liquid state and slow down as the transition is approached from the liquid side, becoming frozen in the solid. To quantify this we have defined the correlation time τm as the time when the correlation has decreased below 1% of its initial value. This correlation time shown in Fig. 4B diverges at the liquid–solid phase transition, suggesting that shear stress autocorrelations can provide a robust metric for the transition. Our work therefore suggests that TFM combined with mechanical inference can provide a tool for the measurement of tissue rheology. Correlating such measurements with cell shape data will provide a stringent test for our theory. Finally, in the liquid state we define an effective viscosity η𝑒𝑓𝑓, using a Green–Kubo-type relation by integrating the correlation function over the duration of the correlation time (41, 42),η𝑒𝑓𝑓=ATkbT𝑒𝑓𝑓∫0τmCss(τ)dτ,[13]where we have used the ideal gas effective temperature kBT𝑒𝑓𝑓=v02/(2μDr) so that η𝑒𝑓𝑓 has dimensions of a shear viscosity in two dimensions. Of course the Green–Kubo relation is based on the existence of a fluctuation–dissipation theorem, which does not hold in active systems. For small values of the persistence time τr, however, the orientational noise in the SPV becomes uncorrelated in time and can be mapped onto thermal noise at an effective temperature T𝑒𝑓𝑓. In this limit we expect η𝑒𝑓𝑓 to indeed play the role of a shear viscosity. Remarkably, the effective viscosity shown in Fig. 4 C and D grows as the tissue approaches the solid state from the liquid side and diverges at the transition. The effective viscosity quickly approaches zero deep in the liquid phase, suggesting that the system behaves as a gas of uncorrelated cells in this region.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Time autocorrelation of interaction shear stress and effective tissue viscosity. (A) Time autocorrelation function of the mean interaction shear stress for P0=3.45 at various v0. (B and D) Heat maps of the correlation time τm (B) and of the effective viscosity (D). The red crosses denote the liquid–solid phase transition boundary determined by q=3.813. The white region corresponds to regions where correlation time τm and effective viscosity η𝑒𝑓𝑓 diverge. Note that these regions largely coincide with the solid regime. (C) The effective viscosity as a function of P0 at various v0. The dashed lines correspond to the critical P0 where the liquid–solid phase transition occurs. All results are for 100 cells, Dr=1, and T=40,000 with periodic boundary condition.