Multiscale dynamics of a heterotypic cancer cell population within a fibrous extracellular matrix

TL;DR

This paper develops a multiscale model of heterotypic cancer invasion within a fibrous extracellular matrix, integrating tissue-scale dynamics with cell-scale molecular processes to better understand tumor spread.

Contribution

It introduces a novel multiscale framework that incorporates fiber network dynamics and enzyme-driven invasion in heterotypic tumors, extending previous models.

Findings

Fiber network influences invasion patterns

Heterotypic interactions affect invasion speed

Model provides insights into tumor-matrix interactions

Abstract

Local cancer cell invasion is a complex process involving many cellular and tissue interactions and is an important prerequisite for metastatic spread, the main cause of cancer related deaths. Occurring over many different temporal and spatial scales, the first stage of local invasion is the secretion of matrix-degrading enzymes (MDEs) and the resulting degradation of the extra-cellular matrix (ECM). This process creates space in which the cells can invade and thus enlarge the tumour. As a tumour increases in malignancy, the cancer cells adopt the ability to mutate into secondary cell subpopulations giving rise to a heterogeneous tumour. This new cell subpopulation often carries higher invasive qualities and permits a quicker spread of the tumour. Building upon the recent multiscale modelling framework for cancer invasion within a fibrous ECM introduced in Shuttleworth and Trucu (2019), in this paper we consider the process of local invasion by a heterotypic tumour consisting of two cancer cell populations mixed with a two-phase ECM. To that end, we address the double feedback link between the tissue-scale cancer dynamics and the cell-scale molecular processes through the development of a two-part modelling framework that crucially incorporates the multiscale dynamic redistribution of oriented fibres occurring within a two-phase extra-cellular matrix and combines this with the multiscale leading edge dynamics exploring key matrix-degrading enzymes molecular processes along the tumour interface that drive the movement of the cancer boundary. The modelling framework will be accompanied by computational results that explore the effects of the underlying fibre network on the overall pattern of cancer invasion.

Tables

Table 1. Table 1 : The parameters in Σ Σ \Sigma

Parameter	Value	Description	Reference
$D_1$	$3.5\times {10}^{-4}$	diffusion coeff. for cell population $c_1$	Domschke et al. [2014]
$D_2$	$7\times {10}^{-4}$	diffusion coeff. for cell population $c_2$	Domschke et al. [2014]
$D_m$	${10}^{-3}$	diffusion coeff. for MDEs	Estimated
${\mu }_1$	$0.25$	proliferation coeff. for cell population $c_1$	Domschke et al. [2014]
${\mu }_2$	$0.25$	proliferation coeff. for cell population $c_2$	Domschke et al. [2014]
${\gamma }_1$	2	non-fibrous ECM degradation coeff.	Shuttleworth and Trucu [2019]
${\gamma }_2$	1.5	macroscopic fibre degradation coeff.	Peng et al. [2016]
$\omega$	0-0.02	non-fibrous ECM remodelling coeff.	Domschke et al. [2014]
${\alpha }_1$	$1$	MDE secretion rate of $c_1$	Estimated
${\alpha }_2$	$1.5$	MDE secretion rate of $c_2$	Estimated
${\text{S}}_{max}$		cell-cell adhesion coeff. matrix	Domschke et al. [2014]
${\text{S}}_{cF}$		cell-fibre adhesion coeff.	Estimated
${\text{S}}_{cl}$		cell-matrix adhesion coeff.	Estimated
$R$	0.15	sensing radius	Shuttleworth and Trucu [2019]
$r$	0.0016	width of micro-fibres	Shuttleworth and Trucu [2019]
$f_{\text{max}}$	0.6360	max. micro-density of fibres	Shuttleworth and Trucu [2019]
$p$	0.15-0.2	percentage of non-fibrous ECM	Estimated
$h$	0.03125	macro-scale spatial discretisation size	Trucu et al. [2013]
$ϵ$	0.0625	size of micro-domain $ϵY$	Trucu et al. [2013]
$\delta$	0.03125	size of micro-domain $\delta Y$	Shuttleworth and Trucu [2019]

Equations

$${\theta }_{{}_{f,\delta Y(x)}}(x,t)=\frac{\underset{\delta Y(x)}{\int }f(z,t)(z-x)dz}{\underset{\delta Y(x)}{\int }f(z,t)dz}.$$

$${\theta }_{{}_f}(x,t)=\frac{1}{\lambda (\delta Y(x))}{\int }_{\delta Y(x)}f(z,t)dz\cdot \frac{{\theta }_{{}_{f,\delta Y(x)}}(x,t)}{||{\theta }_{{}_{f,\delta Y(x)}}(x,t)||}$$

$$F(x,t):=||{\theta }_{{}_f}(x,t)|{|}_2.$$

$$\text{u}(x,t)=[\text{c}(x,t{)}^T,F(x,t),l(x,t){]}^T,$$

$$\rho (\text{u}(x,t))={\vartheta }_v(l(x,t)+F(x,t))+{\vartheta }_c\bar{c}(x,t),$$

$$\bar{c}(x,t):=c_1(x,t)+c_2(x,t).$$

$$\frac{\partial c_1}{\partial t}$$

$$\frac{\partial c_2}{\partial t}$$

$$\begin{split}\mathcal{A}_{n}(x,t,\textbf{u}(\cdot,t),\theta_{{}_{f}}(\cdot,t))=\frac{1}{R}\int_{B(0,R)}\mathcal{K}(\parallel\!y\!\parallel_{2})&\big{(}n(y)(\textbf{S}_{{}_{cc}}\textbf{c}(x+y,t)+\textbf{S}_{{}_{cl}}l(x+y,t))\\ &+\hat{n}(y)\ \textbf{S}_{{}_{cF}}F(x+y,t)\big{)}(1-\rho(\textbf{u}))^{+}\end{split}$$

$$\mathcal{K}(r)=\frac{2\pi R^2}{3}\left(1-\frac{r}{2R}\right),$$

$$\mathbf{n}(y):=\left\{\begin{array}[]{l}y/||y||_{2}\quad\text{if}\quad y\in B(0,R)\setminus\{(0,0)\},\\ \left(0,0\right)\quad\textrm{otherwise}\end{array}\right.$$

$${\text{S}}_{{}_{cc}}(x,t):={\text{S}}_{{}_{max}}{e}^{\left(1-\frac{1}{1-(1-l(x,t){)}^2}\right)}.$$

$$\hat{n}(y):=\left\{\begin{array}[]{l}\frac{y+\theta_{f}(x+y)}{||y+\theta_{f}(x+y)||_{2}}\quad\text{if}\ (y+\theta_{f}(x+y))\neq(0,0)\\[5.69046pt] (0,0)\in\mathbb{R}^{2}\quad\quad\text{otherwise}.\end{array}\right.$$

$$M_{c}(x,t):=\left\{\begin{array}[]{l}\frac{\text{exp}\left(\frac{-1}{\kappa^{2}-(1-v(x,t))^{2}}\right)}{\text{exp}\left(\frac{1}{\kappa^{2}}\right)}\cdot H(t-t_{m})\quad\text{if}\ 1-\kappa<v(x,t)<1,\\ 0\quad\quad\quad\textrm{otherwise},\end{array}\right.$$

$$\frac{dF}{dt}$$

$$\frac{dl}{dt}$$

$$\mathcal{F}(x,t):={\mathcal{F}}_1(x,t)+{\mathcal{F}}_2(x,t)$$

$${\mathcal{F}}_1(x,t):=D_1\nabla c_1(x,t)-c_1(x,t){\mathcal{A}}_1(x,t,\text{u}(\cdot ,t),{\theta }_{{}_f}(\cdot ,t)),$$

$${\mathcal{F}}_2(x,t):=D_2\nabla c_2(x,t)-c_2(x,t){\mathcal{A}}_2(x,t,\text{u}(\cdot ,t),{\theta }_{{}_f}(\cdot ,t)).$$

$$\omega (x,t)=\frac{\bar{c}(x,t)}{\bar{c}(x,t)+F(x,t)}.$$

$$r(\delta Y(x),t):=\omega (x,t)\mathcal{F}(x,t)+(1-\omega (x,t)){\theta }_{{}_f}(x,t).$$

$$z^{\ast }:=z+{\nu }_{{}_{\delta Y(x)}}(z,t)$$

$${\nu }_{{}_{\delta Y(x)}}(z,t)=\left(x_{\text{dir}}(z)+r(\delta Y(x),t)\right)\cdot \frac{f(z,t)(f_{\text{max}}-f(z,t))}{f^{\ast }+||r(\delta Y(x))-x_{\text{dir}}(z)|{|}_2}\cdot {\chi }_{{}_{\{f(\cdot ,t)>0\}}}(z)$$

$$||r(\delta Y(x))-x_{\text{dir}}(z)|{|}_2.$$

$$f^{\ast }=\frac{f(z,t)}{f_{\text{max}}}.$$

$$p_{move}:=\max\big{(}0,\frac{f_{\text{max}}-f(z^{*},t)}{f_{\text{max}}}\big{)}$$

$$\displaystyle\begin{split}1.\quad&g_{\epsilon Y}(y,\tau)=\frac{\int\limits_{\textbf{B}(z,\gamma)\cap\Omega(t_{0})}\alpha_{1}c_{1}(x,t_{0}+\tau)+\alpha_{2}c_{2}(x,t_{0}+\tau)\ dx}{\lambda(\textbf{B}(y,\gamma)\cap\Omega(t_{0}))},\quad y\in\epsilon Y\cap\Omega(t_{0}),\\[19.91684pt] 2.\quad&g_{\epsilon Y}(y,\tau)=0,\quad y\in\epsilon Y\setminus\big{(}\Omega(t_{0})+\{y\in Y|\ ||z||_{2}<\gamma\}),\end{split}$$

$$\frac{\partial m}{\partial \tau }=D_m\Delta m+g_{ϵY}(z,\tau ),z\in ϵY,\tau \in [0,\Delta t].$$

$$c_1(x,0)$$

$$c_2(x,0)$$

$$l(x,0)=min\{0.5,1-c_1(x,0)\},$$

$$l(x,0)=\text{min}\left\{h(x_1,x_2),1-c(x,0)\right\},$$

Full text

Building upon the recent multiscale modelling framework for cancer invasion within a fibrous ECM introduced in Shuttleworth and Trucu [2019], in this paper we consider the process of local invasion by a heterotypic tumour consisting of two cancer cell populations mixed with a two-phase ECM. To that end, we address the double feedback link between the tissue-scale cancer dynamics and the cell-scale molecular processes through the development of a two-part modelling framework that crucially incorporates the multiscale dynamic redistribution of oriented fibres occurring within a two-phase extra-cellular matrix and combines this with the multiscale leading edge dynamics exploring key matrix-degrading enzymes molecular processes along the tumour interface that drive the movement of the cancer boundary. The modelling framework will be accompanied by computational results that explore the effects of the underlying fibre network on the overall pattern of cancer invasion.

††journal: Journal of Theoretical Biology

1 Introduction

Regarded as one of the hallmarks of cancer [Hanahan and Weinberg, 2000], cancer cell invasion of tissue is a highly complex process that occurs over many spatial and temporal levels. The invasion of surrounding tissues is a key process for tumour progression and plays a crucial role during the formation of metastases. One of the first steps of tumour invasion is the production and secretion of matrix-degrading enzymes (MDEs) by the cancer cells. These enzymes degrade the surrounding extracellular matrix (ECM) enabling this way an immediate spatial progression of the cancer cells population into neighbouring tissues via both random motility and cell-adhesion mediated migration as well as enhanced proliferation, and, eventually, via a cascade of invasion events (including angiogenesis) leading to the spread at distant sites in the body (metastasis). In this context, as the prognosis for patients is still poor, with limited treatment options (such as chemotherapy and radiation) and metastases being currently attributed to $90\%$ of cancer related deaths [Sporn, 1996], a deeper understanding of the processes that occur during local invasion is key for the improvement and future development of treatment strategies.

The local invasion of cancer is, in the first instance, stimulated by the secretion of MDEs. There are several classes of MDEs, such as matrix metalloproteinases (MMPs) and the urokinase-type plasminogen activator (uPA), produced by cancer cells and these enzymes degrade and reshape the structure of the ECM. Through full or partial degradation of the surrounding matrix the cells create free space in which they can invade, further advancing tumour progression. MMPs are substrate specific and can be either membrane bound (MT-MMPs), or can move freely around the cells, each working on different components of the ECM ensuring a significant breakdown of the matrix [Brinckerhoff and Matrisian, 2002, Parsons et al., 1997].

The invasive abilities of cancer can be strengthened in many ways, namely, enhanced proliferation, increased migrative capability and adaptive cellular adhesion properties. Cell-cell and cell-matrix adhesion are both key players during cancer invasion and play important roles in tumour progression [Behrens et al., 1989, Kawanishi et al., 1995, Todd et al., 2016, Berrier and Yamada, 2007]. Any perturbations to either factor will contribute to a change in tumour morphology and the direction of migration of the tumour. A decrease in cell-cell adhesion allows the cells to detach from the primary tumour and invade further into the matrix, and coupling this with cell-matrix adhesion increases leads to a notable escalation of tumours’ invasive capabilities [Cavallaro and Christofori, 2001].

In healthy cells, cell-cell adhesion is mediated by a family of specific molecules on the cell surface known as cell adhesion molecules (CAMs). Adhesion is dependent on the cell-cell signalling pathways that are formed due to the interactions between the distribution of calcium-sensing receptors and $\text{Ca}^{2+}$ ions within the ECM. Essential for cellular adhesion is the family of transmembrane glycoproteins cadherins. These are calcium dependent adhesion molecules that interact with intra-cellular proteins, known as catenins. In particular, the subfamily, E-cadherins, are responsible for binding with these catenins, typically $\beta$-catenin, forming an E-cadherin/catenin complex. Any alteration to the function of $\beta$-catenin will result in the loss of ability of the E-cadherin to initiate cell-cell adhesion [Wijnhoven et al., 2000]. The direct correlation between this calcium-based cell signalling mechanism and the regulation of E-cadherin and $\beta$-catenin was first discovered in colon carcinoma [Bhagavathula et al., 2007].

On the other hand, cell-matrix adhesion is regulated by the subfamily of calcium independent CAMs, known as integrins, that enable the cells to bind to different components of the ECM. Integrins link the cytoskeleton inside a cell to the ECM outside and this is aided by their ability to attach to a wide variety of ligands [Humphries et al., 2006]. Additionally, integrins can bind to actin proteins within the cytoskeleton, aiding in cell migration by creating a leading and trailing edge to the cell, resulting in persistent migration in one direction [Delon and Brown, 2007, Moissoglu and Schwartz, 2006].

All of these properties of invasive tumour cells have one thing in common, they all require or interact with the ECM. The extracellular matrix is comprised from a variety of proteins including collagen and elastin, glycoproteins such as fibronectin and laminins as well as a large array of other molecules, MDEs, fibroblasts, etc. It is however the complex network of ECM fibres, such as collagen and fibronectin fibres, that provides not only support for cells, but it acts as a platform through which cells can communicate. This feature is of particular use to cancer cells, which use the matrix as a means of invasion. Through a combination of degradation and cell-matrix adhesion, the ECM is manipulated and exploited to further advance their spatial progression. A main fibrous component that provides ECM with a scaffolding structure and integrity is collagen, which is the most abundant protein in the human body, with collagen type I, II and III making up around $90\%$ of the overall collagen present [Lodish et al., 2000].

Moving away from the scaffolding structure of the ECM, there are other fibrous ECM components that provide functional qualities within the matrix. One such component is the glycoprotein fibronectin. Fibronectin contributes to cell migration, growth and proliferation, ensuring the normal functionality of healthy cells. It also plays a crucial role in cell adhesion, having the ability to anchor cells to collagen and other components of the ECM.

Biological experiments have revealed that increased stromal collagen density promotes tumour formation and results in tumours exhibiting a more invasive phenotype [Provenzano et al., 2008]. In addition, it has been demonstrated that local invasion is further accelerated by collagen reorganisation [Provenzano et al., 2006], and this behaviour is significantly increased in regions of high collagen density. The realignment of fibronectin fibrils has also been associated with increased local tumour invasion [Erdogan et al., 2017] enabling a smooth invasion of the cancer cells. On the other hand, collagen type I has been shown to down-regulate E-cadherin gene expression in pancreatic cancer cell lines which leads to a reduction in cell-cell adhesion and increased proliferation and cell migration [Menke et al., 2001]. This motivated us to explore the migratory behaviour of the overall cancer cell population by assuming in this work the presence of a secondary cancer cell sub-population exhibiting a decrease in cell-cell adhesion and an increase in migratory activity arising in places of high matrix density.

The past two decades or so have witnessed a vast of interest in the modelling of cancer invasion, see, for example, [Perumpanani et al., 1998, Anderson et al., 2000, Chaplain et al., 2006, Hillen, 2006, Chauviere et al., 2007, Szymańska et al., 2009, Andasari et al., 2011, Ramis-Conde et al., 2008a, Chaplain et al., 2011, Scianna and Preziosi, 2012]. With biological experiments advancing, there is an increasing need for more extensive modelling of the processes involved in cancer invasion. Biological and mathematical models of both in vivo and in vitro experiments have given us a deeper insight into many processes involved during tumour invasion. Great focus has been placed on modelling the effects of cell-cell and cell-matrix adhesion [Painter et al., 2010, Armstrong et al., 2006, Anderson, 2005, Turner and Sherratt, 2002, Gerisch and Chaplain, 2008, Domschke et al., 2014, Bitsouni et al., 2017]. On the other hand, recent works such as [Chauviere et al., 2007, Painter, 2008, Hillen et al., 2010, Schluter et al., 2012, Hillen et al., 2013, Engwer et al., 2015] have highlighted the vital importance that the composition of the ECM has on the overall invasion of cancer. Finally, the multi-scale nature of cancer invasion has received special attention over the past decade [Ramis-Conde et al., 2008b, Trucu et al., 2013, Peng et al., 2016, Shuttleworth and Trucu, 2018], with significant advancements towards two-scale approaches appropriately linking the spatio-temporal dynamics occurring at different scales being proposed in Trucu et al. [2013], Shuttleworth and Trucu [2019].

In this paper we extend the model developed in [Shuttleworth and Trucu, 2019] to investigate the invasiveness of a heterogeneous tumour in a multi-phase fibrous tissue environment by introducing a secondary cancer cell subpopulation. This will examine the invasive abilities of a malignant tumour with two cancer cell sub-populations under the presence of a two-component ECM consisting of both an oriented fibre phase distribution (such as collagen and fibronectin), and a soluble phase, that accounts for all other ECM non-fibrous components (such as Ca2+ ions, laminin, and other soluble matrix constituents). To that end, in the context of multiscale moving boundary approach introduced in Trucu et al. [2013] (capturing the key influence of the underlying tumour invasive edge two-scale MDEs proteolytic activity), we will consider here the dynamics of the two cancer cell sub-populations progressing within the surrounding fibrous ECM and explore its multiscale evolution and interaction with the ECM fibre phase in the presence of both homogeneous and heterogeneous non-fibre ECM soluble phases.

The paper is structured as follows. In Section 2 we will detail the distribution of oriented macroscopic (tissue-scale) ECM fibres vector field induced by the microscale (cell-scale) mass distribution of micro-fibres alongside the multiscale interaction that arises between the fibre ECM phase and the two cancer cell populations (both in terms of dynamic micro-scale fibres rearrangements by the cancer cells and the fibre impact within the macro-scale cancer evolution). In Section 3 we will present the numerical approaches used and initial conditions for the computations and in Section 4 we present the simulation results. Finally, we will conclude with Section 5 where we discuss the implications of our results as well as future work.

2 The mathematical model

Here we will build upon the key aspects of the two-part multiscale model introduced and developed in [Shuttleworth and Trucu, 2019] that utilises the two-scale moving boundary framework first proposed in [Trucu et al., 2013] and introduces a two-component ECM in which the fibre phase plays a central role in the two-scale dynamic redistribution of microscopic fibres. As we will briefly detail below, this complex dynamics will be captured by two interconnected multiscale systems that share the same macro-scale cancer dynamics at the tissue-scale, whilst having their own distinct micro-scale dynamics occurring at cell-scale that are linked to the macro-dynamics through two double feedback loops, as illustrated in Fig. 3.

During the invasion process a tumour can become increasingly malignant, whereby the primary cancer cell population acquires the ability to mutate, giving rise to a secondary subpopulation of cancer cells that exhibits more aggressive invasion qualities, including: faster random motility, increased proliferation, and changes in cell-cell and cell-matrix adhesion properties enabling an acceleration of local cancer invasion. Thus, the presence of this secondary cancer cell subpopulation has implications for cancer dynamics at both macro- and micro- scales, and to address all these, we will devote the following sections to develop and adapt the modelling approach introduced initially in [Shuttleworth and Trucu, 2019] to the new context created by the two cancer cell sub-populations.

Using here the same terminology as in [Shuttleworth and Trucu, 2019], let us denote the support of the invading tumour region by $\Omega(t)$, and assume this evolves in the maximal reference tissue cube $Y\in\mathbb{R}^{N}$ with $N=2,3$, centred at the origin of the space. At any spatio-temporal node $(x,t)\in\Omega(t)\times[0,T]$ we consider the tumour to be a mixture of cancer cells $c_{n}(x,t)$, $n=1,2$, with their combined vector denoted $\textbf{c}(x,t)=[c_{1}(x,t),c_{2}(x,t)]$, integrated within a multiphase distribution of ECM, $v(x,t)$, whose components will be defined in the next section.

2.1 Multiscale fibre structure and their dynamic contribution in tumour progression

Adopting the approach and terminology introduced in [Shuttleworth and Trucu, 2019], to capture the dynamics of the two sub-populations of cancer cells in the presence of a two component heterogeneous extra-cellular matrix (consisting of a fibre and non-fibre phase), we proceed as follows. As derived in [Shuttleworth and Trucu, 2019] and briefly outlined below, at any macroscale point $x\in\Omega(t)$, the ECM-fibre phase is represented through a macro-scale vector field $\theta_{{}_{f}}(x,t)$. This vector field captures and represents at tissue-scale not only the amount of fibres distributed at $(x,t)$ but also their naturally arising macroscopic fibres orientation induced by the revolving barycentral orientation $\theta_{{}_{f,\delta Y(x)}}(x,t)$ generated by the microscopic mass distribution of microfibres $f(\cdot,t)$ within the micro-domain $\delta Y(x):=\delta Y+x$. An example of a mass distribution of ECM micro-fibres $f(z,t)$, $z\in\delta Y(x)$, is illustrated in Fig. 1 below and is defined in A.

In brief, while referring the reader to its full derivation presented in Shuttleworth and Trucu [2019], the naturally generated revolving barycentral orientation $\theta_{{}_{f,\delta Y(x)}}(x,t)$ associated with $\delta Y(x)$ is given by the Bochner-mean-value of the position vectors function $\delta Y(x)\ni z\mapsto z-x\in\mathbb{R}^{N}$ with respect to the density measure $f(x,t)\lambda(\cdot)$, where $\lambda(\cdot)$ is the usual Lebesgue measure (see Yosida [1980]), and so this is expressed mathematically as:

[TABLE]

Following on, at any spatio-temporal node $(x,t)$, the macroscopic fibre oriented is defined as

[TABLE]

where $\lambda(\cdot)$ is the usual Lebesgue measure. The macroscopic mean-value fibre representation at any $(x,t)$ is then given by the Euclidean magnitude of $\theta_{{}_{f}}(x,t)$, namely,

[TABLE]

Finally, concerning the non-fibre soluble phase of the ECM, we consider that this include all the other non-fibre components of the ECM, i.e., elastin, laminins, fibroblasts, etc.. Thus, denoting the spatial distribution of the non-fibre ECM phase distribution $l(x,t)$, the total ECM distributed at any spatio-temporal node $(x,t)$ is therefore given by $v(x,t)=l(x,t)+F(x,t)$.

2.2 Macro-scale dynamics

To explore the dynamics of a heterotypic cancer cell population with two sub-populations consisting of a primary and a mutated cell population, within the modelling framework, let’s denote the global macro-scale tumour vector distribution by

[TABLE]

with the tumour’s volume fraction of occupied space being given as

[TABLE]

where $\vartheta_{v}$ and $\vartheta_{c}$ represent the fraction of physical space occupied by the entire ECM (including both the fibre and non-fibre phase) and $\bar{c}(x,t)$ denotes the combined cell populations $c_{1}$ and $c_{2}$, this being given by

[TABLE]

Therefore, adopting the same approach as in Andasari et al. [2011] and Shuttleworth and Trucu [2018], focusing first on the cancer cell population, the dynamics of the two cell sub-populations are much similar. Under the presence of a logistic proliferation law, per unit time, the spatial movement of the primary tumour cells $c_{1}(x,t)$ is described by a combination of local Brownian movement (approximated here by diffusion) and cell-adhesion. They experience a loss of cells through mutation towards a second, more aggressive population $c_{2}(x,t)$. Similar to population $c_{1}(x,t)$, once mutations have begun, per unit time, sub-population $c_{2}(x,t)$ exercises spatial movement through a local Brownian movement and a perturbed cell-adhesion under the presence of logistic proliferation law. Hence, mathematically, the dynamics of these cell populations can be represented as

[TABLE]

where: $D_{n}$ and $\mu_{n}$, $n=1,2$ are the non-negative diffusion and proliferation coefficients of sub-populations $c_{1}(t,x)$ and $c_{2}(t,x)$ respectively, $M_{c}$ describes the conversion from population $c_{1}(t,x)$ to $c_{2}(x,t)$ and finally, the non-local flux $\mathcal{A}_{n}(x,t,\textbf{u}(t,\cdot),\theta_{f}(\cdot,t))$ accounts for the cellular adhesion processes which directly influence the spatial movement of the tumour cell population $c_{n}$, $n=1,2$. Embracing the modelling concept proposed in Shuttleworth and Trucu [2019], the cells will not only interact with other cells, i.e., cell-cell adhesion, but also with the surrounding multi-phase ECM, that in this instance constitutes of both cell-fibre and cell-ECM-non-fibre adhesion. Hence, within a sensing radius $R$ at time $t$, the non-local adhesive flux can be expressed as:

[TABLE]

Whilst we consider the adhesive activities of the cells to become less influential as the distance $r$ from $x$ increases, and account for this through the radial kernel $\mathcal{K}(\cdot)$ defined in Shuttleworth and Trucu [2019] as

[TABLE]

we explore the strength of the adhesion bonds created between the cancer cells distributed at $x$ and the cells or the non-fibre phase of the ECM distributed at $y$ in the direction of the unit normal

[TABLE]

and let us denote the cell-cell and cell-non-fibre-ECM adhesive strengths by $\textbf{S}_{{}_{cc}}$ and $\textbf{S}_{{}_{cl}}$, respectively. Furthermore, proceeding as in [Shuttleworth and Trucu, 2019], we consider $\textbf{S}_{{}_{cl}}$ to be constant whilst the coefficient representing cell-cell adhesion $\textbf{S}_{{}_{cc}}$ is monotonically dependent on the level of extracellular Ca2+ ions enabling strong adhesive bonds between cells [Hofer et al., 2000, Gu et al., 2014]. Hence we assume that cell-cell adhesion is dependent on the density of the underlying non-fibre ECM phase, ranging from [math] to a Ca2+-saturation level denoted $S_{{}_{max}}$, and is taken as

[TABLE]

The final term in (6) describes the adhesive behaviour between the cancer cells and the fibres distributed on $\mathbf{B}(x,R)$. Within this term we account for the role of the fibres in two ways. On one hand, we account for the influence of the macroscopic fibres magnitude $F(x+y,t)$ has over the adhesion that the cells distributed at the spatial location $x$ exercise with adhesion strength $\textbf{S}_{{}_{cF}}$ upon the fibres distributed at $x+y$. On the other hand, their full macroscopic orientation $\theta_{f}(x,t)$ (induced by their micro-scale mass distribution of microfibres), biases the cell-matrix adhesion in the direction of the resultant vector

[TABLE]

The last important aspect that we consider here concerns mutations that enable cells from primary tumour cell subpopulation $c_{1}$ to undergo genetic conversions and become secondary cancer cells $c_{2}$, process that is explored here through the mutation term $M_{c}(\textbf{u},t)$, which captures the direct correlation between the presence of significant ECM levels, as detailed in the following. Indeed, as cellular adhesion is controlled by the binding of the calcium dependent adhesion molecules, E-cadherins, with the intra-cellular proteins, catenins, it is significant to note that collagens present in the ECM (in particular collagen type I) down-regulate the gene expression of E-cadherins, resulting this way in a loss of cell-cell adhesion [Menke et al., 2001]. This loss of adhesion is accompanied by an increase in proliferation and migratory activity, and so, to reflect this behaviour, the mutation term $M_{c}(\textbf{u},t)$ is dependent on the underlying ECM density levels and is taken to be

[TABLE]

where $\kappa$ is a certain level of ECM beyond which mutations can occur and $H(\cdot)$ is the usual Heaviside function, with $t_{m}$ being the time at which mutations begin.

Finally, we describe the dynamics of the ECM, considering the individual constituents of the matrix, namely the fibre and non-fibre component. Both constituents of the matrix are simply degraded by the cancer cells, and so their macroscopic dynamics can be mathematically written as

[TABLE]

where $\gamma_{1}$ and $\gamma_{2}$ are the degradation rates of the fibre and non-fibre components respectively, and $\omega$ describes the rate of non-fibre ECM remodelling. Furthermore, the matrix remodelling, which is important both in the development and in progression of cancer, contributing to processes such as metastasis and tumour cell invasion [Cox and Erler, 2011], is controlled here by the volume filling factor $(1-\rho(\textbf{u}))$.

2.3 Microscopic fibre rearrangement

As explored and modelled in Shuttleworth and Trucu [2019], during their invasion, the cancer cells push the surrounding fibres in accordance with the emerging cell-flux direction and rearrange their micro-fibre mass distribution managing this way to reorient the macro-scale fibres. Indeed, in addition to their macroscopic degradation described in (12), the fibres go through a microscopic rearrangement process induced by the macro-dynamics of the cancer cells. Specifically, at time $t$ and at any spatial location $x\in Y$, the cancer cells will realign the micro-fibres through a microscopic rearrangement process in each micro-domain $\delta Y(x)$ that is induced by the combined macro-scale spatial flux of both cancer cell sub-populations

[TABLE]

where

[TABLE]

The combined flux acts upon the micro-scale distribution $f(z,t)$, $\forall z\in\delta Y(x)$ in accordance to the magnitude that the total mass of cancer cells has relative to the combined mass of cells and fibres at $(x,t)$, which is given by the weight

[TABLE]

At the same time, the spatial flux of cancer cells $\mathcal{F}(x,t)$ is balanced in a weighted manner by the orientation $\theta_{{}_{f}}(x,t)$ of the existing distribution of fibres at $(x,t)$ that is appropriately magnified by a weight that accounts for the magnitude of fibres versus the combine mass of cells and fibres at $(x,t)$ and is given by $(1-\omega(x,t))$. As a consequence, the microscale distribution of micro-fibres $f(z,t)$, $\forall z\in\delta Y(x)$ is therefore acted upon uniformly by the resultant force given by the following macro-scale vector-valued function

[TABLE]

As detailed in Shuttleworth and Trucu [2019], under the uniform incidence of the resultant force $r(\delta Y(x),t)$ upon the mass of micro-fibres distributed at any $z\in\delta Y(x)$, an on-the-fly spatial microscopic rearrangement of this micro-fibres mass takes place. Specifically, under the influence of $r(\delta Y(x),t)$, an appropriate level of micro-fibres mass $f(z,t)$ will undergo a spatial transport towards a new position

[TABLE]

where the relocation direction and magnitude is given by

[TABLE]

Here we have $x_{\text{dir}}(z)=\overrightarrow{x\,z}$ representing the barycentric position vector pointing to $z$ in $\delta Y(x)$, which enables us also the quantification of the position defect of $z$ with respect to $r(\delta Y(x),t)$ namely

[TABLE]

which affects the spatial relocation of micro-fibres mass. Furthermore, $f_{\text{max}}$ represents a level of micro-fibres that can be distributed at $z\in\delta Y(x)$, with the micro-fibres mass relocation in the direction of $\left(x_{\text{dir}}(z)+r(\delta Y(x),t)\right)$ being enabled provided that their level is below $f_{\text{max}}$. Finally, alongside the position defect of $z$ with respect to $r(\delta Y(x),t)$, another aspect that affects the micro-fibres relocation is the level of micro-fibres distributed at location $z$, which is accounted for in (15) through the micro-fibres saturation fraction

[TABLE]

Thus, a micro-fibres mass transport from $z$ to the location $z^{*}$ is exercised provided that micro-fibres level at $z$ are not at their maximum level $f_{\text{max}}$, while lower levels of micro-fibres saturations at $z$ together with a better position alignment given by a smaller position defect lead to a relocation of the micro-fibres mass in direction $\left(x_{\text{dir}}(z)+r(\delta Y(x),t)\right)$ at a greater distance, resulting in reaching a position $z^{*}$ that is further away from $z$. Finally, this micro-fibres transport is also regulated by the capacity available at the new position $z^{*}$, which is explored here through a movement probability

[TABLE]

that enables only an amount of $p_{move}f(z,t)$ of micro-fibres to be transported to position $z^{*}$, while the rest of $(1-p_{move})f(z,t)$ remain at $z$.

2.4 The multiscale moving boundary approach

Let us now briefly revisit the novel multiscale moving boundary framework initially introduced in [Trucu et al., 2013] and later evolved to consider two cancer cell subpopulations in [Shuttleworth and Trucu, 2018], which explore the cell-scale proteolytic activity of MDEs along the invasive edge of the tumour that is non-locally induced by the tissue-scale cancer cell population dynamics and that degrades the peritumoural ECM determining this way the direction and associated displacement for tumour boundary progression. As detailed in [Trucu et al., 2013], the link between the tumour macro-dynamics (4), (5), (12) and (13) and the proteolytic enzyme micro-dynamics is captured via a double feedback loop that is realised by a top-down link describing the source of MDEs induced at cell-scale by the spatial distribution of cancer cells at tissue-scale, and a bottom-up link describing the translation of the resulting boundary relocation to the tissue-scale.

Top-down link

As previously discussed, cancer invasion is a multiscale process in which the MDEs secreted by tumour cells from the outer proliferating ring undergo a cell-scale spatial transport in the neighbourhood of the invasive edge of the tumour and are responsible for the degradation of the surrounding ECM. It is this breakdown of peritumoural ECM that enables the tumour opportunities to expand and proceed with its local invasion. Continuing with the terminology of the framework introduced in Trucu et al. [2013], during a time interval $[t_{0},t_{0}+\Delta t]$, the MDEs micro-dynamics is explored on the invasive leading edge of the tumour $\partial\Omega(t)$ enclosed by a complete cover of $\epsilon$-size half-way shifted overlapping micro-domains $\{\epsilon Y\}_{\epsilon Y\in P(t)}$. Furthermore, the specific topological requirements detailed in full in Trucu et al. [2013] that enable the construction of the covering bundle of microdomains $\{\epsilon Y\}_{\epsilon Y\in P(t)}$, allow us to capture the cell-scale MDEs activity in a neighbourhood of $\partial\Omega(t)$, exploring this molecular dynamics in both the overlapping inner regions $\epsilon Y\cap\Omega(t)$ and the peritumoural outside regions $\epsilon Y\cap\setminus\Omega(t)$ where the MDEs get transported and degrade the ECM.

At each $z\in\epsilon Y\cap\Omega(t)$ a source of MDEs is induced as a collective contribution of both cancer cell sub-populations that arrive during their dynamics in the outer proliferating rim within a maximal distance $\gamma>0$ with respect to $z$ (given by the maximal thickness of the outer proliferating rim), and so this can be mathematically given by

[TABLE]

where $\mathbf{B}(y,\gamma):=\{\xi\in Y\,|\,\parallel y-\xi\parallel_{{}_{\infty}}\leq\gamma\}$; $\alpha_{i}$, $i=1,2$, are the MDE secretion rates for the two cancer cell sub-populations, $\lambda(\cdot)$ is the standard Lebesgue measure on $\mathbb{R}^{2}$; and $\gamma$ is a small distance between the zero source level from outside $\Omega(t_{0})$ and the non-zero source levels on $\epsilon Y\cap\Omega(t_{0})$ where via Urysohn Lemma [Yosida, 1980] we ensure a continuous transition. Once secreted by the cancer cells, the MDEs molecular distribution denoted here by $m(y,\tau)$ exercise a cell-scale cross interface spatio-temporal transport which in the simplest context of a generic MDEs is considered as being given by a local diffusion. Thus, given the topological properties of the covering bundle of overlapping micro-domains $\{\epsilon Y\}_{\epsilon Y\in P(t)}$ detailed in Trucu et al. [2013], this enable us to decompose this cross-interface diffusion process (into a corresponding bundle of micro-processes) and to explore this per each individual $\epsilon Y$, where this micro-dynamics is given as

[TABLE]

Bottom-up link

During the micro-dynamics, the MDEs present in the peritumoural region interact with the ECM distribution captured within each $\epsilon Y$ micro-domain. The local degradation of ECM is dependent on the advancing spatial distribution of MDEs in $\epsilon Y\backslash\Omega(t)$ secreted by the cancer cells. This pattern of degradation gives rise to a movement direction, $\eta_{\epsilon Y}$, and a displacement magnitude, $\xi_{\epsilon Y}$ (detailed in Trucu et al. [2013]), that ultimately dictate the movement of the boundary midpoint $x_{\epsilon Y}^{*}$ to a new position $\widetilde{x^{*}_{\epsilon Y}}$. This process is the catalyst behind the expansion of the macroscopic tumour boundary. Thus, the bottom-up link of the model between the molecular activities of MDEs on the tumour invasive edge and the macroscopic boundary $\Omega(t_{0})$, is realised through the macro-scale boundary movement effectuated by the micro-dynamics of the proteolytic enzymes, resulting in an expanded tumour domain, $\Omega(t_{0}+\Delta t),$ on which the invasion process continues.

2.5 Summary of model

As stated previously, the two-part multiscale model is comprised of two interconnected multiscale systems that share the same macro-scale whilst having their own distinct micro-scale dynamics, as schematically summarised in Figure 3. The macro-scale dynamics governs the spatial distribution of cancer cells and both the fibre and non-fibre components of the ECM. The first multiscale system controls the dynamic redistribution of micro-fibres, weighted according to the cancer cell and macroscopic fibre distributions and triggers a micro-scale spatial rearrangement of micro-fibres via the combined spatial flux of the two cancer cell sub-populations. The second multiscale system reacts to a source of MDEs induced on the boundary by the spatial distribution of cancer cells. The microscopic distribution of MDEs instigates degradation of the peritumoural ECM and it is from this pattern of degradation that the position of the tumour boundary is changed. The microscopic change in the boundary is translated back to the tissue-level and the macro-dynamics continue leading to a newly expanded tumour region, where the cancer invasion process continues its dynamics.

3 Numerical approaches and initial conditions for computations

Building on the multiscale moving boundary computational framework introduced in Trucu et al. [2013] combined with its very recent novel extension introduced in Shuttleworth and Trucu [2019] to address the multiscale fibres dynamics in the bulk of the tumour, we developed a new computational approach to explore the complex multiscale evolution of a heterotypic tumour with two cancer cell sub-populations within a fibrous environment. We consider a uniform spatial mesh of size $h=0.03125$ to solve the macroscale computations on the expanding tumour domain, whilst using an off-grid approach for the calculation of the macroscopic adhesion term that decomposes the cell-sensing region to approximate the adhesive flux at each spatio-temporal node. To complete the tumour macro-dynamics we use a novel predictor-corrector scheme developed in Shuttleworth and Trucu [2019] that accounts for the complexity of the cancer dynamics. Further, to obtain the microscopic boundary relocation, we explore the top-down and bottom-up links as well as solve the micro-dynamics via finite element approximation as detailed in Trucu et al. [2013]. The simulations of the model equations were all performed on MATLAB.

In the following sections we consider the initial cancer cell population $c_{1}$ to occupy the region $\Omega(0)=\mathbf{B}((2,2),0.5)$ positioned in the centre of $Y$, Figure 4, whilst the distribution of cell population $c_{2}$ is initially zero, i.e.,

[TABLE]

where $\psi_{\gamma}$ is the standard mollifier detailed in C that acts within a radius $\gamma<<\frac{\Delta x}{3}$ from $\partial\mathbf{B}((2,2),0.5-\gamma)$ to smooth out the characteristic function $\chi_{{}_{\mathbf{B}((2,2),0.5-\gamma)}}$. Furthermore, for the non-fibre ECM phase, we consider both a homogeneous and heterogeneous initial conditions, with the homogeneous initial conditions illustrated in Figure 5(a) and given by

[TABLE]

and heterogeneous initial condition shown in Figure 5(b) and given by

[TABLE]

where

[TABLE]

these being previous used also in Shuttleworth and Trucu [2019].

For the initial distribution of the ECM fibre phase, we assume a distribution of five pre-assigned micro-fibre distributions of the five different micro-fibre patterns generated along the family of path-networks $\{P^{1}_{i}\}_{i=1..5}$ and assigned randomly onto $\delta Y(x):=x+\delta Y$ as detailed in A). To represent a homogeneous distribution of fibres, we calibrate the maximal height of the microfibres to be uniform across all micro-domains, resulting in a macroscopic fibre distribution $F(\cdot,t)$ that represents a percentage $p$ of the mean density of the non-fibrous homogeneous ECM phase (20), illustrated in Figure 6(a)-(b). On the other hand, a heterogeneous distribution of fibres is achieved by allowing the maximal height of the microfibres to be determined by a heterogeneous pattern, for example, the initial condition for a heterogeneous ECM non-fibre phase (21), Figure 6(c)-(d). We set the maximal height of the micro-fibres in each $\delta Y(x)$, centred at $x$, to correspond to the distribution $l(x,0)$ for all $x\in Y$, resulting in the macroscopic fibre magnitude spatial distribution $F(\cdot,t)$ to represent again a percentage $p$ of the heterogeneous non-fibre ECM phase.

4 Computational simulation results

4.1 Homogeneous non-fibre ECM component

In the first instance, the initial condition for the non-fibre ECM component will be taken as the homogeneous distribution, namely as $l(x,0)=min\{0.5,1-c_{1}(x,0)\}$. The initial condition for both cancer cell populations is given in (18), (19), and the fibres take an initial homogeneous macroscopic distribution of $15\%$ the non-fibre phase, with the combined ECM density $v(x,0)=l(x,0)+F(x,0)$. Using the parameter set $\Sigma$ from D, we show the computational results at the final time stage for the evolution of both individual cancer cell populations in subfigures 7(a), 7(b), the non-fibre ECM phase and fibre ECM phase in subfigures 7(c) and 7(d) respectively, presented with the vector field of oriented fibres coarsened four times, 7(e), and a 3D representation of fibres 7(f).

Figure 7 displays computations at final stage $75\Delta t$, during which time mutations from population $c_{1}$ to $c_{2}$ have begun, starting at stage $5\Delta t$. Previously it was observed that in the presence of a homogeneous ECM, a tumour was limited to a symmetric pattern of invasion [Shuttleworth and Trucu, 2018], however with the addition of an initially homogeneous fibre distribution we witness a different pattern of invasion in 7(c). The primary cell population $c_{1}$, 7(a), has been manipulated and reshaped by the underlying fibre network. Although the fibre phase of the ECM begins as a macroscopically homogeneous distribution, once fibre rearrangement occurs, higher density regions of fibres are formed, 7(d), and these areas become increasing influential on the direction of invasion. Since mutations have begun, we see the emergence of cells in population $c_{2}$, 7(b). The pattern of cells in population $c_{2}$ is correspondent to the fibre orientations displayed in 7(e), where the general direction of fibres are aiming roughly towards the origin $(0,0)$. Thus, the distribution of cells present in Figure 7(b) are more prominent in this direction. This behaviour is attributable to population $c_{2}$ retaining a higher affinity for cell-fibre adhesion than population $c_{1}$. Throughout the simulations, the ECM is degraded by the cancer cells, most notable by the large low density region in 7(c), with the highest levels of degradation relative to the regions of highest cell distribution.

It is clear from the computational results presented in Figure 7 that the tumour boundary has expanded in the general macroscopic fibre direction and strayed away from a symmetric route of invasion. Hence we conclude from these computations and comparisons with previous work, [Shuttleworth and Trucu, 2018], that the underlying fibre network strongly influences the direction of invasion.

To examine the effects of an underlying fibre distribution, we increase the initial homogeneous fibre distribution to $20\%$ of the non-fibre ECM phase. Continuing with the same initial conditions for the non-fibre ECM phase and cancer distributions, we present the computations at final stage $70\Delta t$, Figure 8. Comparing with Figure 7, cells of population $c_{1}$ are exhibiting limited migration within the tumour boundary due to their high self-adhesion rate, 8(a), whilst cells of population $c_{2}$, 8(b), have migrated in a direction consistent with the macroscopic fibre orientations in 8(e). The mutated cell population $c_{2}$ congregates around the cells of population $c_{1}$, this time showing more regions of high distribution, suggesting the increased fibre distribution enables faster migration of cells in population $c_{2}$. The entire ECM (both the fibre and non-fibre phase) presented in subfigures 8(c)-(d) has been degraded in the centre of the domain where the cancer cells are most dense.

In line with the previous simulations, the micro-fibres have been pushed outwards towards the boundary, however this time the macroscopic representation of the fibres appears as a region of very low fibre density surrounded by a mottled pattern of fibre densities, Figure 8(d). This behaviour occurs because of two reasons. Firstly, the tumour boundary is progressing faster than the fibres are being rearranged, and because the rearrangement of fibres is restricted to neighbouring domains only, the fibre distributions will never be found to build up explicitly at the rapidly expanding tumour boundary. Secondly, as the degradation of fibres is kept within the bounds of the tumour and dependent on the cancer cell distribution, out-with the dense main body of cancer cells, we see little degradation of the fibres. The increased size of spread of the tumour indicates that the ratio of the underlying fibre network to the non-fibre ECM phase plays a key role in the success of local tumour invasion.

We can investigate a similar scenario to that presented previously whereby we now explore tumour invasion on an initially heterogeneous fibre network. The initial condition for the macroscopic fibre distributions will be taken as $15\%$ of the heterogeneous distribution defined in (21) whilst the non-fibre ECM phase will be kept as a homogeneous distribution. Figure 9 displays computations at final stage $75\Delta t$. Comparing to Figure 7, the boundary of the tumour has formed lobules, specifically towards low density regions of the fibre network, most notable in subfigure 9(d). A different pattern of mutations from population $c_{1}$ to $c_{2}$ is adopted due to the heterogeneity of the fibre network and thus high distribution regions of cells are formed in 9(b) that differ to that when in the presence of a homogeneous fibre phase.

The macroscopic mean value of fibres has had a substantial impact on the overall invasion of cancer. The boundary of the tumour has grown many lobules, first reaching to the lower density areas of fibres in 9(d) before engulfing the higher density regions as invasion progresses. The high distributions of cells in population $c_{2}$ are reminiscent of Figure 7(b) in that they are building according to the general orientation of the fibre network in 9(e). Under the presence of a heterogenous fibre network, tumour invasion is halted at regions of high macroscopic fibre density, whilst in a homogeneous environment there are no barriers during invasion, evidenced in Figure 7. Here the tumour experiences a lobular expansion as it invades the surrounding ECM whilst the high macroscopic fibre distributions coupled with the higher cell-fibre adhesion rate that population $c_{2}$ carries allows the cells more migratory freedom within the tumour boundary.

We consider the same initial conditions, but we now introduce matrix remodelling to the non-fibre ECM phase, namely, increasing $\omega$ in (13) from [math] to $0.02$. Initially the computational results are very similar to that of the previous section, however, as tumour invasion progresses to stage $75\Delta t$, Figure 10, both the boundary of the tumour and the main body of cells looks very different. Although the density of tumour cells is not as high as in Figure 9, there is an increase in cell distribution within the boundary of the tumour, particularly in cell population $c_{2}$, 10(b), with the entire boundary of the tumour containing a visible amount of cells. This behaviour is as a consequence of matrix remodelling only occurring under the presence of cancer cells, hence, the increased density of matrix within the boundary of the tumour gives more opportunity for cell adhesion and opens more pathways in which the cells can invade. Due to the heterogeneity of the fibre ECM phase, the tumour again exhibits a lobular route of invasion in 10(d), first enveloping the low density regions where space is available for the cells to freely invade and consequently overrunning the high density regions in their path.

4.2 Heterogeneous non-fibre ECM component

Investigating the effects of different initial fibre distributions present in a homogeneous non-fibre ECM domain is effective in showing the influence of fibres on tumour invasion, however, is it crucial to examine these varying fibre distributions under a heterogeneous non-fibre ECM phase. We express heterogeneity of the non-fibre ECM phase using the initial condition (21) and we assume an initial $15\%$ homogeneous fibre distribution. Although the cell-non-fibre-ECM adhesion coefficients are kept low, we see the tumour of the boundary developing a lobular pattern, witnessed in subfigure 11(c), encroaching on the lower density regions of matrix. By first advancing on the low density regions of the matrix, the cancer cells proceed by engulfing the high density regions in their path and thus the tumour region becomes larger. The tumour continues to progress in this manner with small islands appearing over the regions of high matrix density. These islands arise when the cancer cells have failed to fully overrun areas within the ECM. This behaviour can occur when there are not yet enough cells to fill the area, or because the matrix density is simply too high and the cells must first degrade part of the matrix in order to make space available. Population $c_{2}$ is sparse in its spread 11(b), with only a few regions of cells visible, having gathered into pockets of the matrix where density of both fibre and non-fibre ECM is low, subfigures 11(c)-(d). Comparing to the simulations in Figure 9, where initially the non-fibre ECM phase is homogeneous and the fibre phase is heterogeneous, we see a much larger overall spread of the tumour when the non-fibre ECM phase is taken as an initially heterogeneous distribution. This comparison implies that the non-fibre ECM phase plays a key role in the local invasion of cancer and is largely responsible for the progression of the tumour.

Continuing to explore invasion in a heterogeneous non-fibre ECM phase, we increase the initial homogeneous fibre distribution to $20\%$ of the non-fibre phase. Figure 12 give simulations at the final stage $75\Delta t$. The tumour boundary is larger and more advanced when compared with the simulations in Figure 11, which have a lower initial fibre distribution. The increased fibres distributions within the ECM presents the cells with more opportunities for adhesion and as a result the cells spread at an increased rate into the surrounding matrix. Cell population $c_{2}$ has progressed outwards towards areas of the leading edge and formed small masses of cells in 12(b). The direction of migration is similar to that in Figure 8, implying that even though they are subject to a rearrangement process, the initial fibre micro-domains influence the direction of migration of cells. It should also be noted that the density of cells in population $c_{2}$ is lower than that to previous simulations, this occurs as the cells have increased migration due to the increased fibre distributions, thus the cells do not build up in any one space, Figure 12(b). Cells of population $c_{1}$ are surrounded by the more aggressive cells of population $c_{2}$, hence their invasion of the tissue is restricted and they have created denser regions of cells through proliferation, Figure 12(a).

We next consider here a fully heterogeneous ECM, with the fibre phase taking $15\%$ of the heterogeneous non-fibre phase, defined in (21). Figure 13 shows simulations at final stage $70\Delta t$. Cell population $c_{1}$ continues to be confined to the centre of the domain showing little outward movement whilst exhibiting the highest distribution of cells noted throughout all previous simulations. This increase in distribution can be described by the heterogeneous nature of the matrix; as the patches of low density ECM become more prominent, the cells have less areas to adhere and are therefore restricted to proliferation in these areas only, forming larger masses of cells. In contrast, cell population $c_{2}$ has established patches of high cell distributions following the general direction of the fibres, Figure 13(e), and we again see small islands within the boundary forming over the lower density regions of ECM, Figure 13(c), similar to the simulations in Figure 11. The tumour boundary is visibly more lobular and consistent in its pattern by following the structure of the ECM, attributable to the heterogeneity of the entire matrix. In line with previous simulations, the tumour is advancing on the higher density regions of matrix first, hence the fingering behaviour of the boundary towards these regions, Figure 13(c)-(d).

4.3 New family of pre-determined micro-fibre domains

Currently, we have explored the local invasion of a tumour under the presence of randomly placed micro-domains, allocated from a family of five pre-determined micro-domains. Now we consider the impact changing two of these micro-domains to contain a different pattern of micro-fibres (i.e., those given by the second family of paths $\{P^{2}_{i}\}_{i=1..5}$ defined in A) will have on the overall invasion of cancer. Using the same initial condition as before for the cell populations (18), (19) and a homogeneous initial condition for the non-fibre ECM phase along with a $15\%$ fibre phase, we obtain the computational results in Figure 14 at final stage $75\Delta t$. Comparing with the results using the previous family of micro-domains, Figure 7, the simulations appear similar, with only slight differences in the shape of the tumour boundary and the composition of the cells within each subpopulation. Cells in population $c_{2}$ have formed an increased number of high density regions in 14(b), particularly at the lower part of the tumour boundary. The degradation of ECM is consistent to that of the previous results, showing regions of high degradation where the cells are of highest distribution, observed in subfigures 14(c)-(d).

To complete our comparison of the pre-determined fibre micro-domains, we explore cancer invasion within a heterogeneous $15\%$ fibre component and a homogeneous non-fibre ECM component. Figure 15 displays simulations at final stage $75\Delta t$, and when compared to Figure 9 there are definite changes in the migration of cells in population $c_{2}$, Figure 15(b), with pools of high cell distributions forming in different regions to that of Figure 9(b). Here, the tumour boundary is remaining closer to the main body of cells, with degradation of the ECM occurring explicitly within the boundary, subfigures 15(c)-(d). The boundary of the tumour has expanded in a pattern consistent with the heterogeneity of the fibre phase 15(d). Population $c_{1}$ has developed in a similar fashion to that in Figure 9, however cells in population $c_{2}$ have formed regions of high cell distributions around the entire boundary 15(b). From these simulations we cannot strictly conclude which family and arrangement of fibre micro-domains gives the more aggressive spread of cancer, however we can deduce that the initial placement of fibre micro-domains has a strong influence on the overall invasion of cancer, owing to the initial orientation of the fibres.

5 Discussion

We have presented a novel multiscale moving boundary model with dynamic fibre redistribution, first introduced in [Shuttleworth and Trucu, 2019], which now includes a heterotypic cancer cell population, with two sub-populations, $c_{1}$ and $c_{2}$. This model was developed to explore both the random disposition of microscopic fibre distributions and the inclusion of a heterogeneous cancer cell population and their combined effects on the overall invasion pattern during tumour growth and local invasion of tissue. Paying attention in the first instance to the heterogeneous cell population, we introduced a secondary cell subpopulation, mutated from the primary cell population, which holds a lower cell-cell adhesion rate and a higher affinity for cell-fibre adhesion. Mutations from the primary to secondary subpopulation were irreversible, explored through a mutation rate dependent on the underlying ECM density. Focussing next on the fibre component of the ECM, we explored several scenarios; a homogeneous and heterogeneous macroscopic distribution of fibres, varying the initial ratio of macroscopic fibre distributions in relation to the non-fibre ECM phase, and finally the exploration of randomly allocated fibre micro-domains. The multiscale dynamics of the underlying fibre network were explored and modelled within a two-part multiscale model, where the micro-scale dynamics are connected to the macro-scale through a double feedback loop. We considered the macroscopic representations of the micro-fibres and their resulting effects on the macro-scale dynamics, in particular cellular adhesion, whilst at the same time allowing the distribution of both cancer cell subpopulations to cause macroscopic fibre degradation. In addition, through their spatial flux the cancer cells were able to influence the direction of microscopic fibre rearrangement. Finally, we include to this model the previously developed multiscale moving boundary framework developed in [Shuttleworth and Trucu, 2018] that considers the effects of a heterogeneous cancer cell population and its combined contribution to the leading edge proteolytic MDE dynamics. Thus, the model proposed here combines two multi-scale systems that contribute to the same tissue- (macro-) scale dynamics whilst having separate cell- (micro-) scale processes that are simultaneously connected via two double feedback loops. As in [Shuttleworth and Trucu, 2019], the first multi-scale model gives the cell-scale fibre rearrangement process within the tumour region, and the second multi-scale system describes the MDE proteolytic activity within a cell-scale neighbourhood of the tumour boundary.

To accommodate a heterogeneous cancer cell population, we adapt the macro-scale dynamics introduced in [Shuttleworth and Trucu, 2019] to include a similar set of dynamics for a second cell subpopulation which has been mutated from the primary cell population. In doing so, we adjust the term for cellular adhesion to include both cell subpopulations allowing for different parameters for each population [Domschke et al., 2014], implementing a weaker self-adhesion coefficient and a higher rate of cell-fibre adhesion for the mutated population to promote migration and thus promote the overall invasive capabilities of the tumour. The heterogeneous cancer cell population is then incorporated within the micro-scale dynamics of both multi-scale systems. Firstly, it has influence in determining the source term for MDEs within each cell-scale neighbourhood, with the mutated population inducing a higher amount of MDEs than the primary cell population. Secondly, the distribution and spatial flux of cancer cells that determines the direction of micro-fibre redistribution is now taken as the combined cancer cell distributions and the addition of the spatial fluxes of both populations. Finally, we use the combined cell distributions to contribute equally to the degradation of both the fibre and non-fibre components of the ECM.

Comparing with results from the previous framework developed in [Shuttleworth and Trucu, 2019] that describes the local invasion of a single population of tumour cells in the presence of a macroscopically homogeneous fibre network, we witness some similarities with our simulations. In the presence of a heterogeneous component of the ECM (fibre or non-fibre phase) we observe a fingering, lobular pattern of invasion, where the cells first flood the low density cavities and proceed to engulf the higher density regions on their route of invasion. In the presence of a homogeneous non-fibre ECM phase and an initially homogeneous fibre phase, we see the boundary of the tumour stray from the symmetric invasion witnessed in [Shuttleworth and Trucu, 2018], implying that the underlying fibre network indeed plays a role in tumour invasion. Increasing the initial macroscopic fibre density by only $5\%$ from $15\%$ to $20\%$ of the non-fibre ECM phase results in a very fast growing tumour that almost takes over the entire domain. This solidifies the reasoning that the fibre network plays a key role in the local invasion of cancer. Additionally, we have explored this model with two sets of randomly placed micro-fibre domains. Comparing the simulations of each family of micro-domains, particularly under the presence of both components of the ECM being initially homogeneous, the pattern of the subpopulations of cancer cells is affected. The initial orientation of the fibres has influenced the migration of the cells, most obviously observed within population $c_{2}$, where the cell-fibre adhesion rate is higher. As expected, this concludes that the orientation of the fibres is key during tumour growth and development, particularly at the initial stage of invasion.

Looking forward, this modelling framework allows for the opportunity to fully incorporate both multi-scale models; incorporating the micro-scale fibre network within the MDE micro-scale neighbourhoods by investigating the cellular effects of MDEs on the fibres, i.e., the degradation/slicing of the fibres by the matrix metallo-proteinases (MMPs) such as the freely moving MMP-2 and the membrane-bound MT1-MMP molecules. The exploration of TACS (tumour associated collagen signatures) could provide important information of the severity of the tumour. Biological experiments have shown that TACS-3, which is characterised by bundles of straightened and aligned collagen fibres that are orientated perpendicular to the tumour cells provide the poorest prognosis for patients [Provenzano et al., 2008], whilst tumours in a TACS-2 environment, where the angle between cells and the fibres are between $0\degree-30\degree$, have been shown to provide an environment in which tumour progression is slower [Conklin et al., 2011]. Further work will be focussed on these areas, with the hope that as both the model and biological experiments develop, we can obtain a more realistic representation of fibres in vivo and incorporate these into the model to gain a deeper understanding of the role of fibres within the ECM and during the local invasion of cancer.

Appendix A Microscopic fibre domains

For the fibres initial conditions, on the micro-domains $\delta Y(x)$, we consider two families of five distinctive micro-fibres patterns, $\{P^{1}_{i}\}_{i\in J}$ and $\{P^{2}_{i}\}_{i=1..5}$, which are defined by the union of paths $P^{l}_{i}=\bigcup\limits_{j=1..5}h^{l}_{i,j}$, $l=1,2$, which are given as follows.

For the first family of fibre paths $P^{1}$, we have:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For the second family of fibre paths $P^{2}$, we have:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For each of these families of fibre paths $P^{l}$, $l=1,2$, as described in Shuttleworth and Trucu [2019], the micro-scale fibres pattern within each micro-domain $\delta Y(x)$ is given as

[TABLE]

where $\{\psi_{h^{l}_{i,j}}\}_{i,j=1..5}$ are smooth compact support functions of the form

[TABLE]

Here $r>0$ is the width of the micro-fibres and $C_{h^{l}_{i,j}}$ are constants that determine the maximum height of $\psi_{h^{l}_{i,j}}$ along the smooth paths $\{h^{l}_{i,j}\}_{i,j=1..5}$ in $\delta Y(x)$. Finally, $\psi_{\gamma}$ is the standard mollifier defined in C, with $\gamma_{{}_{0}}=h/16$.

Finally, the initial spatial configuration of the pattern of macroscopic ECM fibre phase is selected according to a randomly generated matrix of labels $A=(a_{i,j})_{i,j=1..n}$ corresponding to the entire $n\times n$ grid discretising $Y$, in which the entries $a_{i,j}$ are allocated values randomly selected from the set of configuration labels $\{1,2,3,4,5\}$ that will dictate the choice of micro-fibres pattern among those described above that will be assigned to the micro-domains $\delta Y(j\Delta x,i\Delta y)$, for all $i,j=1..n$.

Appendix B Adhesion matrices

The cell adhesion matrices use in the manuscript were considered as follows:

[TABLE]

Appendix C The mollifier $\psi_{\gamma}$

The standard mollifier $\psi_{\gamma}:\mathbb{R}^{N}\to\mathbb{R}_{+}$ (which was used also in [Trucu et al., 2013, Shuttleworth and Trucu, 2019]) is defined as usual, namely

[TABLE]

where $\psi$ is the smooth compact support function given by

[TABLE]

Appendix D Table for the parameter set $\Sigma$

Here we present a table for the parameter set $\Sigma$.

Acknowledgment

RS and DT would like to acknowledge the support received through the EPSRC DTA Grant EP/M508019/1 on the project: Multiscale modelling of cancer invasion: the role of matrix-degrading enzymes and cell-adhesion in tumour progression.