Jamming and force distribution in growing epithelial tissue

We investigate morphologies of proliferating cellular tissues using a newly developed numerical simulation model for mechanical cell division and migration in 2D. The model is applied to a bimodal mixture consisting of stiff cells with a low growth potential and soft cells with a high growth potential; cancer cells are typically considered to be softer than healthy cells. In an even mixture, the soft cells develop into a tissue matrix and the stiff cells into a dendrite-like network structure. When soft cells are placed inside a tissue consisting of stiff cells (to model cancer growth), the soft cells develop to a fast growing tumor-like structure that gradually evacuates the stiff cell matrix. The model also demonstrates 1) how soft cells orient themselves in the direction of the largest effective stiffness as predicted by the theory of Bischofs and Schwarz (Proc. Natl. Acad. Sci U.S.A., 100, 9274--9279 (2003) and 2) that the orientation and force generation continue a few cell rows behind the soft-stiff interface. With increasing inter-cell friction, tumor growth slows down and cell death occurs. The contact force distribution between cells is demonstrated to be highly sensitive to cell type mixtures and cell-cell interactions, which indicates that local mechanical forces can be useful as a regulator of tissue formation. The results shed new light on established experimental data.


I. INTRODUCTION
Morphology and dynamics of proliferating cells are fundamental issues in cellular development [1][2][3][4][5][6]. They are controlled by a number of factors, but from the physical point of view, morphology is tightly coupled to intercellular force transmission, see e.g. Refs. [6][7][8]. Mechanical forces have been shown to be important in tissue healing after damage [9] and cancer development, and it has been suggested that tumor growth may even become arrested by inter-cellular mechanical forces [10][11][12]. In addition to uniform structures, a plethora of structures with various mechanisms and division modes have been suggested but the issue remains largely unresolved [13][14][15][16].
Among the many complications in investigating force transmission are that at their embryonic state, cells may not yet have developed junctions and may display more fluid-like behavior, and that cell-cell adhesion depends on the cell type and type of adhesion (focal or nonfocal) [7,15,17,18]. Junctions are crucial in cell-to-cell stress transmission [7,8,19,20] but it is, however, challenging to probe the individual junctions experimentally. In addition, the different cellular level signaling mechanism can be coupled depending on local properties and environment, that is, chemical signaling based on molecular interactions and mechanotransduction may depend on each other, see e.g. Ref. [21].
From a coarse-grained point of view, i.e., ignoring chemical details and treating cells as elastic objects, cel-lular systems can be seen as disperse grand canonical soft colloidal systems under evolving pressure. Several studies have tried to capture aspects of growing soft matter systems [22][23][24][25] but even in simple systems many fundamental questions remain open including the precise nature of colloidal phase diagrams when colloids are soft with size dispersity [26], and structure selection via selfassembly [27]. Cellular systems are even more complex since they exhibit additional behaviors such as growth and division, have varying mechanical properties (e.g. elasticity and cell-cell adhesion) and their responses to external stimuli may be sensitive to the local environment.
Dimensionality has an important role in regulation of intra-and inter-cellular forces at different levels, see e.g. Refs. [28][29][30][31][32]. Systems such as epithelial tissues and Drosophila wing discs, are inherently two dimensional which gives them distinct morphological properties due to the nature of cellular packing, and transmission of and response to forces [5,15,33]. In addition, jamming can be very strong in two dimensions; understanding the effects of stiffness, density and inter-membrane friction is crucial for being able to determine how jamming emerges in cellular systems [34,35]. Besides being important in understanding the mechanisms of cell movement under pressure [8], such situations have been proposed to be important in tumor growth [10,11] -cancer cells are often softer than healthy cells [36] although the opposite has also been reported [37]. Our main focus is on the above effects in systems consisting of hard cells in a soft matrix and vice versa. We use our previously introduced  1. Illustration of the model and the forces (Eq. 1). As the cells grow and deform, they reach the division threshold. At that time, a division plane is identified based on the long-axis (Hertwig's) or some other selected rule. As the cell divides, new beads must be added to preserve the topology of the cells. Once division is completed, two child-cells are produced. Colors indicate the different forces and the new beads that need to be created at the time of division. two dimensional computational model called EpiCell2D (Epithelial Cell 2D) [38]. The model is summarized below in Models and Methods and described in full detail in Ref. [38]. Cell stiffness, its measurements and connection to cancer metastasis have been recently reviewed by Luo et al. [39].
Two dimensional cellular systems have been studied with a number of computational methods including offlattice vertex models [14,[40][41][42][43][44][45][46][47][48][49], and Voronoi tessellation or Delaunay triangulation based models [50][51][52][53][54]. Such models approximate cell membranes as edges or planes and they do not include inter-cellular forces. The immersed boundary method of Rejniak et al. [55] which explicitly models inter-membrane interactions is better suited to problems that require inter-cellular forces. Other methods include lattice-based methods such as the Cellular Potts Model [56][57][58][59]. Although versatile, they describe cellular interactions with scalar energy terms, making it impossible to study forces between cells unless they are amended. Phase-field modeling, based on defining a free energy functional and an order parameter (or several parameters) and is a relatively new approach to cell division and migration [60][61][62][63].
The models that describe cells as discrete entities with individual properties can be further classified as either off-lattice [64] or lattice-based [65] agent models. For example, the former include the models of Rejniak et al. [55,66] and our current model, and the latter cellular automata and Potts models. In off-lattice models cells are typically deformable and forces or energies are used in determining their behavior whereas in lattice-based models updates are done based on pre-determined rules rather than forces and the cellular shapes or topologies are typically fixed. Both approaches have their advantages and caveats. A comprehensive review of agentbased cellular models is provided by Van Liedekerke et al. [67] and a review of mechanobiological aspects and morphology during the whole cell cycle is provided by Clark and Paluch [68].

II. MODEL AND METHODS
We employ the two-dimensional EpiCell2D model to study tissue morphologies. Full details, derivation, and parameter mapping with all parameters (force field) are provided in Ref. [38] and not repeated here. Below, we summarize the salient features and provide the new additions. In the context of the models discussed above, EpiCell2D can be classified as an off-lattice agent-based model; cells are well-defined entities, have individual properties, they can deform, divide and forces are used to update their positions and motions.
In EpiCell2D the cell membrane is discretized as beads connected by elastic bonds of stiffness K spr i to form a closed loop, see Fig. 1. In our simulations, the cells had 76 beads ( Fig. 1 has only eight for clarity). The model has been previously tested [38] using 10-100 beads per cell and it was found that the number of beads has no effect on results if N ≥ 40. The model parameters (details in Ref. [38]) were mapped using a cell diameter of 10 µm, mass of 10 −12 kg [69] and Young's modulus from mitotic HeLa cells [70]. The orientation of the division plane was chosen randomly such that it passed through the center of mass of the cell.
In EpiCell2D, the force field is defined as a sum of intra-cellular, inter-cellular, and cell-medium interactions, see Fig. 1. The intra-cellular terms include the internal pressure ( F P ) and spring forces ( F spr ) which provide a mean-field description for the components that give cells their integrity (cell cortex contractility), that is, the cell membrane and cytoskeleton. The inter-cellular terms define the interactions between the cells. The first of them ( F rep ) is repulsion to prevent cells from penetrating each other and the second term ( F adh ) describes cell-cell adhesion. In real cells adhesion occurs, e.g., due to lipids and different adhesive proteins depending on the cell type. For a recent discussion on the physical aspects of adhesion, see, e.g., Schulter et al. [71] and van Helvert et al. [18]. The final term is medium-cell interaction (friction). With these, the force field can be given in the following general form (1) The functional forms of the terms and all the parameters (14 in total) are given in Mkrtchyan et al. [38] and are not repeated here.
In EpiCell2D cellular growth is controlled by a growth pressure. This is motivated by the fact that cells have mechanisms to control their internal hydrostatic pressure particularly before division [70]. Division is triggered by a threshold in cell area (above which cells divide). We would also like to note that different criteria can be used. As pointed out but Streichan et al. [72] at least area and growth rate are possible criteria for triggering cell division.
Importantly, EpiCell2D allows the aggregate topology (the polygonal distribution) to vary spontaneously [38]; note that cellular topology is fixed, that is, the number of nodes per cell must remain constant upon division, Fig. 1.
One important issue in modeling cells is cell death. The role of cell death in tissue dynamics has been discussed in numerous publications, a recent perspective is given by Green [73]. Different approaches have been taken to include cell death in computational models. In some models that is done as a probability for a cell to disappear [74] or via a tunable cell cycle rate [75]. In the current study, cell death occurs due to local stresses. All of these approaches can be justified as cell death is a very complex and multifaceted phenomenon. Our choice is based the observations of Chen et al. [76] and Streichan et al. [72] that mechanical constraints, in particular cell shape and local stress, are critical factors and determinants for cell cycle and death. Similar mechanisms have also been proposed and analyzed by Shraiman [77].
Previously, we focused on model development, parameter mapping and verification against experiments using only one type of cells [38]. Due to the large parameter space (14 in total), it is not feasible, however, to attempt to map a phase diagram. As a new direction, we extend EpiCell2D for simulations of different cell types using two simple approaches with parameters based on experiments: 1. changing cell stiffness, and 2. changing the cell-cell friction between cell membranes; in EpiCell2D cell membrane and cytoskeleton are treated as a coarse-grained single object.
Modification 1) allows for simulations of different cell types. As mentioned above, cancer cells are typically softer than the matrix cells and softness, or higher malleability, is typically associated with the invasiveness of cancer cells [36,39]. This has recently been challenged by Nguyen [78] et al. who measured Young's modulus of pancreatic cancer cells using different cell lines and found the stiffer (than the matrix cells) cells to be more invasive than the softer cancer cells. Whether this is purely mechanical or due to simultaneously occurring biological processes remains unclear; simulations using EpiCell2D indicate that stiffer cells migrate easier following in the wake of a leader while softer cells collectively evacuate stiffer cells due to aggressive growth. This is in excellent agreement with the findings of Trepat et al. who showed that collective effects are essential cell migration [9]. We will return to this issue in Results as well as in Discussion.
Here, we use two types of cells: (i) Stiff cells with a low growth potential with stiffness K spr 1 = 4 µN/µm. The low growth potential means that the cell membrane is so stiff that the applied pressure is barely enough to grow the cell to a size above the division threshold. (ii) Soft cells with stiffness K spr 2 = 1 µN/µm. These cells have a high growth potential which means that cell membrane stiffness is low and the cell area can easily grow beyond the division size; due to the lower elastic modulus, the growth rate of a soft cell is higher even if the internal pressure is the same as for a stiff cell. Note also that in EpiCell2D the plasma membrane and cytoskeleton are treated as a coarse-grained single object referred to as the membrane.
Modification 2) allows for comparisons of systems of cells with different inter-membrane friction coefficients (term F adh in Eq. 1 and Fig. 1). The importance of cellcell friction in mechanotransduction has recently been reviewed by Angelini et al. [79].
In Eq. 1 inter-membrane friction is modeled via term −µ v ij , where µ is the friction coefficient and v ij is the relative velocity between two membranes tangential to the cell that the bead i belongs to. We compare systems with µ = 0.0 µg/s, that is, cells do not interact very much with their neighbors, and strongly interacting cells with µ = 200.0 µg/s. In the simulations, both open and closed boundaries were used.

III. RESULTS
We focus on systems with two cell types, Stiff and Soft, present simultaneously. Soft cells represent tumor cells based on the fact that cancer cells tend to be softer [36], Young's modulus for cancer cells is typically ≈ 0.5 kPa and whereas for normal cells it often ≈1.0−2.0 kPa (although variations are large, see, e.g., Ref. [18]). In real cancer cells this effect can be enhanced by a lesser number (or lack) of adhesion proteins. However, although cancer cells are generally assumed to grow faster than healthy cells, measurements are not trivial as different metrics can be used [80]. An additional complication is the fact that in healthy tissue growth is regulated whereas cancer cells typically lack such regulation. For a review of properties of cancer cells, see, e.g., Hanahan and Weinberg [81] and for biomechanical aspects Fritsch et al. [82].  The initial setups were created by placing equal proportions of Stiff (red) and Soft (blue) cells randomly, see Figure 2a; the simulations were repeated several times and the results did not depend on the initial conditions. Figure 2b shows the tissue structure at the end of one of the simulations. The soft cancer cells have invaded the space while the stiff cells (red) have been compressed into dendrite-like structures. Another distinct feature is that the cells interpenetrate in the regions marked with light purple and arrows in Fig. 2b. Although this behavior may seem as an artifact, it occurs in diverse systems as has been shown by Eisenhoffer et al. for canine, human and zebrafish epithelial cells [83] and discussed at length by Guillot and Lecuit [15] (see Fig. 2 in Ref. [15]). Figure 2c shows the average contact forces between cells. The white dots show the centers of masses of the stiff cells. The contact force peaks correlate highly with the locations of the stiff cells indicating that the Soft cells overwhelm the Stiff cells as the tissue grows.
These interpenetrating cells is where cell death occurs. Although the dead cells are not physically removed (to keep the simulation code faster), they collapse and occupy only minimal space as defined by the non-overlap of the surface nodes.
Next, we examine if the collapse of stiff cells can be mitigated by making their interactions stronger. This is done by changing the magnitude of inter-membrane friction µ. Since cells need space to grow, they need to slide past each other into empty regions. Higher friction, however, induces jamming and thus reaching the division threshold takes a much longer time. The softer cells will also need to counteract this effect to be able to grow.   Figure 3 shows a similar simulation setup as before, except with two different values of friction µ. Figure 3a shows the initial conditions, and Figs. 3b and c, the final states for µ = 0.0 µg/s and µ = 200.0 µg/s, respectively. The tissues in Fig. 3c grow with open boundaries. In both cases, the simulations were run for a time corresponding to 10 division cycles. As expected, at low inter-membrane friction, there are more cells at the end of the simulation indicating faster growth. The high inter-membrane friction system is more porous with slower growth. Physically, the friction-less system (Fig. 3b) Fig. 5a. The color map shows the strain, that is, relative change in local membrane length as indicated by the color bar (1.0 equals to 100% strain). Soft cells orient themselves in the direction of the largest effective stiffness (the soft-stiff boundary). Since strain is given as stress/Young's modulus, it is evident that the soft cells at the boundary layer manifest force dipoles as predicted by Bischofs and Schwarz [84].
very early stages of development when junctions have not yet developed and the latter (Fig. 3c) to a case when cell adhesion molecules are present.
To investigate further, we analyzed cellular areas (Fig. 4a) and the forces acting on the cells (Fig. 4b) at confluence. Both distributions display lower total number of cells in the case of the high friction tissue. The peak in area distributions is just below 100 µm 2 , which is due to the threshold division area (A div = 100 µm 2 ). Some of the cell areas have grown past this limit as division occurs only at discrete time intervals.
The area distribution for µ = 0.0 µg/s (Fig. 4a) shows a small peak at A ≈ 20 µm 2 due to the higher number of collapsed cells in the low friction system. The large-area peak represents the soft-cell majority with approximately Gaussian shape. This is consistent with the results from simulations of non-dividing soft colloids [22]. For µ = 200.0 µg/s the distribution develops only a single peak.
The distribution of contact forces can perhaps be best understood by picturing two different phases of tissue: 1) A granular phase in which cell density and tissue structure resemble a granular material near jamming. Then the contact forces have the characteristics of granular force-chains with an exponential force distribution at the high-force. 2) A densely packed phase in which cells have overcome jamming, cell density approaches space-filling, forces have been equilibrated and approach a Gaussian. These aspects are discussed in detail below.
All contact forces would relax to toward zero if proliferation stopped and forces were measured after a long enough time. During growth and at low cell-cell friction, the force distribution in Fig. 4b is fairly close to a Gaussian (without negative values). There is however a hint of two separate peaks; when friction is increased, the distribution deviates more from a Gaussian, and its overall structure approaches an exponential function indicating that at zero friction the tissue is closer to the dense-packed phase, while at higher friction the tissue is closer to the granular phase. At high friction, the secondary peak in force distribution develops more clearly. Approximately exponential tails have also been reported in the cell cycle experiments of Trepat et al. [9]. We will return to this issue in Discussion.
Cells at and above the peaks at larger forces arise as a result of remaining mechanical frustrations in particular at stiff/soft cell boundaries. These frustrations are, as expected, more pronounced at high friction, and more relaxed at low friction. Using tissue with only one type of cells present, the high-force tails of the distributions vanish entirely as there are no such boundaries.
As the final case, we study a cluster of soft cells surrounded by stiff cells to model tumor growth in a healthy tissue and with closed boundary conditions to be able to study densely packed tissues. Figures 5a,b show morphologies for µ = 0.0 µg/s and at µ = 200.0 µg/ at confluence. Figure 5 shows that the softer (blue) cells introduced into matrices of stiffer cells grow faster when inter-membrane friction is low; weaker cell-cell interactions provide conditions for easier growth. This suggests that inter-cellular interactions can be an indicator of how well epithelial tissue can damp the growth of rogue cells that have a higher growth potential. The high friction system contains a few pores that have formed due to jamming. Cell sizes are roughly equal within the tumor and inside the matrix but along the tumor boundary, however, the matrix cells are compressed and the tumor cells are enlarged. Further analysis also showed that contact forces are lowest far from the tumour, slightly elevated inside, while largest forces are scattered along the tumor boundary. The matrix cells at the boundary are compressed (i.e. they have a negative membrane strain), while the tumour cells are elongtaed (i.e. they have a positive mebranes strain) near the boundary, Fig. 6. This effect is more pronounced in the low friction case.
To quantify the above, Figs. 7a,b show the contact force distributions. Independent of friction, the distributions are Gaussian at small forces with an exponential tail. The Gaussian part suggests that the force distribution is uniform. The exponential tail, on the other hand, is a signature of disorder and jamming induced by mechanical frustrations. Here, it results from large differences in the cell-cell contact forces along and near the tumor boundary.
Distributions have also been measured for soft colloidal systems under compression. It is well established that the distribution has an exponential tail in the vicinity of the jamming transition, see e.g., [22,85,86] and reference therein. Experiments by Jose et al. for 3dimensional packings of soft colloids also show that the distribution well above the jamming transition becomes Gaussian [86]. As Fig. 7 shows, the exponential tail is present at both zero and high friction. The fact that the cells grow also means that their volumes are not conserved (in contrast to typical colloids). This is also the case for the cells that are being pushed and compressed What is clearly different here is the distribution at low forces: The exponential is preceded by a Gaussian distribution. Gaussian peak has been observed in simulations of soft colloids in two dimensions with zero friction [22].
In contrast, in the three dimensional experiments of Jose et al. the low force part of the distribution remained almost flat except well above jamming transition. In addition, van Eerd et al. have reported faster than exponential decay from their high accuracy Monte Carlo simulations [87] although the deviations can be very hard to detect without high accuracy sampling methods.

IV. DISCUSSION
Using proliferating Madin-Darby canine kidney cells, Trepat et al. [9] studied collective migration of cells and came to the conclusion that long-range traction forces drive collective migration. They measured traction distributions and found that in all cases the distributions had exponential tails. Here, we measured the cell-cell contact forces (Figs. 7) and also found exponential tails. Similarly to Trepat et al. [9], the peak of the distribution appears to be Gaussian. This indicates that our model is capable of capturing the behavior of real systems, and it confirms that mechanical aspects can be modeled using a such a coarse-grained approach.
Fritsch et al. [82] point out the conundrum that although cancers cells are softer than healthy ones, tumors appear as hard lumps. In particular, Fritsch et al. [82] write: "At first sight, cell softening is contradictory to the observation that tumours are rigid masses -a notion borne out by the fact that breast tumours are often felt as lumps. Moreover, this apparent softness of tumours would hinder their invasiveness." As Fig. 2 and 5 show, our model is capable of demonstrating the effect of invasive soft tissue and when growth starts within a matrix, it will become balanced by the pressure from the surrounding tissues (Fig. 5), that is, homeostatic pressure.
In separate studies, Basan et al. [88] and Podewitz et al. [74] investigated the role of homeostatic pressure using experiments, numerical models and mesoscale simulations in which cells are represented by two point-particles with potentials for growth, adhesion and exclusion. Constant temperature and pressure ensemble were used for the molecular dynamics simulations [74]. Due to the representation of cells with central potentials, the approach is not capable of modeling cell-shape changes to any significant degree. As their main conclusion based on their model and experiments, they showed that negative homeostatic pressure (due to compression in the bulk) is both possible and stable. They also showed that homeostatic pressure increases linearly with increasing cell compressibility. Our results agree with this and also demonstrate that the situation is complex. For example, there is destruction of healthy cells at the boundary as a result of concentrated pressure from the cancer cells. These results indicate why cancer cells may be softer.
Perhaps the most interesting comparison is to the theoretical work of Bischofs and Schwarz [84]. In particular, they demonstrated that cells prefer to orient themselves in the direction of the largest effective stiffness. This orientational preference is illustrated in Figures 2b and 5a by blue ovals (solid line). The softer blue cells are facing a stiff environment and as a result, they elongate and organize with their long axis perpendicular to the interface. Figure 6 shows the local strain which clearly manifests that force dipoles are present at the boundary. This is exactly what is predicted by the theory of Bischofs and Schwarz for clamped (stiff) boundaries [84]. In the interior, the cell shapes are more isotropic as is also the case with a free boundary, see the oval with the dotted blue line in Figure 2b; the outermost cells are not in contact with the boundary yet and hence do not feel it. This is also clearly visible in the simulations with open boundaries, Figure 2. In addition, Figures 2b and 5a show that when going away from a boundary where the softer cells have elongated, this order may persist for a few cell lines before disappearing depending on the local environment. Thus, forces are generated away from the immediate boundary. This is in excellent agreement with the experiments of Trepat et al. [9] who argued that force generation is a collective long-range effect rather than requiring 'leader cells'. The figure also shows that when the interface between the hard and soft cells has a more complex shape than a straight line, the cellular shapes and their orientations become more complex. This is also visible in the case of high cell-cell friction.
As already discussed in Introduction, different vertex, Potts and Voronoi-type models have been used to model cell division. Typically, the focus has been on glass-like behavior in (dense) tissues and topological transitions such as the neighbor exchange T1 transition [43,44] and in some cases such models have been combined even with dynamics [44,75]. Since the current model gives spontaneously rise to polygonal distributions (as already reported earlier [38]), it would be possible to use it to study topological changes. That is, however, beyond the scope of the current paper.
Using a vertex model, Bi et al. [89] demonstrated a new rigidity transition at constant density (confluence). As the key parameter, the model has a (cellular) perimeter to area ratio which determines deformability of cells. This parameter has a critical value distinguishing between rigid and fluid-like tissues. The same effect can be observed in our model. Cells with high internal pressure form a rigid tissue, while deflated cells easily change their perimeter area allowing them to deform and flow past each other. Cell friction also plays a signigicant role in such behaviour. Using cellular area/rigidity ratio is not an optimal control parameter here, but as Fig. 7 shows, the amount of rigidity is governed by cell-cell friction.

V. CONCLUSIONS
We use the EpiCell2D model to study systems of cells of two different types in 2D. Cell populations are differentiated by their membrane/cortex stiffness. We showed that this simple difference is enough, provided internal pressure is identical for both, to favor soft cell growth. Even if a few soft cells are surrounded by stiff cells, it is enough for the softer cell to grow rapidly. This effect can be mitigated by a higher interaction strength between cells. Force distributions show similarities to non-proliferating soft colloidal systems. Although not studied in detail here, the model allows for tuning the cell-cell friction, an issue that has recently been raised by Vinuth and Sastry for shear jamming [90].
From the modeling perspective, the EpiCell2D approach appears to be very versatile. As discussed in Introduction and elsewhere in the text, there are several different models around each with their own strengths. Epi-Cell2D is, however, able to capture a very wide range of phenomena without additional modifications. As we have already shown earlier [38], it can reproduce the polygonal cell distributions and mitotic indices observed experimentally in epithelial systems, and as shown here, force distributions, and cellular response to different boundaries agree well with experiments and theory.

MK thanks the Discovery and Canada Research Chairs
Programs of the Natural Sciences and Engineering Research Council of Canada (NSERC) for financial support.