Amoeboid-mesenchymal migration plasticity promotes invasion only in complex heterogeneous microenvironments

During tissue invasion individual tumor cells exhibit two interconvertible migration modes, namely mesenchymal and amoeboid migration. The cellular microenvironment triggers the switch between both modes, thereby allowing adaptation to dynamic conditions. It is, however, unclear if this amoeboid-mesenchymal migration plasticity contributes to a more effective tumor invasion. We address this question with a mathematical model, where the amoeboid-mesenchymal migration plasticity is regulated in response to local extracellular matrix resistance. Our numerical analysis reveals that extracellular matrix structure and presence of a chemotactic gradient are key determinants of the model behavior. Only in complex microenvironments, if the extracellular matrix is highly heterogeneous and a chemotactic gradient directs migration, the amoeboid-mesenchymal migration plasticity allows a more widespread invasion compared to the non-switching amoeboid and mesenchymal modes. Importantly, these specific conditions are characteristic for in vivo tumor invasion. Thus, our study suggests that in vitro systems aiming at unraveling the underlying molecular mechanisms of tumor invasion should take into account the complexity of the microenvironment by considering the combined effects of structural heterogeneities and chemical gradients on cell migration.

individual cell motility has been studied using both experimental 7,10,12 and theoretical [13][14][15][16][17][18] approaches. However, it remains unclear how the adaptation responses of amoeboid and mesenchymal migration modes contribute to the tumor invasion process. In particular, it is not known if and how amoeboid-mesenchymal plasticity allows a more effective invasion compared to the non-adaptive amoeboid or mesenchymal modes. So far, only the impact of interactions between non-switching moving cells and the ECM on tumor invasion has been studied 4,6,19 . Hecht et al. 6 considered an agent-based model for fixed proportions of amoeboid and mesenchymal cells migrating in a maze-like environment with directional cue, and showed that an amoeboid cell population can benefit from the presence of a small number of path creating mesenchymal cells.
In this work, we develop a mathematical model to study the consequences of amoeboid-mesenchymal migration plasticity on tumor invasion in a switching population of cells that adopt their migration mode in response to the ECM conditions. We are interested in the question under what environmental conditions the plasticity of amoeboid and mesenchymal migration modes provides an advantage for tumor invasion. We formulate and analyze a cellular automaton model, where cells are able to switch between a slow mesenchymal migration mode with ECM degradation, and a fast amoeboid migration mode without ECM degradation, depending on the biomechanical resistance imposed by the ECM. With computer simulations of the mathematical model, we compare the invasion behavior of a population of cells where each cell is allowed to switch between ameoboid and mesenchymal migration modes with the behavior of a non-switching cell population, under different environmental conditions. In particular, we distinguished spatially homogeneous and heterogeneous ECM resistance conditions, ranging from weakly structured to highly structured, in combination with different responses towards a chemotactic gradient which directs cell migration. Our numerical analysis reveals that ECM structure and presence of a chemotactic gradient are key determinants of the model behavior. The amoeboid-mesenchymal migration plasticity leads to a more widespread invasion compared to the non-switching situation only in complex microenvironments, more specifically, if the ECM is strongly heterogeneous and a chemotactic gradient directs migration. The latter conditions are characteristic for in vivo tumor invasion. This suggests that experimental studies on tumor invasion should represent this complexity of the microenvironment.

Methods
The model. We develop a mathematical model to study the effects of amoeboid-mesenchymal migration plasticity on tumor invasion. To determine the specific impact of migration plasticity of individual cells on overall cell population invasion dynamics, we coarse-grain to a cell-based model, namely a probabilistic cellular automaton (CA), which is analyzed at the population level. Probabilistic cellular automata are a class of spatially and temporally discrete mathematical models which allow to (i) model cell-cell and cell-ECM interactions, as well as cell migration, and (ii) to analyze emergent behavior at the cell population level [20][21][22][23][24][25][26] .
We consider the ECM as a physical barrier which imposes a resistance against the moving cell body. A widely studied parameter which mechanically impedes cell movement is the ECM network density. Other physical properties of the ECM, such as porosity, as well as biomechanical properties like ECM stiffness, have been observed to either enable or restrict cell migration. Importantly, the different ECM properties are not independent but rather connected 2 . Thus, for instance, the density of fibrillar ECM is positively interconnected with stiffness and inversely proportional to pore size, such that alterations of either property impact the overall ECM structure 12 . In view of this, we incorporate ECM properties like density and stiffness into a lumped parameter called "resistance".
We distinguish two migration modes, namely amoeboid-like (A) and mesenchymal-like (M) migration. We assume that amoeboid-like cells migrate in a protease-and integrin-independent way, driven by Rho/ ROCK-mediated contractility. In contrast, mesenchymal-like cells migrate protease-and integrin-dependently, where Rho/ROCK signaling is inhibited. Thus, the amoeboid and mesenchymal migration phenotypes in our model represent two ends in a broad spectrum of individual cell migration modes 11 .
We assume that cells change their migration phenotype in dependence of local ECM resistance. This assumption is supported by several studies. In particular, it has been reported that under soft ECM conditions, cells are round-shaped and their migration is independent of protease activity but rather driven by Rho/ROCK-mediated contractility. In contrast, under higher ECM stiffness conditions, Rho/ROCK signaling is inhibited, cells show elongated morphology and their migration relies on protease activity 4 . In another recent study it has been shown that in the absence of focal adhesions and under conditions of confinement, mesenchymal cells can spontaneously switch to a fast amoeboid migration phenotype 8 . Accordingly, we assume in our model that amoeboid-like cells move relatively fast and are not able to degrade ECM, while mesenchymal-like cells move slower and locally degrade the ECM. We further assume that cell migration speed is highest for low ECM resistance but decreases with increasing ECM resistance. This assumption is based on findings by Wolf et al. 12 migration speed gradually diminishes as a function of decreasing pore size (which is associated with increasing ECM resistance) with a steeper slope when cells migrate in an amoeboid fashion. An optimal migration speed has been observed during migration through tissue of sufficient to high porosity (which corresponds to low ECM resistance). In our coarse-grained model, we do not explicitly account for cell shape changes or Rho/ROCK-and integrin-dependent signaling pathways.
Furthermore, we assume that the direction of cell movement is determined by a constant chemotactic gradient. Cells respond to this external gradient with a certain intensity, independent of their migration phenotype.
Our CA model is described on a two-dimensional regular lattice ⊂ S Z 2 . Each lattice site is in one of a set of states (η, μ) ∈ {0, 1, 2} × [0,1], with η denoting the cell state value and μ the ECM state value. We interpret the cell state value 0 as an unoccupied lattice site, cell state value 1 as an A-cell and cell state value 2 as an M-cell. The ECM state value represents the physical resistance of the ECM.
The time evolution of our model is defined by the following rules: • Phenotypic switch. A cell changes its phenotype depending on the local ECM resistance μ. In particular, an A-cell changes its migratory phenotype with rate αμ and an M-cell with rate β(1-μ), respectively, see Fig. 1(a). • ECM degradation. M-cells locally degrade the ECM with constant rate δ>0.
• Cell migration. A-cells migrate with rate λ A (μ), M-cells with rate λ M (μ), where both λ A (μ) and λ M (μ) depend on the local ECM resistance μ, see Fig. 1(b). The direction of movement is determined by a chemotactic gradient. Cells respond to this gradient with an intensity κ > 0. If κ = 0, cells perform independent random walks, which is equivalent to the case when no gradient field is present.
At each discrete model time step, one cell is selected at random and updated according to the rules (R1)-(R3). On average, each cell is updated once per n consecutive time steps if there are n cells on the lattice. The sequence of n time steps corresponds to one Monte Carlo Step, which is the time unit of our system.
Cell migration, rule (R3), is implemented as an exclusion process 27 , which means that cells move to neighboring nodes if the target node is empty, otherwise the movement is aborted. Thus, the model accounts for steric interactions between cells due to volume exclusion. A detailed description of the model update rules (R1)-(R3) can be found in Supplementary Section 1. A schematic illustration of the update rules is depicted in Fig. 2.
In each simulation, we track individual cells and measure the migration distance (d) from the initial position for each individual cell. At the population level, the model observable characterizing invasion dynamics is the migration distance from the initial position averaged over the entire cell population (d p ). In addition, we measure the maximum migration distance (d max ) of the cells in the population in each simulation. We consider different simulations to compare (d p )-values (respective (d max )-values) between populations of cells where each  cell is allowed to switch between ameoboid and mesenchymal migration modes and those of non-switching cell populations (Δd p , respective Δd max ). Calculation details are given in Supplementary Section 2.
Simulation study. The CA model is simulated on a rectangular lattice with dimensions S 1 × S 2 . The vertical length of the lattice is given by the r 1 coordinate, ≤ ≤ r S 1 1 1 and the horizontal length of the lattice is given by the r 2 coordinate, ≤ ≤ r S 1 2 2 . For the simulations we chose initial and boundary conditions such that we can analyze cell migration using one-dimensional, vertically averaged, cell density profiles. Periodic boundary conditions are applied along the horizontal boundaries (r 1 = 1 and r 1 = S 1 ), and reflecting boundary conditions on the vertical boundaries (r 2 = 1 and r 2 = S 2 ). The initial conditions of the model are illustrated in Fig. 3. At the initial time, 50% M-cells and 50% A-cells are placed at random along the left border at r 2 = 1, as depicted in Fig. 3(a). Notice that the results do not depend on the initial fractions of M-and A-cells but on the phenotypic switch ratio α/β (data shown in Supplementary Fig. S2). Figure 3(b) shows the chemotactic gradient field, which controls the direction of cell movement. The chemotactic gradient is chosen such that the one-dimensional cell density profile advances along the S 2 -axis. In particular, the gradient is specified by a linear function which describes an attractive signal from the target location r 2 = S 2 , see Fig. 3(b). We consider different chemotactic responsiveness parameters κ ≥ 0, reflecting a range from undirected cell movement (κ = 0) to a high response to the presence of a chemotactic gradient (κ ≫ 0). For the ECM resistance, different initial scenarios are implemented. Homogeneous ECM is modeled by a constant ECM resistance, as illustrated in Fig. 3(c). Heterogeneous ECM is modeled by a sinusoidal function  (1). The color code indicates the chemotactic gradient concentration function, where blue corresponds to low, red is related to high concentrations. The arrows illustrate the local gradient direction. (c-e) Initial ECM resistance distribution of differently structured ECM: (c) homogeneous, moderate ECM resistance distribution, modeled by a constant value μ(r 1 ,r 2 ) = 0.5; (d) heterogeneous, weakly structured ECM resistance distribution, defined by equation (2) with heterogeneity parameter θ = 0.1; (e) heterogeneous, highly structured ECM resistance distribution, defined by equation (2) with heterogeneity parameter θ = 0.5. The gray scale in (c-e) indicates the level of matrix resistance on cell migration: low (white) to high (black).
where ξ~ (0, 1) is a fixed uniformly distributed random variable and θ ∈ [0,1] is a heterogeneity parameter. On average, the ECM resistance is μ(r 1 ,r 2 ) ≈ 0.5. The heterogeneity parameter controls the level of ECM structure, ranging from completely random (θ = 0) to completely coherent structures (θ = 1). An example of a less structured ECM is shown in Fig. 3(d), a highly structured ECM is illustrated in Fig. 3(e). We have also tested other functional representations of the ECM resistance distribution. The results we describe below remain valid for other functional forms of μ, as long as μ is not a constant but modeled by a spatial function μ: S 1 × S 2 → [0, 1] that depends on the spatial position (r 1 , r 2 ) and on a random component, where the latter adds local irregularities to an otherwise regular spatial structure. See Supplementary Fig. S1 for further examples.
The A-and M-cell migration rates are modeled by sigmoidal functions where C A and C A are constants, chosen such that the migration speed of A-cells is higher than that of M-cells for low to moderate ECM resistance, see Fig. 1(b). For μ = 0.5, the slopes of the migration rate functions are steepest. The sigmoidal shape of the migration rate functions accounts for three levels of matrix resistance on cell migration: low, middle and high. Cell migration is unhindered for low ECM resistance, linearly decreases for medium ECM resistance and is impeded for high ECM resistance. The focus of our study is on the impact of the switch parameters α and β on the migration behavior of the cell population. We refer to a switching cell population if cells are allowed to change their phenotype (α, β > 0), and to a non-switching cell population if cells do not adapt their phenotype (α = β = 0). For switching populations, we change the phenotypic switch parameters α = 1/10, …, 1, β = 0, β = 1 and α = 1, β = 1/10,…,1, such that the phenotypic switch ratio is α/β ∈ {1/10, 2/10, …, 1, 2, …, 10}. A low switch ratio means that cells predominantly switch to A-type, while a high switch ratio leads to a preferred switch to M-type cells. We measure how far the switching population migrates from the start position. Notice that this migration distance does not depend on the particular choice of α and β as long as the switch ratio α/β remains fixed (simulations not shown). In order to identify novel behaviors under the influence of the phenotypic switching, we compare our results with the non-switching situation α = β = 0, that is when A-and M-cells do not change their phenotype. For the non-switching population, we simulate the same number of cells as in the switching scenario, with changing ratio of A-and M-cells. The fraction of M-cells in a non-switching population is denoted by γ. Notice that the non-switching situation reflects a modeling scenario similar to the one studied by Hecht et al. 6 .
In our simulations, we compare the migration behavior of switching and non-switching cell populations under different environmental conditions: homogeneous or heterogeneous ECM structure in combination with different chemotactic responsiveness. A summary of all possible simulation scenarios is given in Table 1. In each scenario, we compare the average migration distance d p and the maximum migration distance d max of switching and non-switching populations of size n = 50. For the switching population, we change the phenotypic switch Table 1. Overview of simulation scenarios. We study different environmental conditions: homogeneous or heterogeneous ECM structure (left column) in combination with different chemotactic responsiveness (top row). In each scenario, we measure the average migration distance dp of switching and non-switching populations depending on the switch ratio α/β. , while keeping the A-cell migration rate fixed (c A = 1). A small ratio c M /c A indicates a large difference between the two migration rates, if the ratio c M /c A is close to 1, M-cells migrate nearly as fast as A-cells. The ECM degradation rate is set to δ = 0.1 ensuring partial but not complete degradation. An overview of the model parameters is given in Table 2.
As a final step, we study if cooperative effects in a switching population contribute to an efficient migration of the entire cell population. To this end, we compare the migration behavior of a single cell moving alone with that of an individual cell moving in a cell population of size n = 500. In particular, we first simulate n repetitions of the single cell scenario, each time measuring the migration distance d of the single cell, and at the end determine the maximum migration distance d. Secondly, we simulate a cell population of size n, using the same model parameters as in the single-cell scenario. We measure the migration distance d of each individual cell of the population and determine the maximum migration distance d Both simulations are repeated 100 times to account for stochastic fluctuations.

Results
Advantage of non-switching behavior under homogeneous ECM conditions. We first study if phenotypic switching provides an advantage in terms of cell population migration distance under homogeneous ECM conditions. Exemplary cell trajectories are shown in Supplementary Figs S3 and S4. Here, we compare the average and maximum migration distances of switching and non-switching populations under homogeneous ECM conditions, varying the level of μ(r 1 , r 2 ) = constant initially, in combination with a high chemotactic responsiveness. Figure 4 illustrates simulation results under different homogeneous ECM resistance conditions:  Table 2. Overview of all model parameters.
homogeneous, low ECM resistance ( Fig. 4(a)) and homogeneous, high ECM resistance ( Fig. 4(b)). The figure shows the average migration distance d p of the switching cell population as a function of the phenotypic switch ratio α/β. Also shown is the average migration distance d p of non-switching populations with different M-cell fraction γ ∈ {0, 0.3, 0.7, 1}, which remains constant as a function of α/β. In Fig. 4(a), for homogeneous, low ECM resistance, we observe that the average migration distance d p of the non-switching population depends monotonously on the M-cell fraction γ: the higher γ, the smaller the average migration distance d p . Similarly, the average migration distance d p of the switching population is highest for low switch ratio α/β, which means that cells predominantly switch to A-type. As the phenotypic switch ratio α/β increases, that is fewer A-cells and more M-cells are present, the average migration distance d p of the switching population decreases as the phenotypic switch ratio α/β increases. In contrast, under homogeneous, high ECM conditions, the average migration distance d p of a non-switching population increases monotonously with increasing M-cell fraction γ, see Fig. 4(b). Likewise, for the switching population we observe that increasing the switch ratio α/β, which leads to a preferred switch to M-type cells, increases the average migration distance d p . Under both homogeneous ECM conditions, low and high, we find that the possibility to switch the phenotype does not constitute a benefit in terms of migration distance of the entire cell population. Additional simulation data with the maximum migration distance d max , shown in Supplementary Fig. S5, confirms our finding. The critical ECM resistance value above which the greatest migration distance is no longer observed for the pure A-cell population but for the pure M-cell population, is approximately μ(r 1 , r 2 ) = 0.5. The observed behavior can be understood by considering the cell migration characteristics. On the one hand, for low to moderate, homogeneous ECM resistance, cell migration is not hindered and A-cells migrate with higher rate than M-cells. Thus, the higher the fraction of A-cells is, the higher the migration distance d p of the entire cell population. On the other hand, sufficiently high homogeneous ECM resistance presents a barrier for A-cells and only M-cells are able to create space to migrate through the ECM. Therefore, a high M-cells fraction increases the migration distance d p .
Migration plasticity can be advantageous under heterogeneous ECM conditions. In a next step, we determine if a heterogeneous ECM structure provides an environmental condition under which the switching population of cells migrates the greatest distance from the initial position. To this end, we compare the average and maximum migration distance of switching and non-switching populations under heterogeneous ECM conditions (μ as defined in equation (2)), varying the levels of ECM heterogeneity θ, in combination with a high chemotactic responsiveness. Under heterogeneous, weakly structured ECM conditions, we observe that the average and the maximum migration distance is obtained by the non-switching population. In contrast, under heterogeneous, highly structured ECM conditions, we find that the average and the maximum migration distance of the switching population may be lower or higher than of the non-switching population, depending on the migration rate ratio c M /c A . In the following, we present simulation results about the average migration distance d p . Qualitatively similar results obtained by considering the maximum migration distance d max are illustrated in Supplementary Fig. S6. Figure 5(a) shows the average migration distance d p of switching and non-switching populations as a function of the switch ratio α/β under heterogeneous, highly structured ECM conditions and a large  Supplementary  Fig. S6. Simulation parameters are κ = 1, δ= . 0 1. The red arrow indicates the difference Δd p between the switching population with maximum average migration distance d p with respect to varied α/β ratio and the non-switching population with maximum d p with respect to varied γ fraction, which is the observable analyzed in Fig. 7. migration rate ratio (c M /c A ≈ 1). We observe that the higher the portion of M-cells, the higher the average migration distance d p of the non-switching population. Similarly, for the switching population, the higher the switch ratio α/β, that is cells are predominantly M-type, the higher the average migration distance d p . The non-switching cell population with γ = 1 (pure M-cell population) migrates the greatest distance from the initial position. The situation is different if the migration rate ratio is small (c M /c A ≫ 1): while the average migration distance d p of the non-switching population is higher, the more M-cells are present, the maximum average migration distance d p of the switching population is obtained for α/β = 1, that is when both cell types are equally present, see Fig. 5(b). Under the latter conditions, that is the ECM is heterogeneous, highly structured and the difference between Aand M-cell migration speeds is large, we observe that the switching behavior provides an advantage in terms of cell population migration distance. Figure 6 illustrates exemplary individual cell trajectories for switching and non-switching populations under heterogeneous, highly structured ECM conditions with small migration rate ratio (c M /c A ≪ 1). Figure 6(a) shows that A-cells are locally impeded by high ECM resistance regions and move around such regions. In contrast, M-cells do not make a detour and move through highly structured areas, see Fig. 6(b). The combined benefits of fast A-cell migration in low ECM resistance regions and the ability of M-cells to locally degrade the ECM, lead to a high average and maximum migration distance for switching populations, see Fig. 6(c). The advantage of the switching behavior is no longer observed if the migration rate ratio is large (c M /c A ≈ 1), as shown in the trajectory representation in Supplementary Fig. S7. Notice that the trajectory figures point out that the variance of the average migration distance d p within the simulations is quite large (see also Supplementary Fig. S6). Such a large "within-simulation variances" is explained by the stochastic effects of the switching mechanism and environmental conditions. The "between-simulation variances", however, are almost negligible (coefficient of variation below 5%, see Fig. 5). This can be expected since there are 50 cells observed in each simulation such that the law of large numbers applies. Figure 7 shows the observed dependency of the average migration distance d p on the level of ECM heterogeneity and the cell migration rate ratio. In the phase diagram, the difference between the switching population with maximum migration distance d p with respect to varied switch ratio α/β and the non-switching population with maximum dp is illustrated. We observe that the switching behavior is favorable if the ECM is highly structured and the migration rate ratio is small. Otherwise the maximum average migration distance dp is reached by the non-switching population. We conclude that, provided the difference between A-and M-cells is sufficiently large enough, a high chemotactic responsiveness together with a pronounced physical resistance imposed by the ECM, represents an environmental situation under which phenotypic adaptation of cell migration modes is beneficial. We hypothesize that cells are faced with a trade-off between high ECM resistance and chemotactic responsiveness: the higher the chemotactic responsiveness, the stronger the attempt of cells to move in a directional manner Figure 8. An individual cell within a switching population can migrate further compared to a situation where only a single switching cell is moving. CA simulations of (i) a switching population of 500 cells, (ii) a scenario with only a single switching cell repeated 500 times, (iii) a non-switching population of 500 cells and (iv) a scenario with only a single non-switching cell repeated 500 times. The figure shows the maximum migration distance d, which refers to the maximum distance the individual cells within a population migrate from the initial position (i, iii), or to the maximum distance the single cell migrates among the 500 repetitions (ii, iv), respectively. Each simulation, (i)-(iv), is performed for 200 Monte Carlo steps and repeated 100 times to account for stochastic fluctuations. The ECM condition and model parameters are chosen as in Fig. 5, with phenotypic switch ratio α/β = 1. towards a target location. The higher the ECM resistance, the harder it is for the cells to move forward and the locally favorable direction may not be the overall best choice. While rather fast A-cell movement is impeded by high ECM resistance, slower moving M-cells create their own short paths. Thus, if cells in the population are allowed to adapt their phenotype to the local environmental conditions, they are able to follow the shortest path towards the target location and to overcome local ECM constraints.
Chemotactic responsiveness influences the efficiency of cell movement. We now compare the migration behavior of switching and non-switching populations under different chemotactic responsiveness.
To determine if the response towards the chemotactic gradient alters the simulation results described in the previous sections, we repeat the simulation study where we explored homogeneous and heterogeneous ECM conditions, this time however changing the chemotactic responsiveness. The result of this sensitivity analysis is that under homogeneous ECM conditions, as well as under heterogeneous, weakly structured ECM conditions, the migration distance of the non-switching population is higher than of the switching population, independent of the chemotactic responsiveness, as shown in Supplementary Figs S8 and S9. Under heterogeneous, highly structured ECM conditions, we find that the advantage of phenotypic switching in terms of cell population migration distance depends on the chemotactic responsiveness, such that the advantage of the switching behavior becomes more pronounced the higher the chemotactic responsiveness. In particular, if the ECM is highly structured and the migration rate ratio c M /c A is small, the difference between the migration distance of switching and non-switching populations increases with increasing chemotactic responsiveness, see Supplementary Fig. S10.

Cooperative effects increase the migration distance.
To study if cooperative effects in a switching population influence the efficiency of cell movement, we measure the maximum migration distance d of the switching population for different population sizes. In particular, the maximum distance d of a single switching cell migrates among 500 simulation repetitions is compared to the migration distance d of the farthest individual cell in a switching cell population of size n = 500, as described in the Methods section. Figure 8 illustrates simulation results under heterogeneous, highly structured ECM conditions. We observe that the maximum migration distanced is higher for a switching cell population than for a single cell which adapts its migratory phenotype. For comparison, we repeat the simulation procedure for the non-switching situation, with results also shown in Fig. 8. We find that the maximum migration distance d of the non-switching population is lower than the switching situation, and lowest in the case of a single non-switching cell. We conclude that cooperative effects are present in a switching population. We hypothesize that M-cells benefit A-cells by creating space in an otherwise impenetrable ECM. Our argument is in agreement with the findings of Hecht et al. 6 , who showed that in a non-switching-population, paths created by mesenchymal cells are used by amoeboid cells. In fact, our results show that cooperative effects are present in a switching population and allow the cells to migrate longer distances compared to the non-switching situation.

Discussion
We developed a cellular automaton model to study the impact of amoeboid-mesenchymal migration plasticity on tumor invasion. The switch between the two migration modes was assumed to depend on the local microenvironment through the biomechanical resistance imposed by the ECM. Cells in the model are able to switch between a slow mesenchymal migration mode with ECM degradation, and a fast amoeboid migration mode without ECM degradation, depending on ECM resistance. With computer simulations of the mathematical model we analyzed invasion dynamics, characterized by the migration distance of the cell population, under different environmental conditions. In particular, we distinguished spatially homogeneous and heterogeneous ECM resistance conditions, ranging from weakly structured to highly structured, in combination with different responses towards a chemotactic gradient which directs cell migration. Provided that the difference between amoeboid and mesenchymal migration speeds is sufficiently large, as observed in experiments 8 , we found that an amoeboid-mesenchymal migration plasticity provides an advantage in terms of migration distance if the spatial structure of ECM resistance is strongly heterogeneous and if migration is directed. In our model, we only distinguish two migration modes, which we term amoeboid and mesenchymal. The mesenchymal mode is characterized by slow migration speed and ECM degradation, while cells in amoeboid mode move fast and do not degrade the ECM. It is clear that this is a strong simplification of biological reality. Cells may display features of both migration modes, in particular concerning parameters that are not explicitly incorporated in our model, such as Rho/ROCK signaling and cell shape changes. However, our coarse-grained model captures the extreme cases of plasticity behavior. Therefore, we expect that our predictions remain valid also for intermediate situations.
In this paper, we assumed that physical and biomechanical properties of the ECM, such as density and stiffness, can be combined in a single parameter called ECM resistance. We are aware that this simplification does not account for the structural and molecular complexity of the ECM. However, we coarse-grain the extra-and intracellular details into a simplified model to provide a basic understanding of emerging plasticity effects. Our model is a first step to analyze consequences of migration plasticity on tumor invasion.
The chemical gradient which directs cell migration in the model is assumed to be independent of the ECM resistance. It is known that chemotaxis is important for tumor invasion 1 . Other types of directed cell migration have been observed in tumor cells, such as haptotaxis (movement along gradients of ECM ligands) or durotaxis (gradient of ECM stiffness). Our choice to consider ECM-independent chemotaxis serves as a first attempt to incorporate the effects of directed migration into an amoeboid-mesenchymal plasticity model for tumor invasion. We expect that our results hold for other than chemotactic gradients.
In the model, we assumed that amoeboid and mesenchymal cell migration speeds decrease monotonically with increasing ECM resistance. We based our assumption on findings of Wolf et al. 12 , where it was shown, by using breast cancer and fibrosarcoma cell lines, that cell migration speed does decrease monotonically with decreasing ECM pore size but does not depend on ECM stiffness. Both pore size and ECM stiffness contribute to the biomechanical ECM resistance against cell migration. However, experimental findings on the relation between cell migration speed and the ECM environment are controversial. For example, experiments with brain tumor cell lines have shown that an increase of ECM stiffness induces an increased cell migration speed 28 . A possible explanation for different experimental findings might be the type of tumor studied. Furthermore, it is well-known that while ECM porosity and stiffness can be independently controlled under in vitro conditions, in vivo, stiffening of the ECM is highly correlated with alterations of the ECM fiber structure 29 . Clarifying the particular contributions of different ECM properties to cell migration is primarily an experimental task.
In our model, we indirectly account for energetic costs of migration plasticity. Different cellular tasks as migration and matrix degradation consume energy and we assume a trade-off between the energies invested in these tasks. Consequently, cells in mesenchymal migration mode which degrade the ECM move slower than amoeboid cells which do not spend energy for proteolysis. This assumed difference in migration speeds is supported by several biological findings 8,30 . Our results on the existence of cooperative effects in switching cell populations support previous findings of Hecht et al. 6 . They showed that an amoeboid cell population can benefit from the presence of a small number of path creating mesenchymal cells in a model where metabolic cost for ECM degradation is modeled directly as trade-off with migration speed.
To our knowledge, our theoretical study is the first that highlights the effects of plasticity and heterogeneity in cell migration modes on tumor dissemination. Previous models studied the role of different types of cellular heterogeneity and plasticity for tumor growth, dissemination, and invasion 21,31-33 .
Our model is a cell-based, discrete cellular automaton model which we have analyzed with computer simulations. In addition, we have derived an approximate macroscopic description of the automaton dynamics (see Supplementary Section 3). It is given by a system of coupled nonlinear drift-diffusion differential equations. This system offers the potential to relate our discrete model to a continuum model and exploit the analytical tools available for the differential equation model. However, this analysis is beyond the scope of the current study.
Previous theoretical works have studied amoeboid and mesenchymal cell migration under different environmental conditions 6,18 . Here, for the first time, the impact of plasticity of amoeboid and mesenchymal migration modes on tumor invasion is investigated. Contrary to our initial expectation, amoeboid-mesenchymal migration plasticity is not advantageous per se. Interestingly, it is only with the combined influence of ECM heterogeneity and chemotactic responsiveness that the switching behavior provides an advantage in terms of migration distance. So far, experimental studies of amoeboid and mesenchymal migration either focus on ECM structure or on varying chemotactic responsiveness, but not on the combined influence of both environmental factors. There are currently no experimental findings we can compare our results with, although it is well-known that in real tumors the ECM is characterized by structural heterogeneity 19,34 and the presence of chemotactic gradients at the same time 1,35 . Our results emphasize an important aspect of future experimental studies on tumor invasion. In particular, our study suggests that in vitro systems aiming at unraveling the underlying molecular mechanisms of tumor invasion should take into account the complexity of the microenvironment by considering the combined effects of structural heterogeneities and chemical gradients on cell migration.