Mechanical heterogeneity in tissues promotes rigidity and controls cellular invasion

TL;DR

This study reveals how mechanical heterogeneity in tissues influences their rigidity and cellular invasion, showing that a small fraction of rigid cells can significantly enhance tissue stiffness and affect tumor invasion dynamics.

Contribution

It introduces a vertex-based model to analyze cell stiffness heterogeneity and uncovers a scaling law for tissue rigidity controlled by the fraction of rigid cells.

Findings

Tissue shear modulus scales with the fraction of rigid cells.

Rigidity onset occurs at a low fraction of rigid cells (0.21).

Heterogeneity can create an intermediate solid state that impedes invasion.

Abstract

We study the influence of cell-level mechanical heterogeneity in epithelial tissues using a vertex-based model. Heterogeneity in single cell stiffness is introduced as a quenched random variable in the preferred shape index($p_0$) for each cell. We uncovered a crossover scaling for the tissue shear modulus, suggesting that tissue collective rigidity is controlled by a single parameter $f_r$, which accounts for the fraction of rigid cells. Interestingly, the rigidity onset occurs at $f_r=0.21$, far below the contact percolation threshold of rigid cells. Due to the separation of rigidity and contact percolations, heterogeneity can enhance tissue rigidity and gives rise to an intermediate solid state. The influence of heterogeneity on tumor invasion dynamics is also investigated. There is an overall impedance of invasion as the tissue becomes more rigid. Invasion can also occur in the intermediate heterogeneous solid state and is characterized by significant spatial-temporal intermittency.

Equations

$$\epsilon =\sum _{i=1}^N\left[{\kappa }_A(a_i-1{)}^2+(p_i-p_0^i{)}^2\right].$$

$${\tau }_m=(p_i-p_0^i)+(p_j-p_0^j),$$

$$G=\sigma g_{\pm }\left(\frac{\sigma }{|\mu -{\mu }^{\ast }|}\right).$$

$$f_r={\int }_{-\infty }^{{\mu }^{\ast }}{\mathcal{F}}_{\mu ,\sigma }(p_0)d{p}_0,$$

$$\frac{d{\mathbf{r}}_{\nu }}{dt}=\{\begin{array}{ll}-\frac{\partial ϵ}{\partial {\mathbf{r}}_{\nu }}+v_0\mathbf{n}& \text{vertices of invading cell},\\ -\frac{\partial ϵ}{\partial {\mathbf{r}}_{\nu }}& \text{for other vertices.}\end{array}$$

Full text

Heterogeneity amongst cells in a tumor has been recognized as one of the hallmarks of cancerFidler (1978); Marusyk and Polyak (2010); Hanahan and Weinberg (2011); Gay et al. (2016). This so-called intra-tumor heterogeneity is thought to facilitate metastasisFidler (1978); Brock et al. (2009); Shin et al. (2014) by allowing a cellular community the flexibility and efficiency to adapt to new environmentsAckermann (2015). It is also largely responsible for therapeutic resistance, leading to the failure to predict the progression patterns of metastasisFidler (1978); Marusyk and Polyak (2010); Bhatia et al. (2012); Claudi et al. (2014). Cellular differences within a tumor result from an interplay of both genetic and extrinsic influencesZhang and Austin (2012); Anderson et al. (2010); Snuderl et al. (2011); Notta et al. (2011); Marusyk et al. (2012). Whereas fitness or genotypic heterogeneity has been studied extensively, the role of mechanical heterogeneity in a tumor or cellular collective is still not well understood.

Recent experimental evidenceTian et al. (2015); Yallapu et al. (2015); Fritsch et al. (2010); Morawetz et al. (2017) have revealed that tumor cells exhibit a broad distribution of various biomechanical properties. These include intra-tumor heterogeneity in cell stiffnessesPlodinec et al. (2012); Fujii et al. (2019); Guo et al. (2014); Fritsch et al. (2010); Ciasca et al. (2016), stresses, and cell-cell interactionsGuo et al. (2014); Lee et al. (2012). However, there is no consensus on how these heterogeneities affect the mechanical behavior at the tissue level. For example, while biophysical techniquesAlibert et al. (2017) such as AFMCiasca et al. (2016) and optical stretcherMorawetz et al. (2017) show that individual cancer cells are softer than healthy cells, there is an apparent paradox with measurements at the tissue level, which show that tumors are more rigid than healthy tissuesButcher et al. (2009); Sinkus et al. (2000).

Similarly, from a theoretical and modeling perspective, there is a lack of understanding of how heterogeneity exhibited by single cells influence the mechanical state of the tissue. One large class of vertex-based computational modelsNagai and Honda (2001); Farhadifar et al. (2007); Fletcher et al. (2014) has been used to study mechanics of homogeneous epithelial tissues which can exhibit a solid to fluid transitionStaple et al. (2010); Bi et al. (2015); Kim and Hilgenfeldt (2015). However, previous works often treat heterogeneity as a biological noise and therefore unable to model many of the salient features in tumors.

In this work, we explicitly study the effect of heterogeneity on the mechanics of the confluent cell monolayer using the vertex model . Our results show that heterogeneity always acts to enhance the rigidity of a tissue. The mode of enhancement by heterogeneity depends on the spatial distribution of rigid vs. soft cells, which is directly tuned by a single universal parameter $f_{r}$. This results in distinct mechanical regimes that arise from the mismatch of two percolation processes: (1) rigidity percolation of mechanical tensions in a tissue and (2) contact percolation of rigid cells. We also connect tissue rigidity to cellular migration using model for the invasion of a single metastatic cell in a heterotypic microenvironment of the tissue.

Model We focus on the mechanics of a 2D confluent epithelium. Similar to previous vertex model implementations, the tissue is governed by the energy functionNagai and Honda (2001); Farhadifar et al. (2007); Staple et al. (2010); Li and Sun (2015); Fletcher et al. (2014); Kim and Hilgenfeldt (2015) $E=\sum_{i=1}^{N}\left[K_{A}(A_{i}-A_{0}^{i})^{2}+K_{P}(P_{i}-P_{0}^{i})^{2}\right]$, where cell areas $\{A_{i}\}$ and perimeters $\{P_{i}\}$ are functions of the position of vertices $\{{\bf r}_{i}\}$. $K_{A}$,$K_{P}$ are the area and perimeter elasticities. The term quadratic in cell area $A_{i}$ results from resistance to volume changesFarhadifar et al. (2007); Staple et al. (2010); Fletcher et al. (2014). Changes to cell perimeters are directly related to the deformation of the acto-myosin cortexHufnagel et al. (2007); Zehnder et al. (2015). The term $K_{P}P_{i}^{2}$ corresponds to the energy cost of deforming the cortex. The linear term, $-2K_{P}P_{0}P_{i}$, represents the effective line tension by cell-$i$ which gives rise to a ‘preferred perimeter’ $P_{0}$. The value of $P_{0}$ can be modified by an interplay of cell-cell adhesion and cortical tensionFarhadifar et al. (2007); Staple et al. (2010); Bi et al. (2015). Here we assume the individual preferred cell area $A_{0}$ does not vary from cell-to-cell and is set to the average area per cell (i.e. $A_{0}=\bar{A}$). The tissue energy can be non-dimensionalized by choosing $K_{P}\bar{A}$ as the energy unit and $\sqrt{\bar{A}}$ as the length unit

[TABLE]

where $a_{i}=A_{i}/\bar{A}$ and $p_{i}=P_{i}/\sqrt{\bar{A}}$ are the rescaled area and perimeter of the $i^{th}$ cell. $\kappa_{A}=K_{A}\bar{A}/K_{P}$ and $p^{i}_{0}=P^{i}_{0}/\sqrt{\bar{A}}$ is the preferred cell shape indexBi et al. (2015).

In this model, the cell stiffness is determined by the tension $\tau_{m}$ on cell-cell junctions (edges)Hutson et al. (2003); Rauzi et al. (2008); Brodland et al. (2010); Chiou et al. (2012); Ishihara and Sugimura (2012); Ishihara et al. (2013); Brodland et al. (2014), which in turn is directly tuned by the preferred cell shape index $p_{0}$Yan and Bi (2019); Yang et al. (2017),

[TABLE]

where $p_{i}$ and $p_{j}$ are the actual perimeters of cells $i,j$ adjacent to the edge $m$. To capture the experimentally observed variations in single cell stiffnesses and in cell-cell interactionsPlodinec et al. (2012); Fujii et al. (2019); Guo et al. (2014); Fritsch et al. (2010); Ciasca et al. (2016), we introduce variations in the preferred shape indices $p_{0}$. The majority of this work uses a Gaussian distributed set of $p_{0}$’s with mean $\mu$ and standard deviation $\sigma$. We also performed analysis using uniform distributions with the same ($\mu,\sigma$)(SI Section I).

To initialize the simulation, Voronoi cellsBi et al. (2016) are used to provide a set of initial vertex positions, $\{\bm{r}_{i}\}$. Then each cell is assigned a value of $p^{i}_{0}$ drawn randomly from a Gaussian distribution with mean $\mu$ and standard deviation $\sigma$. The set of $\{p_{0}^{i}\}$ then remains as quenched variables throughout the simulation. A combination of FIREBitzek et al. (2006) and conjugate-gradientByrd et al. (1995) algorithms is used to minimize Eq. (1). The network topology is updated using T1-movesWeaire and Hutzler (2001); Bi et al. (2014, 2015) during minimization. This combined algorithm produces final stable states where the net residual force on vertices is less than $10^{-10}$. A wide range of parameters $\mu=3.6-3.95$, $\sigma=0-0.3$ and $\kappa_{A}$ as well a range of system sizes $N=36-900$ are simulated.

Cellular heterogeneity enhances rigidity- We first focus on the case where area and perimeter elasticities are matched in their strengths ($\kappa_{A}=1$). In Fig. 1(a), the shear modulus $G$ is plotted as function of $\mu$ and $\sigma$. At $\sigma=0$, $G$ vanishes at $\mu\approx 3.828$. This recapitulates the rigidity transition in the absence of heterogeneity Farhadifar et al. (2007); Bi et al. (2015); Yan and Bi (2019). Increasing $\sigma$ at fixed $\mu$ always enhances the rigidity, which can occur in two different ways: (i) for a rigid tissue when $\sigma=0$, $G$ always increases with $\sigma$, and (ii) for a fluid tissue, $G=0$ when $0<\sigma\leq\sigma^{*}$ but becomes rigid when $\sigma\geq\sigma^{*}$. The threshold $\sigma^{*}$ where rigidity emerges depends on $\mu$. To better understand the dependence of $G$ on $(\mu,\sigma)$, we first hypothesized that $\sigma$ provides an overall scale for the shear modulus and plot a rescaled $G/\sigma$ as function of $\mu$ for various $\sigma$ (Fig. 1(a, inset)). This yields an intriguing result, where all curves in Fig. 1(a, inset) intersect at a common point $\mu^{*}\approx 3.8$. The point $\mu^{*}$ serves as a cross-over point between two distinct regimes. When $\mu<\mu^{*}$, $G/\sigma$ decreases as $\sigma$ is increased. For $\mu\geq\mu^{*}$ the behavior is flipped and $G/\sigma$ increases with $\sigma$. The intersection of all curves also suggests that closer to $\mu^{*}$, the shear modulus should converge to a scaling of $G\propto\sigma$. Based on these observations, we hypothesized that the behavior of $G$ below and above $\mu^{*}$ may be described by a universal scaling relation

[TABLE]

We re-plot all data based on the ansatz (Eq. 3) in Fig. 1(b). An excellent scaling collapse is obtained, which also allowed us to pinpoint the location of the crossover transition $\mu^{*}=3.812$. The two distinct branches ($+$, $-$) of the scaling function $g$ correspond to $\mu>\mu^{*}$ and $\mu<\mu^{*}$, respectively. The two branches meet at $\mu=\mu^{*}$, where $G$ scales linearly with $\sigma$.

What does this universal scaling in the vicinity of a cross-over transition point $\mu^{*}$ inform us about the nature of rigidity? At the crossover $\mu=\mu^{*}$ , there are $50\%$ of cells with $p_{0}^{i}<\mu^{*}$ regardless of the value of $\sigma$. Below the crossover ($\mu<\mu^{*}$), cells with $p_{0}^{i}<\mu^{*}$ exceed $50\%$, while above the crossover ($\mu>\mu^{*}$) this fraction is below $50\%$. This suggests that the fraction of cells with $p_{0}^{i}<\mu^{*}$ plays an important role. Therefore we define this fraction cells as

[TABLE]

where $\mathcal{F}$ is the probability distribution function and for a Gaussian, $f_{r}=(1/2)\text{erfc}[(\mu-\mu^{*})/(\sqrt{2}\sigma)]$. In Fig. 1(b-inset), a re-plot all data in terms of $G/\sigma$ vs. $f_{r}$ collapses to a common curve. Interestingly, the numerical value $p_{0}=3.812$ was previously found to be the threshold for rigidity in tissuesBi et al. (2015) with homogeneous $p_{0}$. This universal scaling behavior supports the idea that $p_{0}^{i}<3.812$ can serve as a single-cell criterion for rigidity and their fraction $f_{r}$ determines tissue mechanics. We therefore refer to $f_{r}$ as the fraction of rigid cells.

Heterogeneity drives two different percolations- The emergence of rigidity when a certain population of cells exceeds a fraction is reminiscent of a percolation processFeng and Sen (1984); Jacobs and Thorpe (1995, 1996); Tang and Thorpe (1988); Ostojic et al. (2006). To test if $f_{r}$ drives a percolation transition, we first analyze the spatial organization of rigid cells. Snapshots of connected rigid cells are shown in Fig. 2(a) and the probability for rigid cells to form a system-spanning contact network is plotted in Fig. 2(c). Finite-size scaling (Fig. 2(c, inset)) shows a contact percolation transition at $f_{r}^{c}=0.48$ with exponent $\nu=1.62\pm 0.01$. However, according to Fig. 1(b,inset) and Fig. 3(b), the tissue becomes rigid at a much lower $f_{r}\approx 0.21$, suggesting that contact percolation of rigid cells is not necessary for rigidity. We then hypothesized that the percolation of the tension network is the necessary condition for rigidity rather than rigid cell contacts. Fig. 2(d) shows tension percolation occurring at $f_{r}^{*}=0.21\pm 0.01$, coinciding with the rigidity onset. Snapshots of tensions (Fig. 2(b)) highlight an interesting behavior: e.g. at $f_{r}=0.24$, the rigid-cell contact network does not span the system, but the tensions form a system-spanning structure. Therefore the presence of just a few unconnected rigid cells can induce a much more spatially extended tension network. The distribution of $\tau$ also evolves as function of $f_{r}$ in an unconventional way (Fig. 3(a)). For values of $f_{r}>0.21$ corresponding to rigid states, $p(\tau)$ vanishes as $\tau\to 0$ and appears symmetrical about their mean. At lower $f_{r}$, $p(\tau)$ becomes significantly more skewed and heavy-tailed (SI FIG.S2) and also develops an excess number of $\tau\approx 0$ edges. Interestingly, the average tension for the tissue $\langle\tau\rangle$ does not vanish (Fig. 3(b)) even below $f_{r}^{*}$, and does not show drastic changes at $f_{r}^{*}$. In contrast, $G/\sigma$ undergoes a transition $f_{r}^{*}$, coinciding with the tension network percolation. Remarkably, two states with slight differences in their tensions can differ by several orders of magnitudes in $G$(Fig.S3), depending on whether the tension network is percolated.

Next we use a simple meanfield approach to estimate the rigidity transition threshold $f_{r}^{*}=0.21$. Since our simulated tissues share a similar connectivity and topology with random Voronoi tessellations, we recall the bond percolation threshold ($\sim 0.66$) on such networksShante and Kirkpatrick (1971); Hsu and Huang (1999); Becker and Ziff (2009). Therefore if rigidity requires the percolation of edges with finite tensions, the onset of rigidity should correspond to when there are exactly $66\%$ cell edges with finite tensions. For each parameter set $(\mu,\sigma)$ we define the fraction of finite tensions to be $N_{\tau}=\int_{\tau_{0}}^{\infty}p(\tau)d\tau$, where the threshold $\tau_{0}=10^{-5}$ is chosen to define “zero” tensions (coinciding with the noise floor in the numerical calculations). In Fig. 3(c) we show that the point at which $N_{\tau}(\mu,\sigma)=0.66$ closely coincides with the rigidity transition, i.e. $f_{r}^{*}=0.21$. Taken together, these results suggest that a rigidity percolation occurs at $f_{r}^{*}$.

The mechanics of the tissue is summarized in a phase diagram (Fig. 3(c)). Solid and fluid-like states are separated by $f_{r}^{*}$ which corresponds to a rigidity percolation of the edge tensions in the network. Contact percolation of rigid cells occur deeper within the solid phase at $f_{r}=0.48$. The mismatch between contact and rigidity percolations can be observed only in the presence of heterogeneity since at $\sigma=0$, $f_{r}$ is either [math] or $1$. We call the intermediate region between the two percolation transitions the Heterogeneous solid state. Rigidity transition obtained through the meanfield method is shown as the red line in Fig. 3(c).

Cell area elasticity $\kappa_{A}$ weakens the effect of heterogeneity- Next we focus on the Heterogeneous solid state ($0.21<f_{r}<0.48)$. The regime is characterized by spatial heterogeneity (Fig. 2(b)) where isolated cells or isolated clusters exist along side tissue-spanning percolated networks. Mechanical force balance must hold at every vertexYan and Bi (2019) for a solid tissue. However, for an isolated tension cluster, edge tensions alone cannot guarantee force balance, such as a vertex that sits on the boundary of this cluster. The missing component required for full force balance is the intracellular pressure force that acts perpendicularly on an edge given byYang et al. (2017); Yan and Bi (2019) $\Pi_{m}=\kappa_{A}(a_{i}-a_{j}).$ Here $a_{i}$ and $a_{j}$ are the areas of cell $i$ and $j$ adjacent to the edge $m$. Interestingly, since the pressure depends on the value of the cell area elasticity $\kappa_{A}$, the stability of isolated tension clusters must also depend on $\kappa_{A}$. To see the effect of cell area elasticity, we repeated the numerical calculations at various $\kappa_{A}$ values. The dependence of $G/\sigma$ on $f_{r}$ is plotted for various $\kappa_{A}$ in Fig. 4(a). At $\kappa_{A}>0$, $G/\sigma$ suffers a slight dip at $f_{r}\sim 0.48$, but stays finite all the way until $f_{r}$ below $f_{r}^{*}$. In general, as long as $\kappa_{A}$ is finite, the rigidity transition always occurs at $f_{r}^{*}$ and does not shift. The case of $\kappa_{A}=0$ is a singular limit where the rigidity transition shifts suddenly to $f_{r}=0.48$ to coincide with the contact percolation transition. This confirms our hypothesis that when cells cannot exert pressure forces, the tension network can only support mechanical rigidity when rigid cells physically come in contact. In contrast, at $\kappa_{A}>0$, stable tensions can be induced between rigid cells that are separated by a distance. We summarize these results with in Fig. 4(b) and also incorporate the different scaling regimes for the solid phase when $\kappa_{A}>0$. We are able to differentiate two types of solids: (i)when $\kappa_{A}>0,\ 0.21<f_{r}<0.48$, rigidity is strongly enhanced by heterogeneity (lower-branch in Fig. 1(b)) and (ii) when $\kappa_{A}>0,\ f_{r}>0.48$, there is only weak enhancement of rigidity (upper-branch in Fig. 1(b)). When $\kappa_{A}=0$, heterogeneity has no effect.

Heterogeneity and cellular invasion - Having established the static mechanical properties, we next focus on the effect of the heterotypic microenvironment on cell migration. We use a dynamic vertex modelStaddon et al. (2018); Das et al. (2019) to simulate the motion of a single invading cell in the tissue. The invading cell has a propulsive force $v_{0}$ along a polarity vector $\bf{n}$, which undergoes random rotational diffusionFily and Marchetti (2012); Bi et al. (2016) at a slow rate. This mimics the directional motility of a metastatic cell under the influence of strong chemotactic signalsClark and Vignjevic (2015). In the tissue each vertex $v$ evolves according to the overdamped equation of motion

[TABLE]

Eq. 5 is simulated with a fixed value $v_{0}=0.4$ and at various values of $f_{r}$. The mean total displacement of the invading cell is plotted as function of $f_{r}$ in Fig. 5a, the behavior near $f_{r}=0$ and $f_{r}=1$ recapitulate results in the absence of heterogeneity Bi et al. (2016), where a cell is either moving ($f_{r}=0$) or stuck ($f_{r}=1$). However with heterogeneity tissues in the range of $0<f_{r}<1$ become accessible and cells moving through must interact with rigid as well as non-rigid neighbor cells along its path of invasion. This results in a highly intermittent migration dynamics for the invading cell. To quantify intermittency, we plot the kurtosis in the displacement of the invading cellDas et al. (2016) as function of $f_{r}$. The displacement of the invading cells also show a burst-like dynamics when the tissue is in the heterogeneous solid state (Fig. 5(b) and SI Video 1 https://youtu.be/GTIcr0Wc-NQ).

Discussion -

In cancer research it has long been recognized that tumors are stiffer than the surrounding tissue. Yet tumors often contain softer cells and softer cells facilitate invasionOswald et al. (2017). This causes an apparent paradox: how can soft cells give rise to a tumor that is rigid at the tissue level? Our work suggests that tissue mechanics is not just determined by the average softness of the cells, but on the statistical fluctuations in single cell properties. The heterogeneity-driven rigidification can provide a possible resolution of this paradox. As a simple example, we note that in the phase diagram (Fig. 3(c)), it is possible to transform from a fluid to a solid state by increasing $\mu$ (i.e. cells become softer overall) but at the same time increasing $\sigma$. Furthermore, our predictions show that the tissue mechanics is controlled by a single fraction $f_{r}$. This is consistent with recent experimental findingGamboa Castro et al. (2016) that the fraction of mesenchymal cells can serve as a control parameter in describing jamming properties and that increasing mesenchymal fraction leads to the increase in motility. Finally, there is also experimental evidence that heterogeneity can drive more intermittent dynamics in cell migrationLee et al. (2012). The burst-like intermittency in the heterogeneous solid state is highly reminiscent of the pulsating cancer cell migration recently observed in epithelial monolayersLee et al. (2012).