Tumour heterogeneity promotes collective invasion and cancer metastatic dissemination

Heterogeneity within tumour cell populations is commonly observed in most cancers. However, its impact on metastatic dissemination, one of the primary determinants of the disease prognosis, remains poorly understood. Working with a simplified numerical model of tumour spheroids, we investigated the impact of mechanical heterogeneity on the onset of tumour invasion into surrounding tissues. Our work establishes a positive link between tumour heterogeneity and metastatic dissemination, and recapitulates a number of invasion patterns identified in vivo, such as multicellular finger-like protrusions. Two complementary mechanisms are at play in heterogeneous tumours. A small proportion of stronger cells are able to initiate and lead the escape of cells, while collective effects in the bulk of the tumour provide the coordination required to sustain the invasive process through multicellular streaming. This suggests that the multicellular dynamics observed during metastasis is a generic feature of mechanically heterogeneous cell populations and might rely on a limited and generic set of attributes.


Introduction
Metastasis, the process during which cancer cells migrate away from a primary tumour and disseminate in other organs, accounts for more than 90% of cancer fatalities [1]. At the cellular level, malignant cell behaviour has been associated with an accumulation of gene mutations and a malfunction of key regulatory signalling pathways [2,3] such as TGF-β or Twist [4,5]. Moreover, following earlier pioneering works [6], recent analyses have established that most solid tumours are not homogeneous but a mosaic of various cell populations with different genetic and phenotypic traits [7]. Owing to the variability in these factors across different pathologies, patients and even tumours, no unique mechanistic pathway leading to the induction and progression of metastasis has been yet characterized [8,9].
By contrast, recent advances in imaging techniques at the tissue scale [10][11][12][13], along with histopathological studies on patients' biopsies [13], draw a surprisingly unified picture of cancer invasion in terms of cell migration behaviours. In most cancers, despite tremendous variations in gene expression and biochemical environment, tumour cells are able to invade collectively by maintaining intercellular coordination [13][14][15][16][17][18], generating multicellular structures such as cell sheets, strands or clusters remaining cohesive and polarized [14]. The accumulated knowledge about collective cell migration [19] and its quantitative biophysical models [20,21] is therefore likely to inform the question of tumour invasion at a system level. Recent results have also highlighted a link between tumour heterogeneity and metastatic dissemination through collective invasion. Experimental evidence suggests that non-invasive or poorly invasive tumour cells can be driven to invade by invasion-competent cancer cells [17,18,[22][23][24] or stromal fibroblasts [15,16,25] through a 'leader/follower' invasion mechanism [13,26].
In this article, we aim to address the question of the impact of tumour heterogeneity on malignancy. We use a minimal model for collective cell migration in order to identify both the individual and synergistic effects of cell coordination and cell heterogeneity on tumour invasion.

Modelling approach
Our model involves a Cellular Potts algorithm [27] including a self-propelled term to account for the active motion of cells [28][29][30]. It builds on the computational framework introduced previously [30] to study the collective dynamics of cell populations invading two-dimensional resistive environments, most often at the planar interface between two different tissues, along tracks created by blood vessels, myofibres or nerve bundles [12,31].
This approach is appealing, as the Cellular Potts algorithm has already been used to model a variety of tissue behaviours associated with cancer invasion [32][33][34] or heterogeneous collective migration during angiogenesis [35,36]. Furthermore, our model [30] has been already fully characterized in the context of homogeneous motile cell populations, and has made numerous predictions verified in experiments [37][38][39]. Details of the computational method and numerical values of the parameters used in in silico experiments are provided in the electronic supplementary material.
The key physical parameters of the model are the cell size (or area in two dimensions), its compressibility (low compared with other forces at play), the surface energy of the cell interface which controls the shear stiffness, tissue cohesion and some aspects of energy dissipation. In addition, each motile cell i exerts a motile force μ i n i acting on a planar substrate, where n i (t) defines the direction of its front-rear polarity, and μ i its active motile force magnitude. The polarization vector n i (t) also evolves over time and is set along the direction of the mean previous displacements of the cell in the time interval [t − τ , t]. It represents the time-scale at which cell polarity responds to external stimuli and controls the persistence length of the cell trajectory in the absence of other cells [28,30]. Simulating the evolution of the system relies on a stochastic Monte-Carlo algorithm (cf. electronic supplementary material), which defines the arbitrary time unit used in this paper, the Monte Carlo Step (MCS). The associated noise or stochasticity level is also a parameter, although it has been found not to play a critical role in the range of surface energy and motile force explored here [30].
Earlier numerical work on migratory cell populations [30] has evidenced two behavioural transitions in cell dynamics (cf. electronic supplementary material, table S1). Keeping all parameters but the motile force μ constant, the system exhibits: (i) a critical force μ s needed for a single cell to migrate through a non-motile cell population; (ii) in homogeneous cell populations, a transition from static epithelial behaviour to collective streaming, at a critical value of the motile force parameter μ c < μ s ; and (iii) a regime of motile forces μ c < μ < μ s associated with large spatial and temporal correlations of collective cell behaviours. The latter is of particular interest when considering a minimal tumour model, composed of a two-dimensional spheroid of motile cells surrounded by cohesive tissues that mechanically resist invasion as often observed in vivo [12].
As introduced previously [30], homogeneous tumour spheroids here consist of bulk tumour cells of motile force μ b surrounded by a population of non-motile cells, all other parameters of the model being kept uniform in order to avoid biasing cell migration towards surrounding tissues. Cells displayed behaviours ranging from no invasion, μ b < μ c , to widespread single cell dispersal for μ b > μ s (cf. figure 1a; electronic supplementary material, movies S1, S2, S3, S4 and S5). The regime where μ c < μ b < μ s is remarkable. Invasion of tumour cells, although slow, is observed and takes the form of finger-like collective patterns (cf. figure 1a bottom-right quadrant; electronic supplementary material, movie S2), as described in many in vivo studies [10,11,13]. This indicates that although no individual cell is able to invade in this regime, collective effects occasionally lead to larger mechanical forces sufficient to enable invasion. While this might capture a fundamental mechanism involved in the early stages of collective invasion, finger formation, in these in silico experiments, is a rare event and might not account for highly malignant situations observed in vivo. Several reasons may explain this limited amount of collective invasion. The strength of collective invasion [17,40] tends for instance to be reinforced in vivo by complementary biochemical processes such as enzymatic tissue remodelling [40], a fact already reproduced in silico with some success [29]. The analysis presented below demonstrates how mechanical heterogeneities may also dramatically increase the dynamics of cancer invasion. The homogeneous case is further analysed to develop a control dataset from which the effects due to a small sub-population of mechanically stronger cells can be assessed. The model will specifically focus on the onset of invasion, a regime where cell motility is the key driver [22]. We therefore assume that cell growth and proliferation effects can be ignored.

Invasion from homogeneous tumours
In all numerical experiments, we selected an initial tumour radius r 0 of the order of 10 cell diameters, corresponding to a population of around 300 cells (figure 1a), which is a typical size observed in similar in vivo [22] and in vitro experiments [17]. The extent of invasion is characterized by the length Λ b (μ b , t) corresponding to the sum of the distances from the initial tumour boundary to each cancer cell i that is outside the initial boundary at time t [30]. This cumulative quantity encompasses both the spatial extent of invasion and the number of invading cells. Λ b varies almost linearly with time for the first thousand MCS, as shown in figure 1b in the context of a homogeneous population of motile cells. The capacity of bulk tumour cells to invade surrounding tissues can, therefore, be characterized as a rate of invasion at short time scales, Γ b (μ b ) = Λ b (μ b , t)/ t. The evolution of Γ b as a function of the cells' motile force μ b is represented in figure 1c, for t = 1200 MCS. The sensitivity and convexity of the invasion rate to the motile force of the tumour cells is a preliminary indication that heterogeneities in the motile properties of the cancer cells could have a significant effect on invasion rates and patterns.

Invasion from heterogeneous tumours
Heterogeneity is introduced in the system in the form of a small number, N f , of cells positioned at the boundary of the tumour (cf. figure 2a), with a larger motile force, μ f . This configuration is commonly found in vivo [17] and provides a reproducible starting point to study the initiation of the invasion process. We analysed the behaviour of the system for N f = 6, 12 or 24, i.e. about 2%, 4% or 8% of the total motile cell population respectively, proportions again consistent with experimental observations [15][16][17][18]. As for the homogeneous tumour case, the invasive behaviour of bulk tumour cells in the presence of these particular cells is monitored by computing the invasion length Λ f (μ b , μ f , N f , t) and the corresponding invasion rate The effect of heterogeneity on the invasion efficiency is measured as the percentage increase in the invasion rate per stronger cell: 2b). As expected, the presence of stronger cells brings no significant enhancement to the invasion efficiency for μ b ≥ μ s , where the bulk of the cancer cells are able to invade on their own. Furthermore, when both bulk and stronger tumour cells generate a motile force below the threshold for single cell invasion, there is only a marginal effect on invasion (cf. electronic supplementary material, movie S6).  By contrast, when the bulk of the tumour is in the regime of collective invasion (μ c ≤ μ b ≤ μ s ), the addition of a small number of stronger cells (μ f ≥ μ s ) has a dramatic effect. Each of them significantly increases the invasion rate (20% ≤ Ω ≤ 45%), with an optimum for μ b ≈ 0.150 and μ f ≈ 0.350.

Morphological analysis
To understand the origin of this optimum, we analysed the morphology of the tumour at the onset of invasion. Figure 2a illustrates the typical spatial arrangements between stronger and bulk tumour cells in protrusive fingers, in which stronger cells are commonly seen initiating and leading fingers (cf. electronic supplementary material, movie S7). To quantify these morphological features, a graph-based method is used to identify the contact network of cells leaving the tumour. Fingers are defined as groups connected to the tumour body that are more than two cells away from the tumour boundary. The tip is defined as the cell of the finger which is topologically the furthest away from the original tumour boundary. The length of a finger is the length of the shortest path between its tip and the original tumour boundary, measured in number of cells. Figure 3 shows the proportion of fingers led by a stronger cell, pooling all data from experiments with 6, 12 and 24 stronger cells. In the regime of collective invasion (μ c ≤ μ b ≤ μ s ), up to 70% of the fingers can be led by stronger cells, far more than what would be expected from their density at the tumour boundary (approx. 25-30%) in average. Stronger cells, therefore, promote invasion by acting as leader cells.
The trend observed in figure 3 along the μ b axis is expected. When μ b > μ s , all cells can invade and fingers are equally likely to be led by any cell type; consequently, the proportion of fingers led by a stronger cell decreases. Similarly, when μ b < μ c , bulk cells cannot migrate on their own and any of the rare emerging invading structures would be pulled by stronger cells. The fact that the proportion of fingers led by a stronger cell decreases with μ f for large values of μ f is more surprising (cf. region between points 5 and 6 in figure 3). Why would an increase in the motility of the stronger cells reduce their ability to lead fingers? Figure 4a shows that the proportion of stronger cells leaving the tumour by a specific time (approx. 4000 MCS) increases monotonically with μ f , up to more than 60% for μ f = 0.5. This is associated with an increase in the number of fingers pulled per stronger cells (figure 4b). In the range 0.175 < μ f < 0.350, this number matches the proportion of stronger cells leaving the tumour. However, the number of fingers pulled per stronger cell plateaus for μ f > 0.350, suggesting a change in the interaction between bulk   Figure 4d shows the mean contact duration between two bulk cells, and between a bulk cell and a stronger cell, as a function of μ f . These data show that, as μ f increases, there is a growing mismatch between these interaction times; a high motile force prevents the stronger cells from maintaining sustained contacts with bulk cells. This implies that, for large μ f , stronger cells detach from bulk cells too quickly to generate long fingers, leaving short fingers behind that are unable to invade further into the tissue (cf. electronic supplementary material, movie S8). These results shed light on the invasion rate data presented in figure 2b as the invasion rate is expected to scale with both the number of fingers and their length. Indeed, figure 2b can be approximatively replicated by multiplying the increase in number of fingers per stronger cell and relative increase in finger length (cf. electronic supplementary material, figures S4, S5 and S6). Therefore, the optimum invasion is a signature of the distinct requirements for (i) initiating finger like protrusions and (ii) sustaining their penetration into surrounding tissues.
Another observation from the model is the clustering of the invading cell population due to fingers breaking up from the primary tumour. These tumour cell clusters have a typical size of two to six cells and can be heterogeneous (cf. electronic supplementary material, figure S7). Although agreement with quantitative experimental data might be coincidental, it is worth noting that both the size and the degree of heterogeneity are close to the typical values measured for circulating tumour cell clusters in vivo [41][42][43].

Conclusion
By combining two key observations of in vivo cancer metastasis, collective cell behaviour and tumour phenotypic heterogeneity, the minimal model presented here demonstrates that mechanical heterogeneity has a dramatic impact on cancer invasion dynamics. The origin of such effect is not associated with a growth driven instability [44]. We anticipate, however, that growth and cell division would enhance collective invasion behaviours as has been reported [12]. In our model, the invasion rate is maximized when the population contains cells with complementary traits, combining a few invasion competent cells with a bulk population able to collectively migrate and follow leading cells away from the main tumour. Such a qualitative interpretation of the results is consistent with the fact that various forms of tumour heterogeneity lead to similar invasion patterns; stromal fibroblasts [15,16], genetic or phenotypic variability [17,18] or partial epithelial-to-mesenchymal transition [45,46] within the tumour are all able to trigger collective invasion patterns led by highly motile cells. These similarities can be explained by the fact that generic physical interactions are sufficient to reproduce leader/follower dynamics in invading finger-like multicellular protrusions. Although we considered heterogeneities in cell motile forces, we anticipate that differences in other cellular attributes such as proteolytic activity [40], adhesion [46] or cell stiffness [47] are likely to lead to similar invasion patterns when combined with collective effects in the cell population.
A general approach to the physics of heterogeneities in collective systems is still lacking. This work demonstrates that developing such a framework would represent an important contribution to understanding a number of biological processes both in animal development and pathologies.