Mathematical Modelling of Glioblastomas Invasion within the Brain: A 3D Multi-Scale Moving-Boundary Approach

Abstract

Brain-related experiments are limited by nature, and so biological insights are often limited or absent. This is particularly problematic in the context of brain cancers, which have very poor survival rates. To generate and test new biological hypotheses, researchers have started using mathematical models that can simulate tumour evolution. However, most of these models focus on single-scale 2D cell dynamics, and cannot capture the complex multi-scale tumour invasion patterns in 3D brains. A particular role in these invasion patterns is likely played by the distribution of micro-fibres. To investigate the explicit role of brain micro-fibres in 3D invading tumours, in this study, we extended a previously introduced 2D multi-scale moving-boundary framework to take into account 3D multi-scale tumour dynamics. T1 weighted and DTI scans are used as initial conditions for our model, and to parametrise the diffusion tensor. Numerical results show that including an anisotropic diffusion term may lead in some cases (for specific micro-fibre distributions) to significant changes in tumour morphology, while in other cases, it has no effect. This may be caused by the underlying brain structure and its microscopic fibre representation, which seems to influence cancer-invasion patterns through the underlying cell-adhesion process that overshadows the diffusion process.

1. Introduction

Glioblastoma multiforme is a highly invasive and aggressive type of brain tumour, typically with poor patient prognosis [1,2,3,4,5,6,7] (with a median survival rate of less than one year [8]). These tumours arise from abnormal glial cells located in the central nervous system, and shortly after their appearance, they invade the surrounding tissues in a heterogeneous fashion. This heterogeneous invasion pattern leads to tumours whose outer edges are difficult or impossible to determine with current imaging technologies, including for instance magnetic resonance imaging (MRI) and diffusion tensor imaging (DTI), both of which measure the diffusion of water molecules and enable the study of brain structures.

Due to the limited experimental approaches that one can use to study the brain, researchers have started using mathematical models to provide certain biological insights that otherwise would be difficult to obtain experimentally. Such models can help predict how tumours grow for specific patients, aiding clinicians in decision-making, or they can help test and provide new hypotheses about potential anti-tumour treatments. The mathematical modelling of tumours has seen significant advances in recent decades, which have broadened our understanding of tumour dynamics and how cells interact with their environment [9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]. Although the majority of these models do not restrict themselves to a specific tumour type and rather focus on general tumours, there are some models that focus on the evolution of gliomas within the brain [34,35,36,37,38,39,40,41,42]. Recently, models have also started to incorporate the structure of the brain, by including MRI and DTI scans [35,36,37,43,44,45,46]. Even though these images are generated in 3D, most of these models are simulating tumour growth in 2D and only a few of them are 3D models [43,47,48]. Moreover, the majority of published models have focused on tumour progression only on one spatio-temporal scale. However, tumour progression is characterised by various biological processes occurring on different scales, and thus their effects on the overall tumour dynamics cannot be neglected. Hence, recent efforts have been made to establish new multi-scale frameworks for tumour progression [26,27,28,29,30,35,36,37,49], which were able to capture some of these underlying multi-scale biological processes usually involving the extracellular matrix (ECM).

In this paper, we extended the general 2D multi-scale moving-boundary modelling framework introduced in [20,26] to capture the invasion of glioblastomas within a 3D fibrous brain environment. To this end, we incorporated the information provided by both the T1 weighted and DTI scans into our multi-scale framework and used the resulting model to numerically simulate the growth of 3D gliomas within the brain. We focused on a few cases showing tumour growth in different regions in the brain, with different distributions of grey/white matter densities, which leads to different tumour invasion patterns.

The model presented in this paper is a first to advance the multi-scale moving boundary modelling of 3D brain gliomas in two key ways compared to previous ones, as can be seen for instance in [35,36,37]. On the one hand, using this modelling framework, we explicitly address the micro-fibres and their rearrangement dynamics caused by the macro-scale movement of the cancer cell population (a process of great importance as the micro-scale structural changes of the underlying tissue can significantly influence the cell population motility). On the other hand, we also model another vital proteolytic process that occurs along the invasive edge of the tumour, allowing us to explicitly model tumour morphological changes and how the tumour boundary progresses over a time interval.

This paper is organised as follows. First, we formulate our extended multi-scale moving boundary framework in Section 2. Then, following a brief description of the numerical methods, we present the computational simulation results in Section 3. Finally, in Section 4, we summarise and discuss these results.

2. Multi-Scale Modelling of the Tumour Dynamics

To model the evolution of glioblastomas within a three-dimensional brain, we employ a multi-scale moving boundary model that was initially introduced in [20] and later expanded in several other works [26,27,28,29,30,49]. To account for the brain’s structure, we aim to use 3D T1 weighted and DTI scans that ultimately influence the migration of the cancer cells as well as affect both micro-scale dynamics. Hence, here we aim to explore the impact of the brain structure on the interlinked macro-scale and micro-scale tumour dynamics.

2.1. Macro-Scale Dynamics

Since in this work we extend the 2-dimensional (2D) modelling framework introduced in [20,26], we begin by briefly describing some of the key features of this framework and by giving a few useful notations. First, we denote by $\Omega \left(t\right)$ the expanding 3-dimensional (3D) tumour region that progresses over the time interval $[0,T]$ within a maximal tissue cube $Y\subset {\mathbb{R}}^N$ with $N=3$, i.e., $\Omega \left(t\right)\subset Y, \forall t\in [0,T]$; as can also be seen in Figure 1.

Then, at any macro-scale spatio-temporal point $(x,t)\in Y\times [0,T]$, we consider a cancer cell population $c(x,t)$ that is placed within and interacts with a two-phase ECM: the non-fibre $l(x,t)$ and fibre $F(x,t)$ ECM phases [26,27,28,29,30]. On the one hand, the fibre ECM phase accounts for all major fibrous proteins such as collagen and fibronectin, whose micro-scale distribution induces the spatial orientation of ECM fibres. Hence, the macro-scale spatio-temporal distribution of the ECM fibres is represented by an oriented vector field ${\theta }_f(x,t)$ that describes their spatial bias, as well as by $F(x,t):=\Vert {\theta }_f(x,t)\Vert$, which denotes the amount of fibres at a macro-scale point $(x,t)$ [26,27,28,29,30]. On the other hand, in the non-fibre ECM phase, we bundle together every other ECM constituent such as non-fibrous proteins (for example amyloid fibrils), enzymes, polysaccharides and extracellular $C{a}^{2+}$ ions [26,27,28,29,30]. Furthermore, in this new modelling study, we incorporate the structure of the brain by extracting data from DTI and T1 weighted brain scans, and then use these data to parametrise the model. Specifically, we denote by ${\mathbb{D}}_{Water}\left(x\right)$ the water diffusion tensor that is induced by the DTI scan. Furthermore, we denote by $w\left(x\right)$ the white matter density and by $g\left(x\right)$ the grey matter density, both of which are extracted from the T1 weighted image. Finally, to facilitate the description of the model and to make the mathematical notations more compact, we denote by $\mathbf{u}$ the global tumour vector which at each $(x,t)$ is given by

$\mathbf{u}:={(c(x,t),l(x,t),F(x,t))}^{⊺}.$

With this notation, the total space occupied by the macroscopic tissue is denoted by $\rho \left(\mathbf{u}\right)$ and is defined as

$\rho \left(\mathbf{u}\right)=\rho \left(c\right(x,t),l(x,t),F(x,t\left)\right):=c(x,t)+l(x,t)+F(x,t).$

2.1.1. Cancer Cell Population Dynamics

To describe the spatio-temporal evolution of the cancer cell population $c(x,t)$, we first assume a logistic-type growth with rate $\mu$ [26,27,28,29,30,50,51,52]. For the movement of this cell population, we use the structure of the brain by taking into account both the T1 weighted and DTI scans (from the IXI Dataset [53]), as well as the various adhesion-mediated processes [54,55,56,57,58,59]. Hence, the spatio-temporal dynamics of the cancer cell population is described by

(1)${\displaystyle \frac{\partial c}{\partial t}}=\underset{\begin{array}{c}\mathrm{fully} \mathrm{anisotropic}\\ \mathrm{diffusion}\end{array}}{\underbrace{\nabla \nabla :\left[{\mathbb{D}}_T\left(x\right)c\right]}}-\underset{\mathrm{adhesion} \mathrm{processes}}{\underbrace{\nabla ·\left[c{\mathcal{A}}_c(x,t,\mathbf{u},{\theta }_f)\right]}}+\underset{\begin{array}{c}\mathrm{logistic}-\mathrm{type}\\ \mathrm{proliferation}\end{array}}{\underbrace{\mu c{\left[1-\rho \left(\mathbf{u}\right)\right]}^{+}}}.$

Here, the operator $\nabla \nabla :$ denotes the full second order derivative [35], i.e., it is defined as

$\nabla \nabla :\left[{\mathbb{D}}_T\left(x\right)c\right]:=\sum _{i,j=1}^N{\displaystyle \frac{\partial }{\partial x_i}}{\displaystyle \frac{\partial }{\partial x_j}}\left({\mathbb{D}}_{i,j}c\right),\hspace{1em}N=3,$

with ${\mathbb{D}}_{i,j}$ denoting the components of the tumour diffusion tensor ${\mathbb{D}}_T$. Since classical diffusion models with a constant coefficient cannot capture any directional cues, as those provided by the DTI data, in Equation (1), we use a tensor model (involving a fully anisotropic diffusion term) that is able to incorporate the anisotropic nature of the cancer cell movement. These tensor models were proposed in [60,61,62,63] and have been used to mathematically model the gliomas within the brain; see for instance [35,36,37]. The main idea of this approach is to use the measured water diffusivity in the structured, fibrous brain environment characterised by a symmetric water diffusion tensor [64]:

(2)${\mathbb{D}}_{Water}\left(x\right)=\left[\begin{array}{ccc}{d}_{xx}\left(x\right)& d_{xy}\left(x\right)& d_{xz}\left(x\right)\\ d_{xy}\left(x\right)& d_{yy}\left(x\right)& d_{yz}\left(x\right)\\ d_{xz}\left(x\right)& d_{yz}\left(x\right)& d_{zz}\left(x\right)\end{array}\right],$

and appropriately construct a macroscopic diffusion tensor for the cancer cell population. Since this water tensor (2) is naturally symmetric (due to the DTI scans) [64], it can be diagonalised. Denoting its eigenvalues by ${\lambda }_1\left(x\right)\ge \cdots \ge {\lambda }_N\left(x\right)$ and the associated eigenvectors by ${\varphi }_1\left(x\right),\cdots {\varphi }_N\left(x\right)$, we follow [37,65,66] and construct the 3D tumour diffusion tensor as

(3)$\begin{array}{cc}\hfill {\mathbb{D}}_T\left(x\right):=D_c D_{WG}\left(x\right)[& \left(r+(1-r)\left({\displaystyle \frac{cothk\left(x\right)}{k\left(x\right)}}-{\displaystyle \frac{1}{k{\left(x\right)}^2}}\right)\right)I_3\hfill \\ & +(1-r)\left(1-{\displaystyle \frac{3cothk\left(x\right)}{k\left(x\right)}}+{\displaystyle \frac{3}{k{\left(x\right)}^2}}\right){\varphi }_1\left(x\right){\varphi }_1^{⊺}\left(x\right)].\hfill \end{array}$

Here, $D_c>0$ is the diffusion coefficient, $r\in [0,1]$ specifies the degree of isotropic diffusion, $I_3$ denotes the $3\times 3$ identity matrix, and $k\left(x\right)$ is given by

$k\left(x\right):={\mathcal{K}}_{FA}FA\left(x\right),$

with ${\mathcal{K}}_{FA}\ge 0$ being a proportionality constant measuring the sensitivity of the cells to the environments’ direction, and $FA\left(x\right)$ denotes the fractional anisotropy index [67] given by

$FA\left(x\right):=\sqrt{{\displaystyle \frac{{({\lambda }_1-{\lambda }_2)}^2+{({\lambda }_2-{\lambda }_3)}^2+{({\lambda }_1-{\lambda }_3)}^2}{2({\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2)}}}.$

Finally, it is well known that the malignant glioma cells positioned in the white matter exercise quicker motility than those situated in the grey matter [39,68,69,70]. To account for this effect, in (3), we use a regulator term $D_{WG}(·)$ that is given by

(4)$D_{WG}\left(x\right)=a+(1-a)\left(\left(D_G g\left(x\right)+w\left(x\right)\right)\ast {\psi }_{\rho }\right)\left(x\right),$

where $0\le D_G\le 1$ is the grey matter regulator coefficient, ∗ is the convolution operator [71], ${\psi }_{\rho }:=\psi \left(x\right)/{\rho }^N$ denotes the standard mollifier and $g\left(x\right)$ and $w\left(x\right)$ are the grey and white matter densities provided by the T1 weighted image (following an image segmentation process). Finally, $0\le a\le 1$ is a parameter that distinguishes between different cases (see Section 3).

In addition, the movement of the cancer cells is further biased by various adhesion-mediated process [54,55,56,57,58,59]. Due to the increasing evidence that gliomas induce a fibrous environment within the brain [72,73,74,75,76,77,78,79,80,81,82], in (1), we model the overall adhesion process using a non-local flux term that was introduced in [26] (see also [19,27,28,29,30,49,83,84] for similar terms). Specifically, we explore the adhesive interactions of the cancer cells at $x\in \Omega \left(t\right)$ with other neighbouring cancer cells, with the distribution of the non-fibre ECM phase [85,86,87,88] as well as with the oriented fibre ECM phase [89,90], all located within a sensing region $\mathbf{B}(x,R)$ of radius $R>0$. For this, we define the non-local flux term as

(5)$\begin{array}{cc}\hfill {\mathcal{A}}_c(x,t,\mathbf{u},{\theta }_f):=\frac{1}{R}\underset{\mathbf{B}(0,R)}{\int }\mathcal{K}\left(y\right)[& n\left(y\right)\left({\mathbf{S}}_{cc}c(x+y,t)+{\mathbf{S}}_{cl}l(x+y,t)\right)\hfill \\ & +\widehat{n}(y,{\theta }_f(x+y,t)){\mathbf{S}}_{cF}F(x+y,t)]{\left[1-\rho \left(\mathbf{u}\right)\right]}^{+} dy,\hfill \end{array}$

where ${\mathbf{S}}_{cc}$, ${\mathbf{S}}_{cl}$, ${\mathbf{S}}_{cF}>0$ are the strengths of the cell–cell, cell–non-fibre ECM and cell–fibre ECM adhesions, respectively. While we take both ${\mathbf{S}}_{cl}$ and ${\mathbf{S}}_{cF}$ as positive constants, we consider the emergence of strong and stable cell–cell adhesion bonds to be positively correlated with the level of extracellular $C{a}^{+2}$ ions (one of the non-fibre ECM components) [91,92]. Hence, following the approach in [26,27,28,29,30], we describe the cell–cell adhesion strength by

${\mathbf{S}}_{cc}(x,t):={\mathbf{S}}_{min}+({\mathbf{S}}_{max}-{\mathbf{S}}_{min})exp\left[1-{\displaystyle \frac{1}{1-{(1-l(x,t))}^2}}\right],$

where ${\mathbf{S}}_{min}>0$ and ${\mathbf{S}}_{max}>0$ are the minimum and maximum levels of $C{a}^{+2}$ ions. Furthermore, in (5), we denote by $n(·)$ and $\widehat{n}(·,·)$ the unit radial vector and the unit radial vector biased by the oriented ECM fibres [26,27,28,29,30] defined by

$\begin{array}{cc}\hfill n\left(y\right):=& \left\{\begin{array}{cc}{\displaystyle \frac{y}{{\Vert y\Vert }_2}}\hfill & \mathrm{if} y\in \mathbf{B}(0,R)\backslash \left\{0\right\},\hfill \\ 0\hfill & \mathrm{if} y=0,\hfill \end{array}\right.\hfill \\ \hfill \widehat{n}(y,{\theta }_f(x+y)):=& \left\{\begin{array}{cc}{\displaystyle \frac{y+{\theta }_f(x+y,t)}{\Vert y+{\theta }_f{(x+y,t)\Vert }_2}}\hfill & \mathrm{if} y\in \mathbf{B}(0,R)\backslash \left\{0\right\},\hfill \\ 0\hfill & \mathrm{if} y=0,\hfill \end{array}\right.\hfill \end{array}$

respectively (for details on the fibre orientation ${\theta }_f$, see Section 2.2.1). Furthermore, to account for the gradual weakening of all adhesion bonds as we move away from the centre point x within the sensing region $\mathbf{B}(x,t)$ in (5), we use a radially symmetric kernel $\mathcal{K}(·)$ [29,30] given by

$\mathcal{K}\left(y\right)=\psi \left(\frac{y}{R}\right),\hspace{1em}\forall y\in \mathbf{B}(0,R),$

where $\psi (·)$ is the standard mollifier. Finally, in (5), a limiting term ${[1-\rho \left(\mathbf{u}\right)]}^{+}:=max(0,1-\rho \left(\mathbf{u}\right))$ is used to prevent the contribution of overcrowded regions to cell migration [84]. For a schematic of this adhesion process, we refer the reader to Figure 1a.

2.1.2. Two-Phase ECM Macro-Scale Dynamics

In addition to the cancer cell population, the rest of the macro-scale tumour dynamics are described by the two-phase ECM. Here, both fibre and non-fibre ECM phases are assumed to be simply described by a degradation term due to the cancer cell population. Hence, per unit of time, their dynamics is governed by

(6)$\begin{array}{cc}\hfill {\displaystyle \frac{\partial F}{\partial t}}=& -{\beta }_{F}cF,\hfill \\ \hfill {\displaystyle \frac{\partial l}{\partial t}}=& -{\beta }_{l}cl,\hfill \end{array}$

where ${\beta }_F>0$ and ${\beta }_l>0$ are the degradation rates of the fibre and non-fibre ECM phases, respectively.

2.1.3. The Complete Macro-Dynamics

In summary, Equation (1) for cancer cells dynamics and Equation (6) for the two-phase ECM dynamics lead to the following non-dimensional PDE system describing the evolution of a tumour at macro-scale:

(7)$\begin{array}{cc}\hfill {\displaystyle \frac{\partial c}{\partial t}}=& \nabla \nabla :\left[{\mathbb{D}}_T\left(x\right)c\right]-\nabla ·\left[c{\mathcal{A}}_c(x,t,\mathbf{u},{\theta }_f)\right]+\mu c\left[1-\rho \left(\mathbf{u}\right)\right],\hfill \\ \hfill {\displaystyle \frac{\partial F}{\partial t}}=& -{\beta }_{F}cF,\hfill \\ \hfill {\displaystyle \frac{\partial l}{\partial t}}=& -{\beta }_{l}cl.\hfill \end{array}$

To complete the macro-scale model description, we consider zero-flux boundary conditions and appropriate initial conditions (for instance, the ones given in Section 3).

2.2. Micro-Scale Processes and the Double Feedback Loop

Since the cancer invasion process is genuinely a multi-scale phenomenon, several micro-scale processes are closely linked to the macro-scale dynamics [93]. In this work, we consider two of these micro-processes, namely the rearrangement of the ECM fibres micro-constituents [26] and the cell-scale proteolytic processes that occur at the leading edge of the tumour [20]. Here, we briefly outline these micro-processes, in addition to the naturally arising double feedback loop that ultimately connects the micro-scale and the macro-scale.

2.2.1. Two-Scale Representation and Dynamics of Fibres

To represent the oriented fibres on the macro-scale, we follow [26]. There, the authors characterised not only the amount of fibres $F(x,t)$ but also their ability to withstand incoming cell fluxes and forces through their spatial bias. By considering a cell-scale micro-domain $\delta Y\left(x\right):=x+\delta Y$ of appropriate micro-scale size $\delta >0$, both of these characteristics are induced by the microscopic fibre distribution $f(z,t)$, with $z\in \delta Y\left(x\right)$. In fact, both of these are captured though a vector field representation ${\theta }_f(x,t)$ of the ECM micro-fibres [26] that is defined as

(8)${\theta }_f(x,t):={\displaystyle \frac{1}{\lambda \left(\delta Y\right(x\left)\right)}}\underset{\delta Y\left(x\right)}{\int }f(z,t)dz·{\displaystyle \frac{{\theta }_{f,\delta Y\left(x\right)}(x,t)}{\Vert {\theta }_{f,\delta Y\left(x\right)}{(x,t)\Vert }_2}},$

where $\lambda (·)$ is the Lebesgue measure in ${\mathbb{R}}^d$ and ${\theta }_{f,\delta Y\left(x\right)}(·,·)$ is the revolving barycentral orientation given by [26]

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

Hence, the fibres’ ability to withstand forces is naturally defined by this vector field representation (8) and their amount distributed at a macro-scale point $(x,t)$ is given by

$F(x,t):=\Vert {\theta }_f{(x,t)\Vert }_2,$

which is precisely the mean-value of the micro-fibres distributed on $\delta Y\left(x\right)$. Since both of these macro-scale oriented ECM fibre characteristics ($F(x,t)$ and ${\theta }_f(x,t)$) that we use in the macro-scale dynamics (7), genuinely emerge from the micro-scale distribution of the ECM fibres $f(z,t)$, we refer to this link as the fibres bottom-up link. An illustration of a micro-domain $\delta Y\left(x\right)$ and its corresponding macro-scale orientation ${\theta }_f(x,t)$ can be seen in Figure 1b.

On the other hand, there is also a naturally arising link that connects the macro-scale to this micro-scale, namely the fibres top-down link. This connection is initiated by the movement of the cancer cell population that triggers a rearrangement of ECM fibres micro-constituents on each micro-domain $\delta Y\left(x\right)$ (enabled by the secretion of matrix-degrading enzymes that can break down various ECM proteins). Hence, using the fact that the fully anisotropic diffusion term can be rewritten as $\nabla \nabla :\left[{\mathbb{D}}_T\left(x\right)c\right]=\nabla ·[{\mathbb{D}}_T\left(x\right)\nabla c+c\nabla ·{\mathbb{D}}_T\left(x\right)]$, the fibre rearrangement process is kicked off by the cancer cell spatial flux:

(9)${\mathcal{F}}_c(x,t):={\mathbb{D}}_T\left(x\right)\nabla c+c\nabla ·{\mathbb{D}}_T\left(x\right)-c{\mathcal{A}}_c(x,t,\mathbf{u},{\theta }_f),$

which is generated by the tumour macro-dynamics (7). Then, at any spatio-temporal point $(x,t)\in \Omega \left(t\right)\times [0,T]$ this flux (9) becomes naturally balanced in a weighted fashion by the macro-scale ECM fibre orientation ${\theta }_f(·,·)$, resulting in a rearrangement flux [26]:

(10)$r(\delta Y\left(x\right),t):=\omega (x,t){\mathcal{F}}_c(x,t)+(1-\omega (x,t)){\theta }_f(x,t),$

with $\omega (x,t):=c(x,t)/\left(c\right(x,t)+F(x,t\left)\right)$, that acts uniformly upon the micro-fibre distribution on each micro-domain $\delta Y\left(x\right)$. Ultimately, this macro-scale rearrangement vector (10) induces a micro-scale reallocation vector ${\nu }_{\delta Y\left(x\right)}(z,t)$ [26], enabling us to appropriately calculate the new position $z^{\ast }$ of any micro-node $z\in \delta Y\left(x\right)$ as

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

For further details on the micro-fibre rearrangement process, we refer the reader to Appendix B and [26,27,28,29,30].

2.2.2. MDE Micro-Dynamics and Its Links

The second micro-scale process that we take into consideration is the proteolytic molecular process that occurs along the invasive edge of the tumour and is driven by the cancer cells’ ability to secrete several types of matrix-degrading enzymes (MDEs) (for instance, matrix-metalloproteinases) within the proliferating rim [94,95,96,97,98]. Subsequently to the secretion, these MDEs are subject to spatial transport within a cell-scale neighbourhood of the tumour interface and, as a consequence, they degrade the peritumoral ECM, resulting in changes of tumour boundary morphology [93].

To explore such a micro-scale process, we adopt the approach that was first introduced in [20] where the emerging spatio-temporal MDEs micro-dynamics is considered on an appropriate cell-scale neighbourhood of the tumour boundary $\partial \Omega \left(t\right)$. This neighbourhood is represented by a time-dependent bundle of overlapping cubic micro-domains ${\left\{ϵY\right\}}_{ϵY\in \mathcal{P}\left(t\right)}$, with $ϵ>0$ being the size of each micro-domain $ϵY$, which allows us to decompose the overall MDE micro-process, transpiring on ${\bigcup }_{ϵY\in \mathcal{P}\left(t\right)}ϵY$, into a union of proteolytic micro-dynamics occurring on each $ϵY$; see also Figure 1c. Hence, choosing an arbitrary micro-domain $ϵY$ and a macroscopic time instance $t_0\in [0,T]$, we follow the evolution of the MDE micro-dynamics during the time period $[t_0,t_0+\Delta t]$, with appropriately chosen $\Delta t>0$ and within the associated micro-domain $ϵY$. By denoting the spatio-temporal distribution of MDEs by $m(y,\tau )$ at any micro-point $(y,\tau )\in ϵY\times [0,\Delta t]$, we observe that the cancer cell population, located within an appropriately chosen distance ${\gamma }_h>0$ from $y\in ϵY$, induce a source $h(y,\tau )$ of MDEs which can be mathematically described via a non-local term [20]:

(12)$h(y,\tau )=\left\{\begin{array}{cc}{\displaystyle \frac{\underset{\mathbf{B}(y,{\gamma }_h)\cap \Omega \left(t_0\right)}{\int }c(x,t_0+\tau ) dy}{\lambda (\mathbf{B}(y,{\gamma }_h)\cap \Omega \left(t_0\right))}}\hfill & y\in ϵY\cap \Omega \left(t_0\right),\hfill \\ 0\hfill & y\notin ϵY\backslash (\Omega \left(t_0\right){+\{z\in Y | \parallel z\parallel }_2<\rho \left\}\right),\hfill \end{array}\right.$

where $0<\rho <{\gamma }_h$ is a small mollification range and $\mathbf{B}(y,{\gamma }_h)$ denotes the ${\Vert ·\Vert }_{\infty }$ ball of radius ${\gamma }_h$ centred at a micro-node y. Since the calculation of this micro-scale MDE source (12) directly involves the macro-scale cancer cell population $c(·,·)$, we observe a naturally arising MDE top-down link that connects the macro-scale to the MDE micro-scale. In fact, such a source term (12) allows us to describe the spatio-temporal evolution of the MDEs’ micro-scale distribution $m(·,·)$ by [20]

(13)$\begin{array}{cc}\hfill {\displaystyle \frac{\partial m}{\partial \tau }}& =D_m\Delta m+h(y,\tau ),\hfill \\ \hfill m(y,0)& =0,\hfill \\ \hfill {\displaystyle \frac{\partial m}{\partial n}}{|}_{\partial ϵY}& =0,\hfill \end{array}$

where $D_m>0$ is the constant MDEs diffusion coefficient and n denotes the outward normal vector. As it was shown in [20], we can use the solution of the MDEs micro-dynamics (13) to acquire both the movement direction and magnitude of a tumour boundary point $x_{ϵY}^{\ast }$ within the peritumoral area covered by the associated boundary micro-domain $ϵY$. This ultimately causes a boundary movement, and as a consequence, we obtain a new evolved tumour macro-domain $\Omega (t_0+\Delta t)$, the link of which we refer to as the MDE bottom-up link. For the illustration of the boundary movement, we refer the reader to Figure 1c and for further details of the MDE micro dynamics, see Appendix C or [20,26,27,28,29,30].

3. Computational Results: Numerical Simulations in 3D

We start this section by briefly discussing the numerical method that we use to solve the macro-scale dynamics (7), and for details on the numerical approach used for the two micro-scale dynamics (fibres and MDE), we refer the reader to [29,30]. Here, we use the method of lines’ approach to discretise the macro-scale tumour dynamics (7) first in space, and then, for the resulting system of ODEs, we employ a non-local predictor–corrector scheme [26]. In this context, we carry out the spatial discretisation on a uniform grid, where both spatial operators (fully anisotropic diffusion and adhesion) are accurately approximated in a convolution-driven fashion. First, we note that the fully anisotropic diffusion term can be split into two parts:

(14)$\nabla \nabla :\left[{\mathbb{D}}_T\left(x\right)c\right]=\underset{\mathrm{diffusive}}{\underbrace{\nabla ·\left[{\mathbb{D}}_T\left(x\right)\nabla c\right]}}+\underset{\mathrm{advective}}{\underbrace{\nabla ·\left[c\nabla ·{\mathbb{D}}_T\left(x\right)\right]}},$

which enable us to use a combination of two appropriate distinct schemes for an accurate approximation. While for the diffusive part in (14), we use the symmetric finite difference scheme [99,100], for the combination of the advective (14) and adhesion operators (5) (i.e., $\nabla ·\left[c\left({\mathcal{A}}_c(x,t,\mathbf{u},{\theta }_f)+\nabla ·{\mathbb{D}}_T\left(x\right)\right)\right]$), we use the standard first-order upwind finite difference scheme which ensures positivity and helps avoiding spurious oscillations in the solution. Finally, to approximate the adhesion integral ${\mathcal{A}}_c(x,t,\mathbf{u},{\theta }_f)$, we consider an approach similar to [29,30], and use $N_s$ random points located within the sensing region $\mathbf{B}(0,R)$ and sums of discrete-convolutions.

3.1. Initial Conditions

For the numerical simulations presented in this paper, we consider the tissue cube $Y=[0,4]\times [0,4]\times [0,4]$ with the following initial condition for the cancer cells:

$c(x,0)={\displaystyle \frac{1}{2}}exp\left({\displaystyle \frac{-{\Vert x\Vert }_2^2}{0.02}}\right)·{\chi }_{\mathbf{B}\left(\right(2,2,2),0.25)},$

and for the non-fibre ECM phase, the initial condition $l(x,0)$ was acquired by appropriately scaling the T1 weighted image via a normalising constant. Current DTI scans do not provide suitable resolution to determine the underlying micro-fibre distributions, and so here we describe the initial micro-fibre distribution within a micro-domain $\delta Y\left(x\right)$ as follows. When the macro-scale point x that corresponds to the micro-domain $\delta Y\left(x\right)$ is located in the grey matter, then within $\delta Y\left(x\right)$, we randomly draw straight lines until the ratio between the points that belong and the points that do not belong to the collection of lines is approximately $35\%$:$65\%$. On the other hand, when the point x is located within the white matter, we use a set of predefined lines with the same point ratio ($35\%$:$65\%$), ultimately achieving a random orientation within the grey matter and an aligned orientation within the white matter [101]. Finally, the grey matter’s fibre density is assumed to be $1/D_G$ times smaller than the density in the white matter [101]. A schematics of this initial condition for the micro-fibres can be seen in Figure 2. Hence, we also incorporated the information about the white and grey matter tracks provided by the T1 weighted image into our micro-scale fibre distribution.

3.2. Numerical Simulations in 3D

Here, we present the 3D numerical solutions of the multi-scale model described above, for the parameter values listed in Table A1 in Appendix A (any alteration from these values will be stated accordingly). To display the advanced tumours at time $50\Delta t$, we show four panels for each simulation result. In the first three panels, we show the three classical cross-section planes, i.e., the coronal plane (the head of the subject is viewed from behind), the axial plane (the head of the subject is viewed from above) and the sagittal plane (the head of the subject is viewed from the left). In the last panel of each simulation, we show the 3D image of the brain with the embedded tumour alongside the 3D tumour in isolation.

The three figures shown below investigate tumour evolution when the initial tumour starts in different regions of the brain. To present the simulations, we divide each result into four panels: coronal, axial, sagittal and 3D view. Within each coronal, axial and sagittal views, we show the tumour embedded within the brain on the left, the cancer cell density on the top-right and the ECM density on the bottom-right. In the 3D view (the furthest right panel in each result), we show the cross-section of the whole brain with the tumour in the bottom-left corner and in the top-right corner, we show the isolated tumour.

In Figure 3, we present three distinct cases obtained by varying different parameters that appear in the tumour diffusion tensor ${\mathbb{D}}_T\left(x\right)$ defined in (3). In Figure 3a, we assume that the tensor ${\mathbb{D}}_T\left(x\right)$ depends on the white–grey matter and for this purpose we set $r=1$ in (3) and $a=0$ in (4); this results in isotropic tumour diffusion. In Figure 3b, we use the DTI data (i.e., there is no a priori assumption about the preferential direction for cell movement in white matter) and thus we set $a=1$ in (4) (with $r=0.1$, as in Table A1); this results in an anisotropic diffusion that does not explicitly depend on the white–grey matter. In Figure 3c, we use both DTI data and the white–grey matter dependency (i.e., $r=0.1$ and $a=0$), with the baseline parameters from Table A1. Here, it is worth mentioning that even though we do not use the T1 weighted image to obtain the functions $w\left(x\right)$ and $g\left(x\right)$ that appear in $D_{WG}$ (such as $D_{WG}=1$ in Figure 3b, since $a=1$) we still use the T1 weighted image to initialise the micro-scale non-fibre initial density as well as the initial micro-scale fibre distributions as described above.

In all these simulations shown in Figure 3, we place the small initial tumour in the middle-right part of the brain, and we show the results of the three cases at time $50\Delta t$ where we observe significant tumour morphology changes across the three cases. By comparing Figure 3a,b and Figure 3b,c, we see that when we include the white–grey matter dependency function $D_{WG}$ within the tumour diffusion tensor, it leads to a more advanced tumour. On the other hand, comparing Figure 3a with Figure 3c shows that including the DTI data, which creates an anisotropic tumour diffusion term, leads to a slight reduction in tumour spread. Furthermore, in all three cases, we can notice that the advancing tumour tends to mostly follow the white matter tracks and usually avoids the invasion of tissues located in the grey matter. This invasion resulted in the degradation (and rearrangement) of the ECM that we can see in the bottom-right of each panel (coronal, axial and sagittal) which enabled the tumour to further expand into the surrounding tissues.

In Figure 4, we keep the same three cases as in Figure 3, i.e., Figure 4a only white-grey matter dependency, Figure 4b only DTI data and Figure 4c both. However, here we place the initial tumour in the front-right part of the brain and show the results of tumour invasion at the final time $50\Delta t$. Due to the initial position of the tumour, we can see a tumour that is growing away from the skull towards the centre of the brain as well as it is mainly following the white matter. This creates a highly heterogeneous elongated tumour with many branching outgrowths. On the other hand, in Figure 4 we only see slight differences between the three cases. This contradicts the results from Figure 3 and suggests that both the DTI data and white-grey matter dependency may not always be decisive factor in tumour morphology.

Similarly to Figure 3 and Figure 4, in Figure 5 we keep the same three cases (Figure 5a only white–grey matter dependency Figure 5b only DTI data and Figure 5c both) while we place the initial tumour mass in the middle of the brain and present the results at time $50\Delta t$. As a consequence of the initial location, we see a “butterfly”-shaped tumour that branched to both the left and right side of the brain with some asymmetry. Furthermore, as in Figure 4, we can see that all three cases are quite similar, and so the additional information provided by both the DTI data and white–grey matter dependency seems to be unnecessary for this initial condition. However, we must note that the initial conditions (fibre and non-fibre ECM) still use the information provided by the T1 weighted image, and so here, we only investigate the effect of changing the diffusion tensor.

As we mentioned, we see significant differences between the three cases only in Figure 3. This either indicates that the anisotropic diffusion tensor only provides valuable information in certain cases or that the initial micro-fibre density differs from the one that produced the DTI scan, (i.e., the actual distribution). Since we use an artificial micro-fibre structure that does not depend on the DTI scan which also aids the movement of the cancer cell population via the adhesion integral ${\mathcal{A}}_c(·,·,·,·)$ defined in (5), it is possible that in this specific case, the micro-scale fibre distribution introduced a significantly different travelling direction than the DTI data, resulting in discrepancies between the simulations. However, due to the resolution of current DTI scans, it is not possible to construct a unique fibre distribution within a micro-domain $\delta Y\left(x\right)$. Hence, to genuinely capture the underlying brain structures that we can use within a mathematical model, our results suggest that DTI scans with their present resolution may not be sufficient, and one might need to look into either obtaining better resolution DTI scans or combine this with the strength of different technologies such as magnetic resonance elastography. Nonetheless, this exceeds the current scope of this work and requires further investigation.

4. Discussion and Final Remarks

In this study, we further extended the 2D multi-scale moving-boundary framework previously introduced in [20,26], by developing it into 3D and applying it to the study of glioma invasion within the brain. Since experiments are limited within the brain, we focused on incorporating DTI and T1 weighted scans into our framework to provide insights into the structure of the brain, the tumour, and the surrounding tissue.

The original framework developed in [20,26] modelled a generic tumour in a 2D setting, and so to model gliomas within a 3D brain, we extended this modelling approach by considering the structural information provided by both DTI and T1 weighted scans. We used both DTI and T1 weighted scans to construct the tumour diffusion tensor ${\mathbb{D}}_T\left(x\right)$ defined in (3), which resulted in a fully anisotropic diffusion term. While the T1 weighted image can give different diffusion rates based on whether the cancer cells are located in the white or grey matter, the DTI data are used to incorporate the underlying brain structure and to give higher diffusion rates along specific directions based on how the measured water molecules behaved within the brain. The T1 weighted image, which provided the white–grey matter densities, were also used in our initial conditions for both ECM phases. Hence, the initial density of the non-fibre ECM phase was taken as a normalised version of the T1 weighted image, and the initial condition of the micro-fibre distribution and magnitude were also considered to be dependent on the white–grey matter structure. Furthermore, as the available DTI scans lack the adequate resolution required to construct more appropriate micro-fibre distributions, in this work, we considered a simple case where we set the fibre distributions to be either random or oriented based on whether they are positioned in the grey or white matter, respectively.

Therefore, compared to other previous 3D models such as those in [35,36,37], here we can explicitly capture the underlying tissue structure changes via the micro-fibre rearrangement process and trace how the boundary of the tumour changes over time. Moreover, the cancer cells’ movement is not only influenced by the T1 weighted and DTI scans (through an anisotropic diffusion term), but also biased by crucial non-local adhesions such as cell–cell and cell–fibre ECM adhesion, the latter of which is naturally interconnected to the tissue structural changes. Hence, our model presents a novel approach that can incorporate several vital processes for tumour development and have not been available in other previous models. This allowed us to produce exciting and valuable results/outcomes which provide further insights into the overall tumour development process.

To this end, we used this new 3D model to explore the effects of the anisotropic diffusion term for the cancer cell population. Our numerical simulations in Figure 3 show that including an anisotropic diffusion term may lead to significant changes in the overall tumour morphology. However, it seems that these changes depend on the position of the tumour inside the brain, as Figure 4 and Figure 5 do not exhibit changes consistent with the ones observed between the three sub-panels of Figure 3. This may be the result of the underlying brain structure and its microscopic fibre representation, which seems to take a leading role in influencing cancer-invasion patterns through the underlying cell-adhesion process (see Equation (5)), thus overshadowing the diffusion process. More precisely, the simplified fibre representation might not be sufficient for Figure 3, where the initial tumour was positioned in the right-middle part of the brain. However, this fibre representation might be enough for Figure 4 (with the tumour positioned in the front-right of the brain) and for Figure 5 (with the tumour positioned in the middle of the brain), where we did not see significant morphological differences between the three sub-panels considered in each of these figures.

To conclude this study, we mention that further investigation is needed to determine whether these changes in tumour invasion patterns are caused by the lack of directional information on the fibre micro-scale level or an anisotropic diffusive cell motility is necessary to better represent the invasion process. A feasible approach would also be to use a new imaging technology called magnetic resonance elastography, but this is beyond the scope of this current work. Nonetheless, these results are not only interesting and important from a mathematical modelling point of view but from a medical one as well. For instance, we showed how different model parameters change the evolution of the tumour boundary, this being a vital and desired datum in any clinical/medical decision-making process. Such a model with tuned parameters (to specific patients) could help, for example, decide the optimal amount of tissue to resect so that the chance of survival is maximised.

Finally, as our simulations are able to reproduce known tumour patterns of clinically seen growth, future experiments will be refined by prospectively collected MRI data from glioma patients and also incorporate the effects of their radiotherapy and chemotherapy treatments.

Author Contributions

All authors contributed to this work. All authors have read and agreed to the published version of the manuscript.

Funding

S.S., R.E. and D.T. would like to acknowledge the EPSRC DTA EP/R513192/1 award that funded this research.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

MRI	Magnetic resonance imaging
DTI	Diffusion tensor imaging
ECM	Extracellular matrix
MDE	Matrix degrading enzymes
PDE	Partial differential equation

Appendix A. Parameter Values

In Table A1, we summarise the parameter values that were used in the presented numerical simulations.

Table A1

Parameter set used for the numerical simulations.

Variable	Value	Description	Reference
$D_c$	$1.25\times {10}^{-4}$	Diffusion coeff. for the cancer cell population	[37]
$D_G$	$0.25$	Grey matter regulator coefficient	Estimated
r	$0.1$	Degree of randomised turning	[37]
a	0	Model switching parameter	Estimated
${\mathcal{K}}_{FA}$	100	Cell’s sensitivity to the directional information	[37]
${\mathbf{S}}_{max}$	$0.5$	Cell–cell adhesion coeff.	[26]
${\mathbf{S}}_{min}$	$0.01$	Minimum level of cell–cell adhesion	[29]
${\mathbf{S}}_{cl}$	$0.01$	Cell–non-fibre adhesion coeff.	[26]
${\mathbf{S}}_{cF}$	$0.3$	Cell–fibre adhesion coeff.	[19]
$\mu$	$0.25$	Proliferation coeff. for cancer cell population	[19]
${\beta }_F$	$1.5$	Degradation coeff. of the fibre ECM	[30]
${\beta }_l$	$3.0$	Degradation coeff. of the non-fibre ECM	[30]
$\beta$	$0.8$	Optimal tissue environment controller	[20]
R	$0.15$	Sensing radius	[26]
$f_{max}$	$0.636$	Maximum of micro-fibre density at any point	[26]
$\Delta x$	$0.03125$	Macro-scale spatial step-size	[20]
$ϵ$	$0.0625$	Size of a boundary micro-domain $ϵY\left(x\right)$	[20]
$\delta$	$0.03125$	Size of a fibre micro-domain $\delta Y\left(x\right)$	[26]
$N_s$	450	Number of random points used for the	Estimated
		approximation of the adhesion integral ${\mathcal{A}}_c$

Appendix B. Further Details on the Micro-Fibre Rearrangement Process

In Section 2.2.1, we highlighted the fact that the rearrangement of micro-fibre distribution $f(z,t)$ within each $\delta Y\left(x\right)$ is initiated by the macro-scale cell fluxes, resulting in the redistribution of each micro-fibre pixel z to a new position $z^{\ast }$. To calculate this new position $z^{\ast }$, we used the so-called reallocation vector ${\nu }_{\delta Y\left(x\right)}(z,t)$ which takes into account the rearrangement vector $r\left(\delta Y\right(x),t)$, defined in (10), the degree of alignment between $r\left(\delta Y\right(x),t)$ and the barycentral position vector $x_{dir}:=z-x$ and also incorporates the level of fibres at position z. Hence, following [26], we define it as

${\nu }_{\delta Y\left(x\right)}(z,t):=\left[x_{dir}\left(x\right)+r(\delta Y\left(x\right),t)\right]·{\displaystyle \frac{f(z,t)[f_{max}-f(z,t)]}{f^{\ast }+{\Vert r(\delta Y\left(x\right),t)-x_{dir}\left(x\right)\Vert }_2}}·{\chi }_{\left\{f\right(z,t)>0\}},$

where $f_{max}>0$ is the maximum level of fibres, $f^{\ast }:=f(x,t)/f_{max}$ is the saturation level and ${\chi }_{\left\{f\right(z,t)>0\}}$ is the characteristic function of the micro-fibres support. To move the appropriate amount of fibres from position z to the new position $z^{\ast }$, given in (11), we also monitor the available amount of free space at this target position $z^{\ast }$ via a movement probability $p_{move}$ that we define it as

$p_{move}:=max\left({\displaystyle \frac{f_{max}-f(z^{\ast },t)}{f_{max}}},0\right).$

Consequently, we transport a $p_{move}·f(z,t)$ amount of fibres to the new position $z^{\ast }$ and the rest $(1-p_{move})·f(z,t)$ remains at the original position z.

Appendix C. Further Details on the MDE Micro-Scale

Following [20], we briefly detail here the way the MDE micro-dynamics (13) determines the macro-boundary of the progressed tumour domain $\Omega (t_0+\Delta t)$. To that end, on any arbitrary boundary micro-domain $ϵY\in \mathcal{P}\left(t_0\right)$, we consider an appropriate dyadic cubes decomposition ${\left\{{\mathcal{D}}_k\right\}}_{\mathcal{I}}$, and we denote the barycentre of each ${\mathcal{D}}_k$ by $y_k$. Then, a subfamily of small dyadic cubes ${\left\{{\mathcal{D}}_k\right\}}_{{\mathcal{J}}^{\ast }}$ is sub-sampled by selecting only those dyadic cubes that are the furthest away from the boundary point $x_{ϵY}^{\ast }$ while being located outside of the tumour domain $\Omega \left(t_0\right)$ and carrying an above-average mass of MDEs. This enables us to define the associated direction ${\eta }_{ϵY}$ and displacement magnitude ${\xi }_{ϵY}$ of the movement, which are given by

$\begin{array}{cc}\hfill {\eta }_{ϵY\left(x_{ϵY}^{\ast }\right)}& :=x_{ϵY}^{\ast }+\nu \sum _{l\in {\mathcal{J}}^{\ast }}\left(\underset{{\mathcal{D}}_l}{\int }m(y,\tau ) dy\right)\left(y_l-x_{ϵY}^{\ast }\right),\hspace{1em}\nu \in [0,\infty ),\hfill \\ \hfill {\xi }_{ϵY\left(x_{ϵY}^{\ast }\right)}& :=\sum _{l\in {\mathcal{J}}^{\ast }}{\displaystyle \frac{\underset{{\mathcal{D}}_l}{\int }m(y,\tau ) dy}{{\displaystyle \sum _{l\in {\mathcal{J}}^{\ast }}}\underset{{\mathcal{D}}_l}{\int }m(y,\tau ) dy}}|\overrightarrow{x_{ϵY}^{\ast }{y}_l}|.\hfill \end{array}$

Although a movement direction and a displacement magnitude can thus be determined for each boundary point ${\xi }_{ϵY}$, the actual relocation of ${\xi }_{ϵY}$ only occurs if sufficient but not complete ECM degradation will have occurred in the peritumoral region $ϵY\backslash \Omega \left(t_0\right)$. To quantify the amount of ECM degradation, we use a transitional probability that we define by

$q\left(x_{ϵY}^{\ast }\right):={\displaystyle \frac{\underset{ϵY\left(x_{ϵY}^{\ast }\right)\backslash \Omega \left(t_0\right)}{\int }m(y,\tau )dy}{\underset{ϵY\left(x_{ϵY}^{\ast }\right)}{\int }m(y,\tau )dy}},$

Then, the movement of a boundary point is only exercised when the adequate but not complete degradation of the peritumoral ECM occurs, which is characterised by the situation when this transitional probability $q\left(x_{ϵY}^{\ast }\right)$ exceeds a certain tissue threshold $\omega (·,·)$ (as defined in [20]), namely:

$\omega (\beta ,ϵY):=\left\{\begin{array}{cc}sin\left[{\displaystyle \frac{\pi }{2}}\left(1-{\displaystyle \frac{v(x_{ϵY}^{\ast },t_0+\Delta t)}{\beta ·\underset{\xi \in \partial \Omega \left(t_0\right)}{sup}v(\xi ,t_0+\Delta t)}}\right)\right]\hfill & \mathrm{if} {\displaystyle \frac{v(x_{ϵY}^{\ast },t_0+\Delta t)}{\underset{\xi \in \partial \Omega \left(t_0\right)}{sup}v(\xi ,t_0+\Delta t)}}\le \beta ,\hfill \\ sin\left[{\displaystyle \frac{\pi }{2(1-\beta )}}\left({\displaystyle \frac{v(x_{ϵY}^{\ast },t_0+\Delta t)}{\underset{\xi \in \partial \Omega \left(t_0\right)}{sup}v(\xi ,t_0+\Delta t)}}-\beta \right)\right]\hfill & \mathrm{otherwise},\hfill \end{array}\right.$

where $\beta \in (0,1)$ controls the optimal level of ECM for cancer invasion and $v(x,t):=l(x,t)+F(x,t)$.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Schematics of the multi-scale model: (a) Illustration of the sensing region B(x,R), the two vectors y and y+θf(y+x,t) and the overall travelling direction Ac; (b) an example of a fibre micro-domain δY(x) containing fibres (blue lines) that induces an overall orientation θf(x,t) for δY(x); (c) an example of a boundary micro-domain ϵY(x) where the blue volume represents the tumour volume at the current time-step with boundary point x and the red volume represents the evolved tumour at the next time-step with shifted boundary point xϵY∗.

Figure 2. Schematics of the initial condition of the micro-fibres (blue lines) within a micro-domain δY(x) of orientation θf(x,t) located in the (a) white matter and in the (b) grey matter.

Figure 3. 3D computer simulation results: (a) with only white-grey matter dependency (r=1), (b) with only DTI data used (DG=1), (c) with both white-grey matter dependency and DTI data incorporated. To present the simulations, we divide each result into four panels: coronal, axial, sagittal and 3D view. Within each coronal, axial and sagittal views, we show the tumour embedded within the brain on the left, the cancer cell density on the top-right and the ECM density on the bottom-right. In the 3D view (the furthest right panel in each result) we show the cross-section of the whole brain with the tumour on the bottom-left corner and on the top-right corner we show the isolated tumour.

Figure 4. 3D computer simulation results (a) with only white–grey matter dependency (r=1); (b) with only DTI data used (DG=1); and (c) with both white–grey matter dependency and DTI data incorporated. To represent the simulation results, we use the same format as in Figure 3.

Figure 5. 3D computer simulation results: (a) with only white–grey matter dependency (r=1); (b) with only DTI data used (DG=1); (c) with both white–grey matter dependency and DTI data incorporated. To represent the simulation results, we use the same format as in Figure 3.