Correlating Cell Shape and Cellular Stress in Motile Confluent Tissues

TL;DR

This study links cell shape, traction forces, and stress in motile tissues using a model that reveals how mechanical interactions influence tissue rheology and transitions between fluid-like and solid-like states.

Contribution

We developed a generalized inference method connecting traction forces to cellular stresses and characterized tissue rheology and phase transitions in motile confluent tissues.

Findings

Traction-based stresses match those from cell shapes.

Identified a motility-induced swim stress affecting tissue contractility.

Discovered divergence of effective viscosity at the liquid-solid transition.

Abstract

Collective cell migration is a highly regulated process involved in wound healing, cancer metastasis and morphogenesis. Mechanical interactions among cells provide an important regulatory mechanism to coordinate such collective motion. Using a Self-Propelled Voronoi (SPV) model that links cell mechanics to cell shape and cell motility, we formulate a generalized mechanical inference method to obtain the spatio-temporal distribution of cellular stresses from measured traction forces in motile tissues and show that such traction-based stresses match those calculated from instantaneous cell shapes. We additionally use stress information to characterize the rheological properties of the tissue. We identify a motility-induced swim stress that adds to the interaction stress to determine the global contractility or extensibility of epithelia. We further show that the temporal correlation of the interaction shear stress determines an effective viscosity of the tissue that diverges at the liquid-solid transition, suggesting the possibility of extracting rheological information directly from traction data.

Equations

$$E=\sum _{i=1}^N{E}_i,E_i=K_A(A_i-A_0{)}^2+K_P(P_i-P_0{)}^2,$$

$${\partial }_t{\mathbf{r}}_i=\mu {\mathbf{F}}_i+v_0{\hat{\mathbf{n}}}_i,{\partial }_t{\theta }_i=\sqrt{2D_r}{\eta }_i(t),$$

$${\sigma }_{\alpha \beta }^{(i)int}=-{\Pi }_i{\delta }_{\alpha \beta }+\frac{1}{2A_i}\sum _{ab\in i}{T}_{ab}^{\alpha }{l}_{ab}^{\beta },$$

$${\Pi }_i=-\frac{\partial E}{\partial A_i},{\mathbf{T}}_{ab}=\frac{\partial E}{\partial {\mathbf{l}}_{ab}},$$

$${\sigma }_{\alpha \beta }^i={\sigma }_{\alpha \beta }^{(i)int}+{\sigma }_{\alpha \beta }^{(i)swim},$$

$${\sigma }_{\alpha \beta }^{(i)swim}=-\frac{v_0}{\mu A_i}{n}_{\alpha }^i{r}_{\beta }^i$$

$${\sigma }_{\alpha \beta }^{int}=\frac{1}{A_T}\sum _i{A}_i{\sigma }_{\alpha \beta }^{(i)int},{\sigma }_{\alpha \beta }^{swim}=\frac{1}{A_T}\sum _i{A}_i{\sigma }_{\alpha \beta }^{(i)swim}.$$

$${\Pi }_i=-2K_A(A_i-A_0),T_{ab}=2K_P[(P_j-P_0)+(P_k-P_0)],$$

$${\mathbf{F}}_a(\{{\Pi }_i\},\{T_i\})={\mathbf{t}}_a,$$

$${\mathbf{F}}_{grid}(\{{\Pi }_i\},\{T_i\})={\mathbf{t}}_{grid}.$$

$${\sigma }_{n,d}=\frac{1}{2}({\sigma }_{xx}\pm {\sigma }_{yy}),{\sigma }_s=\frac{1}{2}({\sigma }_{xy}+{\sigma }_{yx}).$$

$${\sigma }_n^{int}=\frac{1}{A_T}\sum _i\left[2K_A{A}_i(A_i-A_0)+K_P{P}_i(P_i-P_0)\right].$$

$$C_{ss}(\tau )=⟨{\sigma }_s^{int}(t_0){\sigma }_s^{int}(t_0+\tau ){⟩}_{t_0},$$

$${\eta }_{eff}=\frac{A_T}{k_b{T}_{eff}}{\int }_0^{{\tau }_m}{C}_{ss}(\tau )d\tau ,$$

$${\mathbf{t}}_a={\mathbf{F}}_a,$$

$${\mathbf{F}}_a=-{\nabla }_{a}E=-\left(\frac{\partial E_i}{\partial {\mathbf{r}}_a}+\frac{\partial E_j}{\partial {\mathbf{r}}_a}+\frac{\partial E_k}{\partial {\mathbf{r}}_a}\right).$$

$$\frac{\partial E_i}{\partial {\mathbf{r}}_a}=-{\Pi }_i\frac{\partial A_i}{\partial {\mathbf{r}}_a}+T_i\frac{\partial P_i}{\partial {\mathbf{r}}_a},$$

$$T_i\equiv \frac{\partial E_i}{\partial P_i},{\Pi }_i\equiv -\frac{\partial E_i}{\partial A_i}.$$

$$\frac{\partial A_i}{\partial {\mathbf{r}}_a}=({\hat{\mathbf{n}}}_{ik}{l}_{ac}+{\hat{\mathbf{n}}}_{ij}{l}_{ab})/2,\frac{\partial P_i}{\partial {\mathbf{r}}_a}=-({\hat{\mathbf{l}}}_{ab}+{\hat{\mathbf{l}}}_{ac}),$$

$${\mathbf{t}}_{\mathbf{a}}=T_{ab}{\hat{\mathbf{l}}}_{ab}+T_{ac}{\hat{\mathbf{l}}}_{ac}+T_{ad}{\hat{\mathbf{l}}}_{ad}$$

$$+\frac{1}{2}{\Pi }_{ji}{\hat{\mathbf{n}}}_{ji}{l}_{ab}+\frac{1}{2}{\Pi }_{jk}{\hat{\mathbf{n}}}_{jk}{l}_{ad}+\frac{1}{2}{\Pi }_{ki}{\hat{\mathbf{n}}}_{ki}{l}_{ac},$$

$$T_{ab}=T_i+T_j{T}_{ac}=T_i+T_k{T}_{ad}=T_j+T_k$$

$${\Pi }_{ji}={\Pi }_j-{\Pi }_i{\Pi }_{jk}={\Pi }_j-{\Pi }_k{\Pi }_{ki}={\Pi }_k-{\Pi }_i$$

$${\mathbf{F}}_a(\{{\Pi }_i\},\{T_i\})={\mathbf{t}}_a.$$

$$\mathbf{M}\cdot \mathbf{\varphi }=\mathbf{t},$$

$$\mathbf{\varphi }={\mathbf{M}}^{-1}\cdot \mathbf{t},$$

$$Tr[(\mathbf{M}\cdot \mathbf{\varphi }-\mathbf{t}{)}^T\cdot (\mathbf{M}\cdot \mathbf{\varphi }-\mathbf{t})].$$

$${\mathbf{F}}_{grid}(\{{\Pi }_i\},\{T_i\})={\mathbf{t}}_{grid},$$

$${\mathbf{F}}_{grid}={\int }_{\partial S}\left[T_{ab}{\hat{\mathbf{n}}}_{ab}\delta (l-l_s)+{\Pi }_i{\hat{\mathbf{n}}}_s\right]dl,$$

$${\sigma }_{\alpha \beta }^{int}=-{\Pi }_T{\delta }_{\alpha \beta }+\frac{1}{A_T}\sum _{ab}{T}_{ab}^{\alpha }{l}_{ab}^{\beta },$$

$${\Pi }_T=\sum _i{\Pi }_i\frac{A_i}{A_T}.$$

$${\sigma }_{\alpha \beta }^{int}=\frac{1}{A_T}\sum _i{A}_i{\sigma }_{\alpha \beta }^{(i)int}.$$

Full text

It is now broadly recognized that the transmission of mechanical forces can be as important as genetics and biochemistry in regulating tissue organization in many developmental processes, including embryogenesis, morphogenesis, wound healing and cancer metastasis Lecuit and Lenne (2007); Heisenberg and Bellaiche (2013); Murrell et al. (2015); Etournay et al. (2015); Shraiman (2005); Serra-Picamal et al. (2012); Staple et al. (2010); Guillot and Lecuit (2013); Rauzzi et al. (2008); Trepat et al. (2009); Tambe et al. (2011). In order to make quantitative predictions for large scale cell remodeling in tissues, we must understand their material properties, such as stiffness and viscosity, as well as the forces that build up inside them, characterized by local pressures and stresses. Motivated by experiments highlighting the slow glassy dynamics of dense epithelial tissues, work by us and others has suggested that monolayers of motile cells may form glassy or jammed states, where local cell rearrangements are rare and energetically costly, and that a relatively small change of parameters may trigger a change from an elastic response to a state with fluid-like behavior, where individual cells continuously exchange neighbors Bi et al. (2014, 2015, 2016). This liquid-solid transition is tuned by the interplay of cell-cell adhesion and cortex contractility, as manifest in cellular shape, and by cell motility. This suggestion has been verified experimentally in specific cell types Park et al. (2015), indicating that the paradigm of tissues as active materials may be a useful way of organizing experimental data and classifying large-scale tissue behavior in terms of a few effective parameters.

In contrast, quantifying stresses and pressures in active, motile tissues is largely an open problem. One possible reason is that there are different definitions of stresses and pressures that arise naturally in different experiments and simulations, and it is not immediately clear how they are related to one another or under precisely which conditions each definition applies.

For example, Traction Force Microscopy (TFM) is a powerful tool that probes the dynamic forces cells exert onto soft substrates by measuring the substrate deformations Tambe et al. (2011); Soiné et al. (2014); Tambe et al. (2013); Zimmermann et al. (2014); Trepat et al. (2009). In some experiments, intercellular stresses are extracted from the traction forces using finite element analysis under the condition that traction forces are balanced by cellular interactions. The resulting stress maps reveal a highly dynamical and heterogeneous mechanical landscape, characterized by large spatial and temporal fluctuations in both normal and shear stresses Tambe et al. (2011, 2013); Kim et al. (2013); Trepat et al. (2009); Serra-Picamal et al. (2012). A key assumption used in the TFM approach to infer stresses from tractions is that the cell layer can be described as a continuum linear elastic material Tambe et al. (2013).

A second important set of mechanical inference methods also predict cellular stresses, but in static tissues. These methods rely on advances in imaging techniques that provide a spatially resolved view of tissue development during morphogenesis, with visualization of the cell boundaries of two dimensional cell sheets Ishihara and Sugimura (2012); Nier et al. (2016); Chiou et al. (2012); Noll et al. (2015); Brodland et al. (2010). Assuming mechanical equilibrium, one can then infer the tensions along cell edges and pressures within each cell from the cell configurations. This method provides a spatial distribution of intercellular stresses, but it does not capture temporal stress fluctuations arising from the nonequilibrium nature of the tissue. On the other hand, it directly couples mechanics to morphology and it has been used successfully to characterize cell morphology in the development of the Drosophila wing disk Farhadifar et al. (2007); Chiou et al. (2012), where cellular rearrangements are slow on the time scales of interest. This work has demonstrated that the analysis of cell shapes can provide fundamental insight on the mechanical state of tissues in developmental processes.

A third line of research has focused on the homeostatic pressure that tissues exert on their containers or surroundings. The homeostatic pressure has been proposed as a quantitative measure of the metastatic potential of a tumor Basan et al. (2009); Montel et al. (2011). It is defined as the force per unit area that a confined tissue would exert on a moving piston permeable to fluid, hence it represents an active osmotic cellular pressure. The existing literature has focused on the contribution to the homeostatic pressure from tissue growth due to cell division and death, but in general other active processes, such as cell motility and contractility, will also contribute to the forces exerted by living tissues on confining walls. A related body of work has investigated the pressures generated by motile particles, such as active colloids. In these highly nonequilibrium systems that break time-reversal symmetry, there is a contribution to the total pressure called the swim pressure generated entirely by persistent motility Yang et al. (2014); Takatori et al. (2014); Solon et al. (2015).

So far, there is no unifying theory for these seemingly distinct notions of pressure and stress, or their relationship to material properties. In this paper, we show that a recently proposed Self-Propelled Voronoi (SPV) model of epithelial tissue (Fig. 1) provides a natural framework for unifying these ideas. In the SPV model, cells are described via a polygonal tesselation of the plane. Each cell is endowed with motility with fixed speed and direction randomized by orientational noise and interacts with other cells via forces determined by a shape energy that incorporates the interplay between contractility and adhesion. In a broad class of vertex models Farhadifar et al. (2007) one can develop a self-consistent definition for cellular stress in terms of cell shapes, but this result is only useful in simulations where the energy functional that describes the tissue is precisely known. One of the benefits of the SPV model is that it additionally explicitly accounts for the forces that motile cells exert on the substrate. This allows us to develop a generalized mechanical inference method to infer cellular stresses from traction forces and to show that these match the stresses calculated from instantaneous cell shape, relating for the first time TFM data and mechanical inference techniques in motile tissues. Additionally, our method provides absolute values for junctional tensions and pressure differences. This is in contrast to equilibrium mechanical inference, which only yields relative forces Chiou et al. (2012); Ishihara and Sugimura (2012).

We also wish to know how this definition of cellular stress relates to others in the literature. There are two additive contributions to the mechanical stress that describe the forces transmitted in a material across a bulk plane. The first represents the flux of propulsive forces through a bulk plane carried by particles that move across it. The second describes the flux of interaction forces across a bulk plane. We demonstrate that the generalized mechanical inference measurements probe the latter, which we denote as interaction stresses. The former, which we denote as the tissue swim stress, approximates the contribution from cell motility to the osmotic pressure generated by cells immersed in a momentum conserving solvent on a semipermeable piston, hence to the tissue homeostatic pressure. The tensorial sum of the swim stress and the interaction stress is the total stress. The normal component of the total stress determines whether a tissue will tend to exert extensile or contractile forces on its environment, which is an important consideration in biological processes such as wound healing and cancer tumorogenesis.

An obvious open question, then, is how these stresses vary as a function of material properties. We find that the normal component of the interaction stress is contractile in both the solid and the liquid due to the contractility of the actomyosin cortex, although much more weakly so in the liquid state. In contrast, the normal component of the motility-induced swim stress is always extensile, corresponding to a positive swim pressure, though its magnitude depends on the phase: in a solid the swim pressure is negligible, while in the fluid it can be significant. This can result in a change in sign of the total mean stress: indeed, we find it is always contractile in the solid state but becomes extensile deep in the liquid state when cell motility exceeds actomyosin contractility.

Because the transition from contractile to extensile does not coincide with the fluid to solid transition, it is natural to ask whether the stress displays any signatures of the fluid-solid transition. We develop a definition for the effective viscosity of the tissue that can be extracted from the temporal correlation of the interaction shear stress, and find that it diverges as the tissue transits from the liquid state to the solid state. Importantly, this suggests that TFM combined with mechanical inference can provide rheological information about the tissue.

Results and Discussion

Self-Propelled Voronoi (SPV) Model.

The SPV model describes an epithelium as a network of polygons. Each cell $i$ is endowed a position vector $\bm{r}_{i}$, and cell shape is defined by the Voronoi tessellation of all cell positions (Fig.1). Like for vertex models, tissue forces are obtained from an effective energy functional $E(\{\bm{r}_{i}\})$ for $N$ cells, given by Farhadifar et al. (2007); Staple et al. (2010); Bi et al. (2016); Barton et al. (2016)

[TABLE]

with $A_{i}$ and $P_{i}$ the cross-sectional area and perimeter of the $i$-th cell, respectively. The first term in Eq. (1) arises from incompressibility of the layer in three dimensions and its resistance to height fluctuations, with $A_{0}$ 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 Farhadifar et al. (2007), with $P_{0}$ a preferred perimeter resulting from this competition. We simulate $N$ cells in a square box of area $A_{T}$, with $\bar{A}=A_{T}/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 $\bar{A}=A_{0}=1$ unless otherwise noted, and $K_{A}=K_{P}=1$.

Each Voronoi cell is additionally endowed with a constant self-propulsion speed $v_{0}$ along the direction of polarization $\hat{\bm{n}}_{i}=(\cos\theta_{i},\sin\theta_{i})$ describing cell motility. The dynamics of each Voronoi cell is governed by

[TABLE]

where $\bm{F}_{i}=-\bm{\nabla}_{i}E$ is the force on cell $i$ and $\mu$ is the mobility. The direction of cell polarization is randomized by orientational noise of rate $D_{r}$, with $\langle\eta_{i}(t)\rangle=0$ and $<\eta_{i}(t)\eta_{j}(t^{\prime})>=\delta_{ij}\delta(t-t^{\prime})$. The time scale $\tau_{r}=1/D_{r}$ 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 $D_{0}=v_{0}^{2}/(2D_{r})$ Henkes et al. (2011); Fily and Marchetti (2012); Yang et al. (2014). 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 Bi et al. (2014).

We showed in Ref. Bi et al. (2016) that the SPV model exhibits a transition from a solid-like state to a fluid-like state upon increasing the single-cell motility $v_{0}$, the persistence time $\tau_{r}$, or the cell shape parameter $P_{0}/\sqrt{A_{0}}$ that characterizes the competition between cell-cell adhesion and cortical tension. The phase diagram in the $(P_{0},v_{0})$ plane is reproduced in Fig.1B. The transition is identified by setting the shape index $q=\langle P_{i}/\sqrt{A_{i}}\rangle$ to the value $q=3.813$, where $\langle...\rangle$ denotes the average over all cells. It was shown in Ref. Bi et al. (2016) 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 $A_{T}$, the preferred cell area $A_{0}$ does not affect the interaction forces or cellular shapes. Hence the solid-fluid transition is insensitive to $A_{0}$, as demonstrated numerically and shown analytically in the SI. The preferred area $A_{0}$ 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 vertex model tissue described by the tissue energy Eq. (1), the mechanical state of cell $i$ is characterized by a local stress tensor $\sigma_{\alpha\beta}^{(i)int}$ given by Ishihara and Sugimura (2012); NESTOR-BERGMANN et al. (2016); Chiou et al. (2012)

[TABLE]

where $\Pi_{i}$ is the hydrostatic cellular pressure and $\mathbf{T}_{ab}=T_{ab}\mathbf{\hat{l}}_{ab}$ is the cell-edge tension, with $\mathbf{\hat{l}}_{ab}=\mathbf{l}_{ab}/|\mathbf{l}_{ab}|$. Here we use Roman indices $i,j,k,...$ to label cells and $a,b,c,...$ to label vertices, while Greek indices denote Cartesian components. The summation in Eq. (3a) runs over all edges of cell $i$, $\mathbf{l}_{ab}$ is the edge vector joining vertices $a$ and $b$ when the perimeter of cell $i$ is traversed clockwise, as shown in SI Fig.5A. 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 (see 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

[TABLE]

where, following recent work on active colloids Yang et al. (2014); Takatori et al. (2014),

[TABLE]

describes the flux of propulsive force $\frac{v_{0}}{\mu}\mathbf{n}_{i}$ across a boundary, as 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 $v_{0}$ and vanishes for nonmotile cells. The contribution $\sigma_{\alpha\beta}^{(i)int}$ is still given by Eq. (3a), but depends implicitly on cell motility through the instantaneous values of $\Pi_{i}$ and $T_{ab}$ 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 (see SI), $\sigma_{\alpha\beta}=\sigma_{\alpha\beta}^{int}+\sigma_{\alpha\beta}^{swim}$, with

[TABLE]

In simulations where the energy functional is known, it is simple to directly extract the instantaneous pressures and tensions from cell shapes in order to calculate the interaction contribution $\sigma_{\alpha\beta}^{(i)int}$. The definitions, Eq. (3b), give

[TABLE]

where $ab$ is the cell-cell interface that separates cells $j$ and $k$. Both $A_{i}$ and $P_{i}$ are obtained at every time step of the simulation and implicitly depend on cell motility, as parametrized by speed $v_{0}$ and persistence $\tau_{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 $\mathbf{t}_{a}=-\bm{\nabla}_{a}E$, which balances the interaction force $\mathbf{F}_{a}$. Equilibrium mechanical inference methods express the interaction force $\mathbf{F}_{a}=-\bm{\nabla}_{a}E$ at each vertex in terms of cellular pressures and edge tensions and the measured geometry of the cellular network (see Eqs. (19-21) of SI). Pressures and edge tensions are then obtained by inverting the equations $\mathbf{F}_{a}(\{\Pi_{i}\},\{T_{ab}\})=0$. For a nonequilibrium epithelial layer of motile cells we invert the force balance equations

[TABLE]

where the edge tensions $T_{ab}$ have been written as the sum of cortical tensions $T_{i}$ of the adjacent cells (see Eq. (20) of SI), reducing the number of independent unknowns. The interaction contribution to the local cellular stresses is then again 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 square minimization for the mechanical inference (See 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 utilizes experimentally accessible traction forces averaged over a square grid, with a grid spacing of the order of a cell diameter. Pressures and tensions $(\Pi_{i},T_{i})$ are then calculated by inverting the force balance equation at the center of each grid element,

[TABLE]

An outline of the coarse-graining procedure is given in the SI, with full details to be published elsewhere Yang et al. (2016). 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 Tambe et al. (2013). 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.

Within the framework of the SPV model, we have validated the coarse-grained method by comparing the resulting traction-based stresses to the shape-based stresses computed exactly from the simulations. The two show perfect agreement (Fig. 2C,F).

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 in the solid states. Thus mechanical measurement such as those provided by TFM can be used to characterize the rheological properties of the tissue.

The stress tensor $\sigma_{\alpha\beta}$ is symmetric and has three independent components in $2d$. Both the mean and local stress components are most usefully expressed in terms of normal stress $\sigma_{n}$, shear stress $\sigma_{s}$ and normal stress difference $\sigma_{d}$, with

[TABLE]

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, hence by interaction stresses. In contrast, the swim components of stress and pressure cannot be probed in TFM, but contribute to the pressure $\Pi=-\sigma_{n}=\Pi_{int}+\Pi_{swim}$ that the tissue would exert laterally on a confining piston. As we will see 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 (See SI Text)

[TABLE]

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 $\bar{A}=A_{0}$ and the normal stress comes mainly from the cortical tensions. Equation 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 become extensile deep in the liquid phase. We show in Fig. 2 snapshots of the local interaction normal stress in the solid state (left column) and in the liquid state (right column). In both the solid and the liquid states the interaction normal stress is on average contractile (red), with relatively weaker spatial fluctuations, but much larger mean value in the solid state, 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 SI Fig.6) across the solid-liquid transition. The color map shows that the total normal stress is contractile in the solid state and across the transition line denoted by the black crosses, but changes sign and becomes extensile deep in the liquid state. While the interaction stress is always positive due to cell contractility and consistent with experimental observations Lecuit and Lenne (2007); Bambardekar et al. (2015); Murrell et al. (2015); Heisenberg and Bellaiche (2013); Rauzzi et al. (2010); Kasza et al. (2014), 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 (see SI Fig.6), 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 Yang et al. (2014); Takatori et al. (2014). 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.

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,

[TABLE]

where $\langle...\rangle_{t_{0}}$ denotes the average over the length $t_{0}$ of the simulation, is shown in Fig. 4A for various $v_{0}$ 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, ultimately becoming frozen in the solid state. To quantify this we have defined the correlation time $\tau_{m}$ as the time when the correlation has decreased below one percent 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 $\eta_{\mathit{eff}}$ using a Green-Kubo-type relation by integrating the correlation function over the duration of the correlation time Hansen and McDonald (2005); Levashov et al. (2011),

[TABLE]

where we have used the ideal gas effective temperature $k_{B}T_{\mathit{eff}}=v_{0}^{2}/(2\mu D_{r})$ so that $\eta_{\mathit{eff}}$ has dimensions of a shear viscosity in $2d$. 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 $\tau_{r}$, however, the orientational noise in the SPV becomes uncorrelated in time and can be mapped onto thermal noise at an effective temperature $T_{\mathit{eff}}$. In this limit we expect $\eta_{\mathit{eff}}$ to indeed play the role of a shear viscosity. Remarkably, the effective viscosity shown in Fig.4(C,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.

Discussion and Conclusions

Using the Self-propelled Voronoi Model (SPV) we have formulated a unifying framework for quantifying the contributions from cell shape fluctuations and cell motility to mechanical stresses in an epithelial tissue. Cell shape fluctuations from actomyosin contractility and cell-cell adhesion control the interaction stress, while cell motility determines the swim stress that is generically present in all self-propelled systems Yang et al. (2014); Takatori et al. (2014).

Unlike monolayer stress microscopy that computes interaction stress from traction forces by assuming the tissue to be a continuum linear elastic material Tambe et al. (2013), the traction-based mechanical inference method developed here incorporates spatial and temporal deformations of the tissue due to actomyosin contractility and cell-cell adhesion, and can be generalized to account for cell division and apoptosis. In contrast to equilibrium mechanical inference techniques Chiou et al. (2012); Ishihara and Sugimura (2012), our approach does not require cells to be in or close to static mechanical balance, and it also provides the absolute scale of the junctional tensions and pressure differences. This can for instance be important for testing hypotheses involving mechanosensitive biomolecules. Experimentally, our method provides new ways to extract inter-cellular interaction stresses from existing traction force data and segmented cell images.

The swim stress, on the other hand, cannot be measured using TFM as it represents the flux of propulsive forces across a bulk plane in the tissue. It contributes to the homeostatic pressure at the lateral boundary of the tissue. The sum of the swim stress and the interaction stress approximates the total stress at the tissue boundary, which is generally contractile but can become extensile when the tissue is deep in the liquid state and cell motility exceeds actomyosin contractility. This change in sign may be observable in wound healing assays where the transition from contractile to extensile behavior can result in tissue expansion upon removal of confinement by neighboring tissue. The application of our method to the wound healing geometry will be described in future work.

We have extracted an effective tissue viscosity from the temporal correlation of the interaction shear stress. The correlation time and effective viscosity display a slowing down and arrest at the transition to the solid regime, thus serve as a direct probe of tissue rheology. Moreover, we observed a similar behavior for the temporal correlations of traction forces as demonstrated in the SI. Therefore, our work suggests that TFM measurement Tambe et al. (2011) combined with mechanical inference could provide information on tissue rheology. To our knowledge, this has not been attempted yet on experimental data.

Our work sets the stage for examining the feedback between cell activity and tissue mechanics that is apparent in many tissue-level phenomena. Recent work has shown that mechanical stresses influence cell proliferation in tumor spheroids Montel et al. (2011) and regulate cell growth in the developing Drosophila wing Pan et al. (2016). Regulation of cell motility, as in contact inhibition of locomotion, has been proposed to explain stress patterns during collective cell migration J et al. (2016). TFM has revealed the tendency of cells to move along the direction of minimal shear stress, a phenomenon termed “plithotaxis” Tambe et al. (2011). Our model provides a unifying framework for quantifying the relative roles of various cell properties, such as shape as controlled by contractility and cell-cell adhesion, motility and growth, on the mechanics of the tissue.

We thank Daniel Sussman for valuable discussions. This work was supported by the Simons Foundation through a Targeted Grant Award No. 342354 (MCM & MC) and an Investigator Award No. 446222 (MLM, MC, & MM) in the Mathematical Modeling of Living Systems, by the National Science Foundation through awards DMR-1305184 (MCM & XY), DMR- 1609208 (MCM), DMR-1352184 (MLM & DB), and the IGERT grant DGE-1068780 (MCM & MC), the NIH through R01GM117598-02 (MLM), and by the computational resources provided by Syracuse University and through NSF ACI-1541396. All authors acknowledge support from the Syracuse University Soft Matter Program.

Supporting Information (SI)

Traction-based Mechanical Inference

Force balance in the vertex model

Force balance at each vertex requires

[TABLE]

where $\mathbf{t_{a}}$ is the traction force vertex $a$ exerts on the substrate, which is balanced by the interaction force $\bm{F}_{a}=-\nabla_{a}E$. We consider for concreteness four adjacent vertices $a,b,c,d$ and corresponding cells $i,j,k$, as shown in Fig.5A. Using the additivity of the shape energy, we write

[TABLE]

The right hand side of Eq. (15) can be decomposed in terms of tensions and pressure jumps along and perpendicular to the cell edges, as illustrated in Fig.5A, by writing

[TABLE]

where $T_{i}$ and $\Pi_{i}$ are the cortical tension and hydrostatic pressure for cell $i$, respectively, defined as Chiou et al. (2012)

[TABLE]

The geometric factors in Eq. 16 can be written as

[TABLE]

where $\hat{\bm{n}}_{ik}$ is the unit vector perpendicular to edge $ac$ pointing from cell $i$ to cell $k$, $l_{ac}$ is the length of edge $ac$, and $\hat{\bm{l}}_{ad}$ is the unit vector along edge $ad$ pointing from vertex $a$ to vertex $d$. Combining Eq. (14)-Eq. (18), we obtain

[TABLE]

where

[TABLE]

are edge tensions and pressure jumps along and across cell edges, respectively. Eq. (19)-Eq. (21) form the basis for the traction-based mechanical inference and can be written in a general form as

[TABLE]

Operationally, we write the force balance equation in a matrix form

[TABLE]

where $\bm{M}$ is a $4N\times 2N$ structural matrix incorporating the orientation and length of the edge vectors, $\bm{\phi}$ is a $2N\times 1$ column vector whose elements are the cortical tensions and pressures, and $\bm{t}$ is the $4N\times 1$ column vector of the traction forces at each vertex. The tensions and pressures are obtained by inverting the structural matrix

[TABLE]

which is equivalent to a least square minimization of

[TABLE]

In practice, this is done with the Moore-Penrose pseudoinverse algorithm Chiou et al. (2012).

Force balance on the coarse-grained grid

In real tissues, the traction force acts not only at cell vertices, but also at focal adhesions and throughout the basal area of a cell, and can be measured using TFM on a coarse-grained grid. Therefore, we generalize the traction-based mechanical inference by coarse-graining the tissue on a Cartesian grid and impose force balance on each grid element (Fig.5B)

[TABLE]

where $\mathbf{F}_{grid}$ is the interaction force on each grid element as calculated from tensions and pressures as shown in Fig.5B, and $\mathbf{t}_{grid}$ is the traction force the grid element exerts on the substrate. A grid element is a closed region of the tissue. In a Cartesian grid, it is a square element as highlighted in Fig.5B. The interaction force acts on the grid element through the boundary and is generally expressed as a surface integral, or a line integral in 2D, as

[TABLE]

where tension is a point force concentrated at the intersection $l_{s}$ between the cell edge and the grid element boundary, while pressure induces a normal force in the normal direction $\hat{\bm{n}}_{s}$ of the boundary. Notice that the pressure varies as the boundary crosses different cells. Using cell segmentation data together with Eq. (27) and Eq. (20), the interaction force can be written explicitly in terms of cortical tensions and cell pressures. The traction force the grid element exerts on the substrate can be measured using TFM. By inverting the force balance Eq. (26) using the Moore-Penrose pseudoinverse algorithm, one can obtain tensions and pressures, from which the local stress tensor can be constructed using Eq. (3a) in the main text.

To test the coarse-grained mechanical inference, we generate cell network configurations using the SPV model. The traction force of the grid element is calculated as the sum of tensions (normal forces) on the vertices (cell edges) within the grid element. It is numerically verified that the traction force is exactly balanced by the interaction force for each grid element. We have tested the coarse-grained mechanical inference in our numerical simulation with a grid element size of the order of a single cell size. The inferred stress agrees perfectly with the accurate value obtained from the cell shape using the energy functional of the tissue (Fig.2 in the main text).

This coarse-grained mechanical inference requires knowledge of cell shapes as well as traction forces. The cell shape allows one to express the interaction force in terms of tension and pressure and resolves the stress down to a single cell level. The resolution of the coarse-graining can be adjusted to adapt to the resolution of the traction force microscopy, where a higher TFM resolution will yield a more accurate result for the inferred tensions and pressure differences by over-constraining the system. A systematic study of the method will be published elsewhere Yang et al. (2016).

The traction-based mechanical inference does not require knowledge of the energy functional or the material property of the tissue, and is based on the assumption of force balance and the decomposition of interaction forces into junctional tensions and cell pressures. We expect our method to have general applicability to epithelial tissues where cell shape and traction forces can be measured simultaneously.

Tissue Stress

From cellular tensions and pressures one can also evaluate the mean interaction stress for the whole tissue, given by

[TABLE]

where the summation is over all edges, $A_{T}$ is the area of the tissue, and the pressure conjugated to the tissue area $A_{T}$ is given by

[TABLE]

Combining Eq. (3a) in the main text with Eq. (28)-Eq. (29), we show that the interaction tissue stress is related to the cellular stress as

[TABLE]

Next we show that for the energy functional given by Eq. (1) of the main text, the tissue interaction normal stress can be expressed in terms of cell areas and perimeters. We first focus on the interaction normal stress of cell $i$ defined as

[TABLE]

Using Eq. (3a) in the main text, we rewrite Eq. (31) as

[TABLE]

where we used the relation in Eq. (20) and the summation is over the neighbors of cell $i$. $l_{ij}$ is the length of the edge shared by cell $i$ and cell $j$. Using the energy functional (Eq. (1) in the main text) and the definitions for tension and pressure (Eq. (17)), we have

[TABLE]

Substituting Eq. (33) into Eq. (30), we obtain the interaction normal stress for the tissue in terms of cell areas and perimeters

[TABLE]

where we used the identity

[TABLE]

Eq. (34) expresses the interaction tissue stress in terms of cell shape parameters, providing a prescription for extracting mechanical information directly from cell images. We plot the mean normal interaction stress (Eq. (34)) and the mean normal swim stress (Eq. (6)) as a function of $v_{0}$ and $P_{0}$ in Fig. 6. Notice that the interaction stress is contractile and the swim stress is extensile. The latter vanishes in the solid state and dominates deep in the liquid state, resulting in a change of sign of the total normal stress as shown in Fig. 3 in the main text.

Effect of $A_{0}$ on cell dynamics

For a system with periodic boundary conditions and fixed size, the parameter $A_{0}$ plays no essential role in the cell dynamics. It only renormalizes the pressure of the system by a constant, but does not affect cell shapes, forces on vertices, or the anisotropic part of the stress tensor. To see this, we rewrite the energy of the system using the average cell area $\bar{A}$. Up to a constant offset, the total energy reads:

[TABLE]

Note that the parameter $A_{0}$ only appears in the last term of Eq. (36), which merely offsets the pressure of the system by $2K_{A}A_{0}$. Thus, all forces $\bm{F}_{i}$ and shear stresses are independent of $A_{0}$. The liquid-solid phase transition is unaffected by $A_{0}$ as shown in Fig.7.

Effect of Voronoi constraint on the tissue stress

The SPV model can be mapped to an active vertex model (AVM) with effective constraints imposed by the Voronoi tessellation on the vertex positions. In a vertex model free of Voronoi constraint, the interaction stress tensor $\sigma_{\alpha\beta}^{int}$ as defined by Eq. (6) in the main text represents precisely the stress acting on the periodic boundary. This is not, however, the case for the SPV model, because the creation of cell shapes using a Voronoi tessellation adds constraints to a vertex model. In this section, we show that the results presented in the paper are robust against the Voronoi constraint.

The normal stress as defined in Eq. (10) in the main text is unaffected by the Voronoi constraint, but the shear stress does not correspond to the shear stress at the boundary of the system. This is because an affine shear deformation of all cell positions does not induce an affine shear transformation of all vertex positions. Indeed, the $2N$ vertices would in total provide $4N$ degrees of freedom, but the Voronoi tessellation reduces this to the $2N$ degrees of freedom of the $N$ cell positions.

In the SPV model, the stress experienced by the boundary can be computed by quantifying the total energy change induced by a deformation of the periodic box described by the strain tensor $G_{\alpha\beta}$ while affinely displacing all cell positions and respecting the Voronoi tessellation:

[TABLE]

where the tilde is used to denote the stress at the boundary under Voronoi constraint. The shear component of this tensor is different from that of the stress $\sigma_{\alpha\beta}^{int}$ defined in Eq. (6) in the main text.

To numerically study how different $\tilde{\sigma}_{\alpha\beta}^{int}$ and $\sigma_{\alpha\beta}^{int}$ typically are, we impose a sufficiently small strain on the cell positions and calculate the tissue energy difference before and after the imposed strain. To incorporate the Voronoi constraint, the cell network is constructed based on the cell positions before and after the imposed strain according to the Voronoi tessellation. An expansion along $x$ and $y$ axis is implemented to compute the normal components of the stress $\tilde{\sigma}_{xx}^{int}$ and $\tilde{\sigma}_{yy}^{int}$, using

[TABLE]

respectively. Here, $\bm{r}$ is the cell position and $\epsilon\ll 1$ is the imposed strain. Similarly, simple shear deformations along $x$ and $y$ axes are implemented to compute the off-diagonal components of the stresses $\tilde{\sigma}_{xy}^{int}$ and $\tilde{\sigma}_{yx}^{int}$, using the transformations

[TABLE]

respectively.

We present the results of the stress at the boundary under Voronoi constraint in Fig.8. Comparing with the results in the main text, the effect of the Voronoi tessellation does not affect the behavior of the normal and the shear stress as presented in the main text. This suggests that for practical purposes, one can apply the definitions of stress of the vertex model to approximate the stress in the SPV model. The insensitivity of the behavior of the stress to the geometrical constraint also implies general applicability of the SPV model in describing the mechanics of confluent epithelial tissue in experiment.

Autocorrelation of traction forces

The slowing down of the autocorrelation of the shear stress implies a similar behavior for the autocorrelation of traction forces due to the force balance condition $t_{\alpha}=\partial_{\beta}\sigma_{\alpha\beta}^{int}$. We have shown this explicitly by evaluating the autocorrelation function of the vertex tractions averaged over a small patch of tissue as a function of $v_{0}$ and $P_{0}$, as shown in Fig. 9. In contrast to the shear stress, we cannot average the traction over the entire system, as due to force balance the net traction force is zero. Instead we compute the two-time correlation function of the traction forces averaged over a circular patch of tissue of radius $R$. Such a mean local traction is calculated in two reference frames: the Eulerian frame where the circular patch of radius $R=0.1L$ is fixed in space (Fig. 9A) and the Lagrangian reference frame where cells initially within the patch are tracked in time (Fig. 9D). To reduce the finite size fluctuations, we compute the autocorrelation using $7$ patches at different positions and average over the patches. Experimental traction forces are usually measured in an Eulerian reference frame by TFM. Such measurements could, however, be transformed to a Lagrangian reference frame if cell tracking data are available.

Figure 9 B and E show the autocorrelation of the traction, $C_{tt}(\tau)=\langle t_{x}(t_{0})t_{x}(t_{0}+\tau)\rangle_{t_{0}}$, with $t_{x}$ the $x$ component of the the average traction in the Eulerian and Lagrangian frame, respectively. When calculated in the Lagrangian frame, the correlations show the same dramatic slow down when the solid is approached from the liquid side as found in the autocorrelation of the mean shear stress. In the Eulerian frame, however, a second time scale associated with cells moving in and out of the fixed patch masks the slowing down of the relaxation at long times. The same behavior is obtained when calculating correlations of shear stress averaged over a local patch instead of over the entire system. The slowing down of the relaxation is quantified in Fig. 9 C and F that show heat maps of the correlation time $\tau_{m}$ defined as $C_{tt}(\tau_{m})/C_{tt}(0)=0.02$ in the Eulerian and Lagrangian frames, respectively. The correlation time increases as the system transits from the fluid state to the solid state in both reference frames and diverges deep in the solid state. The behavior is less clear in the Eulerian reference frame, where cell transition across the patch introduces a new time scale that masks the divergence. In summary, this confirms that the divergence of traction and shear relaxation can in principle be seen in TFM data, especially if combined with cell tracking so that correlations can be computed in a Lagrangian frame. Our results suggest a new way of analyzing TFM data not previously attempted to provide rheological information about the tissue.