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

Abstract

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 mg tissue in pancreas, to 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 spatially-expanding 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 cell-cell 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–18] to represent the cells, and a multi-scale continuous framework [19–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 with and time interval . The state of an arbitrary individual cell is described by , where denotes the position of the cell center, the cell radius, the current cell age, 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 denotes the number of neighbouring cancer cells at time t. Here, where is the number of cancer cells at time .

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 the density of the cancer cell population in a square neighbourhood is given bywhere is the Lebesgue measure in , is the usual indicator function, and describes the spatial region occupied by the body of an individual cell within the neighbourhood that is given by a ball of an appropriately chosen radius (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 and , 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 .

Finally, for compact notation we denote by the total space occupied at position x, i.e.,where u represents the global three-dimensional tumour vector

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 that progresses at the same rate as time, and cell maturation age 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 ( of the whole cell cycle), S (), G2 (), and M (), 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 of a growing cell is increasing linearly until it reaches the size of the mature cell . 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., ), it will divide into two daughter cells that are placed symmetrically around the mother cell’s nucleus (t) within a distance equal to the half of the cell maximal radius [16, 17]:where ϕ is an angle randomly chosen from , is the position of the mother cell, and is the maximal cell radius. The initial ages of the two new cells are set to zero (), and their respective maturation ages and are inherited from the mother cell maturation age with a small noise term to avoid cell cycle synchronisation between the cells [17, 26]:where is a small number drawn from a uniform distribution and is the division age of the mother cell. Finally, both initial radii of the daughter cells are set to (i.e., ) 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 in the specified neighbourhood of radius :where again denotes the indicator function and is the set of all cancer cells that are close to the cell , i.e.,

When the number of neighbouring cells reaches a specified threshold (i.e., ), the cell (t) is 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 and (with ), the repulsive force between them is given by:where 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 . Taking into consideration now the whole cancer cell population, the overall repulsive force that acts on the cell is given by the sum of all repulsive forces between the cell and any other overlapping cell , with . 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 and resting length of 2 $R_{\text{max}}$[17]. Hence, such forces are activated between cells separated by the distance larger than , but not exceeding . This prevents from over-activation of the cell-cell adhesion between cells that are located far off each other. Therefore, considering two arbitrary cells and (with ), 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 is given by

Continuous Cell-ECM Adhesive Interactions

Besides the cell-cell interactions described above, of particular importance, are the cell-ECM adhesions [28–31] that we explore here through a cell-non-fibre ECM adhesion [32–35] as well as a cell-fibre ECM adhesion [36, 37]. In the existing literature [19–23, 38–42], this type of interaction is usually modelled by a non-local adhesion integral with a sensing region 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 aswhere R represents the maximum range within which a cell can establish adhesive bonds with the surrounding ECM constituents, i.e., within and and are assumed to be the constant cell-non-fibre ECM and cell-fibre ECM adhesion strengths, respectively. Furthermore, in Eq. 5 is the unit radial vectorand similarly is the unit radial vector biased by the orientation of the fibres, i.e., it is given bywhere 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 within the sensing region , we simply use a radially symmetric kernel which is given bywhere is the standard mollifier defined as

Finally, in Eq. 5 ensures that any overcrowded space within the sensing region does not influence the overall cell-ECM adhesion.

For completeness, here we also briefly discuss the numerical approach for the cell-ECM adhesion integral given in Eq. 5. To this end, we adopt the partitioning of the sensing region from Shuttleworth and Trucu [19], and so we divide into annulus sectors shown in Figures 1C,D.

FIGURE 1

Thus, as shown in Figure 1C, the sensing region is split into s annuli and then each of these annuli is further split into uniformly distributed sectors, i.e., the number of annulus sectors is given by

Further, let us now denote these annulus sectors by , with , and their corresponding off-grid barycentres by , as shown in Figure 1D. Then, by using simple bi-linear shape functions, on each annulus sector 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 whose value is given by a linear combination of the mean-values. Finally, these mean values are used to obtain the numerical approximation for 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 depends on each of the forces discussed above (repulsion, cell-cell adhesion and cell-ECM adhesion). To this end, the spatial dynamics of a cell 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 is assumed to be proportional to the velocity of the cell and so the dynamics is given by [16–18].where η is the damping coefficient and , and are defined in Eqs 3–5, respectively. To solve the spatial dynamics Eq. 6 of the cell , we use the classical forward Euler method. We illustrate the overall travelling direction of a cell in Figure 2.

FIGURE 2

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 (, and defined in Eqs 3–5, respectively) on the macro-scale. For this, let us use the average of the collective contribution of all cell generated forces within a neighbourhood of the macro-node x. This leads to the following continuous cancer cell flux:where again is the Lebesgue measure in , is the ball used also in Eq. 1, 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 . 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 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 (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 and spatial bias . Considering a micro-scale distribution of the micro-fibres on a cell-scale micro-domain of appropriate micro-scale size , we observe that both of these two characteristics of the fibres ECM can be captured through a vector field representation 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 at micro-scale (on ), 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 bywhere is the Bochner-mean-value [49] of the vector-valued function with respect to the measure 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 on , this being given by

In this context, the second characteristic of the fibres ECM, i.e., its amount is precisely given by the Euclidean norm of the vector field , namelywhich in fact represents the mean-value of the micro-fibre distribution on the micro-domain . Since both of these macro-scale fibre ECM representations are derived by using the micro-scale distribution of the micro-fibres , 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 that was defined in Eq. 7. Then, such flux that acts uniformly on naturally gets regulated in a weighted manner by the orientation of the fibres , resulting in a micro-fibres rearrangement vectorwhere the weight is the appropriate mass fraction of the cancer cells and fibres ECM given by

Ultimately, the new position of each micro-fibre is calculated by using a reallocation vector which accounts not only for the rearrangement vector in Eq. 10 but also upon the degree of alignment between and the barycentral position vector as well as it utilises the level of micro-fibres at position z. Hence, this reallocation vector is given bywhere is the maximum possible level of fibres at any micro-point z, is the fibres saturation level and is the characteristic function of the micro-fibres support. Therefore, by using Eq. 11 the new position is given byand the amount of fibres that is reallocated from z to is determined by a movement probabilitywhich tracks the free space available at the new position . Thus, the level of fibres that gets moved to is given by and the rest of the fibres, i.e., remain at the current position z. In Figure 1B we illustrate the reallocation of a micro-fibre using a typical example of the rearrangement and reallocation vectors defined in Eqs. 10, 11, respectively.

Finally, since this rearrangement process that redistributes the fibres ECM micro constituents is initiated by the macro-scale cancer cell flux , the macro and micro-scales are connected via a top-down link which is also illustrated in Figure 1A.

For a diagrammatic description of the full hybrid model, with the exact stages where the agent-based model of the cells connects with the continuous model of the ECM, please see Figure 15 in Appendix B.

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 (that will be defined for each simulations) of the computational domain , 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–23]: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 (as well as throughout the rest of this Section) we present the non-fibre ECM initial conditions only on the domain of . Moreover, for both Figures 3A,B (as well as for all other figures in this Section), the units of the x and y axes are in ‘mm’. Finally to highlight the phase of each individual cell, we use five colours. Specifically, we colour the cells that are in the G1 phase white, those in the S phase yellow, those in the G2 phase orange, those in M phase red, and the ones that are both growth-arrested (due to contact inhibition) and in the M phase are coloured black in all simulations.

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, or , on tumour shape and invasion pattern. We also investigate the impact of different fibres to non-fibres 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., fibres), for two different sensing radii: in Figure 4, and in Figure 5.

FIGURE 4

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 (, ), and investigate again the effect of changes in fibre ECM density (i.e., 10–30–% fibres), for two different sensing radii: in Figure 6, and 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 [sub-panels (A)-(C)] or by aligned fibres [sub-panels (A’)-(C’)].

FIGURE 6

FIGURE 7

First, we see that cell-ECM interactions can lead to irregular-shaped tumour aggregations as well as elongated tumour aggregations. The random fibres are associated with stationary cell aggregations [sub-panels (A)-(C)], while the aligned fibres are associated with moving cell aggregations [sub-panels (A’)–(C’)].

The increase in sensing region (from in Figure 6 to 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 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 (), which impacts the overall shape of tumour aggregations. Let us now decrease the cell-cell adhesion strength and consider (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 which, as before, tends to favorise more round-shaped tumour aggregations.

FIGURE 8

FIGURE 9

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 (while keeping ). In Figures 10, 11 we see that lowering cell-fibre adhesion leads to a highly-irregular tumour shape, especially for the case of fibres to non-fibres ECM ratios.

FIGURE 10

FIGURE 11

Consider now an even lower cell-ECM adhesion strength, , and at the same time a lower cell repulsion . 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.

FIGURE 12

FIGURE 13

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–C′, namely: , , , , , and fibres and non-fibres ECM ratios. In Figure 14 we present six scenarios corresponding to different regions of random and aligned ECM fibres. In sub-panels (A)–(C) we plot the distribution of cancer cells and the ECM degradation pattern at a single time snapshot [i.e., (A) T = 5283 min, (B) T = 4601 min, (C) T = 9518 min]. In panels (D)–(F) we plot the distribution of cancer cells and the ECM degradation patterns at three different times (T = 9500 min, T = 15000 min, T = 20000 min), to illustrate the time-evolution of the system.

FIGURE 14

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 spatio-temporal 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 , cell-cell repulsion stiffness ), and parameters associated with cell-ECM interactions (i.e., cell-fibre adhesion coefficient , and cell-non-fibre adhesion coefficient ), 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 ) and Table 2 (for the case ). By comparing the 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 ), can be further applied to understand various other fibrotic diseases characterised by ECM remodelling [57].

References

• 1.

  WalkerCMojaresEdel Río HernándezA. Role of Extracellular Matrix in Development and Cancer Progression. Int J Mol Sci (2018). 19:3028. 10.3390/ijms19103028

  • CrossRef
  • Google Scholar

• 2.

  XuSXuHWangWLiSLiHLiTet alThe Role of Collagen in Cancer: from Bench to Bedside. J Transl Med (2019). 17:309. 10.1186/s12967-019-2058-1

  • CrossRef
  • Google Scholar

• 3.

  PadhiANainAS. ECM in Differentiation: A Review of Matrix Structure, Composition and Mechanical Properties. Ann Biomed Eng (2019). 48:1071–89. 10.1007/s10439-019-02337-7

  • CrossRef
  • Google Scholar

• 4.

  LinJShiYMenYWangXYeJZhangC. Mechanical Roles in Formation of Oriented Collagen Fibers. Tissue Eng B: Rev (2020). 26:116–28. 10.1089/ten.teb.2019.0243

  • CrossRef
  • Google Scholar

• 5.

  YasuiTTohnoYArakiT. Determination of Collagen Fiber Orientation in Human Tissue by Use of Polarisation Measurement of Molecular Second-Harmonic-Generation Light. Appl Opt (2004). 43:2861–7. 10.1364/ao.43.002861

  • CrossRef
  • Google Scholar

• 6.

  SudaKAbeKKanedaK. Changes in the Orientation of Collagen Fibers on the Superficial Layer of the Mouse Tibial Bone after Denervation. Scanning Electron Microscopic Observations. Arch Histology Cytol (1999). 62:231–5. 10.1679/aohc.62.231

  • CrossRef
  • Google Scholar

• 7.

  HolmesDFGilpinCJBaldockCZieseUKosterAJKadlerKE. Corneal Collagen Fibril Structure in Three Dimensions: Structural Insights into Fibril Assembly, Mechanical Properties, and Tissue Organization. Proc Natl Acad Sci (2001). 98:7307–12. 10.1073/pnas.111150598

  • CrossRef
  • Google Scholar

• 8.

  WeberKT. Cardiac Interstitium in Health and Disease: the Fibrillar Collagen Network. J Am Coll Cardiol (1989). 13:1637–52. 10.1016/0735-1097(89)90360-4

  • CrossRef
  • Google Scholar

• 9.

  LinJPanSZhengWHuangZ. Polarization-resolved Second-Harmonic Generation Imaging for Liver Fibrosis Assessment without Labeling. Appl Phys Lett (2013). 103:173701. 10.1063/1.4826516

  • CrossRef
  • Google Scholar

• 10.

  ConklinMWEickhoffJCRichingKMPehlkeCAEliceiriKWProvenzanoPPet alAligned Collagen Is a Prognostic Signature for Survival in Human Breast Carcinoma. Am J Pathol (2011). 178:1221–32. 10.1016/j.ajpath.2010.11.076

  • CrossRef
  • Google Scholar

• 11.

  MalandrinoAMakMKammRDMoeendarbaryE. Complex Mechanics of the Heterogeneous Extracellular Matrix in Cancer. Extreme Mech Lett (2018). 21:25–34. 10.1016/j.eml.2018.02.003

  • CrossRef
  • Google Scholar

• 12.

  FriedlPAlexanderS. Cancer Invasion and the Microenvironment: Plasticity and Reciprocity. Cell (2011). 147:992–1009. 10.1016/j.cell.2011.11.016

  • CrossRef
  • Google Scholar

• 13.

  MaXSchickelMEStevensonMDSarang-SieminskiALGoochKJGhadialiSNet alFibers in the Extracellular Matrix Enable Long-Range Stress Transmission between Cells. Biophysical J (2013). 104:1410–8. 10.1016/j.bpj.2013.02.017

  • CrossRef
  • Google Scholar

• 14.

  SchallerGMeyer-HermannM. Continuum versus Discrete Model: a Comparison for Multicellular Tumour Spheroids. Phil Trans R Soc A (2006). 364:1443–64. 10.1098/rsta.2006.1780

  • CrossRef
  • Google Scholar

• 15.

  ChamseddineIMRejniakKA. Hybrid Modeling Frameworks of Tumor Development and Treatment. Wires Syst Biol Med (2020). 12:e1461. 10.1002/wsbm.1461

  • CrossRef
  • Google Scholar

• 16.

  KimMReedDRejniakKA. The Formation of Tight Tumor Clusters Affects the Efficacy of Cell Cycle Inhibitors: A Hybrid Model Study. J Theor Biol (2014). 352:31–50. 10.1016/j.jtbi.2014.02.027

  • CrossRef
  • Google Scholar

• 17.

  KarolakAPoonjaSRejniakKA. Morphophenotypic Classification of Tumor Organoids as an Indicator of Drug Exposure and Penetration Potential. Plos Comput Biol (2019). 15:e1007214. 10.1371/journal.pcbi.1007214

  • CrossRef
  • Google Scholar

• 18.

  BerrouetCDorilasNRejniakKATuncerN. Comparison of Drug Inhibitory Effects (IC₅₀) in Monolayer and Spheroid Cultures. Bull Math Biol (2020). 82:68. 10.1007/s11538-020-00746-7

  • CrossRef
  • Google Scholar

• 19.

  ShuttleworthRTrucuD. Multiscale Modelling of Fibres Dynamics and Cell Adhesion within Moving Boundary Cancer Invasion. Bull Math Biol (2019). 81:2176–219. 10.1007/s11538-019-00598-w

  • CrossRef
  • Google Scholar

• 20.

  ShuttleworthRTrucuD. Cell-scale Degradation of Peritumoural Extracellular Matrix Fibre Network and its Role within Tissue-Scale Cancer Invasion. Bull Math Biol (2020). 82:65. 10.1007/s11538-020-00732-z

  • CrossRef
  • Google Scholar

• 21.

  ShuttleworthRTrucuD. Multiscale Dynamics of a Heterotypic Cancer Cell Population within a Fibrous Extracellular Matrix. J Theor Biol (2020). 486:110040. 10.1016/j.jtbi.2019.110040

  • CrossRef
  • Google Scholar

• 22.

  SuvegesSEftimieRTrucuD. Directionality of Macrophages Movement in Tumour Invasion: A Multiscale Moving-Boundary Approach. Bull Math Biol (2020). 82. 10.1007/s11538-020-00819-7

  • CrossRef
  • Google Scholar

• 23.

  SuvegesSEftimieRTrucuD. Re-polarisation of Macrophages within a Multi-Scale Moving Boundary Tumour Invasion Model. arXiv:2103.03384 (2021). 1–50.

  • Google Scholar

• 24.

  AlbertsBJohnsonALewisJRaffMRobertsKWalterP. Molecular Biology of the Cell. New York: Garland Science (2002).

  • Google Scholar

• 25.

  TysonJJNovakB. Temporal Organization of the Cell Cycle. Curr Biol (2008). 18:R759–R768. 10.1016/j.cub.2008.07.001

  • CrossRef
  • Google Scholar

• 26.

  Pérez-VelázquezJRejniakKA. Drug-Induced Resistance in Micrometastases: Analysis of Spatio-Temporal Cell Lineages. Front Physiol (2020). 11:319. 10.3389/fphys.2020.00319

  • CrossRef
  • Google Scholar

• 27.

  CooperG. The Cell : A Molecular Approach. Washington, D.C. Sunderland, Mass: ASM Press Sinauer Associates (2000).

  • Google Scholar

• 28.

  HudaSWeigelinBWolfKTretiakovKVPolevKWilkGet alLévy-like Movement Patterns of Metastatic Cancer Cells Revealed in Microfabricated Systems and Implicated In Vivo. Nat Commun (2018). 9:4539. 10.1038/s41467-018-06563-w

  • CrossRef
  • Google Scholar

• 29.

  PetrieRJDoyleADYamadaKM. Random versus Directionally Persistent Cell Migration. Nat Rev Mol Cel Biol (2009). 10:538–49. 10.1038/nrm2729

  • CrossRef
  • Google Scholar

• 30.

  WeigerMCVedhamVStueltenCHShouKHerreraMSatoMet alReal-time Motion Analysis Reveals Cell Directionality as an Indicator of Breast Cancer Progression. PLOS ONE (2013). 8:e58859–12. 10.1371/journal.pone.0058859

  • CrossRef
  • Google Scholar

• 31.

  WuP-HGiriASunSXWirtzD. Three-dimensional Cell Migration Does Not Follow a Random Walk. Proc Natl Acad Sci (2014). 111:3949–54. 10.1073/pnas.1318967111

  • CrossRef
  • Google Scholar

• 32.

  GhoshSSalotSSenguptaSNavalkarAGhoshDJacobRet alp53 Amyloid Formation Leading to its Loss of Function: Implications in Cancer Pathogenesis. Cell Death Differ (2017). 24:1784–98. 10.1038/cdd.2017.105

  • CrossRef
  • Google Scholar

• 33.

  GrasSL. Surface- and Solution-Based Assembly of Amyloid Fibrils for Biomedical and Nanotechnology Applications. In: KoopmansRJ, editor. Engineering Aspects of Self-Organizing Materials. Advances in Chemical Engineering, 35. Academic Press (2009). p. 161–209. 10.1016/S0065-2377(08)00206-8

  • CrossRef
  • Google Scholar

• 34.

  GrasSLTicklerAKSquiresAMDevlinGLHortonMADobsonCMet alFunctionalised Amyloid Fibrils for Roles in Cell Adhesion. Biomaterials (2008). 29:1553–62. 10.1016/j.biomaterials.2007.11.028

  • CrossRef
  • Google Scholar

• 35.

  JacobRSGeorgeESinghPKSalotSAnoopAJhaNNet alCell Adhesion on Amyloid Fibrils Lacking Integrin Recognition Motif. J Biol Chem (2016). 291:5278–98. 10.1074/jbc.m115.678177

  • CrossRef
  • Google Scholar

• 36.

  WolfKAlexanderSSchachtVCoussensLMvon AndrianUHvan RheenenJet alCollagen-based Cell Migration Models In Vitro and In Vivo. Semin Cel Develop Biol (2009). 20:931–41. 10.1016/j.semcdb.2009.08.005

  • CrossRef
  • Google Scholar

• 37.

  WolfKFriedlP. Extracellular Matrix Determinants of Proteolytic and Non-proteolytic Cell Migration. Trends Cel Biol (2011). 21:736–44. 10.1016/j.tcb.2011.09.006

  • CrossRef
  • Google Scholar

• 38.

  ArmstrongNJPainterKJSherrattJA. A Continuum Approach to Modelling Cell-Cell Adhesion. J Theor Biol (2006). 243:98–113. 10.1016/j.jtbi.2006.05.030

  • CrossRef
  • Google Scholar

• 39.

  GerischAChaplainMAJ. Mathematical Modelling of Cancer Cell Invasion of Tissue: Local and Non-local Models and the Effect of Adhesion. J Theor Biol (2008). 250:684–704. 10.1016/j.jtbi.2007.10.026

  • CrossRef
  • Google Scholar

• 40.

  DomschkePTrucuDGerischAA. J. ChaplainM. Mathematical Modelling of Cancer Invasion: Implications of Cell Adhesion Variability for Tumour Infiltrative Growth Patterns. J Theor Biol (2014). 361:41–60. 10.1016/j.jtbi.2014.07.010

  • CrossRef
  • Google Scholar

• 41.

  DomschkePTrucuDGerischAChaplainMAJ (2017). Structured Models of Cell Migration Incorporating Molecular Binding Processes. J. Math. Biol.75, 1517–1561. 10.1007/s00285-017-1120-y

  • CrossRef
  • Google Scholar

• 42.

  PengLTrucuDLinPThompsonAChaplainMAJ. A Multiscale Mathematical Model of Tumour Invasive Growth. Bull Math Biol (2017). 79:389–429. 10.1007/s11538-016-0237-2

  • CrossRef
  • Google Scholar

• 43.

  NagaseHVisseRMurphyG. Structure and Function of Matrix Metalloproteinases and TIMPs. Cardiovasc Res (2006). 69:562–73. 10.1016/j.cardiores.2005.12.002

  • CrossRef
  • Google Scholar

• 44.

  VisseRNagaseH. Matrix Metalloproteinases and Tissue Inhibitors of Metalloproteinases. Circ Res (2003). 92:827–39. 10.1161/01.RES.0000070112.80711.3D

  • CrossRef
  • Google Scholar

• 45.

  HanahanDWeinbergRA. Hallmarks of Cancer: The Next Generation. Cell (2011). 144:646–74. 10.1016/j.cell.2011.02.013

  • CrossRef
  • Google Scholar

• 46.

  WeinbergRA. The Biology of Cancer. New York: Garland Science (2006).

  • Google Scholar

• 47.

  StixBKähneTSlettenKRaynesJRoessnerARöckenC. Proteolysis of AA Amyloid Fibril Proteins by Matrix Metalloproteinases-1, -2, and -3. Am J Pathol (2001). 159:561–70. 10.1016/s0002-9440(10)61727-0

  • CrossRef
  • Google Scholar

• 48.

  LiaoM-CVan NostrandWE. Degradation of Soluble and Fibrillar Amyloid β-Protein by Matrix Metalloproteinase (MT1-MMP) In Vitro. Biochemistry (2010). 49:1127–36. 10.1021/bi901994d

  • CrossRef
  • Google Scholar

• 49.

  YosidaKFunctional Analysis. 6th ed.Springer-Verlag (1980).

  • Google Scholar

• 50.

  MurphyDACourtneidgeSA. The 'ins' and 'outs' of Podosomes and Invadopodia: Characteristics, Formation and Function. Nat Rev Mol Cel Biol. (2012). 12:413–26. 10.1038/nrm3141

  • CrossRef
  • Google Scholar

• 51.

  KlemkeRL. Trespassing Cancer Cells: 'fingerprinting' Invasive Protrusions Reveals Metastatic Culprits. Curr Opin Cel Biol (2012). 24:662–9. 10.1016/j.ceb.2012.08.005

  • CrossRef
  • Google Scholar

• 52.

  ProvenzanoPPInmanDREliceiriKWTrierSMKeelyPJ. Contact Guidance Mediated Three-Dimensional Cell Migration Is Regulated by Rho/ROCK-dependent Matrix Reorganization. Biophysical J (2008). 95:5374–84. 10.1529/biophysj.108.133116

  • CrossRef
  • Google Scholar

• 53.

  HanWChenSYuanWFanQTianJWangXet alOriented Collagen Fibers Direct Tumor Cell Intravasation. Proc Natl Acad Sci USA (2016). 113:11208–13. 10.1073/pnas.1610347113

  • CrossRef
  • Google Scholar

• 54.

  FriedlP. Rethinking Research into Metastasis. eLIFE (2019). 8:e53511. 10.7554/eLife.53511

  • CrossRef
  • Google Scholar

• 55.

  UchugonovaAZhaoMWeinigelMZhangYBouvetMHoffmanRMet alMultiphoton Tomography Visualizes Collagen Fibers in the Tumor Microenvironment that Maintain Cancer-Cell anchorage and Shape. J Cel Biochem. (2013). 114:99–102. 10.1002/jcb.24305

  • CrossRef
  • Google Scholar

• 56.

  SunYYouSTuHSpillmanDRChaneyEJMarjanovicMet alIntraoperative Visualization of the Tumor Microenvironment and Quantification of Extracellular Vesicles by Label-free Nonlinear Imaging. Sci Adv (2018). 4:eaau5603. 10.1126/sciadv.aau5603

  • CrossRef
  • Google Scholar

• 57.

  CoxTRErlerJT. Remodeling and Homeostasis of the Extracellular Matrix: Implications for Fibrotic Diseases and Cancer. Dis Models Mech (2011). 4:165–78. 10.1242/dmm.004077

  • CrossRef
  • Google Scholar

• 58.

  ShashniBAriyasuSTakedaRSuzukiTShiinaSAkimotoKet alSize-based Differentiation of Cancer and normal Cells by a Particle Size Analyzer Assisted by a Cell-Recognition Pc Software. Biol Pharm Bull (2018). 41:487–503. 10.1248/bpb.b17-00776

  • CrossRef
  • Google Scholar

• 59.

  NCI. NCI-60 Human Tumor Cell Lines Screen [internet]. on line (2015). Available from: https://dtp.cancer.gov/discovery˙_development/̵nci- ̵60/.

  • Google Scholar

• 60.

  BaumgartnerWDrenckhahnD. Transmembrane Cooperative Linkage in Cellular Adhesion. Eur J Cel Biol (2002). 81:161–8. 10.1078/0171-9335-00233

  • CrossRef
  • Google Scholar

• 61.

  RejniakKADillonRH. A Single Cell-Based Model of the Ductal Tumour Microarchitecture. Comput Math Methods Med (2007). 8:51–69. 10.1080/17486700701303143

  • CrossRef
  • Google Scholar

• 62.

  ShimolinaLEIzquierdoMALópez-DuarteIBullJAShirmanovaMVKlapshinaLGet alImaging Tumor Microscopic Viscosity In Vivo Using Molecular Rotors. Sci Rep (2017). 7:41097. 10.1038/srep41097

  • CrossRef
  • Google Scholar