Mathematical modelling of glioblastomas invasion within the brain: a 3D multi-scale moving-boundary approach

Brain-related experiments are limited by nature, and so biological insights are often restricted 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 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 explicitly the role of brain micro-fibres in the 3D invading tumours, in this study we extend 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.

the T1 weighted and DTI scans into our multi-scale framework and use the resulting model 48 to simulate numerically the growth of 3D gliomas within the brain. We focus on a few cases 49 showing tumour growth in different regions in the brain, with different distributions of 50 grey/white matter densities, which leads to different tumour invasion patterns. 51 The paper is organised as follows. First, we formulate our extended multi-scale moving 52 boundary framework in Section 2. Then, following a brief description of the numerical 53 methods, we present the computational simulation results in Section 3. Finally, in Section 4 54 we summarise and discuss these results. 55 56 To model the evolution of glioblastomas within a 3-dimensional brain, we employ a 57 multi-scale moving boundary model that was initially introduced in [20] and later expanded 58 in several other works [26][27][28][29][30]49]. To account for the brain's structure, we aim to use 3D T1 59 weighted and DTI scans that ultimately influence the migration of the cancer cells as well 60 as affect both micro-scale dynamics. Hence, here we aim to explore the impact of the brain 61 structure on the interlinked macro-scale and micro-scale tumour dynamics. 62 63 Since in this work we extend the 2-dimensional (2D) modelling framework introduced 64 in [20,26], we begin by describing briefly some of the key features of this framework and 65 by giving a few useful notations. First, we denote by Ω(t) the expanding 3-dimensional 66 (3D) tumour region that progresses over the time interval [0, T] within a maximal tissue 67 cube Y ⊂ R N with N = 3, i.e., Ω(t) ⊂ Y, ∀t ∈ [0, T]; see also Figure 1. Then at any macro- 68 scale spatio-temporal point (x, t) ∈ Y × [0, T] we consider a cancer cell population c(x, t) 69 embedded within a two-phase ECM, consisting of the non-fibre l(x, t) and fibre F(x, t) 70 ECM phases [26][27][28][29][30]. On the one hand, the fibre ECM phase accounts for all major fibrous 71 proteins such as collagen and fibronectin, whose micro-scale distribution induces the spatial 72 orientation of the ECM fibres. Hence, the macro-scale spatio-temporal distribution of the 73 ECM fibres is represented by an oriented vector field θ f (x, t) that describes their spatial 74 bias, as well as by F(x, t) := θ f (x, t) which denote 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 76 every other ECM constituent such as non-fibrous proteins (for example amyloid fibrils), 77 enzymes, polysaccharides and extracellular Ca 2+ ions [26][27][28][29][30]. Furthermore, in this new 78 modelling study we incorporate the structure of the brain by extracting data from DTI and 79 T1 weighted brain scans, and then using this data to parametrise the model. Specifically, 80 we denote by D Water (x) the water diffusion tensor that is induced by the DTI scan. Also, 81 we denote by w(x) the white matter density and by g(x) the grey matter density, both of 82 which are extracted from the T1 weighted image. Finally, for compact writing, we denote 83 by u the global tumour vector described as

Macro-Scale Dynamics
and by ρ(u) the total space occupied that is defined as 2.1.1. Cancer cell population dynamics 86 To describe the spatio-temporal evolution of the cancer cell population c(x, t), we 87 first assume a logistic-type growth with rate µ [26][27][28][29][30][50][51][52]. For the movement of this 88 cell population, we use the structure of the brain by taking into account both the T1 89 weighted and DTI scans (from the IXI Dataset [53]), as well as the various adhesion mediated 90 processes [54][55][56][57][58][59]. Hence, the spatio-temporal dynamics of the cancer cell population is 91 described by Here, the operator ∇∇ : denotes the full second order derivative [35], i.e., it is defined as with D i,j denoting the components of the tumour diffusion tensor D T . Since classical 94 diffusion models with constant coefficient cannot capture any directional cues, as those 95 provided by the DTI data, in Eq. (1) we use a tensor model (involving a fully anisotropic 96 diffusion term) that is able to incorporate the anisotropic nature of the cancer cell movement. 97 These tensor models were proposed in [60][61][62][63] and have been used to mathematically model 98 the gliomas within the brain; see for instance [35][36][37]. The main idea of this approach 99 is to use the measured water diffusivity in the structured, fibrous brain environment 100 characterised by a symmetric water diffusion tensor and appropriately construct a macroscopic diffusion tensor for the cancer cell population.

102
Since this water tensor (2) is assumed to be symmetric, it can be diagonalised. Denoting its 103 eigenvalues by λ 1 (x) ≥ · · · ≥ λ N (x) and the associated eigenvectors by φ 1 (x), · · · φ N (x), 104 we follow [37,64,65] and construct the 3D tumour diffusion tensor as Here D c > 0 is the diffusion coefficient, r ∈ [0, 1] specifies the degree of isotropic diffusion, . Finally, it is well known that the malignant glioma cells positioned in the white matter 110 exercise quicker motility than those situated in the grey matter [39,[67][68][69]. To account for 111 this effect, in (3), we use a regulator term D WG (·) that is given by where 0 ≤ D G ≤ 1 is the grey matter regulator coefficient, * is the convolution operator [70], 113 ψ ρ := ψ(x)/ρ N denotes the standard mollifier and g(x) and w(x) are the grey and white 114 matter densities provided by the T1 weighted image (following an image segmentation 115 process). Finally, 0 ≤ a ≤ 1 is a model switching parameter distinguishes between different 116 cases (see Section 3).

117
In addition, the movement of the cancer cells is further biased by various adhesion 118 mediated process [54][55][56][57][58][59]. Due to the increasing evidence that gliomas induce a fibrous 119 environment within the brain [71][72][73][74][75][76][77][78][79][80][81], in (1) we model the overall adhesion process using 120 a non-local flux term that was introduced in [26] (see also [19,[27][28][29][30]49,82,83] for similar 121 terms). Specifically, we explore the adhesive interactions of the cancer cells at x ∈ Ω(t) with 122 other neighbouring cancer cells, with the distribution of the non-fibre ECM phase [84][85][86][87] 123 as well as with the oriented fibre ECM phase [88,89], all located within a sensing region 124 B(x, R) of radius R > 0. For this, we define the non-local flux term as where S cc , S cl , S cF > 0 are the strengths of the cell-cell, cell-non-fibre ECM and cell-fibre 126 ECM adhesions, respectively. While we take both S cl and S cF as positive constants, we 127 consider the emergence of strong and stable cell-cell adhesion bonds to be positively 128 correlated with the level of extracellular Ca +2 ions (one of the non-fibre ECM component) 129 [90,91]. Hence, following the approach in [26][27][28][29][30], we describe the cell-cell adhesion strength where S min > 0 and S max > 0 are the minimum and maximum levels of Ca +2 ions. Further-132 more, in (5) we denote by n(·) and n(·, ·) the unit radial vector and the unit radial vector 133 biased by the oriented ECM fibres [26][27][28][29][30] respectively (for details on the fibre orientation θ f see Section 2.2.1). Also, to account for the 135 gradual weakening of all adhesion bonds as we move away from the centre point x within 136 the sensing region B(x, t) in (5), we use a radially symmetric kernel K(·) [29,30] given by where ψ(·) is the standard mollifier. Finally, in (5) a limiting term [1 − ρ(u)] + := max(0, 1 − 138 ρ(u)) is used to prevent the contribution of overcrowded regions to cell migration [83]. For 139 a schematic of this adhesion process, we refer the reader to Figure 1(a). In addition to the cancer cell population, the rest of the macro-scale tumour dynamics 142 are described by the two-phase ECM. Here, both fibres and non-fibres ECM phases are 143 assumed to be simply described by a degradation term due to the cancer cell population.

144
Hence, per unit time, their dynamics is governed by where β F > 0 and β l > 0 are the degradation rates of the fibre and non-fibre ECM phases, 146 respectively.
To complete the macro-scale model description, we consider zero-flux boundary conditions 152 and appropriate initial conditions (for instance the ones given in Section 3). 153

154
Since the cancer invasion process is genuinely a multi-scale phenomenon, several 155 micro-scale processes are closely linked to the macro-scale dynamics [92]. In this work, 156 we consider two of these micro-processes, namely the rearrangement of the ECM fibres 157 micro-constituents [26] and the cell-scale proteolytic processes that occur at the leading 158 edge of the tumour [20]. Here we briefly outline these micro-processes, in addition to the 159 naturally arising double feedback loop that ultimately connects the micro-scale and the 160 macro-scale. To represent the oriented fibres on the macro-scale, we follow [26]. There, the authors them are captured though a vector field representation θ f (x, t) of the ECM micro-fibres [26] 168 that is defined as: where λ(·) is the Lebesgue measure in R d and θ f ,δY(x) (·, ·) is the revolving barycentral 170 orientation given by [26] 171 Hence, the fibres' ability to withstand forces is naturally defined by this vector field repre-172 sentation (8) and their amount distributed at a macro-scale point (x, t) is given by which is precisely the mean-value of the micro-fibres distributed on δY(x). Since both of 174 these macro-scale oriented ECM fibre characteristics (F(x, t) and θ f (x, t)) that we use in 175 the macro-scale dynamics (7) , the fibre rearrangement process is kicked off by the cancer cell spatial flux which is generated by the tumour macro-dynamics (7). Then, at any spatio-temporal point (9) gets naturally balanced in a weighted fashion by the 188 macro-scale ECM fibre orientation θ f (·, ·), resulting in a rearrangement flux [26] 189 with ω(x, t) := c(x, t)/(c(x, t) + F(x, t)), that acts uniformly upon the micro-fibre distri- For further details on the micro-fibre rearrangement process, we refer the reader to  pendix B and [26][27][28][29][30]. The second micro-scale process that we take into consideration is the proteolytic 197 molecular process that occurs along the invasive edge of the tumour and is driven by the 198 cancer cells' ability to secrete several types of matrix-degrading enzymes (MDEs) (for instance, 199 matrix-metalloproteinases) within the proliferating rim [93][94][95][96][97]. Subsequent to the secretion, 200 these MDEs are subject to spatial transport within a cell-scale neighbourhood of the tumour 201 interface and, as a consequence, they degrade the peritumoral ECM, resulting in changes of 202 tumour boundary morphology [92].

203
To explore such a micro-scale process, we adopt the approach that was first introduced where 0 < ρ < γ h is a small mollification range and B(y, γ h ) denotes the · ∞ ball of 217 radius γ h centred at a micro-node y. Since the calculation of this micro-scale MDE source where D m > 0 is the constant MDEs diffusion coefficient and n denotes the outward normal

231
We start this section by briefly discussing the numerical method that we use to solve 232 the macro-scale dynamics (7), and for details on the numerical approach used for the two 233 micro-scale dynamics (fibres and MDE) we refer the reader to [29,30]. Here, we use the 234 method of lines approach to discretise the macro-scale tumour dynamics (7) first in space, 235 and then, for the resulting system of ODEs, we employ a non-local predictor-corrector 236 scheme [26]. In this context, we carry out the spatial discretisation on a uniform grid, 237 where both spatial operators (fully anisotropic diffusion and adhesion) are accurately 238 approximated in a convolution-driven fashion. First, we note that the fully anisotropic 239 diffusion term can be split into two parts which enables us to use a combination of two appropriate distinct schemes for an accurate 241 approximation. While for the diffusive part in (14), we use the symmetric finite difference 242 scheme [98,99], for the combination of the advective (14) and adhesion operators (5) (i.e.,

243
∇ · c A c (x, t, u, θ f ) + ∇ · D T (x) ) we use the standard first-order upwind finite difference 244 scheme which ensures positivity and helps avoiding spurious oscillations in the solution.

245
Finally, to approximate the adhesion integral A c (x, t, u, θ f ), we consider an approach similar 246 to [29,30], and use N s random points located within the sensing region B(0, R) and sums of 247 discrete-convolutions. orientation within the grey matter and an aligned orientation within the white matter [100].

261
Finally, the grey matter's fibre density is assumed to be 1/D G times smaller than the density 262 in the white matter [100]. A schematics of this initial condition for the micro-fibres can be 263 seen in Figure 2. Hence, we also incorporate the information about the white and grey 264 matter tracks provided by the T1 weighted image into our micro-scale fibre distribution.

266
Here, we present the 3D numerical solutions of the multi-scale model described above, 267 for the parameter values listed in Table 1 in Appendix A (any alteration from these values 268 will be stated accordingly). To display the advanced tumours at time 50∆t, we show four 269 panels for each simulation results. In the first three panels we show the three classical  In Figure 3 we present three distinct cases obtained by varying different parameters 277 that appear in the tumour diffusion tensor D T (x) defined in (3). In Figure 3 (a) we assume 278 that the tensor D T (x) depends on the white-grey matter and for that purpose we set r = 1  in (3) and a = 0 in (4); this results in isotropic tumour diffusion. In Figure 3 (b) we use 280 the DTI data (i.e., there is no a-priory assumption about the preferential direction for cell 281 movement in white matter) and thus we set a = 1 in (4) (with r = 0.1, as in Table 1); this 282 results in an anisotropic diffusion that does not depend explicitly on the white-grey matter. Similarly to Figure 3 and Figure 4, in Figure 5 we keep the same three cases ( Figure 5 312 (a) only white-grey matter dependency, Figure 5 (b) only DTI data and Figure 5 (c) both) 313 while we place the initial tumour mass in the middle of the brain and present the results 314 at time 50∆t. As a consequence of the initial location, we see a "butterfly" shaped tumour 315 that branched to both the left and right side of the brain with some asymmetry. Also, as in 316 Figure 4 we can see that all three cases are quite similar, and so the additional information 317 provided by both the DTI data and white-grey matter dependency seems to be unnecessary 318 for this initial condition. However, we must note that the initial conditions (fibre and 319 non-fibre ECM) still uses the information provided by the T1 weighted image, and so here, 320 we only investigate the effect of changing the diffusion tensor.

321
As we mentioned, we see significant differences between the three cases only in Figure   322 3. This either indicates that the anisotropic diffusion tensor provides valuable information 323 only in certain cases or that the initial micro-fibre density differs from the one that produced 324 the DTI scan (i.e., the actual distribution). Since we use an artificial micro-fibre structure that 325 does not depend on the DTI scan which also aid the movement of the cancer cell population

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

342
The original framework developed in [20,26] modelled a generic tumour in a 2D 343 setting, and so to model gliomas within a 3D brain, we extended this modelling approach 344 by considering the structural information provided by both DTI and T1 weighted scans. 345 We used both DTI and T1 weighted scans to construct the tumour diffusion tensor D T (x) 346 defined in (3), which resulted in a fully anisotropic diffusion term. While the T1 weighted 347 image can give different diffusion rates based on whether the cancer cells are located in 348 the white or grey matter, the DTI data is used to incorporate the underlying brain structure 349 and to give higher diffusion rates along specific directions based on how the measured 350 water molecules behaved within the brain. The T1 weighted image, which provided the 351 white-grey matter densities, were also used in our initial conditions for both ECM phases.

352
Hence, the initial density of the non-fibre ECM phase was taken as a normalised version 353 of the T1 weighted image, and the initial condition of the micro-fibre distribution and 354 magnitude were also considered to be dependent on the white-grey matter structure.

355
Furthermore, as the available DTI scans lack the adequate resolution required to construct 356 more appropriate micro-fibre distributions, in this work we considered a simple case where 357 we set the fibre distributions to be either random or oriented based on whether they are 358 positioned in the grey or white matter, respectively.

359
We used this new 3D model to explore the effects of the anisotropic diffusion term for 360 the cancer cell population. Our numerical simulations in Figure 3 showed that including an 361 anisotropic diffusion term may lead to significant changes in the overall tumour morphology.

362
However, it seems that these changes depend on the position of the tumour inside the brain,    Consequently, we transport p move · f (z, t) amount of fibres to the new position z * and the 407 rest (1 − p move ) · f (z, t) remains at the original position z.