“Patchiness” in mechanical stiffness across a tumor as an early-stage marker for malignancy

Abstract

Mechanical phenotyping of tumors, either at an individual cell level or tumor cell population level is gaining traction as a diagnostic tool. However, the extent of diagnostic and prognostic information that can be gained through these measurements is still unclear. In this work, we focus on the heterogeneity in mechanical properties of cells obtained from a single source such as a tissue or tumor as a potential novel biomarker. We believe that this heterogeneity is a conventionally overlooked source of information in mechanical phenotyping data. We use mechanics-based in-silico models of cell-cell interactions and cell population dynamics within 3D environments to probe how heterogeneity in cell mechanics drives tissue and tumor dynamics. Our simulations show that the initial heterogeneity in the mechanical properties of individual cells and the arrangement of these heterogenous sub-populations within the environment can dictate overall cell population dynamics and cause a shift towards the growth of malignant cell phenotypes within healthy tissue environments. The overall heterogeneity in the cellular mechanotype and their spatial distributions is quantified by a “patchiness” index, which is the ratio of the global to local heterogeneity in cell populations. We observe that there exists a threshold value of the patchiness index beyond which an overall healthy population of cells will show a steady shift towards a more malignant phenotype. Based on these results, we propose that the “patchiness” of a tumor or tissue sample, can be an early indicator for malignant transformation and cancer occurrence in benign tumors or healthy tissues. Additionally, we suggest that tissue patchiness, measured either by biochemical or biophysical markers, can become an important metric in predicting tissue health and disease likelihood just as landscape patchiness is an important metric in ecology.

Introduction

Cancer arises from malfunctioning cells where a deregulation of normal signaling pathways leads to the acquisition of hallmark features such as chronic proliferation and evasion of apoptosis, which disrupts normal tissue structure and function [1,2,3,4]. While the majority of focus in literature has been on how changes to biological signaling pathways can aid cancerous behavior, less attention has been drawn to the effects of altering mechanical signaling pathways [2,3,4]. Mechanical signals, generated via cell-cell and cell-matrix interactions and transcribed by corresponding changes in intra-cellular tension, cell size, cell shape and nucleus architecture, alter the localization, activation, and interaction of various protein effectors within the cytosol and the nucleus [5,6,7,8,9]. Modifications to these mechanical signaling pathways can disrupt a variety of cellular processes such as cell division, cell death, and cell differentiation via mechanotransduction - a process by which cells sense their physical environment and relay mechanical cues into biochemical signals [10, 11]. Malignant transformation has been associated with changes to both extra-cellular and intra-cellular mechanical properties, with studies suggesting that these altered physical properties are significant in promoting tumor progression [2, 4, 12,13,14]. One such manifestation of these altered mechanical properties in malignant cells is that they are significantly softer than their healthy counterparts, an observation consistent across many cancer types [4, 15,16,17]. This decrease in cell stiffness has been linked to certain cancerous features like uncontrollable proliferation, evasion of apoptosis, and increase in motility [3, 16, 18].

Studies have suggested that the stiffness of cells can grade its metastatic potential, where highly invasive malignant cells are on an average (at the population level) softer than less invasive cells [19, 20]. Since no two individuals are perfectly identical, cell stiffness measurements, irrespective of stiffness measurement techniques, produce a distribution of cell stiffness values for any given cell population. This distribution arises from phenotypic heterogeneity across individual cells and can be persistent over several cell generations [21,22,23] (Fig. 1). While the heterogeneity in cell mechanical properties (such as cell stiffness or adhesion) has been noticed and reported previously, it is the mean values of these distributions that is used to differentiate between the two cell types, and potentially grade the malignancy and metastatic potential of a particular cell population [14, 22,23,24]. We posit that, the distribution in cell mechanical properties across a population of cells comprising a tumor or even a normal tissue can provide significant insight into how the population will evolve over time and potentially lead to malignancy and metastasis (Fig. 1). This position is based on the fact that even healthy cell populations have broad distributions in cell mechanical properties indicating the presence of at least a small number of cells that have cancer-type phenotypes at the extremes of these distributions [19, 21,22,23, 25, 26]. Additionally, recent work on micron-scale, in-situ mechanical characterization of tumors has shown that tumors are made up of regions of soft and stiff cells [27], and the definitions and organization of these regions correlates with the aggressiveness of tumors [28]. We build on this position with the help of computational simulations of cell-cell interactions within a 3D tissue environment, observing the evolution of cell population with heterogenous mechanical properties over time. We focus mainly on the cell stiffness as a key mechanical property where individual cells show heterogeneity, and specifically look for the growth of cells with lower stiffness (softer cells) since this is a cell trait which is strongly associated with cancer [21, 22, 29].

Fig. 1

(a) Mechanical phenotypic diversity in a healthy tissue, forming a spatially homogeneous tissue (b) Mechanically similar cells are in close proximity with one another, forming a spatially heterogeneous tissue, (c) Based on spatial arrangement and mechanical variation within the tissue, some cells are able to replicate at a faster rate, overtaking the tissue and migrating through the basement membrane

Methods

Mechanics-based mathematical and computational models have been a common tool used to isolate and study the sole effect that specific tumor mechanical properties have on cancer progression [2]. Many multi-scale models exist, all intended to answer specific questions regarding cancer mechanics and tumor progression [3, 30,31,32]. We have previously developed one such model to understand how mechanically distinct populations of stiff and soft cells can interpenetrate [33, 34] or promote and drive tumor growth within stiff environments [3, 34]. The focus in these prior works was on the interactions between two mechanically distinct phenotypes of cells. Here, we build on these models to incorporate a heterogeneous population of cells with stiffness values distributed around a predetermined mean corresponding to a healthy cell’s stiffness.

In the model, cells are described as viscoelastic shells, characterized by a dense actin cortex and a liquid core, able to compete for space while interacting with other cells in the tissue [3] (Fig. 2). The position of each cell is defined by a single point while the cell shape and its neighbors are obtained by Voronoi Tessellations about these individual cell points. The polyhedral cell area, volume and the interface between neighboring cells are used to compute the mechanical energy stored within each cell at any given time using Eq. 1 (3,33)-

$${U_i} = \frac{1}{2}k\frac{{{{\left( {{A_i} - {A_0}} \right)}^2}}}{{{A_0}}} - \frac{1}{2}\sum\nolimits_n {{A_j}{\gamma _j}{\sigma _i} - \lambda {\rm{log}}\frac{{{V_i}}}{{{V_0}}}}$$

(1)

Fig. 2

(a) Individual cells are modeled as spherical in free solution, (b) cells interact with one another by forming a flat interface with neighboring cells, (c) a 2D representation of cells interacting with one another in tissue, where cells are modeled as Voronoi polyhedrons, (d) a 3D representation of a cluster of cells modeled as Voronoi polyhedrons, (e) a tissue comprised entirely of cells

The equation consists of three terms which estimate the energy of the stretched actin shell, the energy released by formation of inter-cellular bonds and the work done by the osmotic pressure inside the cell in changing its volume from ${V}_{0}$ in free solution to ${V}_{i}$ inside the tissue. $k$ is the stiffness of the actin cortex, ${A}_{i}$ is the area of the cell, and ${A}_{0}$ is the area of the cell had it been spherical, n is the number of neighbors, and ${A}_{j}$, ${\sigma }_{i}$, and ${\gamma }_{j}$ are the area of the interface, the bond density, and the bond energy between the cell and its j^th neighbor, respectively. λ is calculated as P₀V₀, where P₀ and V₀ are the osmotic pressure and volume of a cell free in solution. Cellular rearrangements, obtained by displacing the individual cell points, are dependent on the overall mechanical energy of the system (sum of mechanical energies of individual cells). Rearrangements that lower the total mechanical energy of the system ($\sum {U}_{i}$) are always accepted, while rearrangements that increase the energy of the system are accepted based on the probability of acceptance given by Eq. 2, following the Monte-Carlo Metropolis algorithm [3].

$$ {P}_{accept}=\left[\text{e}\text{x}\text{p}\left(\frac{-\left({U}_{current}-{U}_{prev}\right)}{{\left({k}_{B}T\right)}_{eff}}\right)\right]$$

(2)

The term U_current refers to the tissue’s total energy in the current configuration, U_prev refers to the total energy in the previous configuration and (k_bT)_eff refers to the internal energy of the cells, analogous to the work a cell can do via filopodial protrusions [3, 35].

Under physical stretch either via cell-cell interactions or external forces, cells show a higher likelihood of entering the S-phase of the cell division cycle, triggered by either stretch activated membrane channels, localization of key transcription factors or the restructuring of the chromatin within the nucleus [36,37,38]. Based on these observations, cell fate (death or division) in the model system is a stochastic function of cell stretch, derived empirically based on cell death and division likelihoods quantified in [39, 40], and given by Eqs. 3 and 4 (3).

$$ {p}_{death}={p}_{0}^{({A}_{i}-{A}_{sp})/({A}_{m}-{A}_{sp})}$$

(3)

$$ {p}_{div}={\left(1+\frac{(1-{p}_{0})}{{p}_{0}}{\left[\frac{1-{p}_{0}}{100-{p}_{0}}\right]}^{\raisebox{1ex}{$\left({A}_{i}-{A}_{m}\right)$}\!\left/ \!\raisebox{-1ex}{$\left({A}_{m}-{A}_{sp}\right)$}\right.}\right)}^{-1}$$

(4)

${A}_{i}$ is the area of cell i, which is combination of both cell shape and size in 3D and provides a good estimate of the cell-stretch. ${A}_{sp}$ is the area of the cell were it perfectly spherical, and ${A}_{m}$ is the area of an average cell (mean of ${A}_{i}$) in a homeostatic, healthy tissue. Cell death and division events drastically increase the system energy and are thus interspersed between a large number (∼ 70,000) of cellular rearrangements that drive the tissue back to a lower energy state as dictated by ${\left({k}_{B}T\right)}_{eff}$.

All of the parameter values are restricted to biologically relevant values based on either experimentally determined observations or prior theoretical estimations. For example, experimentally measured cell elastic modulus values (∼ between 100 and 1000 Pa) [21, 41] and actin cortex thickness values (∼ 1 μm) [42, 43] are used to estimate the cell cortex stiffness, k. Similarly, cell-cell adhesion strengths are based on the strength of individual cadherin bond strengths, ${\gamma }_{j}$, and the average density of cadherins on cell surfaces, ${\sigma }_{i}$. The choice of ${\left({k}_{B}T\right)}_{eff}$, accounts for the active motility of cells in tissues and its value is based on the average work a cell can do (force x displacement) while migrating using cellular protrusions. The parameter values used in the above calculations are summarized in supplementary table S1 along with sources from where they are derived.

The initial tissue/tumor configuration for a homogeneous tissue environment is obtained by randomly dispersing all the cell defining points into a 3D space (Fig. 2e). These location points act as the centers of their encapsulating, polyhedral cells. A log-normal distribution of cell stiffness values with a fixed mean of 500 Pa and differing amounts of variance is specified for each cell to model mechanical heterogeneity in the tissue [22, 29, 41]. Either high, moderate, or low amounts of mechanical heterogeneity is considered, modeled by changing the mode of the log-normal distribution. The primary results presented here use log-normal distributions for cell stiffness values to relate most closely with experimental measurements [22, 26]. The mode of the log-normal distributions for the different populations are set to either 465 Pa, 455 Pa, 445 Pa, 435 Pa or 425 Pa (corresponding to a standard deviation of 115 Pa, 130 Pa, 145 Pa, 160 or 175 Pa respectively).

In addition to modeling the mechanical variation, the effects of the spatial arrangement of cells within the healthy tissue are also tested, as cell-to-cell interactions may further impact tumor incidence [3] (Fig. 3). To setup up the initial conditions for a spatially heterogeneous tissue environment, a small number of cells (8, 27, or 64) are first randomly distributed into the 3D space and their mechanical properties are randomly chosen from the overall population distribution. The rest of the tissue is then sequentially seeded by placing cells with randomly assigned mechanical stiffness values in proximate locations to existing mechanically similar cells. This leads to clustering of mechanically similar cells while maintaining overall population heterogeneity. By altering the number of seed cells first introduced, we can alter the sizes of the clusters in the tissue (Fig. 3). The size of the resulting clusters is dependent on the number of initial seed cells, where fewer seed cells amount to larger, but fewer clusters of cells with similar mechanical properties within the healthy tissue environment. For each spatial arrangement and stiffness distribution, a complementary non-clustered tissue is simulated as a control. A large number of cell rearrangements are performed immediately after seeding all the cells (∼ 70,000 iterations), without any cell death or division, in order to reach an initial low energy configuration for the model tissue system.

Fig. 3

Initial healthy cell stiffness distributions applied to tissue with various spatial arrangements (a) Log-normal cell stiffness distribution with a mean of 500 Pa and mode 465 Pa (standard deviation of 115 Pa) applied to a spatially homogeneous tissue (b) Log-normal distribution cell stiffness distribution with a mean of 500 Pa and mode 425 Pa (standard deviation of 175 Pa) applied to a spatially non-homogeneous tissue, where clusters are formed between mechanically similar cells (64 seed cells around which other cells are sequentially populated)

This simulation adequately mimics the environment of a dense 3D tissue or tumor where cells are packed together in packing fractions close to 1. Cellular rearrangements drive the tissue to a steady energetic state. When a tissue is in an energetically steady state, cell fate decisions are made based on the stochastic Eqs. 3 and 4, where cell stretch and area provide a driving signal for cell death and division [39]. Empirically, perfectly spherical cells have a low probability of cell division and a high probability of cell death [44], while more stretched out cells with larger areas have a low probability of death and a high probability of division [45]. By using this probabilistic method of cell death and division, an important aspect of cell fate decisions – mechanosensing, is integrated within the model while preserving the homeostatic nature of the tissue. When the cell divides, the resulting two daughter cells will either retain the same mechanical stiffness as its parent cell or vary slightly, where mechanical phenotype is strongly conserved [46,47,48]. We assume a ± 5% (max) variation in daughter cell stiffness compared to the parent cell. After cell fate decisions are made, cellular rearrangements again drive the tissue to a steady state energy value with the same Monte-Carlo Metropolis thermalization process. Cell rearrangement takes place over a timescale on the order of days, whereas the cell death and division processes occur within a couple of hours [49]. This separation of timescales allows for the cellular rearrangement and cell death and division events to be simulated sequentially.

The model is implemented in MATLAB. We simulate an 80 μm x 80 μm x 80 μm cubic space, containing 512 cells (average cell volume of ∼ 1000 µm³) with periodic boundary conditions. Here we assume that healthy cell population has an average cell stiffness value of 500 Pa (starting mean stiffness of all cells in the tissue). Individual cells within this population with a stiffness value lower than 400 Pa are considered tumor-like, based on experimental results that determine that cancer cells are found to be at least 20% softer than healthy cells [19, 24, 33]. The cell shape and sizes are defined by Voronoi Polyhedra generated about the cell points and their nearest neighbors. The cell volume, surface area and shape are obtained by analyzing the convex polyhedron formed by the vertices of the Voronoi polyhedron. A single simulation cycle involves ∼ 70,000 cell rearrangement iterations where a few (7 here) randomly selected cells are moved a small distance (exponentially distributed around a mean of 0.5 μm) in each iteration. The Voronoi Tessellation is reperformed and new cell volume and shape parameters are obtained. The mechanical energy of each cell affected by the move is recalculated based on Eq. 1 and the new total energy of the tissue (defined as the sum of mechanical energies of each cell) is obtained. The new arrangement of cells is either accepted or rejected based on the Metropolis algorithm as described above (Eq. 2). The rearrangement iterations are continued until the average energy of the tissue (averaged over 1000 iterations) reaches a steady state (usually less than ∼ 70,000 iterations, supplementary figure S1). Once the tissue reaches a steady state where further cellular arrangements fluctuate the energy around a low value dictated by the effective temperature, ${\left({k}_{B}T\right)}_{eff}$, of the system, each cell can either divide, die or do neither based on probabilities calculated using Eqs. 3 and 4. The probability of doing neither is calculated as (1-(${p}_{div}+{p}_{death}$)). Death and division events push the tissue systems out of the low energy state steady, and the cellular rearrangement iterations are reinitiated to drive the tissue back to a new steady state${\left({}_{}\right)}_{}$.

The changes in the cell populations are tracked over 20 cellular rearrangement and death/division cycles, to determine if there are any shifts within the population and if there is an increase in the tumor-like cell population. For analysis purposes, we assume that a drop of 15 Pa in the mean cell stiffness within the tissue environment over the 20 death and division cycles implies a malignant transformation of the tissue. This drop in mean tissue stiffness can be a result of an increase in the softer cell population, a decrease in the stiffer cell population, or a combination of both. We can track the number of healthy (stiffness > 400 Pa) and cancer cells (stiffness < 400 Pa) within the population. The choices of tissue size, cell numbers, and simulation cycles are primarily restricted by the computational cost of this model. However, we believe the overall insight from this model system can be applied to larger collection of cells in tissues.

To quantify global vs. local heterogeneity within the tissue environment, we borrow the concept of patchiness from environmental ecology [50,51,52,53]. The patchiness index is calculated as the ratio of global to local heterogeneity using Eqs. 5,

$$ c= \frac{-{\Sigma }{P}_{i}.ln{P}_{i}/\text{ln}{N}_{total}}{mean(-{\Sigma }{P}_{i}.ln{P}_{i}/\text{ln}{n}_{cluster})}$$

(5)

where, P_i is the fraction of the population within the sample that is similar (in this case has similar cell stiffness) and N (N_total or n_cluster) is the population size of the sample. The numerator is the global heterogeneity calculated for the total population, while the denominator is the local heterogeneity calculated for a small region of the system. To calculate the P_i, the cell stiffness values are combined into 10 Pa increment bins, and P_i is calculated as the ratio of the number of cells in bin i to the total number of cells. The same process is used for both local and global populations. The local population heterogeneity is obtained using a moving window approach, where a cubic window of predetermined dimensions is sampled, the heterogeneity of the cell population within that window is calculated, and the window is then shifted by specified distance to obtain a new local sample population. The mean of the local heterogeneities thus calculated is used in Eq. 5 to obtain the patchiness index. Here we use window sizes of 30 μm x 30 μm x 30 μm, 40 μm x 40 μm x 40 μm, and 50 μm x 50 μm x 50 μm and shift the window by ∼ 1 cell length (10 μm) each time, along each of the three dimensions.

Results

As mentioned above, we use a decrease in mean cell stiffness by 15 Pa over 20 simulation cycles as the benchmark value when considering tumor occurrence and malignant transformation. This value is determined to be significant based on the constant drop of the mean value of the distribution, as well as the increase of tumor-like cells in the tissue after every cell fate decision. If the simulated tissue were to undergo more cell death and division cycles, one would expect larger drops in the mean value as well as an increased number of tumor-like cells in the tissue (supplementary figure S2). In the instances when tumor incidence is observed as defined by a 15 Pa drop in mean cell stiffness, it is often accompanied by an increase in tumor-like cells and a drop in the number of stiffer cells (see supplementary figure S3). However, there are also instances when the number of tumor-like cells in the tissue remains constant and the number of stiffer cells drops, which we believe is due to the inherent stochasticity in our system (see supplementary figure S4). In rare cases, the mean value of the distribution increases by 15 Pa and is accompanied by increases in the number of stiffer cells and a drop in the number of tumor-like cells (see supplementary figure S5). Although this is infrequent, it may be beneficial to understand the interplay between the mechanical variation and spatial clustering, and how it can sustain the growth of stiffer cells.

Overall, we simulate tissue systems with 5 different initial cell stiffness distributions (log-normal distributions with mean 500 Pa and modes − 465 Pa, 455 Pa, 445 Pa, 435 Pa or 425 Pa). For each of these 5 tissue systems, we further simulate four different levels of clustering within the tissue (high – 8 seed cells, medium – 27 seed cells, low – 64 seed cells and no clustering). For each of these 20 scenarios (listed in supplementary table S2), we simulate 10–12 unique instances of cell population evolution using our model, running each system for 20 cell death/division and reorganization cycles. For each of the main 20 scenarios, we count the runs where malignant transformation is observed as defined by a decrease in the mean stiffness by 15 Pa and divide it by the total number of simulations run for that scenario to find the probability of malignant transformation. We observe that both the mechanical variance and spatial arrangement collectively play influential roles in increasing the likelihood of malignant transformation (Fig. 4). Individually, these factors are not significant enough in inducing tumorigenesis but are instrumental in doing so together. As both the mechanical variation and clustering of cells with similar mechanical stiffness increases in an initial healthy tissue, the probability for tumor occurrence increases (see Fig. 4).

Fig. 4

The probability of tumor occurrence is dependent on increases in the mechanical variance and the spatial clustering between cells with similar mechanical properties

The need for spatial clustering of cells with similar mechanical properties for tumorigenesis reveals the dependency of interactions between tumor-like cells that is needed to promote their increased proliferation. This may be due to the enlarged surface area of the tumor-like cell clusters within the tissue, which gives them the ability to resist the force imposed by their surrounding stiff cells [3]. This can be related to what has been seen experimentally, where tumor cells are able to withstand their stiff environments as they grow and multiply, driven by an increased homeostatic pressure for these cells [18, 54].

We combine the effect of overall variance (heterogeneity) in cell mechanotype and the local clustering of mechanically similar cells using the patchiness index for the tissue. We find that independent of the cubic window size used to estimate the patchiness index (3 cell lengths, 4 cell lengths or 5 cell lengths along each dimension), there is a significant likelihood for tissues with patchiness index greater than 0.85 to show a decrease in mean cell stiffness and thus a transition towards a malignant mechanotype (Fig. 5, supplementary figure S6). In these figures, the stars denote p < 0.05 on a single population ttest, indicating a non-zero change in mean cell stiffness for tissues with that patchiness index. We note here that we are limited in our analysis to these small window sizes because of the small size of our simulation domain (8 cell lengths along each dimension).

Fig. 5

Change in mean stiffness of the cells against c value (ratio of local to global diversity). Window size of 4 × 4 × 4 cell lengths (40 μm x 40 μm x 40 μm) used to calculate the patchiness index using Eq. 5. Each open red circle represents the change in mean cell stiffness over 20 cell death and division cycles for a single simulation, plotted against the patchiness index of the starting initial tissue configuration for that simulation run. The error bars are standard errors of the mean. Stars above the bar plots indicate p values < 0.05 for a single population ttest

Discussion

Here we have used a purely theoretical setup to study the evolution of a mechanically diverse population of cancer cells within a densely packed 3D tissue environment. We have ignored the presence of nutrient gradients, extracellular matrix, or significantly different cell types, concentrating only epithelial or near epithelial cell types. Cell-cell interactions are assumed to be purely mechanical and cell fate decisions are also driven primarily by mechanosensitive feedback based on empirical experimental observations. While all these amount to a gross simplification of a complex tissue system, such simplified models have previously been employed to understand and describe processes related to tissue dynamics [55, 56], tissue growth [57], tumor incidence [3], tumor growth and metastatic invasion [58,59,60]. Additionally, the small length scale of the tissue system simulated in this case allows for the assumption of a dense collection of tightly packed cells without the inclusion of an extracellular matrix, nutrient gradients, or vasculature. While additional complexity can be added to the model, we leave these for future extensions of the model.

As mentioned above, this study focuses on mechanical heterogeneity at the individual cell level and patchiness at lengths scales of a few cells (∼ 40 μm). This is relevant to the length scales of recent tissue mechanotype characterization studies [27, 61] as well as tracking the effect of stiffness changes on individual cell fates in tissue and tumors. The origin of mechanical heterogeneity in cell populations at this scale can be attributed to random genetic mutations, gene methylation states, or phenotypic changes arising from local environmental cues such as extra-cellular matrix properties and nutrient availability. However, it does not answer questions about macroscale heterogeneities in tissues and factors that contribute to bulk tissue stiffness measurements such as recruitment of various other cell types such as fibroblasts, endothelial cells and vasculature, and immune cells to the tumor environment, and consequent modifications to the stromal tissue surrounding tumors.

Within the confines of the limitations described above, we use our model to understand how heterogeneity in the mechanical properties of tumor cells, specifically, cell stiffness, as observed and quantified by experimental tools, may be used to predict the occurrence of tumor growth and a malignant shift in a seemingly healthy population of cells. We find that heterogeneity at the cell population level, but homogeneity at a local level, quantified by the patchiness of cell distributions within the tissue system, predicts an increasing likelihood of tumor occurrence and an overall shift towards malignant populations within an initially health tissue. It is interesting to note that overall high heterogeneity, by itself is not detrimental to tissue fate, so long as the population remains well mixed. However, local clustering of similar populations, which may occur due to an underlying heterogeneity in the extra-cellular environment, nutrient availability, stochastic cell division events, or motility-dependent segregation of cells, tips the balance to drive malignant transformation in highly heterogenous tissue systems [33, 34, 62]. These observations are indeed limited to tissue level heterogeneity at the short length scales of 10 to 100 s of µm and do not consider greater tumor environment heterogeneity associated with varying cell types accumulating within the tumor stroma as well as degrees of vascularization and nutrient perfusion. However, since the focus of this discussion is on the transition from a healthy or pre-malignant state to a malignant state, we posit that the commonly tumor-associated macroscale heterogeneities will not have yet accumulated within the tissue environment. Instead, the heterogeneity at the tissue cell population level, still at homeostasis with its environment, holds more significance in driving malignant transformation.

From an ecology of cancer perspective, these results augment several ideas proposed before [63, 64]. For example, a high diversity in the cell population defined by a high Evo-index, or the presence of a suitable environmental niche defined by a high Eco-index have been proposed as ecological indices predictive of tumor progression, disease outcome as well as treatment efficiency [65]. The model we have presented looks at the contribution of both the overall heterogeneity of the native tissue system as well as the formation of favorable local niches. Here the local niches are dictated by cell-cell interactions, rather than cell-matrix or cell-environment interactions, which is along the lines of cell fate control in localized ecological niches via cell-cell interactions as proposed by Adler and Gordon [66]. However, additional control based on extracellular factors such as nutrient and oxygen gradients, presence of cytotoxic and genotoxic agents, and direct interactions between cell and the extra-cellular matrix can be integrated into future modeling efforts along with the present cell-cell mechanosensitive interactions [67, 68]. Another key aspect of both ecology and tumors, that is homeostatic competition [69, 70], is also captured in the present model. Homeostatic competition implies a steady state population density at the global population level driven by resource competition or cooperation between the various species inhabiting an ecological landscape. Maintenance of homeostasis at the tissue level is a key condition built into the healthy tissue model described here. The low diversity tissue maintains homeostasis through competition between the cell types for available space. Even a high diversity, well mixed tissue system maintains homeostasis, where the extreme cancer-like cell types are kept in check by their more normal neighbors. However, when sufficient cancer-like cells cluster together, their mutual interactions promote cell division and suppress cell death, allowing them to outcompete their neighbors for the available space and begin to dominate the tissue. These observations also align with parallels between dormancy, the Allee effect and growth lag in ecological and cancer communities [71].

Analysis of the spatial heterogeneity of tumor environments has further highlighted parallels between ecological landscapes and tumors. The tumor environment not only contains a high diversity in the profile of its native tissue cells, but also diversity in the cells of the surrounding stroma, specific immune cells associated with the tumor and the overall architecture and organization of these cell types within a fibrous extracellular matrix [72]. Based on these large-scale tumor observations, a few parametric and non-parametric measures of spatial heterogeneity have been proposed to provide diagnostic and prognostic insights about the tumor [73,74,75]. These measures have shown success, especially when applied to advanced tumors, but not as much when predicting the likelihood of benign neoplasms turning carcinogenic, or low-grade tumors. The question of why and when benign neoplasms, which are much more common across all forms of life, including mammals, progress to aggressive cancers that grow, and spread is an important question in oncology [76]. The transition from a homogeneous to a heterogeneous population (which is inevitable in most biological systems), and a transition from a mixed to a segregated population with regions of low heterogeneity (which is most likely the rate limiting step), as shown by our results, might hold the key to answering this question. The patchiness index used here may also help draw better parallels between cancer and ecological systems with high number of neoplasms but no catastrophic events that unbalance the ecosystem.

We have restricted our study to the heterogeneity in the cell mechanotype defined by cell stiffness, since this is one of the key mechanical properties known to distinguish cancerous vs. normal cell populations and can also help grade cancer cell populations based on their aggressiveness [19, 20, 77, 78]. Additionally, there are tools being developed that can measure and map cellular mechanical stiffness in situ within tissue samples [79, 80]. However, other cellular mechanical properties such as cell adhesion strength, cell contractility and cell nuclear deformability also show similar heterogeneity levels within tissue and tumor cell populations [81,82,83]. Both cell adhesion strength and cell contractility have similar effects to that of cell stiffness on cell organization, shape, and size within vertex-based cell-cell interaction models [84], and we believe will lead to similar observations as those presented. There are few models currently incorporating cell nuclear stiffness when considering cell-cell interactions and tissue dynamics, and this is something that needs to be worked on in the future. Lastly, these mechanical differences between cells arise from biochemical differences, likely driven by differences in gene expression profiles within a cell population. Thus, cell phenotypic differences may directly correlate to cell mechanotype differences, and the analysis of tissue patchiness may be extended to any spatial characterization of cells within a tissue environment. Indeed, with invention of accurate single cell gene expression analysis and tracking tools, as well as large scale spatial phenotypical and mechanical profiling of cells in-situ [85, 86], it might be possible to experimentally characterize tissues based on their patchiness in both mechanical and biochemical properties. Based on the results presented here, we propose that tissue patchiness can provide key diagnostic and prognostic insights for malignant growth in health tissues and benign tumors.

Heterogeneity is a norm in biology and cannot be done away with. This is true across scales and species. On the other hand, limiting population de-mixing and clustering into local niches may be a potential strategy that can be employed to avoid tipping of the scales and allowing one sub-species to dominate. Patchiness may be avoided by limiting factors that promote proliferation of only certain sub-populations, or factors that limit the mobility of cells within certain regions of the tissue. Biological events and biochemical or biomechanical factors that potentially increase or decrease patchiness within normal tissues or tumors, and their relation to actual tumor growth and malignancy needs to be further investigated and discussed. However, we strongly believe that the idea of patchiness might be just as applicable to cancer ecology as it is to environmental ecology and population dynamics.

Conclusion

Here we present a model that focuses on cell interactions and mechanoreciprocity as drivers of tissue dynamics. We use this model to understand how global and local heterogeneity in cellular mechanics across a dense population of cells may lead to an incidence of malignant tumor growth defined by an increase in the population of cancer-like soft cells and decrease in the population of normal, stiff cells. Based on our results, we find that tissue patchiness as defined by the ratio of global to local heterogeneity may be an excellent metric to predict malignant transformation in healthy tissues or benign tumors. While limited by the purely theoretical nature of this study and a sole focus on cell mechanics, the potential for such a metric is extremely appealing for early diagnosis and intervention in cancer patients as well as prognosis in patients with benign tumors. Additionally, the model is not only limited to studying the effects of cellular mechanical properties on tumorigenesis in a healthy tissue, but it can also serve in studying a broader scope of tissue mechanics in processes like wound healing and aging. Our model also provides a versatile platform that can be built upon, with the ability to study the effects other mechanical properties that have been linked to cancerous behavior or potential interventional strategies focused on manipulating the ecological landscape of tissues and tumors.