Mechanical influences on in silico tumor evolution

Experimental insight and conceptual understanding of tumor growth are steadily growing and leading to new therapeutic interventions. Experiments and clinical studies are able to link single-cell properties to macroscopic tumor attributes. The development of cellular subpopulations in heterogeneous tumors can be understood as an evolutionary system with different cell types competing over both space and nutrients. However, to predict the growth trajectory and development of a tumor, fitness and trade-offs of cell properties in the context of the surroundings are required and often inaccessible. The optimum of the evolutionary trajectory provides a target for intervention, but can mostly not be identified. We propose that the optimal value of cellular properties is influenced by the tumor surrounding. Computational multiscale-modeling of tissue enables the observation of the trajectory of each cell while modeling the tumor surrounding. We model a 3D spheroid tumor and the fitness of individual cells and the evolutionary behavior of the tumor are quantified and linked to global parameters. Cell–cell adhesion and cell motility are two important mechanical properties for cell development and used as free parameters. Mechanical properties alone are able to drive the tumor towards low adhesion.We implement a dynamically changing nutrient surrounding representing the fluctuating blood-supply through blood vessel collapse and angiogenesis. We find that the evolutionary speed depends on the frequency of the fluctuations. We identify a frequency domain in which the evolutionary speed is significantly increased over a tumor with constant nutrient supply. The findings suggest that mechanically-induced fluctuations can accelerate tumor evolution. Author summary Limited space and nutrients together with competing cell types drive an evolutionary process inside tumors. This process selects for the fittest cell types and optimizes the growing behavior for its local surroundings. An expanding tumor exerts mechanical forces on its cells and its surroundings, leading to a fluctuating nutrient supply through collapsing blood vessels. Here, we observe the influence of a dynamically changing surrounding on the evolutionary behavior of heterogeneous tumors in a high-resolution computational model. We find that the evolutionary speed depends on the frequency of the fluctuations and a fitness advantage of low-adhesion cells.

The emergence and development of tumors in humans still presents a significant 10 challenge to medicine and is one of the leading causes of death. Knowledge about the 11 development of tumors is continuously expanding. Especially, the large variety of tumor 12 types, differences between patients and tumors consisting of many types of tumor cells 13 (termed tumor heterogeneity) present challenges. Here, a tumor in surrounding solid 14 tissue and the development of tumor properties is observed over time in a 15 computational model. The development of a tumor can be described as an evolutionary 16 system in which cell types present species competing over resources [1,2]. The evolution 17 of a tumor is driven by mutations. These mutations change the properties (e.g., motility 18 or cell-cell adhesion) of cells, which then can provide advantages or disadvantages for 19 further growth. Competition over resources and space leads to the selection of the most 20 advantageous tumor cell types, whose subpopulation will overtake the rest. However, 21 cells can not optimize all parameters at the same time, several trade-offs in cellular 22 properties in cancer have been described and studied [3,4]. We assume a limited energy 23 budget for each cell and introduce a trade-off of proliferation against motility for cell 24 types. 25 While the properties and parameters of a cell are accessible in experiments, the 26 fitness of a cell is very hard to determine experimentally. The fitness determines the 27 ability of a cell to reproduce in its surroundings and multiply, therefore cells with high 28 fitness are selected for in evolution. To determine the fitness of cells experimentally, the 29 lineage of cells has to be tracked over time, making it necessary to track each individual 30 cell over a long time, which is currently not possible in in vivo tumors. Models of the 31 fitness landscape based on genetic changes and driver mutations have been 32 introduced [5], yet phenotypic changes in cells can also drive mutation [6,7]. Naturally, 33 an ensemble of cells with different cell types evolves towards the cell type of highest  The mechanical interactions of the tumor with its environment and between the 40 tumor cells influence the trajectory of the tumor [8]. Cells with low cell-cell adhesion 41 mechanically sort to the outside of tumors, providing a higher nutrient supply and 42 therefore evolutionary advantages. The nutrition levels on the surface of tumors are 43 higher than in the center of the tumor, leading to an induced advantage of low adhesion 44 cells. Therefore, we will test this prediction and observe whether those cells have an 45 evolutionary advantage in an evolving heterogeneous tumor. Cell motility describes the 46 active movement of cells, we predict that random cellular movements should provide no 47 significant evolutionary advantage in a constant surrounding. Phenotypical changes 48 March 6, 2021 2/16 towards higher motility cells and the epithelial to mesenchymal transition are associated 49 with tumor invasion [9][10][11]. In a dynamically changing nutrient environment, cells with 50 higher motility can dynamically occupy the most advantageous positions. Enhanced 51 motility is known to be a driver of the formation of metastases [12], therefore a dynamic 52 nutrient field could induce metastases.

53
In this work, we focus on two variables in cellular properties, namely cell-cell 54 adhesion and motility. Increased motility is a landmark for the development of 55 metastases [12]. We hypothesize that mechanical and geometric constraints alone are 56 sufficient to drive the evolution of a tumor towards high motility.

57
The computational modeling of tumor development and heterogeneity has evolved 58 and was successfully applied, e.g., a cellular automaton model was used to observe the 59 influence of cell dispersal and turnover on tumor heterogeneity [13]. The use of 60 computational models in clinical applications and personalized treatment plans is 61 gaining momentum [14,15]. Work on tumor internal evolution has been recently done 62 by Büscher et al. [16] where they compared adhesion and 'growth strength' in an 63 evolving tumor and found that in some cases a mixture of different cell types emerges as 64 a stable state. The importance of cell-cell adhesion for tumor invasion has been recently 65 highlighted by Ilina et al. [17]. In Ref. [18] the authors showed with a mathematical 66 model using the evolutionary game theory that tumor heterogeneity and the optimal 67 cell phenotype depends on the microenvironment and position within the tumor. The

68
'go-or-grow' hypothesis describes reversible phenotype changes between motile and 69 proliferative phenotypes in cancer. This hypothesis was proposed and tested with 70 positive [19] and negative correlations [20] for different cancer cell lines. Trade-offs 71 between two or more cell properties have been studied by Gallaher et al. [4]. Here, we 72 are not focusing explicitly on the 'go-or-grow' hypothesis, since we assume each cell to 73 retain its phenotype throughout its lifetime and changes of cell properties can only 74 occur at mutation events during cell division.

75
Cell adhesion and cellular motility have been recently shown to strongly influence 76 the invasion behavior of breast cancer and to differentiate between solid-like, fluid-like, 77 and gas-like behavior [17]. Confinement by the extracellular matrix and cellular motility 78 were recently investigated systematically and a phase space of tumor invasion modes 79 was proposed [21]. blood vessels [22]. This interplay of formation and collapse of blood vessels can lead to a 87 fluctuating availability of nutrients for the cells in the tumor.

88
During the development of tumors, single cell effects are of major importance, with 89 mutations initially occuring in a single cell. Stochasitcity and rare events are 90 non-negligible during tissue development. As we showed in [23], single cell effects are 91 important for the patterning of tissue. Cell migration is facilitated by many concurrent 92 processes and represents a multiscale process, therefore requiring multiscale 93 modeling [24]. To capture the complexity of tumor growth, computational models with 94 single cell resolution enable the incorporation of single cell behavior that would be 95 difficult to capture using continuous ODE models. Since evolution is a stochastic process 96 requiring many iterations to find a stochastically meaningful result and the requirement 97 of large-scale tissues makes the use of high-performance computing necessary.

98
The recently developed framework CellsInSilico [25] is used for the simulations, it is 99 based on the cellular Potts model (cpm) [26] that has been established for the simulation of tumor growth [27]. Through its parallelization and optimization for 101 supercomputers, the framework enables large-scale three-dimensional simulations and 102 high numbers of simulations, which is necessary to generate sufficient statistics to study 103 the proposed problems. 104 We implement a discrete evolution of two independent parameters, namely cell-cell 105 adhesion and cell motility. Cells can alter and change their parameters stepwise, with a 106 small probability at each division. Parameters are changed incrementally, which implies 107 a continuous evolution and neglects possible mutations with larger effects. For cell 108 motility, we introduce a trade-off on division rates since we assume cells have a limited 109 energy contingent. We simulate a spheroid tumor in an external gradient of nutrients. 110 The tumor is surrounded by a population of non-mutating host tissue. We observe the 111 evolution of the tumor composition along both free parameters. Different mechanisms 112 that couple cell division rates to nutrient availability are compared. We compare the 113 effect of temporally changing nutrient availability.

Constant environment 116
First, we investigate the dynamics of the system without the dependency of cell division 117 and death on nutrient availability. Division and death are determined by the respective 118 rates and age and volume constraints. A spheroid tumor develops from an initial tumor 119 seed, consisting of ≈ 3700 tumor cells (cf. Fig. 1) in a surrounding of nondividing and available space. In the emerging steady state, the tumor size remains constant with cells 125 dying and dividing at equal rates. The cell division rate is higher than the death rate 126 (cf. Fig. 2 a-c) and the absolute number of cells is geometrically constrained, thereby the 127 effective division rate is limited by the death rate.

128
Observing the statistical occurrence of cell events in relation to the radial distance to 129 the tumor center, we find that cell deaths are located in the tumor center while cell 130 divisions are located at the tumor margin (see Fig. 2 a). This leads to an inward 131 movement of cells (as seen by the stripes in Fig. 1). There is a buildup of pressure 132 inside the tumor, this is facilitated by the inwards movement of the cells and the central 133 potential. This behavior is well studied in experiments and computational models [28]. 134 The simulation is started with a single cell type in the center of the phenotype space. 135 Through mutations during tumor growth, more cell types are introduced which leads to 136 a distribution in the phenotype space around the initial cell type. With the progressing 137 simulation, the distribution moves in the phase space towards an evolutionary favorable 138 position. The center of mass, spread, and movement speed of this distribution are   Macroscopically, the introduction of a nutrient dependency affects the evolution 160 speed of the tumor (see Fig. 2 d). TBD decreases the evolutionary speed of the system 161 while LRD accelerates the evolution. Enhanced directionality is visible in the different 162 sizes of the spread in the phenotype space. A smaller spread is achieved through a more 163 pronounced directionality of evolution. TBD increases the spread while LRD lowers it. 164 This increase of the evolutionary speed of LRD over TBD can be explained by the larger 165 regions, in which cell divisions and cell deaths occur. In these critical regions (tagged 166 red in Fig 2 a-c, right), the competition between different cell types is most pronounced 167 since cells that can stay in this area or escape outside will survive, while cells that are 168 pushed to the inside will die. Despite the different evolutionary speeds and spreads, the 169 evolution is highly directional towards low adhesion for both mechanisms and no The mechanism of nutrient dependency with LRD shows a larger selective pressure 178 on the tumor and therefore leads to a higher speed of evolution. Furthermore, a linear 179 dependency is biologically more reasonable, since cells do not binarily up-or 180 down-regulate cell division in most cases, but adapt continuously [29]. Hence, LRD is 181 used in the subsequent manuscript if not explicitly stated otherwise. Values are averaged between 1500 − 3000 kMCS. In the first group (left), the simulations are performed with constant nutrient availability and a varying radius of the dip in nutrient availability (r dip ). In the following groups, the nutrient availability is dynamic, and the dip moves on a circle or a line in the x,y-plane with radius A and period T . In the dynamic cases is r dip = 90. The green dashed line indicates the value of the reference simulation with constant surroundings. b The average center of mass in the conformation space on the 'motility-1/division-rate' axis. Cases are distinguished and data is collected as in a).

Gradient steepness 183
Next, we observe how the gradient steepness which is linear to the gradient extent 184 influences evolution. Evolutionary speed increases with decreasing steepness of the 185 nutrient gradient(see Fig. 3 a, left). The speed is increased because the extent of the  Fig. 6 c). The tumor center is never 229 affected by the nutrient drop, and therefore always provides optimal conditions. This  Fig. 4 b). This is reasonable since the 239 nutrient availability is directly coupled to the position. Looking at the dynamic case, 240 the nutrient availability first decreases, but then increases. This is due to the fact that 241 cells divide in a region in which other cells die, cell death is introduced by a shortage of 242 nutrients. After division, cells move inward towards the tumor center and the nutrient 243 drop continues to move. Here, preferably cells survive that quickly exit the nutrient 244 drop, or divide at the 'rising edge' of the moving gradient. We hypothesise that this is 245 the main driver of the change in evolutionary optimum on the motility scale since 246 fast-moving cells have a statistical advantage in this case.

247
We define the lifetime of a cell as the time between the last division and cell death. 248 The mean lifetime and the extent of the lifetime distribution, significantly increases 249 with the introduction of a dynamic surrounding that enables a 'save spot' (cf. SI 250 Fig.12 a). Here, the nutrient dip moves around the tumor in a circle that is larger than 251 the tumor. The nutrient availability in the tumor center is therefore always optimal, 252 providing zero evolutionary pressure and therefore a 'natural reserve' on the population 253 scale. We observe the development of the mean lifetime of the tumor cells over time, see 254 SI. Fig. 6 e. Overall, an increase in the mean lifetime is visible for both dynamic and We present a computational model of a spheroid tumor in surrounding tissue. Mutation 277 of cells is enabled by a change of phenotype at cell divisions. Two parameters can be 278 changed during a mutation, cell-cell adhesion and the motility of the cell. We introduce 279 a division rate trade-off for motility. The system is allowed to evolve freely and the 280 tumor composition is tracked in parameter space over time. 281 We find that the mechanical and geometrical properties of the system are sufficient 282 to drive the ensemble towards low-adhesion cell types. This mechanical effect in tissue 283 evolution has been described in [16]. Mechanical properties alone drive proliferation at 284 the tumor edge and cell death in the center. 285 We introduce a dependency of cell divisions and deaths on nutrient availability, 286 which is linearly decreasing towards the tumor center. Using a linear dependence on 287 nutrient availability for proliferation and inverse linear dependence for cell death leads 288 to a higher evolutionary speed than a threshold-based dependency.

289
In in vivo tumors, solid stresses through tissue displacement are built up that are 290 able to compress and block blood vessels [22]. This can lead to fluctuating and nutrient 291 availability in tumors. We investigate how fluctuating nutrient availability influences 292 tumor evolution. 293 We find that a temporally variable nutrient surrounding introduces a larger life span 294 for cells. Especially a 'save spot' enables much longer lifetimes and a broader 295 evolutionary spread 296 We find a significant dependency of dynamic nutrient surroundings on the 297 evolutionary speed in phenotype space. The speed shows a frequency dependency, with 298 a lower evolutionary speed for fast fluctuations followed by an increase over the constant 299 case for lower frequencies. A critical time scale exists for the fluctuation of nutrient 300 availability that provides a distinct peak in evolution speed, which we find to be 301 between T = 100 kMCS and T = 200 kMCS.

302
The fitness of cells can be determined by lineage tracking and the fitness is linked to 303 the cell types and their parameters. With this, the effect of parameter values of a cell 304 on its fitness can be determined. While a clear preference in fitness is visible for low 305 adhesion cells, no clear change in the preferential direction of evolution can be identified 306 along the motility axis. A trend towards higher motility is visible for large radii of 307 nutrient fluctuations. We predicted that dynamic nutrient availability influences the 308 fitness optimum of tumor evolution, this could not be conclusively be answered and has 309 to be explored further by extending the range of possible motility and introducing 310 different trade-offs.

311
Experimental work on spheroid tumors could provide a verification of the results 312 found here. The single-cell motility of cells grown in different spheroid cultures could be 313 measured and compared [30]. The nutrient surrounding of the growing spheroid culture 314 can be varied from a constant availability to a periodically changing nutrient 315 concentration in the surrounding media. We expect to find cells with higher motility in 316 the latter setup.  The cell-cell adhesion is proportional to the contact area between cells and 326 independent of the duration of the adhesion. The strength is not limited or quantized 327 by focal adhesion but only determined by the adhesion parameter between the cell types 328 and the shared area. To divide, cells need to exceed the age of 2 kMCS and their volumes have to exceed 351 90% of their goal volume. Similarly, cells can divide once they exceed the age of 4 kMCS. 352 There are three different cases for dependency of cell division and death on nutrient 353 availability: Division rates vary from cell-type to cell-type since they are determined by the division 366 rate -motility trade-off.

367
The tissue that surrounds the tumor does not participate in cell death or division.

368
The cells that make up the surroundings, therefore participate in the entire simulation 369 and act as a medium that the tumor cells at the tumor edge can interact with and 370 redistribute forces and pressure.

371
Evolution speed and spread: The speed of evolution is calculated by tracking the 372 center of mass of the distribution of cell types in the phenotype space. The spread is 373 measured by the extent of the distribution.