Modeling the impact of granular embedding media, and pulling versus pushing cells on growing cell clones

In this paper, we explore how potential biomechanical influences on cell cycle entrance and cell migration affect the growth dynamics of cell populations. We consider cell populations growing in free, granular and tissue-like environments using a mathematical single-cell-based model. In a free environment we study the effect of pushing movements triggered by proliferation versus active pulling movements of cells stretching cell–cell contacts on the multi-cellular kinetics and the cell population morphotype. By growing cell clones embedded in agarose gel or cells of another type, one can mimic aspects of embedding tissues. We perform simulation studies of cell clones expanding in an environment of granular objects and of chemically inert cells. In certain parameter ranges, we find the formation of invasive fingers reminiscent of viscous fingering. Since the simulation studies are highly computation-time consuming, we mainly study one-cell-thick monolayers and show that for selected parameter settings the results also hold for multi-cellular spheroids. Finally, we compare our model to the experimentally observed growth dynamics of multi-cellular spheroids in agarose gel.


Introduction
Most in vivo tumors expand into the surrounding tissue. To do so, a growing tumor must overcome mechanical barriers by either exerting mechanical stress on the host tissue (the tumor micro-environment) or secreting matrix-degrading enzymes that modify the local environment by cutting fibers such that the macroscopic mechanical resistance is lowered.
Biomechanically induced interactions are increasingly discovered to play an important role in the growth control of tissues and tumors [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. Biomechanical interactions can serve as morphogenetic regulators [3] and permit cells to compare and adjust their growth rate to surrounding cells [5], help cells to choose their orientation according to their environment [6] or to adjust their growth rate to a local geometric property such as tissue curvature; this has led to the conjecture that sometimes functions follow form [7]. Morphogenetic strain rates, cell shape change and intercalation have been demonstrated to be interlinked [18].
Tumors have been observed to grow slower when their embedding medium is stiffer, and their shape may reflect the geometrical constraints of their environment [1,14]. Chaplain et al [19] studied this situation in a multi-scale, multi-phase model considering tumor cells, the extracellular matrix (ECM) and the host tissue including molecular factors controlling the growth. Momentum balance was mimicked by an approach of the Darcy type, hence assuming the tumor behaves as a porous medium. They find an exponential growth for small tumors saturating at large stress. Helmlinger et al [1] showed the influence of the biomechanical properties of the tumor micro-environment by growing tumor spheroids embedded in agarose gels of different levels of concentration and thus rigidity. Regardless of the host species and the tissue of origin, increased mechanical stress led to significantly decreased maximal tumor spheroid sizes, which for example in human colon carcinoma decreased from a diameter of 400 µm (in 0.5% agarose) to 50 µm (in more rigid 1% agarose). Chen et al [20] studied the saturation of growth within a multi-phase model of a growing tumor in external poroelastic material in spherical coordinates to explain the Helmlinger experiment. The tumor was modelled as a two-phase solid-liquid material similar to a porous material where the volume fraction of the interstitial liquid was assumed to decrease upon compression of the tumor. They found saturation after a certain time, as it was observed in the experiments and linear growth if the poroelastic gel was removed. More recently, another group demonstrated that also the shape of tumor spheroids is dictated by the shape of the solid stress field by using agarose gels and co-embedding fluorescent micro-beads [14]. Further analysis revealed that one reason for this observation was suppression of proliferation and induction of apoptosis in regions of high mechanical stress. A number of hypotheses have been proposed on the question of how and to what extent the physical interaction of tumor and host tissue influences tumor morphology and growth kinetics during cancerogenesis [21,22]. Nevertheless, many aspects of these complex interactions remain to be experimentally elucidated [23].
Comparisons of experiments with mathematical models have shown that the increase of the cell population diameter as well as of the cell proliferation pattern, in both growing monolayers [24,25] and growing multi-cellular spheroids [24][25][26][27], could largely be explained by a biomechanical form of contact inhibition, controlled by a force threshold, a pressure threshold or a deformation threshold, above which cells become quiescent. For example, a careful analysis of EMT6/Ro multi-cellular spheroids revealed that the size of a multi-cellular spheroid is almost unaffected even when the external glucose concentration is varied by a factor of 20 from 0.8 to 16.5 mM, while at the same time the cell population sizes varied significantly. The glucose medium concentration only affected the size of the necrotic core, but not the tumor size. Even the growth kinetics of well-vascularized xenografts of human NIH3T3 cells subcutaneously injected into mice could be explained by the same form of biomechanical growth control [28].
These examples clearly indicate that biomechanical effects have a potentially important role in growth and form.
However, unlike physical particles, biological cells have the capability to change their physical properties by intracellular control processes and thereby modify their phenotype. A prominent example is the epithelial-mesenchymal transition (EMT; see, e.g. [29]) in cancer where the tissue phenotype changes from an epithelium-type phenotype to a mesenchymal phenotype by active intracellular regulation facilitating cell detachment [30][31][32]. Another example is that cells are able to enter blood vessels by a process called intravasation involving cancer-cell-induced src activation, leading to the degradation of contacts among endothelial cells, and thereby decrease the mechanical resistance of the blood vessel walls against cancer cell invasion [33]. Cell birth and death may facilitate cell movement in tissues which on short time scales appear as elastic material [16].
These examples suggest that there is a complex interplay of mechanical and active regulative components during different stages of tumor growth, development and invasion. To find out the contributions of different active versus passive components it is important to examine how far the cancer phenotype can be explained by physical effects alone. Otherwise, morphotypes may be erroneously attributed to active regulatory behavior of cancer cells. Invasive fingers, for example, can also be triggered purely by physical effects as known from viscous fingering [34] or other instabilities known from non-equilibrium physics (see, e.g. [35,36]). Examples are the Mullins-Sekerka instability of a growing crystal in a supercooled melt driven by undercooling at the solid-liquid interface (see, e.g. [37,38]), and an elasticity-driven growth instability as a consequence of an applied stress [39], a buckling instability driven by the competition of cell proliferation and stabilizing effects such as bending or shear (see, e.g. [8,[40][41][42]), or an undulation stability as it may occur at the epithelial stromal 4 interface again driven by cell proliferation ( [43]; see also text below). Different morphotypes have observed bacterial growth which could be explained by different types of instabilities [44].
In this paper, we explore how far physical properties of the embedding medium of a growing cell clone affect the emerging spatial-temporal growth pattern. For this purpose we consider two situations, namely a growing clone of cells in a granular medium and a growing clone in an environment of inert cells. Within our model the main difference between the two situations is that for cells we assume an active component of the movement, the cell micromotility, whereas for the granular objects we do not. We do not take into account the active regulation mechanisms for cells as reported in, for example [30,33]. We mimic each biological cell individually within the framework of 'agent-based models'. Our model is based on an earlier approach where cells have been mimicked as adhesive homogeneous isotropic elastic objects capable of active migration, growth and division of particles. This model has permitted us to explain the growth pattern of monolayers and multi-cellular spheroids [24][25][26][27] and predicted a subsequently validated order principle in the regenerating liver [45].
We consider a growing monolayer as the reference situation. Recent experiments in [46] suggest, for expanding monolayers of MDCK (Madin-Darby canine kidney) cells, an active migration into the direction of the outwards-pointing normal of the cell population border. The cell density as well as the cell traction has been reported to increase from the border to the middle of the monolayer. This suggests a pulling type of movement-perhaps in combination with an increasing contraction of the cell projection area from the border to the center [47]. In a recent paper [45], we could show that active directed cell movement is necessary for explaining the regeneration of a CCl 4 -induced tissue lesion in mouse liver. Passive movement triggered by cell proliferation alone was insufficient. However, in both cases the cell populations were not constrained by the strong resistance of other cells or a capsule.
In contrast, in liver regeneration after partial hepatectomy, the growth of a cell population must occur against the resistance of a capsule enclosing the liver lobes. Hence, growth in such an environment must occur against mechanical resistance.
The above-reported observations [45,46] suggest that both types of movement may exist: (i) pushing-type fronts if cell division leads to local mechanical stress pushing neighbor cells into the direction of minimal mechanical compression. (ii) Pulling-type movement may be triggered by the emergence of a border, resulting in active cell migration into the free zone. For this reason we compare the effect of a biased micro-motility with the existence of an un-biased micro-motility. In [24][25][26][27] we had only considered the case of an un-biased (i.e. uniformly distributed) micro-motility.
We first consider growing monolayers with a free border and briefly discuss cell movement triggered either by cell division ('pushing-type movement') or by an active outwards 'pullingtype movement' of cells.
For growing clones embedded in tissues or tissue-like media, a free cell population border does not exist. Instead, an external medium exerts a pressure on the expanding clone resembling a carcinoma located in an organ (such as e.g. a liver carcinoma) or a population of hepatocytes enclosed by a capsule rather than a tumor expanding in free space or in a soft tissue environment. In such a situation a pulling-like movement of the tumor cells seems unlikely to occur prior to the EMT.
So, in summary, the main purpose of this paper is to explore the growth and invasion pattern of biological cells in a minimum model mimicking cells as homogeneous isotropic elastic, sticky particles capable of self-reproduction and migration.

5
In section 2, we present the basic model and its variants used in this paper. Then we consider pushing-and pulling-type expansion of monolayers in a liquid environment. We then study monolayers expanding in external media of different properties, namely granular objects and cells of different types. We mainly consider simulations in two dimensions noting, however, that in three dimensions the results should be equivalent as long as glucose or oxygen do not become limiting, which seems to hold over a wide range of glucose and oxygen concentrations [24]. We demonstrate this by simulations of multi-cellular spheroids embedded in granular objects in three dimensions.

The model
We model each individual cell by an isotropic, elastic and adhesive object. The model cells are capable of active migration, growth and division and are defined by cell-kinetic, biophysical and cell-biological parameters that can be experimentally measured. Model cells interact with each other, with other cells and with the underlying substrate.

Cell shape and cell-cell interaction
In particular, epithelial-derived tumor cells and epithelial cells may reveal a largely spherical shape in isolation as shown in [45] (see supporting information in that reference). In culture at high cell densities they adopt a polygonal shape.
When a cell gets into touch with another cell it can form an adhesive cell-cell contact. With decreasing distance between cells (e.g. upon compression), the contact area between them, and with it the number of adhesive bonds, increases, resulting in increasing attractive interaction. On the other hand, if cells in isolation are spherical, an increasing contact area is accompanied by an increasing deformation, which results in a repulsive interaction. Hence, the preferred distance of a pair of interacting cells is where adhesive and repulsive forces balance. If adhering cells are pulled from each other they reveal a hysteresis effect: they elongate and stick beyond the distance at which they had come into contact. At the rupture point they still have a finite interaction area. These observations had led the authors of [24] to approximate the interaction forces between cells by the Johnson-Kendall-Roberts (JKR) model. The same model has independently been shown in [48] to apply to pairwise interacting S180 cells using micro-pipette experiments. The JKR model mimics the force between a pair of interacting sticky homogeneous isotropic spherical objects. It directly includes adhesion and relates the contact area to the elastic material properties and the adhesion strength and therefore gives in the limit studied in [48] a proper quantitative description.
where d i j is the distance between the centers of two interacting spheres i and j that is calculated from two implicit equations [49]: where a is the contact radius. The effective radiusR is defined byR −1 = R −1 i + R −1 j , where R i is the radius of cell i. d i j = R i + R j − δ is the distance between the centers of model cell The force F i j as a function of the distance d i j for two interacting spheroidal cells i, j according to the JKR model for three different elastic moduli, E i = E j = 300, 450 and 1000 Pa. The dashed red lines illustrate different cellcycle entrance control conditions that we tested. Rule 1 assumes a compression threshold, assuming that a cell i can only enter the cell cycle if its distance d i j to its neighbor cell j remains larger than d th . This rule introduced in [24] is not considered here. Rule 2 assumes that cells can only enter the cell cycle if the pressure-like measure p i = j F i j /A i j (A i j is the contact area) is below a threshold p th > 0 (the sum is over all neighbors of i; for details see text). Rule 3 assumes that only cells that experience a negative pressure are able to enter the cell cycle. Rule 4 assumes that only if the cell-cell contact to at least one neighbor cell is stretched is it able to enter the cell cycle. Inset: cell division in the program. Cells double their volume in interphase and deform into a dumbbell in mitosis. (ν i = ν j = 0.4, ς m = 10 15 m −2 .) i and cell j, where δ = δ i + δ j is the sum of the deformations of each cell (upon compression it is the overlap of the two spheres) along the axis linking the centers of these cells.Ẽ i j is the composite Young's modulus defined byẼ −1 We approximatê γ ≈ ς m W s . m is the density of cell-cell adhesion molecules in the cell membrane, W s the energy of a single bond. If the density of cell adhesion molecules would differ in two cell i and j in contact, then m = min( i , j ). In our simulations we assume the same cell adhesion molecule density for each cell. Equation (2) has to be solved implicitly to obtain a(F JKR i j ). The value of a is then inserted into equation (1) to give δ(a), (1) and (2) are not amenable to an explicit solution ifγ > 0. We fitted the resulting plot by a polynomial of degree three and used this polynomial in the simulations. In this way, were able to avoid a prolongation of the simulation times, which would occur if the JKR model equations were solved numerically during the simulation runs with the multi-cellular growth model, and to keep the advantages of 7 the JKR force outlined above: The force vector denoting the force of cell j on cell i results from Certainly, although the JKR-force model seems to be a plausible approximation for the interaction between cells, it neglects inhomogeneities of the cell material, deviations of the cell shape from a spherical shape prior to deformation and plastic deformations that are expected to occur on longer time scales as a consequence of re-organization processes of the cytoskeleton. Moreover, the force calculated by the JKR model is a two-body force and neglects forces that may emerge from multi-body interactions. However, in [26] we have shown that multicellular growth phenomena in liquid suspension or monoclonal monolayers show only moderate differences if the interaction force model is varied, for example, by using an extended Hertz model or a linear spring-like interaction force between cells. From this we would not expect any qualitative impact from replacing the JKR-force model by a more detailed description of interaction forces even though we consider it as an important challenge in the future, as far as quantitative experiments permit a more detailed measurement of the forces to allow increasingly quantitative model predictions. On the other hand, we consider it useful to take the JKR force instead of a cruder interaction force, as it is parameterized in measurable quantities and is expected to be correct for sufficiently fast deformations.
We here consider perfect one-cell-thick monolayers. They can experimentally be realized in two ways. As we have shown earlier [25], cells seeded on a flat substrate do not detach from the substrate if their growth is contact inhibited and if at the same time the cell-substrate forces are sufficiently strong to ensure that the cells enter the cell cycle rest (G 0 -) phase before they would detach from the substrate. Therefore, contact-inhibited cell lines with sufficiently strong substrate adhesion would form a one-cell-thick monolayer. Another experimental possibility to ensure the formation of a perfect one-cell-thick monolayer might be to constrain the monolayer from above so that cells cannot leave the one-cell-thick layer configuration even if they are not contact-inhibited. We use the JKR force also for the adhesion between cells and the underlying substrate.

Cells and the enclosing medium
The expanding cells can grow and divide as explained below. We consider two types of environmental objects: non-dividing cells or inert granular objects. We label the cells of the expanding clone 'A' and the enclosing objects 'B'. Below we refer to the cells of the expanding clone as 'A-cells' and to objects (cells or granular objects) of the enclosing medium as 'Bobjects'.

Cell growth and division
During G 1 -, S-and G 2 -phase we assume that a cell increases its volume by increasing the radius R in small steps R R until it has doubled its initial 'intrinsic' volume to V DIV = 2V INIT , where V INIT is its volume immediately after cell division (figure 1). Here, the intrinsic volume V i of a model cell i is approximated by the model cell i deforms into a dumb-bell at constant volume. Then it divides into two daughter cells of radius R. We assume that a cell i in G 1 switches into a rest phase G 0 if the measure p tot i , defined by overcomes a threshold value p 0 . p tot i is a pressure-like measure (hereinafter referred to as 'pressure').
Here, F AA i j is the JKR force between A-cell i and a neighbor j of type A, F AB i j is that between A-cell i and a neighbor j of type B, F AS iS is that between A-cell i and the flat substrate, and A i j are the corresponding contact areas calculated with the JKR model. In the first sum of the above relation, j runs over the nearest neighbors of type A, and in the second sum, it runs over the nearest neighbors of type B. The pressure is calculated directly after the division of a cell. If p tot i p 0 , then the cell i enters the cell cycle and passes through the whole cell cycle again. If the forces on cell i were isotropically distributed, acting uniformly on the whole cell surface, and if the adhesion strength were zero, the measure is identical to the definition of hydrostatic pressure on the cell.
The orientation of cell division in the monolayer simulations is assumed to be random uniformly distributed parallel to the two-dimensional (2D) substrate. During the mitosis phase we assume that cells adopt for a short time period a dumb-bell shape. Forces on cells at this stage can lead to a torque. The orientation can change as a consequence of the torque which we take into account by orientation changes. For simplicity we modeled these by energy minimization (the Metropolis algorithm) instead of numerically integrating equations for the torques. Energy minimization provides an alternative to a forced-based single-cell dynamics [26]. Within each time interval t for each cell a rotation trial around three space-fixed axes by angles δβ i with i = 1, 2, 3, δβ i ∈ [0, δβ max ), with δβ max π/2 was performed, using the algorithm of Barker and Watts as explained in [26]. The probability P Rot that the rotation trial is accepted was calculated by P Rot = min (1, e is the total potential energy after the orientation change and V (t) the total potential energy before the orientation change. F T is a reference energy [24]. The energy and force can directly be linked by F i j = −∂ V i j /∂r i . The Metropolis algorithm ensures that orientation changes that lead to a decrease of the energy of the multi-cellular configuration are always accepted, while those which lead to an energy increase are only accepted with probability e − V /F T . For spherical cells we do not consider orientation changes since they do not change the total energy of the multi-cellular configuration. For the monolayers the cell-substrate forces were assumed to be strong enough to prevent cells from moving out of the layer.

Cell migration
First, the modeling framework for monolayers is developed, where cells can crawl on a 2D substrate. After this, we discuss the difference between monolayers and multi-cellular spheroids, which are 3D cell aggregates.
Monolayers. In the absence of chemotactic signals, isolated cells in suspension or culture medium have been observed to perform an active random-walk-like movement in monolayers [50,51]. We assume that the random component of the active cell motion, the 9 cell micro-motility, can be characterized by the cell diffusion constant D of isolated randomly moving cells. We used D 0 ≈ 10 −12 cm 2 s −1 as the reference value. More motile cells are characterized by a larger D. In our model, D is also used to quantify the micro-motility of cells in contact with other cells.
While in mechanical contact with other cells, proliferating cells exert a pressure on their neighbors. The neighboring cells try to escape this pressure by moving against the friction caused by the other neighbor cells and extracellular material, for example ECM [52]. The movement could be partly passive, due to pushing, and partly active [9,46], if cells migrate into the direction in which they escape the mechanical (or morphogen) stimulus. We simulate cell migration as a friction-dominated over-damped motion with a stochastic contribution by a stochastic equation of motion for each cell [26,45,53]. Schienbein et al [50] have demonstrated that the migration of an ensemble of isolated cells in culture subject to an external force can be mimicked by a Fokker-Planck equation for the single-particle distribution function having a drift and a diffusion term.
Dickinson and Tranquillo [54,55] have developed equations of motions for different types of cell movement from cell receptor dynamics and obtained a stochastic dynamics reminiscent of the approach of [50]. More recently, in [46], in monolayers of MDCK cells a net traction of cells perpendicular to the border of the layer was observed, suggesting an active pulling-like movement to increase the spread of the monolayer. Our equation of motion tries to capture the main aspects of these observations.
The equation of motion of cell i of the expanding cell type (denoted by the superscript A) is determined by (if X = A) or (ii) object j of type B (if X = B) or (iii) the substrate S (if X = S). Generally, the friction may be decomposed into a perpendicular and a parallel component, so that Here, u i j = (r j − r i )/|r j − r i | with r i denoting the center of cell i of type A (or if cell i is in mitosis and has a dumb-bell shape, the position of the closest sphere of the dumb-bell). r j denotes the center of object X ∈ {A, B, S} i.e. a cell of type A, an object of type B, or the substrate S. '⊗' denotes the dyadic product. F XY i j denotes the JKR force between cells i and j of type A (X = Y = A), objects i and j of type B (X = Y = B) or between cell i of type A and object j of type B (X = A, Y = B or X = B, Y = A). F AS i is the JKR-force between cell i of type A and the substrate S. The JKR-forces between cells A and the substrate and between objects B and the substrate are assumed to be large enough to prevent detachment and are not explicitly modeled here. We assume the substrate that the cell crawls on to be flat and denote it by a sphere of infinite radius. I is the unity matrix. ( u ⊗ u )v i is the projection of v i on u. This can be immediately seen as, after 10 some re-writing, ). If, for example, X = S, then u ik is the unit vector of the line connecting the center of the sphere of infinite radius representing the substrate k = s and the center of cell i of type A. (u is ⊗ u is )v i then denotes the projection of the velocity component on the connecting line, i.e. perpendicular to the substrate, and (I − u is ⊗ u is )v i is the velocity component parallel to the substrate. Accordingly, γ AS denotes the tangential friction coefficient of a cell of type A with the substrate. Correspondingly, (u is ⊗ u is )v i would contribute to a deformation of a sphere against inner cell friction characterized by γ AS ⊥ . In general, the friction coefficients γ XY ⊥ , γ XY may differ for cells of type A and objects of type B. We assume here that they are the same i.e. γ AA ⊥ = γ B B ⊥ and γ AA = γ BB . We further assume throughout this paper that γ XY ⊥ = 0 and that γ XY ∝ A i j for two interacting objects i, j of type X, Y (XY = AA, AB, BB, AS or BS). As we set γ XY ⊥ = 0, we in the following use the simplified notion γ XY ≡ γ XY || . Assuming that γ XY ∝ A i j assumes that the density of surface adhesion molecules responsible for the strength of the friction force is nearly constant. If the number of surface molecules were constant, then γ XY should be independent of A i j .
F AA i j is the JKR force the cell i of type A experiences from neighbor cell j of type A, F AB i j is the force that A-cell i experiences from a neighbor j of type B, and F AS i j summarizes the force A-cell i experiences from the flat substrate, ECM and a possible plate or mesh preventing movements perpendicular to the substrate. In monolayers, cells are usually attached to a substrate, usually the culture plate. F active,A i is then the active migration force of cell i of type A from crawling on the substrate.
The equation of motion of the external objects (denoted by the superscript B), either external medium cells or granular particles, is mimicked by Equivalently, F BB i j is the JKR force B-cell i experiences from neighbor cell j of type B, and F BS i j is the JKR force of B-cell i by the flat substrate, extracellular material and a constraining plate or mesh.
If the B-type objects are granular particles, F active,B i = 0. For the other cases, the active migration force of cell B is as defined below.
Beysens et al [56] have proposed an analogue of an 'Einstein relation', linking the cell diffusion constant with an effective temperature. Note, however, that the physics for cells and colloidal particles is very different. For cells, the effective temperature parameter is controlled by the cell surface receptor dynamics and not, as for Brownian particles, by collisions with smaller solvent particles. If the cell motion occurs on a flat substrate, then the random micromotility is not isotropic in three dimensions but needs to be split into components parallel and perpendicular to the substrate.
For A-type cells we consider two types of micro-motility which are presented below for a one-cell-thick monolayer.
Case 1: Unbiased micro-motility. We assume that the micro-motility is isotropic within the plane of the substrate. In this case, is the random force component parallel to the substrate. We assume the perpendicular components to be zero. For the parallel components, we assume that F active,A i; in formal analogy with colloidal particles. Moreover, the random force components of different cells as well as between different vector components of the same cell are uncorrelated. However, as discussed above the Einstein relation in the case of cells only represents a formal analogy because unlike colloidal particles cells can control their migratory activity on their own. Hence the amplitude of the force autocorrelation function generally is a function of cell internal processes or states.
Case 2: Biased micro-motility. We assume that the active migration is biased into the direction of minimum stress exerted by A-cells if no embedding medium exists. We found that this assumption was able to explain the regeneration of the necrotic zone that emerges after intoxication by CCl 4 in liver lobules [45]. Moreover, recent monolayer experiments with MDCK cells suggest that cells may actively move towards free spaces and pull other cells behind [46]. We mimic aspects of this type of micro-motility in two alternative ways: 1. By the selection of random force contributions into directions of decreasing local pressure: is the Heaviside function which is one if x 0 and is zero otherwise and is a measure of the stress of A-type cell i by external forces not taking into account the adhesion of cells. The term mimics the case when cells can distinguish neighbor cells of their own type from neighbor cells of a different type and move in order to relax the compression within cells of their own type. If we would consider in the calculation of p A i also forces between A-type cells and B-type objects, then the effect at the interface between A-type cells and the B-type object would become negligibly small.
2. By the addition of a constant force term directed into the local outward-normal of the border of A-cells. In this case, is the random component of the active force obeying the same relations as in case 1. Some authors (see, e.g., [57]) assume that migrating cells try to adopt a certain target speed v 0 . This can be taken into account by replacing in the first term of the equation of motion We will not consider this case here but note that it can be easily included in our concept.

12
For B-type objects we consider two cases: If B-objects are considered to be cells, then we assume that the active term F active,B i (t) is mimicked in the same way as for A-type cells in case 1, i.e. as an unbiased active random movement. If B-objects represent granular particles, we assume that the active force is zero, Multi-cellular spheroids. In the last part of the paper, we consider multi-cellular spheroids, i.e. 3D cell aggregates. The tissue architecture of (3D) multi-cellular spheroids is usually closer to the in vivo situation and shows the formation of extensive ECM (see, e.g. [58]). For example, staining for collagen V in SK-MES-1 cells, a non-small-cell lung cancer cell line shows massive ECM presence [59]. If this 3D network of ECM is stiff enough, it may allow cells to anchor and move as described in [51]. In this case the interpretation of the terms remains the same as that in the monolayer situation. The factor '4' in the autocorrelation function is then replaced by '6' (= 2d; d is the space dimension) and the diffusion is assumed to be isotropic (we use diffusion constant D). However, some multi-cellular spheroids show only a very sparse matrix. In this case, the main component of active migration occurs from forces exerted from a cell on its neighbor cells (i.e. not on the ECM anymore). Those forces are balanced by the respective reaction forces of the neighbors on that cell. Stochastic forces from one cell on its neighbor cells may also emerge from fluctuations of the cells' membranes in which case also momentum is conserved; hence, the random force terms for neighbor cells have to balance each other. Taking this into account would modify the equations of motion slightly but the results would not be expected to show any major changes for the situations presented in this paper. (i) We are here interested in pattern formation on long time scales compared to the characteristic times of the stochastic fluctuation terms representing the cell micro-motility and, as the spheroids are embedded in a dense arrangement of granular objects and not swimming in liquid suspension, cells do not detach and move away as a consequence of the random movement. (ii) For multi-cellular spheroids we consider below pushing-type movement emerging from cell proliferation. Here active directed movements are not considered and small fluctuation forces do not play a role.
That the precise form of the noise terms representing the micro-motility does not play a role is supported by the simulation results at varying micro-motility of the B-objects shown in figures 6(A)-(D), suggesting that the random component in the movement of cells has statistically no effect on the formed pattern for the situations considered in this paper. However, it might be worth noting that the movement of the cells in liquid suspension could be mimicked in the same framework: cell movements due to collisions with fluid particles generate a random Brownian movement component, usually with smaller amplitude than for active random movement, so in this case the equations of motions are the same as explained above but with a smaller diffusion constant.
Overview of simulations. In our simulations below, we consider one-layer-thick monolayers. For A-type cells the initial number at the start of the simulation at t = 0 is always N (t = 0) = 1.
In the first part we investigate the effect of pushing/pulling-type movement on the growth kinetics of a monolayer not embedded in another external material of B-type objects.
In the second part, we consider an expanding monolayer in an embedding material of B-type objects (granular objects or cells). We vary the following parameters: (i) motility (diffusion constant) and mobility (friction constant) of the embedding medium, (ii) density of the embedding objects, (iii) elasticity of the embedding objects, (iv) adhesion strength of the embedding objects, (v) object size of the embedding objects and (vi) cytolysis versus no cytolysis of apoptotic cells of the expanding clone. The reference parameters used in the simulations are given in table 1. In the third part, we consider 3D multi-cellular spheroids in an embedding material of granular objects and compare our simulation results directly with the findings of Helmlinger et al [1] and Galle et al [60]. We from now on always use D to denote the diffusion constant. For monolayers, this corresponds to using the symbol D instead of D .
To quantify our simulation results, we measure the cell population size N (t), and the radius of gyration, where R CM = 1 N i r i is the position of the center of mass. The sums run over the expanding (the growing) clone. If the monolayer is compact with a circular border, then R = √ 2R gyr , where R is the radius of the monolayer. For this reason, we usually display R or the diameter, L = 2R, calculated from the radius of gyration. In addition, we store the full spatial-temporal configuration of cells and objects at different times during the growth process.

Reference situation: growing monolayers without an embedding medium
In [46], the authors found an active migration of cells towards the free edge of an expanding monolayer culture. Hoehme et al [45] found that an active migration of cells towards a druginduced necrotic zone is necessary for explaining regeneration of liver after toxic damage. Here, D A = 3.1D 0 . Insets: if a cell actively moves into the free space around the monolayer, the cell density in the surface region is smaller (upper left inset) than if active migration is purely random (lower right inset). Active migration accelerates monolayer expansion. In the lower half, white cells are proliferating and gray cells are in G 0 phase. Figure 2 shows the results of a biased and an unbiased micro-motility in a growing onecell-thick dense monolayer. The strength of the bias has been varied up to the size at which cells detach.
As B-type objects are not considered in this part, and as F AS iS is approximately constant by construction, the variable part of the pressure is given by The results show an acceleration of growth which can be explained by an increase of the proliferating rim by the biased cell micro-motility. Due to the active movement towards the free edge of the monolayer the cell-cell contacts are stretched. This leads to a reduction of cell compression and thereby facilitates reentrance of cells into the cell cycle. Moreover, the cell density close to the border is smaller for biased than for unbiased micro-motility. We conclude that a bias in the micro-motility promotes the growth of a monolayer. We expect the same result for multi-cellular spheroids. Transferred to the case of a tumor in the soft tissue environment, any morphogen that biases micro-motility or an increase of its amplitude or both will promote the speed of growth of the expanding tumor. Besides this acceleration effect, we observe that biased micro-motility reduces the random motion of cells inside monolayers because cells in almost isotropic multi-cellular environments then have a smaller tendency to move. On the other hand, the micro-motility is larger close to the border.
So far we have used model rule 2 of figure 1 to mimic the cell cycle entrance control. Rule 1 of figure 1 has been tested in [24]. In figure 3, we test different hypotheses for cell cycle entrance, namely that the entrance into the cell cycle is triggered if p tot i < 0, i.e. if normal stress − p tot i on the cell is positive (rule 3 in figure 1), and that cell cycle entrance only occurs if a cell is under stretch (rule 4 in figure 1). Note that for rule 4 it is sufficient if at least a single contact is stretched (i.e. d i j − (R i + R j ) > 0 for at least one j ∈ { j}) even if other contacts are under compression. It is not necessary that the cell is overall stretched into all directions (i.e. We again find that an increase of the active migration pulling force accelerates the growth. Above a certain pulling force (see the curve for 3 nN), cells detach and actively migrate into the environment. This is reflected in the radius of gyration versus time curve by the characteristic bump at about 2-3 days, indicating free migration (arrow in figure 3). However, the condition that cell entrance is possible only if a cell is under tension (positive normal stress) implies smaller growth velocities as for active directed migration forces and a positive pressure threshold (compression; implying negative normal stress) for the cell cycle entrance.
Interestingly, the largest and smallest growth speeds are observed for rule 3. The growth speed is largest for pulling forces of 1 nN or more and smallest in the absence of pulling forces. For rule 4, stochastic force fluctuations are able to trigger the overstretching of a single cell-cell contact for a moderate pulling force, resulting in a larger monolayer expansion speed than for rule 3.
In summary, the growth curve without any embedding medium shows that: (i) a directed force into the direction of the outwards-pointing interface normal leads to acceleration of growth compared to random uniformly directed micro-motility, (ii) a growth control mechanism permitting growth entrance only if a cell experiences a negative force slows down the growth speed for small pulling forces but accelerates it for large pulling forces, and (iii) a growth control mechanism that permits only cell cycle entrance for stretched cells leads to a smaller (larger) growth speed than for negative force control for large (small) pulling forces.

Growing monolayers in an embedding medium
In the following, the term 'embedding material' encompasses both the B-type cells and the B-type non-cellular objects. A distinction is made only if necessary.
Compared to earlier studies (see, e.g. [24,26]), the introduction of an embedding material is an important step towards a modeling of the in vivo situation. We systematically modified the properties of the embedding material by varying the motility, initial density, elasticity, adhesion properties and average cell size of the embedding objects (B-type). We assumed the growing cells and the embedding material to be differentially adhesive with ζ m = ζ AA = 10 15 m −2 , where ζ AA denotes the density of adhesion molecules in the case of an interaction between two cells of  (Continued) increasing outwards-directed active migration force, the expansion velocity increases and for an active directed migration force of about 1 nN it becomes larger than the reference curve without the directed migration component. Beyond about 3 nN the cells detach and move freely. Note that the expansion speed for a directed micro-motility force of 0.8 nN and positive pressure threshold (= negative normal stress; pink curve) is larger than that for a directed micro-motility force of 1 nN if the cell cycle entrance is triggered by a negative pressure (= positive normal stress). The arrow shows the characteristic bump if cells detach. The bump denotes a regime where R gyr ∝ √ t and is characteristic of the free random movement component (see also [63]). type A. For B-objects we assume in most cases that ζ BB = ζ AA so that the strength of adhesion between B-objects is the same as the strength of adhesion between A-cells. Furthermore, unless stated otherwise, we assume that A-cells and B-objects do not adhere, i.e. ζ AB = 0 m −2 . The elastic moduli of A-cells and B-objects are denoted by E A and E B , respectively, and the diffusion constants for A-cells are denoted by D A and for B-cells by D B . Whereas D 0 = 10 −12 cm 2 s −1 denotes the reference diffusion constant, γ 0 denotes the reference effective friction constant for cells with the substrate.
We started our simulations with a single cell of type A embedded in about 10 4 -10 6 particles of the granular medium or model cells of type B. We used the number of embedding objects to control the initial density of the environmental objects. All elements of the embedding material were randomly arranged (isotropic homogeneously distributed) within a circular environment (figure 4) that was enclosed by a circular impermeable wall mimicking a virtual Petri dish. The initial placement of the B-objects on a regular square lattice was found to lead to symmetry artifacts in the emerging growth pattern. The shape of the impermeable wall fed back to the final shape of the growing aggregate (this was also experimentally observed in [1] for cell populations growing in a tube filled with agarose gel, see also figure 15), which is why we chose a circular border. The total area was held constant and was the same in all simulations. During the simulations, the growing population of cells of type A pushed away and compressed the surrounding objects. Figure 4 shows that a directed active migration component according to case 2 does not affect the growth speed of the expanding cell population as long as the directed force component is insufficient to cause detachment of cells from the expanding monolayer cell population.
For this reason, below we only consider pushing-type movement. We also consider the case when cell cycle entrance occurs only if p tot i > 0. First we study growth without apoptosis and then include apoptosis with subsequent cytolysis. Figure 5 shows the growth scenario with the reference parameter set.

Variation of the friction and motility of embedding objects.
As a first step we analyzed the impact of varying the friction and micro-motility of the embedding objects of type B on the morphology of the expanding cell clone A.  Note there is no difference as long as the active force is below the rupture force (blue curve with circle symbols and red curve with triangle symbols). Only if the active force is above the rupture force (green curve with square symbols and orange curve with diamond symbols), is the growth moderately accelerated and the saturation diameter elevated. As A-type cells are in an environment of B-type objects, not all cells are able to detach even if the directed active movement force amplitude F 0 is above the rupture force. Rather, the probability of rupture increases, leading to an increase of growth speed with increasing F 0 . However, the differences in growth speed are far below that of a free monolayer (compare figures 2 and 3). In all simulations in this figure, strictly diffusion (fluctuation) with dissipation (friction) by a (thermal) temperature does not exist even though a formal analogy has been proposed [56]. Variations of the diffusion constant reflecting the micro-motility, an active motion component of the cell, may be caused by, for example, molecular changes in the complex coordinated cytoskeletal actions that are required for active cell movement [64]. Cell-substrate friction may be influenced by modification of the ECM, the liquid medium, or the covering of the dish surface. For objects B of the embedding tissue, either the cell diffusion constant D B or the friction γ BS of cell type B with the substrate is modified.
While we found that a variation of D B has only a negligible impact on the shape of the expanding multi-cellular population (figures 6(A)-(D)), we found that the variation of γ BS leads to major changes in both growth velocity and population morphology (figures 6(E)-(H); At sufficiently high friction between B-objects and medium, the border of the A-clone forms dendrites. With increasing medium friction, the growth velocity of the expanding A-cell clone decreases (compare, for example, lines A and B in figure 11), and the critical wavelength (the size at which a domain splits into two branches so that dendrites form) decreases (figures 6(E)-(H)). Effectively such an increase of medium friction of the host tissue increases the viscosity of the tumor micro-environment, which is reminiscent of a Saffman-Taylor-like instability [65], a viscous fingering.
The additional cell mass generated by cell proliferation generates a pressure gradient moving one 'fluid' (the expanding cell clone) against another 'fluid' (the extracellular objects). If the friction of the external objects with the dish (or a medium between the objects) is too large, they cannot move away with the same speed as the growth front. So cells of the growing clone expand and form fingers into small holes emerging by chance. At the tip of the dendritic fingers the local stress is smaller than that in the interior ( figure 8(E)).
Another argument supporting the view of a viscous fingering instability is that simulation results of dividing isolated cell clones (not embedded in an external medium of B-objects) with the single-cell-based model used in this paper could be mimicked within a continuum description using mass conservation with a pressure-controlled growth term, ∂ t n + ∇ · (nv) = n H ( p 0 − p)/τ in combination with Darcy's law [66], v = −(k/η)∇ p [66]. Here, n(r,t) is the local cell density, H the Heaviside step function, v the velocity of the expanding front at the border of the A-cell clone, τ the effective cell cycle time that depends also on the shape of the cycle time distribution [53], p the local pressure (probably better referred to as normal stress), k the permeability and η the viscosity. As the simulation results with the single-cell-based model in [66] suggested a step function-like cell density profile dropping sharply at the population border, accompanied by a less sharp decrease of the pressure p towards the population border, the state equation n( p) looked to be a step-like function of n with a jump from zero density to the bulk density as p > 0 so n( p) did not need to be considered explicitly. The adhesion between cells of the expanding clone (of cell type A) generates a cohesive force term, which can be related to a sort of surface tension σ . The friction between B-objects and the substrate is larger than between A-cells and the substrate. So the growing cell clone may be interpreted similarly to a radial Hele-Shaw cell in which a (less viscous) fluid is injected into a more viscous fluid, the two fluids being constrained by two parallel plates. This is known to generate a viscous fingering instability (see, e.g. [34]). Further support for this view emerges from an analysis of the branching pattern, which we exemplify for the simulations shown in figures 6(E)-(H) (see figure 7). In figure 7(A), we draw the number of branches with increasing distance to the center of mass, given by R cm = N i=1 r i . At too large distances, approximately R > 50, the number of branches decreases again as then the population of A-cells come close to the border. Up to this value, the mean branch width is approximately constant even for a single realization of the growth process ( figure 7(B)).
Either by dividing the circumference at a distance R, given by U = 2 R by the number of branches using the data in figure 7(A), or using the branch size from figure 7(B), one obtains a measure proportional to the domain size λ crit at which an instability occurs. From a classical Saffman-Taylor instability, one would expect λ crit ∝ 1/[v(η BS − η 0 )] 1/2 where η 0 is the viscosity of a fluid of small viscosity expanding in another fluid with higher viscosity η BS . v is the growth velocity of the expanding fluid usually noted as dR/dt. In order to test this relationship for growing cell populations, we replaced the difference of viscosities by the difference of the friction coefficients of cell-substrate friction for A-cells and B-objects N /dt as a measure of the velocity of a growing front which can be calculated from the data on N(t) presented in figure 6. However, we found that the current data are insufficient to confirm the classical relation for λ crit even though the scenarios in figures 6(A)-(H) and the findings in figure 7 show the same qualitative tendencies as for a Saffman-Taylor instability. If one replaces the observed expansion velocity v by the intrinsic velocity v 0 of the expanding monolayer in the absence of the embedding B-objects, then we find approximately λ crit ∝ 1/(γ BS − γ 0 ) 1/2 . However, further simulations with larger systems and more realizations are necessary for identifying the precise relationships. These are very time consuming as the present simulations already took several weeks' simulation time.
After the instability has occurred, growth occurs almost only in the tip of the dendrites, even though a change in the width of the proliferating rim with increasing γ BS could not be  figure 8(I)). Due to the limitation of growth to the tip of the fingers, the total number of proliferating cells is reduced. Accordingly, the growth velocity of the expanding clone decreases with increasing friction of the embedding B-type objects with the medium, γ BS (figure 11), as this leads to more branches (figures 6(E)-(H)). On the other hand, dendritic populations show a significant (+50%) increase of the radius of gyration at saturation compared to compact cell aggregates of A-type, reflecting that a dendritic structure results in a larger radius of gyration compared to a compact circular cell distribution at equal numbers of cells (compare also the growth curves for R gyr and N in figure 11).

Variation of the host tissue cell density.
In a further step, we varied the initial density of B-type objects by varying the number of B-type objects seeded around the A-type start cell. This is a relevant situation, as both in vitro and in vivo the average cell density varies from cell line to cell line and also depends on environmental factors, for example the local nutrient concentration or cytotoxic substances [67][68][69].
We studied initial cell densities (ρ B ) between 3500 and 7000 cells per mm 2 , which encompasses the cell densities in many cell lines. The reference cell density (ρ 0 ) (see figure 5) was 5250 cells per mm 2 .
An increase of the initial cell density in the host tissue (ρ B = 7000 cells per mm 2 ) led to a significantly reduced expansion velocity of the growing monolayer, an earlier saturation of growth at a smaller diameter and a proliferating rim of smaller width during growth (figures 8(B) and 11).
Additionally, we studied the impact of the host tissue in case ρ B and γ BS were simultaneously varied. We found the influence of the discussed variations of ρ B to be consistent, regarding the expansion velocity (figure 11), saturation size and width of the proliferating rim (figures 8(C) and (D)), with those reported above. The growth of a cell population embedded in host tissue of both increased medium friction (γ BS = 20γ 0 ) and increased initial cell density (ρ B = 7000 cells per mm 2 ) is highly inhibited ( figure 8C).
On the other hand, comparison of figure 8(A) with (B) and of figure 8(D) with (C) indicates that only a moderate effect on the viscous fingering instability is caused by variations of the initial density ( figure 8(D)) even though the critical wavelength λ crit may be slightly smaller at higher densities of B-objects.

Variation of host tissue elasticity.
Another important property of the tumor host tissue is its elasticity. For example, it has been found that tumors and their associated stroma have a lower elasticity (higher Young's modulus) than normal tissue [70] which enables them to more effectively compress and eventually collapse blood and lymphatic vessels in their vicinity [71] even though cancer cells themselves are usually softer than cells of the tissue they originate from [72].
In general, cell elasticity is largely determined by the biophysical properties of the cytoskeleton, a complex and dynamic cross-linked protein network anchored to the inner surface of the fluid lipid bilayer. Experiments suggest that cellular elasticity can be influenced by substances that depolymerize the cytoskeleton [73]. Furthermore, cellular elasticity is as closely coupled to dynamic cell behavior as cell movement, cell shape changes or cell division. However, the highly nonlinear elastic properties of a living cell are difficult to study and a detailed, experimentally validated model remains elusive. Nevertheless, the elasticity of the cytoskeleton can be measured by an increasing number of experimental techniques such as bulk rheology, traction force microscopy or atomic force microscopy [74][75][76].
In the simulations above, Young's modulus of the embedding objects E B was chosen to be the same as Young's modulus E A of the growing monolayer (type A) for which we assumed a Young's modulus of elasticity E A = E 0 = 450 Pa [77,78]. A decrease of Young's modulus from E B = 450 Pa to E B = 300 Pa leads to an increase of the proliferating rim and to a smaller critical wavelength (figure 8(F)), while an increase of Young's modulus of the embedding B-type objects (to E B = 1000 Pa in figure 8(G)) leads to a decrease of the proliferating rim and a larger critical wavelength. Furthermore, a decrease of Young's modulus also leads to increased growth velocity of the tumor accompanied by saturation at an increased population size (figure 11), reflecting that the resistance of the embedding objects is reduced.
If, in addition to the variation of E B , the medium friction of the host tissue is increased (γ BS = 20γ 0 ), we observe an amplification of the described influence (figures 8(H) and (I)). In case elasticity was decreased (E B = 1000 Pa) and the medium friction in the tumor microenvironment was simultaneously increased (γ BB = 20γ 0 ), the growth of the tumor monolayer was found to be highly inhibited ( figure 8(H)). Apparently, in this case the critical pressure threshold above which cells are not able to enter the cell cycle is reached quickly.

Variation of host tissue adhesivity.
We further studied the impact of differential cell-cell adhesion [79,80]. We consider changes of the adhesion strength of cells of the same type as well as between cell type A-and B-type objects. In general, intercellular adhesion is mediated by transmembrane proteins (selectins, cadherins and the immunoglobulin (Ig) superfamily). Except for the last, all of them require Ca 2+ and Mg 2+ . Therefore, many adhesive interactions are Ca 2+ or Mg 2+ dependent [69]. Moreover, in tumor spheroids, a significantly increased adhesion was found for CoCl 2 -induced hypoxia, while ionizing radiation had an inverse effect [81]; hence adhesion is subject to changes.
In recent decades, many techniques have been developed to measure cellular adhesion; for example, Benoit et al [82] used a modified version of atomic force microscopy to measure adhesion forces between living cells and surfaces at the molecular level.
In our model the strength of the adhesive forces was determined by the density ζ m of the corresponding receptors on the cell surface. As introduced before, we have usually chosen ζ AA = ζ BB = 10 15 m −2 [83,84] and ζ AB = 0 m −2 . Adhesion among B-type objects leads to the formation of patches with empty spaces in between if the density of B-type objects ρ BB is not space-filling (figure 9(I), inset). The comparison of simulations without and with adhesion among B-type objects suggests a greater tendency towards finger formation if B-type objects adhere among each other (figures 9(A)-(C)) as the formation of fingers is facilitated in the free spaces between the patches. If the adhesion between A-cells and B-objects is stronger than A-A cell or B-B object adhesion, then inclusions of B-objects are observed in the expanding A-cell population (figure 9(F)). The adhesion between A-cells and B-objects, moreover, reduces the stress at the A-B interface compared to purely repulsive (Hertz) forces as in figure 9(D).

Variation of the average cell size of the host tissue.
In the preceding sections the size (diameter) l A of the growing tumor cells (type A) and the size l B of the elements of the embedding tissue (B-type objects) were assumed to be identical (l B = l A ). However, cell size is known to vary significantly between different cell types and tissues. Most eukaryotic cells have an average size (diameter) of 10-25 µm. However, the average cell size is not constant. It can be influenced by various factors, for example, cytotoxic substances [69] or genes involved in growth control [85]. The cell size of a specific cell line or tissue can, in general, be measured with bright field and laser scanning microscopy. In experiments using granular objects, the size of the objects can easily be varied. For these reasons, in this section we study the influence of the size of embedding B-type objects on the expanding A-type clone by varying the size of the B-type objects covering the range of typical eukaryotic cell sizes from 7 to 30 µm.
We find a decreasing tendency to form an instability with increasing object size (figures 10(A)-(D)). Perhaps, small objects can be more easily re-arranged than large objects.  Furthermore, we found that larger embedding objects resulted in a slightly decreased expansion velocity, whereas smaller embedding objects led to a slightly increased growth velocity (not shown). In both situations, the width of the proliferating rim and the saturation size of the growing population remained unchanged. This can be explained by the larger finger formation if B-type objects are larger: as the A-type population grows mainly at the tip of the fingers, more (and smaller) fingers as found with increasing the size of B-type objects lead to an overall smaller number of proliferating cells, thereby reducing the expansion speed of the A-type cell clone. However, the effect of a stronger finger formation at smaller B-object size is surprising, because for infinitely small objects one might not expect a tendency towards stronger fingering. Further computational studies have been started to study this limit.

Apoptosis with cytolysis of apoptotic cells.
Cytolysis, which is also designated as osmotic lysis, is a degenerative cellular process that involves the destruction of the outer cell membrane. This causes excess water to move into the cell and ultimately leads to the dissolution of the affected cell. Neighbor cells are able to phagocytize a dead cell. A majority of current assays for measuring cytolysis are based on the detection of changes in plasma membrane permeability and either the subsequent leakage of components (e.g. cytoplasmic enzymes) or uptake of dyes that are normally not able to enter the cell (e.g. Trypan blue or propidium iodide).
In our model, we introduced cytolytic processes to validate the robustness of our findings from preceding sections in a more dynamic and regenerative environment.
In order to model cytolysis, we assumed that after a certain time t cl , apoptotic tumor cells are completely dissolved. In comparison with the preceding sections, where the main fraction of cell proliferation took place at the borders of the population, the introduction of cytolysis led to much more complex and highly dynamic proliferation activity patterns (figures 12(A) and (B)). Snapshots shown in figures 12(A) and (B) reveal proliferating regions (colored white) even in the interior of the tumor cell cluster. In this situation of constant cell renewal, no tumor cells remained arrested in G 0 for a long time (as for example those colored red or black in figure 12). The reason for this was that densely packed cells in G 0 became apoptotic and then were removed by cytolysis after t cl . Thereafter surviving tumor cells in the vicinity of recently dissolved cells (black in figures 12(A) and (B), see arrows) were able to migrate into the free space and re-enter the cell cycle because of temporarily lowered cell density and thus lowered pressure (figures 12(E) and (F)). Cytolysis largely destroys (figures 12(C) and (D)) the typical layered proliferation pattern (e.g. figure 5) of growing cell populations.
Furthermore, the constant cell renewal introduced oscillations into the growth kinetics (figure 11, lines L). Importantly, the radius of gyration remains largely the same as for the corresponding simulations without cytolysis ( figure 12, upper part). Moreover, we found no further changes of any of the previously described biomechanical influences of the tumor host tissue on the morphologies and growth dynamics of tumor cell populations. We therefore consider the results presented in previous sections to be robust and especially not depending on cell proliferation to be confined to the border of a multi-cellular population.

3.2.7.
Growth and apoptosis in the expanding clone and host tissue. Finally, we considered the case of embedding cells capable of proliferation and apoptosis as for the expanding clone testing the idea of a homeostatic pressure proposed in [15]. We model this process in several steps. Firstly, the A-type clone was grown in the embedding medium of B-type cells suppressing cell proliferation and death until the A-type clone had a size of about N = 10 000 cells, which was the case at a gyration radius of about 950 µm (figure 13). Then we introduced a relaxation phase (blue curve in figure 13). At about t = 15 days we restarted the simulation with now both the A-type cells and the embedding B-type cells proliferating and undergoing apoptosis. Both, cell cycle entrance and cell apoptosis, were triggered by a pressure threshold. At the first pressure threshold, a cell becomes quiescent, and at the second pressure threshold that was larger than the first threshold, it becomes apoptotic. We consider two cases. Firstly, we assume that the thresholds are pairwise the same for A-and B-type cells ( p

Growing multi-cellular spheroids in an embedding medium
In this section, we consider the growth of multi-cellular spheroids in a B-type object environment. We assume that cells secrete ECM in which they are able to anchor and move. This line of argument is supported by the observation in SK-MES-1 multi-cellular spheroids which show the presence of collagen V in the region of the tumor spheroid. The equations of motion for this case are adapted to three dimensions. Our studies in this section are inspired by the findings of Helmlinger et al [1].
We again observe a fingering instability if the friction coefficient between B-type objects and matrix is much larger than for A-type cells and matrix, γ BS γ AS = γ 0 ( figure 14). We calibrated the initial density and size of the flask to mimic the modified degree of agarose by different initial densities of B-type objects. We made a comparison with the experimental growth curves of Helmlinger et al [1] for LS174T cell populations and of Galle et al [60] for WiDr cells.
(We find that the initial density of B-type objects determines the saturation size, while the size of the flask influences the slope.) The density of B-type objects at saturation is in each case about 6 × 10 5 B-objects per mm. The idea of this study was not to mimic in detail the constitutive relation of agarose but rather to demonstrate that the pressure-mediated cell cycle control mechanism is able to reproduce the kinetics of growing cell populations in media of different mechanical properties. The inset in figure 14 shows that removal of B-objects, resulting in a release of external mechanical stress on the A-cell population, leads to a continuation of growth of the A-cell population as observed experimentally in [1]. Our findings about the tumor growth kinetics are consistent with those found in a multi-phase continuum model of Chen et al [20] already introduced above. However, assuming spherical symmetry the authors studied the time dependence of the radial coordinate. Moreover, also the experimentally observed effect of the shape of the flask on the spatial shape of the growing A-cell population can be reproduced ( figure 15). Growth occurs in the direction which minimizes pressure.

Discussion
As pointed out in [61,62], isolated monolayers show linear growth kinetics after an initial exponential growth phase. In previous work [24], we have shown that the linear growth kinetics of dense monolayers can be explained by cell proliferation as the main trigger for the movement of cells towards the monolayer edge and thereby for expansion of the monolayer. In that model the entrance of a cell into the cell cycle could occur only if that cell was not subject to too strong deformation. Proliferation-triggered movement leads to a pushing-type movement of cells. If the growth and apoptosis pressure thresholds are equal in the expanding A-type cell clone and the B-type embedding cells, the average size of both clones remains equal (black curve), while if the pressure at which B-cells undergo apoptosis is below the pressure at which A-cells undergo apoptosis, then the A-clone outcompetes the embedding B clone (green curve). The pressurecontrolled cell cycle entrance and apoptosis mechanism thereby provides a possible representation of the concept of homeostatic pressure introduced in [15].
However, recent findings in [46] suggest that cell micro-motility may not be uniformly distributed but instead is directed towards the outwards-pointing normal of the monolayer border leading to traction of the cells. Cell proliferation was observed but the trigger of cell proliferation has not been studied. In order to probe key aspects of these observations, we studied in this paper an outwards-oriented active component of micro-motility modeled as a constant drift force on border cells in combination with three different growth control mechanisms: (i) a cell can only enter the cell cycle if the pressure on it is below a critical positive threshold.
(ii) A cell can only enter the cell cycle if the total pressure is negative. (iii) A cell can enter the cell cycle only if it is stretched into the direction of at least one of its neighbor cells. We mimicked different mechanisms and found that the growth characteristics are very robust with regard to the qualitative dynamics, but the growth speed varies depending on the strength of the directed migration force component. Generally, a directed component in the micro-motility speeds up the monolayer expansion as the movement generates space for the dividing cells. In the case of control mechanism ((i); rule 2 in figure 1), even local patches of cells under tension (negative pressure) could be observed if the directed movement was sufficiently large; in the absence of such a force the cells are all under compression, so the pressure is positive. For the control mechanism ((ii); rule 3 in figure 1) the pressure is further reduced as only negative pressure is able to trigger cell cycle entrance. For the control mechanism ((iii); rule 4 in figure 1) the movement is a pulling movement: only if cells are stretched, are cells able to enter the cell cycle, leading to a further reduction of the local pressure. In both (ii) and (iii), the pressure at the border was negative. However, if entrance into the cell cycle occurs only if cells experience tension (case (ii)), or if they are stretched  [1] for LS174T spheroids: 1.0% agarose gel (filled black circles), 0.8% agarose gel (filled red triangles), 0.5% agarose gel (filled blue diamonds) and 0.3% agarose gel (filled green squares). Model simulation referring to LS174T spheroids using ρ 0 = 611 × 10 3 cells mm −3 and a flask diameter of 500 µm (green line) and 207 × 10 3 cells mm −3 (red line) and a flask diameter of 500 µm. Experimental data of Galle et al [60], WiDr spheroids: 5.0% agarose gel (red squares) and 2.5% agarose gel (blue circles). Blue line: model simulation referring to WiDr spheroids, ρ 0 = 225 × 10 3 cells mm −3 with a flask diameter of 250 µm. Inset: experimental data on growth alleviation by Helmlinger et al [1] for 0.7% agarose gel (symbols) and model simulation (line) using ρ 0 = 215 × 10 3 cells mm −3 , a flask diameter of 500 µm and removal of embedding cells at t = 30 days. For all simulations in (C), γ BS = γ AS .
(case (iii)), this slows down the monolayer growth speed compared to a re-entrance as long as cells are not too much compressed (case (i)).
In the second part of the paper, we considered a growing population of cells embedded in a population of either granular objects or cells in order to study the emerging growth pattern. In this case, a directed movement had almost no consequence as long as the directed force component did not lead to rupture of the cell-cell contacts of the expanding cell clone.
We calculated until the growing clone enclosed by other objects had a population size of about 30 000 cells and varied the parameter values of the model over physiologically reasonable ranges. At this population size we could observe a viscous fingering instability which was particularly pronounced if the friction of the B-objects with the substrate was much larger than the friction between A-cells and the substrate. With increasing friction, the critical wavelength decreased. Moreover, the fingering pattern was moderately increased with increasing object size of the external medium as well as for increasing adhesion with between B-type objects. The micro-motility of the external objects has no detectable influence on the instability, and the density of B-type objects seem to only influence the growth speed and the saturation size 6 A higher resolution version can be found at http://msysbio.com/videos njp. of the population but not the fingering pattern. Moderate changes of Young's modulus from 300-1000 Pa mainly reduce the expansion speed of the growing clone. Even apoptosis, triggered if the local pressure overcomes a critical value, almost does not change the growth pattern. Instead, dying cells are quickly replaced by dividing neighbor cells.
We compared the relation for the critical wavelength of the instability with that of the Saffman-Taylor instability. Despite good qualitative agreement, larger systems will have to be studied to determine the quantitative dependence of the critical wavelength on the model parameters. These simulations are currently being performed, but they are very time consuming with one simulation running over several weeks to several months.
Recently, Basan and co-workers [43] considered the undulations at the interface between the epithelial tissue layer and the underlying stroma often observed in skin or mucosa. Epithelial cells flow towards the free surface from where they desquamate. Proliferation may trigger an instability as a consequence of shear between adjacent regions that becomes important as differences in the local epithelium height cause differences in the epithelial flow in adjacent regions. Higher proliferation, a higher local thickness of the proliferating layer and a higher epithelium viscosity each increase shear between adjacent regions, while, for example, a higher stroma viscosity stabilizes the interface. The latter is in contrast to our findings where higher viscosity in the environmental material facilitates the formation of fingers. Lowengrub et al [39] discuss an elastically induced splitting instability in which an elastic particle splits under deviatoric applied stress. They find that the onset and direction of the splitting instability depends on the magnitude of both the misfit strain and applied stress as well as the details of the applied stress, the elastic constants and the initial shape of the particle. The results in figure 10 show that the size of the embedding objects affects the instability. Figure 9 shows that the strength of adhesion also matters. Hence it seems likely that the instability observed goes beyond the classical Saffman-Taylor instability despite the most prominent effect of the friction of A-cells and B-objects with the substrate.
In a next step, we considered birth and death in both the inner expanding clone and the embedding material, hence mimicking a competition process involving birth and death in the inner and the enclosing cell population. If we assume that the pressure threshold for birth in both populations is equal and the pressure threshold for death is equal, then the average population size of both populations remains the same, mimicking tissue homeostasis. If, however, the pressure threshold at which cells of the inner clone are still able to enter the cell cycle is equal to or larger than the pressure at which the embedding cells undergo apoptosis, then the inner cell population outcompetes the outer cell population. This observation is compatible with the recently proposed concept of a homeostatic pressure [15].
Finally, we considered d = 3 dimensions. Firstly, we noted that the fingering instability occurring with increasing friction of B-objects with the substrate can again be observed. Then we tested the impact of the initial B-object density on the growth kinetics and found that the growth kinetics corresponds to that found in [1] and [60] for cell clones growing in agarose gel.
We conclude that mechanical stress or strain-controlled cell cycle entrance mechanisms are able to explain many features of the growth pattern in monolayers and multi-cellular phenotypes both in the absence of a viscous or visco-elastic embedding medium and in its presence. The model predicts a fingering instability in a granular embedding medium which becomes particularly pronounced when the friction between the embedding objects and the medium between the embedding objects (the substrate) is increased.