Collective Cell Migration in a Fibrous Environment: A Hybrid Multiscale Modelling Approach

The specific structure of the extracellular matrix (ECM), and in particular the density and orientation of collagen fibres, plays an important role in the evolution of solid cancers. While many experimental studies discussed the role of ECM in individual and collective cell migration, there are still unanswered questions about the impact of nonlocal cell sensing of other cells on the overall shape of tumour aggregation and its migration type. There are also unanswered questions about the migration and spread of tumour that arises at the boundary between different tissues with different collagen fibre orientations. To address these questions, in this study we develop a hybrid multi-scale model that considers the cells as individual entities and ECM as a continuous field. The numerical simulations obtained through this model match experimental observations, confirming that tumour aggregations are not moving if the ECM fibres are distributed randomly, and they only move when the ECM fibres are highly aligned. Moreover, the stationary tumour aggregations can have circular shapes or irregular shapes (with finger-like protrusions), while the moving tumour aggregations have elongate shapes (resembling to clusters, strands or files). We also show that the cell sensing radius impacts tumour shape only when there is a low ratio of fibre to non-fibre ECM components. Finally, we investigate the impact of different ECM fibre orientations corresponding to different tissues, on the overall tumour invasion of these neighbouring tissues.


INTRODUCTION
A large proportion of cancer research is currently focused on the role of tumour microenvironment on cancer progression, and in this context particular attention has been given to the interactions between the cancer cells and the extracellular matrix (ECM) and its components [1]. The ECM is composed of water, minerals, proteoglycans and fibrous proteins secreted by various cells, and these components appear in specific percentages in each organ in the body, according to the needs of the tissues [1]. Collagen is the most abundant fibrous protein in the ECM, and its importance is not only given by its role in the ECM architecture [1], but also by its role on tumour tissue stiffness, in promoting tumour metastasis and in the regulation of tumour immunity [2]. Different tissues have different levels of collagen [3]: from ≈ 75.4μg/ mg tissue in pancreas, to ≈ 90μg/ mg tissue in lung.
Moreover, the orientation of collagen fibres plays a very important role in the mechanical properties of tissues, as well as in their physiological and biochemical functions [4]. In different tissues these fibres have different orientations addressing different functional needs [4,5]: e.g., they are arranged in parallel in the tendons [5], or in concentric rings in the bone [6], or in orthogonal grids in the cornea [7] or hypertrophied myocardium [8].
Diseases can lead to a loss of orientation of collagen fibres as well, which adopt a more random distribution; this was observed, for example, in liver fibrosis [9], in heart diseases [8], or in denervated bones with reduced mechanical stress [6]. Cancer cells are also involved in the re-orientation of collagen fibres around tumours, by changing the typical tangled and disorganised collagen fibres that can be found within the stroma, towards thickened collagen fibres that are aligned perpendicularly to the boundary of the invading tumour [1,10]. The density and orientation of collagen fibres inside the ECM plays an important role also in the evolution of cancer, as cancer cells not only move along collagen fibres [1,10], but they also remodel the ECM through degradation [e.g., via matrix metalloproteases (MMPs)] and deposition of new matrix proteins [11], which interfere with cell-cell and cell-ECM adhesion, as well as cell polarity [1].
Despite all this wealth of information regarding the role of ECM in cancer progression, there is still a poor understanding of the roles of collagen fibre orientation on the individual and collective migration of cancer cells (and the non-local interactions between cells via these stiffened fibres), how the cells remodel the ECM, and how these complex cell-fibre interactions impact the overall tumour shape and invasion pattern. In a seminal paper, Friedl and Alexander [12] have reviewed the different mechanisms of individual and collective cancer cell migration: from single cells amoeboid and mesenchymal migration, to multicellular amoeboid and mesenchymal streaming, and various types of collective cell migration (clusters, strands, or files). Friedl and Alexander [12] also mentioned the spatially-expanding tumours that undergo growth and thus passively move by pushing the surrounding tissue. In contrast to the cancer cells that actively move in a collective manner along collagen fibres, the spatiallyexpanding tumours can be found within a capsule of ECM formed of aligned collagen fibres that are circularly disposed around this growing tumour. The different modes of migration are unstable and can change upon variations in cell-cell and cell-ECM adhesion [12]. While many studies [12] discuss the various types of invasion of tumour aggregations (as a result of single cell interactions with other cells or with the ECM), it is still not fully understood how cells perceive other cells further away (although it seems that they can mechanically sense and react to the presence of other cells up to 100 µm away [13]), and how this perception can impact the overall tumour shape. Moreover, it is still not fully understood how the various tissue types can impact the migration of tumour cells and tumour aggregations (as tumours can develop at the boundaries of different tissues with different characteristics).
The goal of this study is to investigate migration cell patterns in various tissues with different levels of ECM fibres and different alignment levels, as we vary: 1) cells sensing radius, 2) cell-cell and cell-ECM adhesion strengths, 3) the orientation of ECM fibres and the ratio of fibres to non-fibres ECM components, 4) the structure of the domain, with various tissue patches that have different fibre orientations. To this end, we consider a hybrid multi-scale modelling approach where cells are modeled as discrete entities while the ECM (with its two phases: fibrous and non-fibrous) is continuous. We show that this hybrid model can reproduce a variety of cell migration types (e.g., cluster, strand or files) as well as a variety of tumour shapes: circular, elongated or irregular (with finger-like protrusions). We also show that cell sensing radius has an impact on tumour shape and migration pattern only when the ratio of fibres to non-fibres ECM components is increased from lower to higher ratios. Finally, we show that the evolution of tumours that arise at the boundaries of different tissues with different fibre orientations is influenced by the directionality of these fibres.
The paper is structured as follows. In The Multi-Scale Hybrid Model section we describe the multi-scale hybrid model for cellcell and cell-ECM interactions. In Results section we discuss variety of numerical simulations obtained with this model, as we vary the above mentioned parameters. We conclude in Summary and Discussion section with a summary and a discussion of the importance of these results.

THE MULTI-SCALE HYBRID MODEL
There are two main types of models that are often used to capture the dynamics of tumour development, namely discrete and continuous models. Both have advantages and drawbacks over the other [14] and to minimise these disadvantages, recent extensive efforts have been made to combine these models into hybrid ones [15].
We employ here a hybrid modeling framework [15] that combines the off-lattice agent based model MultiCell-LF [16][17][18] to represent the cells, and a multi-scale continuous framework [19][20][21][22][23] to represent the microenvironment. To facilitate the description of this multi-scale hybrid model, let us first introduce some useful notations from both frameworks. The model is defined within a maximal tissue cube Y ∈ R d with d 2 and time interval [0, T]. The state of an arbitrary individual cell is described by where X i (t) denotes the position of the cell center, C i rad (t) the cell radius, C i age (t) the current cell age, C i mat the cell maturation age (i.e., the age when the cell is ready for division which is assigned at birth and so it is independent of time) and C i neigh (t) denotes the number of neighbouring cancer cells at time t. Here, To connect the discrete model to the continuous one, we first need to generate a cancer cell density by using the individual cell properties. For this, let us first observe that such density can be determined by using the fraction of unit space that is occupied by the cancer cells. Therefore, at any macro-scale spatio-temporal position (x, t) ∈ Y × [0, T] the density of the cancer cell population in a square neighbourhood B(x, c Δx ) is given by where λ(·) is the Lebesgue measure in R d , I (·) (·) is the usual indicator function, and B · 2 (X i (t), C rad i (t)) : {z ∈ Y| X i − z 2 ≤ C rad i } describes the spatial region occupied by the body of an individual cell C i within the neighbourhood B(x, c Δx ) : {z ∈ Y| x − z ∞ ≤ c Δx } that is given by a · ∞ −ball of an appropriately chosen radius c Δx > 0 (which is proportional to the spatial step-size of the discretised computational domain Y).
Besides the discrete cancer cell population, here, also we consider the dynamics of a continuous multi-scale two-phase ECM. Specifically, we consider a fibre ECM phase, accounting for all major fibres (for instance collagen and fibronectin) whose macro-scale density and spatial bias are denoted by F(x, t) and θ f (x, t), respectively (for further details see Fibre Representation and Fibre Rearrangement Process). On the other hand, the ECM also consist of several other macromolecules such as non-fibrous proteins (for example amyloid fibrils), minerals, various enzymes, proteoglycans and polysaccharides. Hence, to capture the ECM as a whole, we bundle the rest of these non-fibres ECM constituents (i.e., everything that is not considered to be fibrous proteins) into the second phase that we refer to as the non-fibre ECM phase and denote it by l(x, t).
Finally, for compact notation we denote by ρ(u) the total space occupied at position x, i.e., where u represents the global three-dimensional tumour vector u : (c(x, t), F(x, t), l(x, t)) u .

The MultiCell-LF Model
For each single cell in the MultiCell-LF (Multi-Cell Lattice-Free) model, several individually-regulated life processes are included, such as cell ageing, cell growth, cell division, cell-cell and cell-ECM interactions, and cell contact inhibition.

The Cell Cycle
The lifespan of each cell is traced with the current cell age C age i (t) that progresses at the same rate as time, and cell maturation age C mat i that is assigned at the cell birth and varies slightly between the cells to avoid synchronization of cell divisions. The cell cycle is divided into the usual four phases [24]: the G1 phase (gap 1) during which the cells are growing in size, the S phase (synthesis) when biological cells replicate their DNA, the G2 phase (gap 2) in which cells complete the growth and replication processes in preparation for the M phase (mitosis) in which cells physically divide into two daughter cells. Following ours and others previous work, the length of the cell cycle is divided as follows: G1 (45% of the whole cell cycle), S (35%), G2 (15%), and M (5%), respectively [16,17,25]. Within the figures, we indicate the phase of an individual cell by different colours, i.e., we use white for the G1 phase, yellow for the S phase, orange for the G2 phase, red for the M phase and black for cells that are both in M phase and contact inhibited.

Cell Growth and Division
The radius C rad i (t) of a growing cell is increasing linearly until it reaches the size of the mature cell R max . The radius increment has been chosen to assure that the process of cell growth is completed before the M phase. Once the cell current age reaches its division age (i.e., C age i (t) C mat i ), it will divide into two daughter cells that are placed symmetrically around the mother cell's nucleus X i (t) within a distance equal to the half of the cell maximal radius [16,17]: where ϕ is an angle randomly chosen from [0, 2π], X i (t) is the position of the mother cell, and R max is the maximal cell radius. The initial ages of the two new cells are set to zero (C age i1 (t) C age i2 (t) 0), and their respective maturation ages C mat i1 and C mat i2 are inherited from the mother cell maturation age with a small noise term ϵ mat to avoid cell cycle synchronisation between the cells [17,26]: where ϵ mat is a small number drawn from a uniform distribution [−0.05, 0.05] and C mat i is the division age of the mother cell. Finally, both initial radii of the daughter cells are set to 0.65R max (i.e., C rad i1 (t) C rad i2 (t) 0.65R max ) and as a result initially the two daughter cells are overlapping Eq. 2. Therefore, right after the cell division, both daughter cells experience repulsive forces between each other, as well as with other nearby cells. These forces will push the cells apart until they reach an equilibrium and are no longer overlapping (see The Hybrid Cell Movement).

Cell Contact Inhibition
Once the whole cell colony grows in size, individual cells may become overcrowded and growth-arrested due to the contact inhibition signals from the neighbouring cells. The overcrowding condition is modelled by counting the number of cells C neigh i (t) in the specified neighbourhood of radius R neigh 4.5R max : where I (·) (·) again denotes the indicator function and SS i is the set of all cancer cells that are close to the cell C i (t), i.e., When the number of neighbouring cells reaches a specified threshold N neigh (i.e., C considered overcrowded and growth-arrested. However, such an arrested cell remains metabolically active [27], and can resume its active cell cycle when the contact inhibition conditions change.
The time spent in the growth-arrested state does count toward the length of the active cell cycle.

The Hybrid Cell Movement
To model cancer cell passive relocation (due to cell-cell interactions) and active migration (due to cell-ECM interactions), we combine both discrete and continuous approaches. While for cell-cell interactions (both repulsive and adhesive forces) we use the agent-based approach, the cell-ECM (both fibre-based and non-fibre adhesions) are modelled by using a continuous technique. Ultimately, these forces collectively influence the direction of motion of each individual cell, leading to complex tumour dynamics and to the emergence of various tumour morphologies.

Discrete Cell-Cell Repulsive Interactions
The repulsive forces between nearby cells are introduced to maintain cell volume and to avoid cell overlapping during its movement or division (see Cell Growth and Division). Following our previous work [17,18], these forces are modelled as linear Hookean springs. Hence, considering two arbitrary but distinct cells C i (t) and C j (t) (with i ≠ j), the repulsive force between them is given by: where F rep > 0 is the constant spring stiffness and the spring resting length is set be equal to the sum of the two cells' radii, i.e., to C rad i (t) + C rad j (t). Taking into consideration now the whole cancer cell population, the overall repulsive force that acts on the cell C i (t) is given by the sum of all repulsive forces between the cell C i (t) and any other overlapping cell C j (t), with j ∈ {1, . . . , M(t)}\{i}. This cumulative force is given by

Discrete Cell-Cell Adhesive Interactions
The adhesive forces are activated between non-overlapping but nearby cells in order to keep the cell cluster compact. These forces are modelled using Hooke's law with a constant spring stiffness F adh > 0 and resting length of 2 R max [17]. Hence, such forces are activated between cells separated by the distance larger than 2R max , but not exceeding 2.25R max . This prevents from over-activation of the cell-cell adhesion between cells that are located far off each other. Therefore, considering two arbitrary cells C i (t) and C j (t) (with i ≠ j), the cell-cell adhesion force acting between them is given by Similarly to the repulsive forces, each cell can be subjected to multiple adhesive forces, thus the cumulative adhesive force which acts on cell C i (t) is given by

Continuous Cell-ECM Adhesive Interactions
Besides the cell-cell interactions described above, of particular importance, are the cell-ECM adhesions [28][29][30][31] that we explore here through a cell-non-fibre ECM adhesion [32][33][34][35] as well as a cell-fibre ECM adhesion [36,37]. In the existing literature [19][20][21][22][23][38][39][40][41][42], this type of interaction is usually modelled by a non-local adhesion integral with a sensing region B(0, R) of radius R. Since in our hybrid model the two phase ECM is modelled in a continuous manner using densities, we can adopt this approach to describe the present cell-ECM interactions. To this end, let us define the cell-ECM adhesion force/velocity for an arbitrary cell C i as where R represents the maximum range within which a cell C i can establish adhesive bonds with the surrounding ECM constituents, i.e., within B(0, R) and S cl > 0 and S cF > 0 are assumed to be the constant cell-non-fibre ECM and cellfibre ECM adhesion strengths, respectively. Furthermore, in Eq. 5 n(·) is the unit radial vector n y : and similarly n(·, ·) is the unit radial vector biased by the orientation of the fibres, i.e., it is given by is the macro-scale orientation of the fibres. For further details about the fibre orientations see Fibre Representation section. Also, in Eq. 5 to account for the gradual weakening of these adhesive bonds as we move away from the centre of the cell X i within the sensing region B(0, R), we simply use a radially symmetric kernel K(·) which is given by where ψ(·) is the standard mollifier defined as Finally, in Eq. 5 [1 − ρ(u)] + : max(0, 1 − ρ(u)) ensures that any overcrowded space within the sensing region B(0, R) does not influence the overall cell-ECM adhesion.
For completeness, here we also briefly discuss the numerical approach for the cell-ECM adhesion integral F ECM i given in Eq. 5. To this end, we adopt the partitioning of the sensing region B(0, R) from Shuttleworth and Trucu [19], and so we divide B(0, R) into N s annulus sectors shown in Figures 1C,D.
Thus, as shown in Figure 1C, the sensing region B(0, R) is split into s annuli and then each of these annuli is further split into 2 m+(k−1) uniformly distributed sectors, i.e., the number of annulus sectors is given by Further, let us now denote these annulus sectors by S ] , with ] 1, . . . , N s , and their corresponding off-grid barycentres by b S ] , as shown in Figure 1D. Then, by using simple bi-linear shape functions, on each annulus sector S ] we are able to approximate the mean-values of the non-fibres ECM, the fibres ECM, as well as their associated macro-scale orientations. Ultimately, this enables us to appropriately construct the integral of the step functions associated with each annulus sector S ] whose value is given by a linear combination of the mean-values. Finally, these mean values are used to obtain the numerical approximation for F ECM i by following the numerical steps described in full in [19,22,23].

The Overall Direction of Cell Movement
Ultimately, the overall movement direction of a cell C i (t) depends on each of the forces discussed above (repulsion, cellcell adhesion and cell-ECM adhesion). To this end, the spatial dynamics of a cell C i (t) is based on Newton's second law of motion where we assume that each cell returns to its equilibrium without oscillations. Further, the overall force acting on a cell C i (t) is assumed to be proportional to the velocity of the cell and so the dynamics is given by [16][17][18].

The Continuous Multi-Scale ECM
Tumour development is a complex and multi-scale phenomena where the macro-scale events are accompanied by several related micro-scale processes. One of the most important of these processes is the re-distribution of the fibres ECM micro-constituents triggered by various cell-generated forces. Hence, here we first describe the macro-scale dynamics of the continuous two-phase ECM and then we detail its micro-scale dynamics.
Similarly to the cancer cell density, to achieve connection between the discrete and continuous parts of the model, we first need to appropriately approximate the generated cell forces (F where again λ(·) is the Lebesgue measure in R d , B(x, c Δx ) is the · ∞ ball used also in Eq. 1, I (·) (·) denotes the indicator function.

The Two-phase ECM
Besides the cancer cells, in this work we also take into account the dynamics of a two-phase ECM to capture not only the fibrous proteins in the fibres ECM, but also every other ECM constituents which are collected into the non-fibre ECM phase. It is well-known that MMPs are responsible for degrading the fibrous ECM proteins [43,44], and since cancer cell do produce several types of MMPs [45,46], here we consider the fibres ECM to be degraded by the cancer cell population at a constant rate β F > 0. However, existing biological evidence also suggest that besides the fibrous proteins, components of the non-fibre ECM phase can also be degraded by several classes of MMPs (for instance amyloid fibrils [47,48]). Hence, we can assume a constant non-fibre ECM degradation rate β l > 0 by the cancer cells. The non-dimensional macro-scale dynamics of both fibre and non-fibre ECM phases are mathematically formalised as Here, we used the macro-scale density of the cancer cell population c(x, t) (see Eq. 1), which was derived using the properties of the individual cells provided by the agent-based part of our hybrid model.

Fibre Representation
Let us now focus our attention to the two macro-scale representations of the fibres ECM, namely its amount F(x, t) and spatial bias θ f (x, t). Considering a micro-scale distribution of the micro-fibres f (z, t) on a cell-scale microdomain δY(x) : x + δY of appropriate micro-scale size δ > 0, we observe that both of these two characteristics of the fibres ECM can be captured through a vector field representation θ f (x, t) of the fibres. This captures the naturally emerging macro-scale orientation of the ECM fibres that is induced by the distribution of micro-fibres constituents f (z, t) at micro-scale (on δY(x)), and represents the spatial bias of the ECM fibres distributed at x to withstand incoming spatial cell fluxes generated by the collective cell migration. Following Shuttleworth and Trucu [19], this vector field representation of the ECM fibre phase is given by where θ x,δY(x) (x, t) is the Bochner-mean-value [49] of the vectorvalued function δY(x) ∋ z1z − x ∈ R d with respect to the measure f (z, t)λ(·) and represents the uniquely emerging revolving barycentral orientation (which was originally derived and justified in detail in Shuttleworth and Trucu [19]) induced by the micro-fibre distribution f (z, t) on δY(x), this being given by In this context, the second characteristic of the fibres ECM, i.e., its amount F(x, t) is precisely given by the Euclidean norm of the vector field θ f (x, t), namely which in fact represents the mean-value of the micro-fibre distribution on the micro-domain δY(x). Since both of these macro-scale fibre ECM representations are derived by using the micro-scale distribution of the micro-fibres f (z, t), there is an emerging link between the two scales which we refer to as the fibre bottom-up link and illustrate it in Figure 1A. Furthermore, the detailed micro-scale rearrangement process is illustrated in full in Figure 1B.

Fibre Rearrangement Process
As the individual cancer cells move and invade, they interact with the micro-fibre structure and thereby push them in the direction of travelling, ultimately resulting in the rearrangement of the micro-fibres. Following Shuttleworth and Trucu [19], here we consider this complex process to be initiated by the spatial flux of the cancer cell population F (x, t) that was defined in Eq. 7. Then, such flux that acts uniformly on δY(x) naturally gets regulated in a weighted manner by the orientation of the fibres θ f (x, t), resulting in a micro-fibres rearrangement vector where the weight ω(x, t) is the appropriate mass fraction of the cancer cells and fibres ECM given by .
Ultimately, the new position z* of each micro-fibre z ∈ δY(x) is calculated by using a reallocation vector ] δY(x) (z, t) which accounts not only for the rearrangement vector r(δY(x), t) in Eq. 10 but also upon the degree of alignment between r(δY(x), t) and the barycentral position vector x dir (z) : z − x as well as it utilises the level of micro-fibres at position z. Hence, this reallocation vector is given by where f max > 0 is the maximum possible level of fibres at any micro-point z, f p : f (z, t)/f max is the fibres saturation level and χ {f (z,t) > 0} is the characteristic function of the micro-fibres support. Therefore, by using Eq. 11 the new position z* is given by and the amount of fibres that is reallocated from z to z* is determined by a movement probability which tracks the free space available at the new position z*. Thus, the level of fibres that gets moved to z* is given by p move · f (z, t) and the rest of the fibres, i.e., (1 − p move ) · f (z, t) remain at the current position z. In Figure 1B we illustrate the reallocation of a micro-fibre using a typical example of the rearrangement r(δY(x), t) and reallocation vectors ] δY(x) (z, t) defined in Eqs. 10, 11, respectively. Finally, since this rearrangement process that redistributes the fibres ECM micro constituents is initiated by the macroscale cancer cell flux F (x, t), the macro and micro-scales are connected via a top-down link which is also illustrated in Figure 1A.

RESULTS
In the context of the above described multi-scale hybrid framework, here we present some numerical simulations of the model, which highlight different migration cell patterns in various tissues with different levels of ECM fibres and different alignment levels. To this end, we start each simulation with a single cancer cell with well defined properties located at a point (x 0 1 , x 0 2 ) (that will be defined for each simulations) of the computational domain Y [−1280mm, 1280mm] × [−1280mm, 1280mm], and any alteration from this will be stated accordingly. Then, this single cancer cell is considered to be embedded within the following (scaled) non-fibre ECM environment [19][20][21][22][23]: l(x, 0) : which can be seen in Figure 3A that is in-silico representation of the heterogeneous distribution of the non-fibre ECM. For illustrative purposes, in Figure 3  To shed some light on the importance of ECM characteristics on cancer invasion, throughout this section we investigate numerically the type of cancer migration and invasion patterns obtained when we consider different cell sensing radii. We note here that the sensing radius is determined not only by the length of cell membrane protrusions called pseudopodia (which have lengths greater than 2 µm [50,51]), but also by the long-range cell sensing due to stress transmission via aligned ECM fibres (which allows cells to sense other cells up to 100 µm away [13]). Throughout this study we investigate the impact of two sensing radii, R 30μm or R 50μm, on tumour shape and invasion pattern. We also investigate the impact of different fibres to nonfibres ECM ratios (i.e., 10: 90%, 20: 80%, or 30: 70%), different ECM fibre distribution (i.e., random, aligned, or mixed), and different cell-cell and cell-ECM interaction strengths.

Simulations Without Cell-ECM Adhesions in a Random Fibrous Environment
We start our numerical investigation of this hybrid multi-scale model by showing some baseline simulations for the case without cell-ECM adhesions, as we vary the fibre ECM density (i.e., 10 − 30% fibres), for two different sensing radii: R 30μm in Figure 4, and R 50μm in Figure 5.
By comparing the numerical results in Figures 4, 5 we can conclude that in the absence of any cell-ECM adhesion, the tumour colonies have an almost circular shape, and they are stationary (i.e., they do not move through the domain). Moreover, the larger cell sensing radius does not seem to have a significant impact on tumour structure.

Simulations With Cell-ECM Adhesions in a Random Fibrous Environment
Next, we consider cell-ECM adhesion forces (S cl 0.01, S cF 1.8), and investigate again the effect of changes in fibre    ECM density (i.e., 10-30-% fibres), for two different sensing radii: R 30μm in Figure 6, and R 50μm in Figure 7. Since now we consider cell-ECM interactions, in these two figures we also show model dynamics when the ECM environment is described by random fibres   The increase in sensing region (from R 30 in Figure 6 to R 50 in Figure 7) leads to an important difference. Even though we have now introduced the cell-ECM adhesion, in a lower fibre ECM density (10-20%) the tumour still kept its more round shape due to the increase of the sensing region. One possible explanation for this is that by increasing the sensing region in a random ECM environment, the different forces that emerge from the fibre orientations, cancel each other out and so the cancer cells are not subject to a strong leading force.
Finally, the increase in cell sensing radius to R 50 in an aligned fibrous environment, leads to a more compact tumour aggregation for the case of 10: 90% fibres: non-fibres ECM ratios. For higher ratios, the tumour is extremely elongated and the sensing radius has almost no impact.

Decreasing the Cell-Cell Adhesion Strength
Previous simulations were performed with a relatively high cell-cell adhesion strength (F adh 0.12), which impacts the overall shape of tumour aggregations. Let us now decrease the cell-cell adhesion strength and consider F adh 0.04 (while keeping all other parameters as before). In Figures 8, 9 we can see that a decrease in cell-cell adhesion causes more irregularly-shaped tumour aggregations. This is counterbalanced a bit by an increase in the sensing radius to R 50 which, as before, tends to favorise more round-shaped tumour aggregations.

Decreasing the Cell-Fibre ECM Adhesion Strength
Since the previous simulations had low cell-cell adhesion but high cell-fibre ECM adhesion, next we investigate tumour invasion  patterns when we lower the cell-fibre adhesion strength to S cF 1.2 (while keeping F adh 0.04). In Figures 10, 11 we see that lowering cell-fibre adhesion leads to a highly-irregular tumour shape, especially for the case of 30% : 70% fibres to non-fibres ECM ratios. Consider now an even lower cell-ECM adhesion strength, S cF 0.3, and at the same time a lower cell repulsion F rep 0.5. In Figures 12, 13 we see that tumour aggregations have recovered their circular shape when the ECM fibres are distributed randomly, and have either circular or slightly elongated structures when the ECM fibres are aligned.
We emphasise that the simulations in sub-panels (B')-(C') in the above figures were ran for shorter times, since for longer times the cells leave the domain.

Different Tissues With Different Fibre Orientations
Finally, we investigate numerically what happens with tumour shape and its invasion pattern when the ECM domain is formed of patches of random and aligned fibres, corresponding to different types of tissues (see the discussion in the Introduction). For this we start with a cluster of cells (e.g., 19 cells, for illustration purposes) instead of just one cell. Moreover, we consider here again the parameter values used for the results in Figures 6C-C9, namely: R 30μm, F rep 0.9, F adh 0.12, S cl 0.01, S cF 1.8, and 30% : 70% fibres and non-fibres ECM ratios. In Figure 14 we present six scenarios corresponding to different regions of random and aligned ECM fibres. The simulations in Figure 14 show that the initial position of the cluster of tumour cells, together with the ECM alignment in that region, influences the direction of migration of these cells. In sub-panels (A)-(B) the cells migrate in the direction of fibre alignment, while in sub-panel (C) the cells migrate along the boundary between two regions with opposite alignment. A similar migration pattern is observed in sub-panel (D). Finally, in sub-panels (E) and (F) we observe that when the original tumour cell cluster is positioned in an area surrounded by different ECM alignment, the tumour tends to stay in that area and to very slowly invade the neighbouring tissues (with different ECM orientations). All these different tumour invasive patterns are consistent with the experimental results in [10,52], which showed that aligned collagen fibres perpendicular to tumour boundary are associated with tumour invasion along those fibres, and the experiments in [52] which showed that aligned collagen fibres parallel to tumour boundary impede invasion.

SUMMARY AND DISCUSSION
The composition and structural characteristics of the extracellular matrix (ECM) are known to vary widely among different tissues [11], and this has a significant impact on the evolution of cancer. Since there is an increasing need for understanding the spatiotemporal changes in ECM and their roles in cancer progression, in this study we introduced a new multi-scale hybrid mathematical model for cell-cell and cell-ECM interactions, and used it to investigate numerically the impact of ECM fibre orientation and fibre density on cancer cell invasion patterns. To   this end, we varied a number of parameters associated with the cell-cell interactions (e.g., cell sensing radius R, cell-cell adhesion stiffness F adh , cell-cell repulsion stiffness F rep ), and parameters associated with cell-ECM interactions (i.e., cell-fibre adhesion coefficient S cF , and cell-non-fibre adhesion coefficient S cl ), as well as we varied the random or aligned distribution of fibres. This enables us to explore their impact on the formation of cell aggregations and their collective dynamics.
Through numerical simulations, we tested a variety of scenarios: from the impact of smaller/larger cell sensing radius on tumour aggregations, to the impact of random vs. aligned fibres on the migration and invasion of cancer cells into nearby tissues, the role of different fibre orientation in the neighbourhood of initial tumour, and the roles of cell-cell and cell-ECM strengths on the shape of solid tumours. We summarise all these numerical results in Table 1 (for the case R 30μm) and Table 2 (for the case R 50μm). By comparing the 3 rd column (Description) with the 10th column (Fibre orientation) it is clear that a random ECM fibre distribution leads to stationary tumour aggregations, while an oriented fibre distribution leads to moving tumour aggregations. These numerical results are consistent with experimental results showing that oriented collagen fibres direct tumour invasion, with the alignment of collagen fibres coinciding with the cell invasive direction [10,53]. Moreover, our numerical results in Figure 14 showed that collagen fibres that are relatively parallel to tumour boundary (and create a capsule-like orientation) can limit tumour invasion. These results are consistent with a recent note by Friedl [54], as well as the experiments in [52], about to the importance of ECM orientation relative to the tumour (see also Figure 1 in [54]).
The results in Figure 14 suggest that if we could know the orientation of collagen fibres surrounding the tumour (which may be spanning different tissues with different ECM structures) we could predict the fast or slow evolution of the solid tumour, as well as potential cell migration directions. Since collagen can be easily visualised in standard histopathology slides through second harmonic generation microscopy [10] or via multiphoton tomography [55], such predictions on the evolution of tumours might help treatment decisions, e.g., by deciding to resect more surrounding tissue if the imaging shows collagen fibres aligned perpendicular to tumour boundary in a certain area, or by deciding to preserve more surrounding tissue if the imaging shows collagen fibres parallel to tumour boundary. Intraoperative visualisation of tumour microenvironment, including collagen alignment in resected tissues [56], could be further combined with mathematical simulations in 2D (and even in 3D, although this is expected to carry a considerably higher computational cost) to improve such treatment decisions.
Finally, our numerical simulations suggest that the sensing radius R (determined by the length of membrane protrusions called pseudopodia [50,51], as well as by the long-range sensing due to stress transmission via aligned ECM fibres [13]) impacts the shape of tumour colonies in an aligned fibrous ECM environment with relatively low ratios of fibres to non-fibres. More precisely, when the ECM ratio of fibres to non-fibres is 10%: 90%, larger R leads to more compact tumor cell colonies (compare with the cluster invasion pattern in [12]). For higher fibres: non-fibres ratios this sensing radius has almost no impact, the tumour being extremely elongated as the cells move quickly along the collagen fibres (compare with the strands and files patterns in [12]).
Overall, this study not only confirmed some of the previous experimental results on the importance of alignement of ECM collagen fibres on tumour invasion [10,52], but also proposed new hypotheses on the biological mechanisms involved in the shape of the tumour colonies (e.g., the roles of sensing radius vs. fibres to non-fibres ratios; the roles of sensing radius vs. cell-cell and cell-ECM mechanical forces). Moreover, this multi-scale hybrid modelling, which incorporates ECM remodelling and fibre rearrangement in the direction of cell travelling (via cell flux F (x, t)), can be further applied to understand various other fibrotic diseases characterised by ECM remodelling [57].

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

APPENDIX A: PARAMETER SET
In Table 3 we summarise the parameter values (and parameter ranges) used throughout the numerical simulations performed in this study.