Tissue evolution: mechanical interplay of adhesion, pressure, and heterogeneity

The evolution of various competing cell types in tissues, and the resulting persistent tissue population, is studied numerically and analytically in a particle-based model of active tissues. Mutations change the properties of cells in various ways, including their mechanical properties. Each mutation results in an advantage or disadvantage to grow in the competition between different cell types. While changes in signaling processes and biochemistry play an important role, we focus on changes in the mechanical properties by studying the result of variation of growth force and adhesive cross-interactions between cell types. For independent mutations of growth force and adhesion strength, the tissue evolves towards cell types with high growth force and low internal adhesion strength, as both increase the homeostatic pressure. Motivated by biological evidence, we postulate a coupling between both parameters, such that an increased growth force comes at the cost of a higher internal adhesion strength or vice versa. This tradeoff controls the evolution of the tissue, ranging from unidirectional evolution to very heterogeneous and dynamic populations. The special case of two competing cell types reveals three distinct parameter regimes: two in which one cell type outcompetes the other, and one in which both cell types coexist in a highly mixed state. Interestingly, a single mutated cell alone suffices to reach the mixed state, while a finite mutation rate affects the results only weakly. Finally, the coupling between changes in growth force and adhesion strength reveals a mechanical explanation for the evolution towards intra-tumor heterogeneity, in which multiple species coexist even under a constant evolutionary pressure.


Introduction
Mutations change the cell fitness and thus its chance to survive and proliferate [1].Advantageous mutations are more likely to persist due to natural selection, which drives the evolution of a tissue towards fitter cells [2].Cancer represents an example of evolution on a short time scale [3].Furthermore, cancer is a multistep process, i.e. several mutations are needed for a tumor in order to develop and become malignant [4].Hence, tumorigenesis might be expected to happen in a serial manner, i.e. a cell acquiring a "beneficial" mutation and taking over the whole tissue.After some time, a daughter cell acquires another mutation and again takes over.Interestingly, however, tumors do not consist of a single cell type, but instead several subpopulations coexist within the same tumor.This is called intra-tumor heterogeneity [5].
Each mutation changes certain biochemical properties of a cell.This ranges from misfunction in the error correction machinery during DNA replication and disruptions in signaling pathways to epigenetic changes in the expression level of certain proteins [1,6,7].All these changes can also affect the mechanical properties of the mutated cell, e.g.mutated cells which express less adhesion proteins might be able to detach from the primary tumor more easily [8], necessary to form metastases. On the other hand, mechanics feeds back onto growth in several ways, e.g.increased apoptosis rate due to mechanical stresses [9,10] or dependence of the growth of tissue spheroids on the properties of the surrounding medium [11][12][13].
It is the mechanical contribution to tissue development that we want to focus on in this work.For mechanically regulated growth, homeostatic pressure plays an important role [14].In the homeostatic state, when apoptosis and division balance each other, a tissue exerts a certain pressure onto its surrounding, the homeostatic pressure P H .The tissue is able to grow as long as the external pressure P is smaller than P H .For the competition between different tissues for space, it has been suggested that the tissue with the higher homeostatic pressure grows at the expense of the weaker tissue.Several theoretical studies employ this concept in order to describe interface propagation between two competing tissues [15][16][17].A metastasis would need to reach a critical size, below which the additional Laplace pressure due to surface tension would cause the metastasis to shrink and disappear [14].However, reduced adhesion between tissues, which increases surface tension, leads to an enhanced growth rate at the interface between them, stabilizing coexistence even for differing homeostatic pressures [18].
In this work, we study the influence of mutations that change the mechanical properties of cells on the competition dynamics, especially the interplay between changes in the adhesive properties and the strength with which a cell pushes onto its surrounding.Particularly interesting is the case where loss of adhesion comes at the cost of lower growth strength.This is motivated by the observed down-regulation of E-cadherin, an adhesion protein in epithelia, in many types of cancer [19].Interestingly, E-cadherin is also involved in signaling processes connected to cell growth [20].We find that in this case several cell types with different mechanical properties can coexist and that the cell type with the highest homeostatic pressure does not necessarily dominate the competition.

Results
Several models have been developed previously in order to study tissue growth [21], in combination with different simulation techniques, including vertex [22,23] and particle-based [24,25] models as well as Cellular Potts models [26,27].We employ the two particle growth (2PG) model of Refs.[18,28,29].A cell is described by two particles which repel each other via a growth force with strength G, unit vector r ij , distance r ij between the two particles and a constant r 0 .Different cells interact via a soft repulsive force F V ij on short distances, maintaining an excluded volume, and a constant attractive force F A ij on intermediate distances, modeling cell-cell adhesion, with with exclusion coefficient f 0 , adhesion strength coefficient f 1 , and cut-off length R PP .A cell divides when the distance between its two particles reaches a size threshold r ct .A new particle is then placed close (randomly within a short distance r d ) to each of the two particles of the divided cell.Each of these pairs then constitutes a new cell.Apoptosis is modeled by removing cells randomly at a constant rate k a .We employ a dissipative particle dynamics-type thermostat, with an effective temperature T , to account for energy dissipation and random fluctuations.We choose the value of T such that cells can escape local minima, but other thermal effects are negligible.Note that all parameters can be set individually for each cell type as well as between different cell types for inter-cell interactions.We only vary the growth-force strength G α and adhesion strength f αβ 1 between cells of the same (α = β) and different (α = β) cell types, respectively, where α and β are cell-type numbers.We report simulation parameters relative to a standard host cell type (see Materials and methods for numerical values), denoted with a dagger, e.g.G † = G/G 0 .Time is measured in terms of the inverse apoptosis rate k a , distance in units of the pair potential cut-off length R PP and stresses in units of G 0 /R 4 PP .Quantities reported in these units are denoted by an asterisk * .All simulations are performed in a cubic box with edge length L = 12 • R PP and periodic boundary conditions in all directions, unless stated otherwise.
Tumor cells even within the same tumor are not all identical, but vary in terms of all kind of attributes, e.g.expression levels of different proteins [30] or their reaction to certain treatments [31].Hence, there is not only a competition between the tumor and the host, but also between cell-subpopulations of the tumor.Different models exist to describe tumor heterogeneity, e.g.cancer stem cells [32] or clonal evolution [33].In the latter case, a tumor originates from a single mutated cell, which can acquire additional mutations over time, yielding additional subpopulations.We model this behaviour by defining a fixed number n of different "genotypes", each having a different growth-force strength G α and adhesion strength f αα 1 .Mutations are implemented by offering each daughter cell after a division event the chance to change its genotype with a certain probability.
In tissues, several adhesion mechanisms exist, serving a variety of different functions to maintain tissue integrity.Between epithelial cells, the strength of cell-cell adhesion is to a large degree regulated by anchoring junctions, e.g.adherens junctions, which connect the actin cytosceletons of neighbouring cells.Adherens junctions are mediated by cadherins, which form homophilic bonds between cells.Thus, the strength of adhesion between cells is limited by the cell expressing less cadherin, or, in terms of our simulation model f αβ 1 = min(f αα 1 , f ββ 1 ).A reduced adhesion strength yields a higher homeostatic pressure [29], which is otherwise dominated by the growth-force strength G.For free parameter evolution, the tissue thus evolves to a strong-growing and low-adhesive genotype (see Fig 1), as predicted by the homeostatic pressure approach [14].
However, E-cadherin also plays a role in signaling processes connected to cell growth, and thus a reduced expression might come at the cost of a lower growth-force strength G, which in turn yields a lower homeostatic pressure.We thus turn our attention to the case where an increase in growth-force strength G α comes at the cost of a higher self-adhesion strength f αα 1 .We assume the relations as with genotype number α in the range [−(n − 1)/2, (n − 1)/2], evolutionary distance After a division event, each daughter cell might mutate into a new genotype with probability p m .If the cell mutates, its genotype number is changed to α mother ± 1 randomly.This yields a mutation rate k m = 2p m k a .
Figure 2 displays results of such simulations for four different cases: only variation of growth-force strength (τ = 0), balanced tradeoff (τ = 1), adhesion strength varied twice as much as growth-force strength (τ = 2) and only variation of adhesion strength (τ → ∞).Without tradeoff (Fig 2a)), the tissue evolves towards the strongest growing genotype or, equivalently, the one with the highest homeostatic pressure.Similarly, for τ → ∞ (Fig 2d)), the system evolves towards the lowest adhesive genotype (again, the one with the highest P H ). We find the most dynamic evolution for a balanced tradeoff (Fig 3 and 2b)).At first, the system evolves to stronger growing and more adhesive genotypes.Over time a noticable fraction of cells evolves also towards weak-growing, less adhesive genotypes.The cell-number fractions φ α = N α /N (with individual and total number of cells, N α and N ), show large fluctuations (see Fig 3b ) and c)), with individual genotypes not being populated at all for certain time periods.Besides this highly dynamic temporal evolution, after an initial time period the system is dominated by genotypes with increased growth force and adhesion strength at all times, with the one at the upper boundary having the highest cell-number fraction for most of the time (see Fig 3a)).This result comes at a surprise, as this is also the genotype with the lowest homeostatic pressure, while the one at the lower boundary, which is basically never populated, has the highest P H .For a higher tradeoff (Fig 2c)), we still find a broad distribution of genotypes, with less adhesive genotypes dominating over the stronger growing ones, i.e. the loss in growth-force strength is overcompensated by a lower adhesion strength.
In order to gain insight into the underlying mechanism of this dynamic evolution, we study the competition between two genotypes and no mutations (p m = 0).Simulations are started from a single mutated cell (with increased/decreased growth force and adhesion strength) in a host tissue at the homeostatic state (we label the mutant with M and the host (wild type) with W).Even in this simplified case, we find one parameter regime in which the mutant is not able to grow, one regime with stable coexistence in a highly mixed state and another regime in which the mutant outcompetes the host.Figure 4 shows the averaged number fractions of the mutant at the steady state.For reduced growth force and adhesion strength (Fig 4a)), the mutant can only grow against the host if its adhesion strength is reduced below a critical f crit 1 .In terms of Eq. ( 4), the value of f crit 1 roughly corresponds to a balanced tradeoff (τ ≈ 1).Already for 1 , the homeostatic pressure of the mutant exceeds the one of the host, i.e. a parameter regime exists in which the mutant is not able to grow, despite of the higher P H .The reverse happens when growth force and adhesion strength are increased.The mutant completely takes over the compartment, although its homeostatic pressure is smaller than that of the host.Again, coexistence is only found when the adhesion strength is increased above f crit 1 .In the coexistence regime, the mutant number fraction scales as ).Altogether, the competition between two genotypes alone yields the same qualitative results as the more complex multi-genotype case discussed before.Still, the question remains how a genotype with lower homeostatic pressure can outcompete a stronger 5/13  for various (reduced) growth-force strengths G M † .Error bars are obtained via block-averaging method (hidden behind markers) [34].Dashed vertical lines indicate the points below which the mutant has a higher homeostatic pressure, solid lines are fits to Eq. 10. b) Same as in a) but for increased growth force and adhesion strengths of the mutant.
genotype.The answer can only lie in the adhesion strength between mutant and host cells.This choice of cross-adhesion strength breaks symmetry, as the stronger adhering genotype has more free space at the interface, which favors divisions [18].
To address this question, we develop a phenomenological model which incorporates pressure-dependent growth as well as interfacial effects, in order to obtain a qualitative explanation of the simulation results.
We start with the expansion of the bulk growth rate k b around the homeostatic pressure, with the pressure response coefficient κ.Due to the high degree of mixing, the number fractions φ M/W and hence the strengh of interfacial effects vary locally.In a mean-field approximation, we take the interfacial effects to be proportional to φ M (1 − φ M ), with individual prefactors ∆k M/W s for each genotype.The time evolution is then given by with the difference in homeostatic pressure ∆P H = P W H − P M H . Addition of Eqs. ( 6) and ( 7) yields the pressure Thus, the pressure is given by the homeostatic pressures of the two genotypes weighted by their number fraction plus an interfacial term.
We discuss this result for the case of reduced growth force and adhesion strength of the mutant.∆k M s might be expected to vanish, as and mutant cells thus would not feel whether neighbouring cells are mutant or host cells.However, in order to grow, a cell needs to impose a strain on its surrounding.Host cells adhere more strongly to each other, thus it is harder for a mutant cell to impose a strain when surrounded by host cells.Hence, ∆k M s is actually negative and the homeostatic pressure of the mutant needs to exceed the host pressure by −∆k M s /κ in order to be able to grow against the host.At this point, φ M 3 becomes positive, as long as ∆k M s + ∆k W s > 0. Host cells can impose a strain more easily when surrounded by mutant cells and, additionally, have more free space than when surrounded by other host cells.Hence, |∆k M s | < ∆k W s and the above mentioned condition is fulfilled.Similarly, coexistence can be found for increased growth force and adhesion strength when ∆P H > −∆k W s /κ.The above mentioned scaling of the mutant number fraction can be obtained by an expansion of ∆P H and ∆k M/W s to linear order in terms of : The zeroth order terms of ∆k M/W s vanish as there are no interfacial effects when the adhesion strength between host and mutant cells is equal to their self-adhesion strength, while ∆P 0 H can be non-zero due to a changed growth-force strength.Indeed, Eq. ( 10) reproduces the simulation data reasonably well (see Fig 4).A discussion of the numerical values of the fitted parameters and additional results can be found in S1 Appendix.
Figure 5 displays similar results as shown in Fig 4, but now as a function of the tradeoff τ in Eq. (4).For τ < 1 the genotype with higher growth-force strength outcompetes the weaker genotype, for 1 < τ < 2 a transition towards the less adhesive genotype occurs, while for even higher values of the tradeoff τ > 2 the less adhesive genotype outcompetes the second genotype.This transition from strongly growing, 7/13 adhesive to weakly growing, less adhesive genotypes is found in the same range of τ as in the competition between many genotypes.Hence, the simplified case of two competing genotypes captures the essential physics to explain the coexistence between many competing genotypes and, additionally, provides a quantitative description.
Next, we turn our attention to the effect of a finite mutation rate on the evolution of the system.Figure 6a) shows the number fraction of the mutant as a function of k m for different combinations of evolutionary distance D α and tradeoff τ , in comparision to the number fraction reached for a single mutation event.As expexted, the number fraction converges towards 1/2 with increasing k m for all combinations.For moderate mutation rates, however, the number fraction largely fluctuates around the same average as of a single mutation event.The single mutation leads to a stable coexistence of the two genotypes -additional mutations quickly relax back to this state.Siginificant deviations occur only if in the steady state of the single mutation event the number fraction of the weaker genotype is close to zero.In that case, the weaker genotype consists only of one or very few small cohesive clusters of cells, because cells of the weaker genotype need to detach from the primary cluster in order to form new clusters, but are likely to die when they do so, as they are only surrounded by cells of the stronger genotype.Hence, the distribution of cells is highly non-homogenous.Compared to the single mutation event, even a small mutation rate leads to the formation of multiple small cluster all over the system, thus increasing the number fraction of the weaker genotype (see Fig 6b ) for comparision in terms of number of clusters and Materials and methods for further discussion).This result explains why at least two genotypes, in addition to the dominating genotype, are populated as well in the cases shown in Fig 2a ) and d).When the number fractions of both genotypes are sufficiently large (for 1 ≤ τ ≤ 2), deviations from the average of a single mutation are still small for the standard mutation probability.Additionally, in the competitions between many genotypes, mutations change the genotype to α ± 1 randomly and not in a preferred direction.Hence, we conclude that the precise value of the mutation probability does not play an important role in the regime where we find a heterogeneous distribution of genotypes, as long as it is reasonably small (k m k a ).
1 < f crit 1 , p s increases linearly with further decreasing adhesion strength.On the other hand, when growth force and adhesion strength are increased, the survival probability first shows a plateau, whose value increases with increasing growth force strength, from which it will probably drop to zero with further increase.Simulations in this regime are difficult, because a mutated cell can easily grow to a few cells, but will hardly reach the number threshold nor completely vanish again.Due to the high self-adhesion strength on the one hand, it becomes hard to detach from the other cells, but on the other hand easy to grow against the host when only few or no other mutant cells are around.This explains the larger error bars at the highest values of the adhesion strength, where the sample size is small.

Discussion
We have shown how intra-tumor heterogeneity, the existence of multiple subpopulations within the same tumor, can arise due to mechanical interactions alone.The simultaneous change of the adhesion and growth-force strength stabilizes the coexistence of multiple subpopulations, in a highly dynamic state.A higher growth-force strength alone, as well as a lower adhesion strength, favor proliferation of a single subpopulation and the evolution of the system to cell types with the highest growth-force strength, or lowest adhesion strength, respectively.A tradeoff between the two, however, yields coexistence between multiple subpopulations of different cell types.Interestingly, the expression of the adhesion protein E-cadherin, which also affects cell growth, has been found to be down-regulated in many real tumors [19].
The simulations also reveal that the homeostatic pressure of a cell type is not 9/13 necessarily the only quantity that determines the result of a competition.Interactions between different cell types, in our model determined by the adhesion between them, can lead to a completely reverse outcome, i.e. a cell type with lower homeostatic pressure can outcompete a stronger one completely.A phenomenological model explains the results on a qualitative level.The evolution of each cell type is governed by mechanically-regulated growth, while mutation rates only play a minor role in the dynamics.
An interesting future aspect to be studied is the influence of open boundaries.A tissue with a negative homeostatic pressure then naturally grows to a spheroid of finite size, with an enhanced rate of division at the surface [29].For competing cell types, this would lead to an interplay between surface and interfacial effects.

Standard (host) tissue and simulation parameters
We define a set of reference simulation parameters, which we refer to as host parameters.Table 1 shows the values in simulation units.In simulations we keep the host W fixed and vary the parameters of the mutant M around the values of the host.

Cluster analysis
As explained in the results section, a constant rate of mutation leads to an enhanced formation of clusters when the weaker genotype is barely able to grow against the stronger genotype and consists of only one or few clusters for a single mutation event.We define a cluster as all cells of the same genotype that are in interaction range to at least one other member of the cluster (DBSCAN clustering algorithm with number of minimal points equal to one). Figure 6b) displays the number of clusters of the weaker genotype in the competitions displayed in Fig 6a ), in comparison to the result of a single mutation event.Indeed, when the number fraction of the weaker genotype is small for a the single mutation event (τ = 1), we find significant deviations even for small mutation rates.In this case, the number of clusters first strongly increases with mutation rate, with roughly a tenfold increase at the peak.For even higher mutation probability, the number of clusters decreases again, due to merging of clusters, finally leading to percolation.

Fig 2 . 1 (
Fig 2. Time evolution of the cell-number fractions φ α of each genotype for tradeoff paramter a) τ = 0, b) τ = 1, c) τ = 2 and d) τ → ∞, d → 0. Simulations start from a host (standard) tissue at homeostasis, with n = 21 genotypes, p m = 0.01 in all and d = 0.025 in a)-c).White space corresponds to times where no cells of the genotype exist.Color is coded on a logarithmic scale.Curves above display homeostatic pressure P α * H (black solid), growth-force strength G α † (red dashed) and self adhesion strength f αα † 1 (green dotted) of the corresponding genotype.

Fig 4 .
Fig 4. a) Average number fraction φ M of the mutant in terms of its adhesion strength f MM † 1

1
Fig 5. a) Average number fractions of the mutant (same simulations as shown in Fig 4a)) as a function of the tradeoff τ of Eq. (4) for different evolutionary distances D α .b) Same as in a) but with results from Fig 4b).Error bars are obtained via block-averaging method (hidden behind markers).

Fig 6 .
Fig 6. a) Average number fraction φ M of the mutant as a function of the mutation rate k * m for different values of the tradeoff τ , for evolutionary distances D α = 0.025 (circles) and D α = −0.025(triangles).Horizontal dotted lines display results for a single mutation event.Vertical black dashed line indicates standard mutation rate.Solid lines are a guide to the eye.b) Average number of clusters of the weaker genotype N c measured in the same competitions in a).Horizontal dotted lines display results for a single mutation event.Error bars are obtained via block-averaging method.

Table 1 .
Simulation parameters and measured properties of the standard (host) tissue .