Post–operative glioblastoma cancer cell distribution in the peritumoural oedema

Glioblastoma multiforme (GBM), the most aggressive primary brain tumour, exhibits low survival rates due to its rapid growth, infiltrates surrounding brain tissue, and is highly resistant to treatment. One major challenge is oedema infiltration, a fluid build-up that provides a path for cancer cells to invade other areas. MRI resolution is insufficient to detect these infiltrating cells, leading to relapses despite chemotherapy and radiotherapy. In this work, we propose a new multiscale mathematical modelling method, to explore the oedema infiltration and predict tumour relapses. To address tumour relapses, we investigated several possible scenarios for the distribution of remaining GBM cells within the oedema after surgery. Furthermore, in this computational modelling investigation on tumour relapse scenarios were investigated assuming the presence of clinically relevant chemo-radio therapy, numerical results suggest that a higher concentration of GBM cells near the surgical cavity edge led to limited spread and slower progression of tumour relapse. Finally, we explore mathematical and computational avenues for reconstructing relevant shapes for the initial distributions of GBM cells within the oedema from available MRI scans. The results obtained show good overlap between our simulation and the patient’s serial MRI scans taken 881 days into the treatment. While still under analytical investigation, this work paves the way for robust reconstruction of tumour relapses from available clinical data.

1 Introduction

Glioblastoma multiforme (GBM) is a devastating and highly invasive brain tumour that presents a significant treatment challenge. Despite the best efforts of medical professionals, the 5–year survival rate for patients with GBM is only 7.2% (1, 2). To improve treatment outcomes, researchers have been exploring new approaches to tackling this aggressive disease. One promising avenue of investigation is the use of mathematical models to simulate tumour evolution and explore potential new treatment strategies (3–7).

GBM is typically treated with surgical resection if possible, followed by chemoradiotherapy. The Stupp protocol is the standard of care treatment regimen that involves a total of 60 grays (abbreviated as Gy, and representing the unit of measurement for absorbed radiation) of radiotherapy delivered in daily doses of 2 Gy over 6 weeks, along with the chemotherapy drug Temozolomide (TMZ). During radiotherapy, patients take 75 mg of TMZ per square meter of body surface area every day for 7 days a week. After radiotherapy is completed, TMZ (adjuvant) is given in 6 cycles of 150–200 mg per square meter for 5 days every 28 days (8).

Treating GBM is a formidable challenge due to several factors. Even after maximal surgical resection and adherence to the Stupp protocol, approximately 90% of patients experience local recurrences (9–12). Another significant challenge is the high infiltration and heterogeneity of GBM, which makes it difficult to identify tumour margins accurately. GBM grows with microscopic finger–like projections that extend beyond what MRI scans (the gold standard for brain tumour imaging) can detect (2). Furthermore, GBM cells invade the brain through the peritumoural oedema (PTE), a condition in which fluid accumulates in the extracellular spaces of brain tissue surrounding the tumour. PTE is formed by tumour cells, reactive astrocytes, and inflammatory cells. The infiltrating GBM cells in the PTE are phenotypically distinct from those isolated from the corresponding mass. Residual GBM cells located at the resection margin are known to proliferate more quickly and be more invasive than GBM cells found in the tumour center (9, 13). Therefore, it is crucial to examine the PTE as this could lead to tumour recurrences (14).

The limited effectiveness of traditional GBM treatments underscores the need for innovative approaches (4, 15). In recent years, mathematical models have emerged as a promising tool for gaining insights into GBM tumour growth and progression (7, 16–21). By incorporating clinical data and biological parameters, mathematical models can provide a more comprehensive understanding of tumour behaviour than traditional experimental techniques alone (22, 23). However, most of these studies, whether 2D or 3D, are limited to simulating tumour growth on one spatio–temporal scale, i.e., focused on modelling tumour growth based mostly on the macro–scale dynamics (24–28). Nevertheless, significant progress has been made in developing multiscale moving boundary modelling and computational frameworks for tumour growth (3, 6, 7, 29–31). As detailed below, the combination of these modelling approaches paves the way for the work discussed here.

In this work we aim to explore the distribution of GBM cells within the oedema. The underlying motivation for this is the understanding of the relationship between the spatial distribution of cancer cells within oedema that remain post–surgery and the likelihood of post–surgical tumour recurrence. This will combine novel mathematical multiscale moving boundary modelling with NHS clinical data assimilation using MRI scans from a single patient with diagnosed GBM. We explore two scenarios: the first utilizes a standard mollifier to describe cell distribution inside the oedema, while the second uses a Gaussian distribution.

This paper presents a multiscale moving boundary model for simulating GBM evolution, incorporating treatment effects and clinical data. After introducing our multiscale modelling for GBM progression, we formulate our tumour relapse hypothesis and outline the mathematical and computational strategy for clinical data inversion (i.e., assimilate MRI images within our modelling to enable tumour recurrence predictions). Details about prospectively collected MRI scans (from GBM patients at Ninewells Hospital) alongside their pre–processing pipeline are also included. The manual tumour segmentation was carried out under the supervision of consultant neurosurgeons, Dr. Kismet Ibrahim (referred here as KHI) and Mr. Mohamed Okasha (referred here as MO). Finally, we describe the multiscale numerical scheme involved in approximating the mathematical model computationally, and present the simulation results as well as discuss future research avenues.

2 Materials and methods

This section details the mathematical model that we developed to simulate the evolution of GBM within a three–dimensional fibrous brain environment. Our framework expands the work of Suveges et al. (7) by incorporating the effects of various treatment modalities such as surgery, chemotherapy, and radiotherapy. Furthermore, we postulate our hypothesis and formulate a minimisation problem. Finally, we leverage clinical data from T1, T1+C, T2, and DTI scans to account for factors like brain structure, tumour location and extent, and oedema.

2.1 Mathematical multiscale model for GBM progression

2.1.1 Macro–scale dynamics

Following the work from Trucu et al. (3) Suveges et al. (7)Shuttleworth and Trucu (32), we denote by Ω(t) the expanding 3–dimensional (3D) tumour region that progresses over the time interval [0,T] within a maximal tissue cube Y ⊂ℝ³. At any macro–scale spatio–temporal point (x,t) ∈ Y ×[0,T], we consider a cancer cell population, denoted by c(x,t), which interacts with a two–phase heterogeneous ECM (consisting of: a non–fibre l(x,t) and fibre F(x,t) ECM phases (7)), while consuming the available nutrients, denoted by σ(x,t), which are present in the environment. The fibre ECM density, F(x,t), accounts for all fibrous proteins such as collagen and fibronectin. On the other hand, the non–fibre ECM density, l(x,t), comprises of non–fibrous proteins (for example, amyloid fibrils), extracellular Ca²⁺ ions, enzymes and polysaccharides (7). Following the methods introduced in Suveges et al. (7), we also incorporate the structure of the brain by extracting data from the modified DTI scan, T1 and T2 brain scans. Finally, we denote by u(x,t) the global tumour vector which embodies the cancer cell population and the fibre and non–fibre ECM components, given by

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

Therefore, the total space occupied by the macroscopic tissue and tumour volume is denoted by ρ(u) and is defined as

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

for all $(x,t)\in \text{Ω}(t)\times [0,T]$.

2.1.1.1 Nutrients: σ(x, t)

As in this study we focus on avascular tumours, the uptake of nutrients that are available in the outside tissue and are absorbed through the outer tumour boundary plays an important role in the overall tumour development. This nutrients absorption is assumed here to occur at the constant rate d_σ > 0 and is enabled in the model through the presence of nutrient Dirichlet boundary condition at the evolving tumour boundary ∂Ω(t). Furthermore, the spatio-temporal nutrient transport is assumed to be in diffusion equilibrium, with an autonomous transport diffusion coefficient ${\mathcal{D}}_{\sigma }=D_{\sigma }/(c+F+p_0)$ that takes account of both the presence of the cancer and ECM fibres distributions as well as the baseline permeability p₀ > 0 (which is here assumed to be a media constant), while D_σ > 0 is a constant standing for the maximal diffusive nutrients transport possible in the tissue. Thus, the nutrients dynamics is mathematically given by:

$\begin{array}{ll}\begin{array}{rrr}\hfill 0=& \hfill \nabla \cdot ({\mathcal{D}}_{\sigma }\nabla )\sigma -d_{\sigma }c\sigma ,& \text{\hspace{1em}on Ω}(t),\forall t\in [0,T],\hfill \\ \hfill \sigma (x,t)=& {\sigma }_{nor},\hfill & \hfill \forall x\in \partial {\text{Ω}}_0(t),\forall t\in [0,T],\end{array}& (1)\end{array}$

where σ_nor is the normal level of nutrients in the outside tissue and is considered to be constant, while ∂Ω₀(t) represents the outside tumour boundary as defined in Appendix 1 in Supplementary Material. Similar to Suveges et al. (33), certain tumour regions become necrotic as soon as the nutrients level σ drop below a critical necrotic threshold denoted σ_n > 0, while σ_p > 0 represents a nutrient for optimal cancer proliferation regime. Hence, we have the following relationship between these three values: σ_nor > σ_p > σ_n.

Further, considering here a simpler context than the one in Suveges et al. (33) by focussing only on two nutrient effects (namely, on cell proliferation and cell death rates), we assume that: (1) very low nutrient levels impede cell proliferation (having no proliferation at all in the necrotic regions); and (2) extremely high nutrient levels cannot increase cell proliferation rate by more than a certain maximal proliferation rate Ψ_p,max > 0 which corresponds to nutrient levels σ ≥ σ_p. Thus, mathematically, these two assumptions are accounted for in the modelling via the following nutrient–dependent proliferation function:

$\begin{array}{ll}{\text{Ψ}}_p(\sigma ):=\{\begin{array}{ll}0,\hfill & \text{if\hspace{0.17em}}\sigma \le {\sigma }_n,\hfill \\ {\text{Ψ}}_{p,max},\hfill & \text{if\hspace{0.17em}}\sigma \ge {\sigma }_p,\hfill \\ \text{Φ}(\sigma ,{\text{Ψ}}_{p,max},0,{\sigma }_p-{\sigma }_n),\hfill & \text{otherwise},\hfill \end{array}& (2)\end{array}$

where $\text{Φ}(\sigma ,.,.,.)$ describes the smooth transition between the two extrema and is defined to be:

$\begin{array}{ll}\begin{array}{l}\text{Φ}(\sigma ,{\text{Φ}}_{max},{\text{Φ}}_{min},{\text{Φ}}_L):=\\ \frac{{\Phi }_{max}-{\Phi }_{min}}{2}\left[cos\left(\frac{\pi (\sigma -{\sigma }_n-{\Phi }_L)}{{\sigma }_p-{\sigma }_n}\right)+1\right]+{\text{Φ}}_{min},\end{array}& (3)\end{array}$

where Φ_L controls the phase shift of the cosine function.

Finally, the effect that the nutrients absence/presence have on cancer cell death is characterised via a function Ψ_d(σ) that is of similar type as the one given in Equation 2. Specifically, here we consider a maximal death rate Ψ_d,max > 0 in necrotic regions, while we assume no death for cancerous cells when the level of nutrients is σ ≥ σ_p. Thus, using again the transition function from Equation 3, the effect over the death rate of cancer cells is mathematically expressed as:

$\begin{array}{ll}{\text{Ψ}}_d(\sigma ):=\{\begin{array}{ll}{\text{Ψ}}_{d,max},\hfill & \text{if\hspace{0.17em}}\sigma \le {\sigma }_n,\hfill \\ 0,\hfill & \text{if\hspace{0.17em}}\sigma \ge {\sigma }_p,\hfill \\ \text{Φ}(\sigma ,{\text{Ψ}}_{d,max},0,0),\hfill & \text{otherwise}.\hfill \end{array}& (4)\end{array}$

2.1.1.2 Cancer cell dynamics: c(x,t)

The spatio–temporal dynamics of the cancer cell population considered in this work accounts for available movement characteristics enabled by T1 and DTI scans (34), based on which the fully anisotropic diffusion tensor, denoted by ${\mathbb{D}}_T$ (7, 35–37). In addition to that, the cell population movement is further biased by adhesion processes, which are mathematically captured through a term denoted by $\mathcal{A}(x,t,\mathbf{\text{u}},{\theta }_f)$ that will be detailed below. Furthermore, we assume a logistic type proliferation law of the form:

$\begin{array}{ll}P(u):=\mu {\text{Ψ}}_p(\sigma )c{(1-\rho (u))}^{+},& (5)\end{array}$

where µ > 0 is the proliferation rate regulated by the available nutrients, represented here by the nutrient proliferation function Ψ_p(σ) given by Equation 2. Additionally, the term (1 − ρ(u))⁺ guarantees that we do not experience cell population overcrowding within the available space.

Further, while it is well known that one of the hallmarks of cancer is resisting death (38), nevertheless, due to the abnormal peritumoural vasculature and the degradation of the ECM, nutrient delivery is reduced inside the tumour, ultimately leading to necrosis (33). Therefore, we assume a death rate d > 0 that is regulated by the cancer cell death function Ψ_d(σ) given by Equation 4. Thus, mathematically the cancer cell death is captured here by the term:

$\begin{array}{ll}Q(u):=d{\text{Ψ}}_d(\sigma )c.& (6)\end{array}$

Finally, the population of cancer cells is being reduced further by the effects of chemotherapy and radiotherapy, which are cross–referenced with the patient’s post treatment MRI scans. Hence, the spatio–temporal cancer population dynamics is given mathematically by the following partial differential equation:

$\begin{array}{ll}\begin{array}{l}\frac{\partial c}{\partial t}=\underset{\text{Diffusion}}{\underbrace{\nabla \nabla :[{\mathbb{D}}_T(x)c]}}-\underset{\text{Adhesion}\hspace{0.17em}\text{interactions}}{\underbrace{\nabla [c\mathcal{A}(x,t,\mathbf{\text{u}},{\theta }_f)]}}+P(\mathbf{\text{u}})-Q(\mathbf{\text{u}})\\ -Radiotherapy(c,t)-Chemotherapy(c,t).\end{array}& (7)\end{array}$

The first term in Equation 7, $\nabla \nabla :[{\mathbb{D}}_T(x)c]$, denotes the full second order anisotropic tumour diffusion, with the 3D diffusion tensor ${\mathbb{D}}_T$ being constructed from DTI scans of the brain (7, 39). Moreover, the second term in Equation 7, ∇[c𝒜(x,t,u,θ_f)], describes adhesion processes that bias the movement of the cell population due to the adhesion bonds that the migratory cells establish with both the surrounding cell and the ECM components. However, for a compact presentation, we defer the description of both $\nabla \nabla :[{\mathbb{D}}_T(x)c]$ and $\nabla [c\mathcal{A}(x,t,u,{\theta }_f)]$ to Appendix 3 in Supplementary Material, where these terms are explained with full details.

The governing equation also accounts for the effects of radiotherapy and chemotherapy. Radiotherapy is administered in multiple sessions scheduled according to five days a week sequence (Monday to Friday) in equal amounts of doses that is captured here mathematically via a subsequence of days {j_m}_k=1…N_radio ⊂ {1,…,N_final} (where {1,…,N_final} represents the entire period of treatment). The intensity of each radiotherapy fraction follows the linear–quadratic model introduced in Bashkirtseva et al. (40) and is delivered here according to an appropriate per-day radiotherapy distribution function $\overline{r}$: {1,…,N_radio}→ (0,∞), given by $\overline{r}$ (j_m) = αD(j_m) + βD(j_m)², where α > 0 and β > 0 are linear and quadratic coefficients of cell damage, and D(·): {1,…,N_radio}→ (0,∞) is the per-day radiation dose level distribution (i.e., indicating the dose administered in each scheduled day). Finally, we account here also for the time–overlapping effect of radiotherapy treatment over each time interval $(T_{i_k}-l,T_{i_k}+d)$ via the asymmetric mollifier-type function ${\psi }_{j_m}^{radio}(t)$, given in Appendix 5, Equation S6 in Supplementary Material, ∀m ∈{1,…,N_radio}, we have that mathematically the radiotherapy treatment delivery and its effect on the tumour is given by

$\begin{array}{ll}Radiotherapy(c(x,t),t):=\sum _{m=1}^{N_{radio}}\overline{r}(j_m){\psi }_{j_m}^{radio}(t)c(x,t).& (8)\end{array}$

Chemotherapy is incorporated based on the Norton–Simon hypothesis (40), which suggests that tumours are more susceptible to treatment when they have grown for a shorter period of time. Following a chemotherapy scheduling given by a selected subsequence of days {i_k}_k=1…N_chemo ⊂ {1,…,N_final}, we deliver N_chemo doses of chemotherapeutic drug, according to the corresponding per–day chemo agent distribution function ρ_g: {1,…,N_chemo}→ {1,1.1,1.5,2,2.4,2.5,2.8}× chemo_dose, with chemo_dose > 0 being the initial chemo dose. For this specific patient, the TMZ dose on day 1 is chemo_dose = 130 milligrams, but for example, on day 79, the dose increases to 265 milligrams of TMZ, which results into a corresponding dosage upscaling factor of 265/130 ≈ 2. The time–overlapping effect of the chemotherapy over the interval $(T_{i_k}-l,T_{i_k}+d)$ is accounted here via a function ${\psi }_{i_k}^{chemo}(t)$, given in Appendix 5, Equation S6 in Supplementary Material, which is similar in shape to the one for radiotherapy. Further, to account for the fractional cell kill impaired by cytotoxic agent, we adopt an Exponential Kill Model given by b(1−e^ζW), where b > 0 represents the relative maximum fractional cell kill, W > 0 stands for the drug concentration, and ζ > 0 describes tumour cells’ sensitivity to the chemo drug. Moreover, the decrease in fractional cell kill as tumour cell population gets closer to its carrying capacity K > 0 (representing the maximum cumulative distribution of cells and ECM supported by an infinitesimal volume of tissue) is described here through a Holling type II functional µK/(K + sc), where µ > 0 is the growth rate, and s > 0 controls the extent of the Norton–Simon effect, i.e., a larger s leads to a steeper decline, effectively amplifying the Norton-Simon effect by significantly reducing cell kill effectiveness when the tumour is close to its capacity. Conversely, a smaller s results in a more gradual decline, making the Norton-Simon effect less pronounced and allowing for potentially higher cell kill even at larger tumour sizes (40). Thus, chemotherapy delivery and its effect on the tumour is given mathematically by:

$\begin{array}{ll}\begin{array}{l}Chemotherapy(c(x,t),t):=\\ \mu b\frac{K}{(K+sc(x,t))}(1-e^{\zeta W})\sum _{k=1}^{N_{chemo}}{\rho }_g(i_k){\psi }_{i_k}^{chemo}(t)c(x,t)\end{array}& (9)\end{array}$

Thus, the governing equation for cancer dynamics finally becomes

$\begin{array}{ll}\begin{array}{l}\frac{\partial c}{\partial t}=\underset{\text{Diffusion}}{\underbrace{\nabla \nabla :[{\mathbb{D}}_T(x)c]}}-\underset{\text{Adhesion}\hspace{0.17em}\text{interactions}}{\underbrace{\nabla [c\mathcal{A}(x,t,\mathbf{\text{u}},{\theta }_f)]}}+P(\mathbf{\text{u}})-Q(\mathbf{\text{u}})\\ -\underset{\text{Radiotherapy}}{\underbrace{{\displaystyle \sum _{m=1}^{N_{radio}}\overline{r}}(j_m){\psi }_{j_m}^{radio}(t)c}}-\underset{\text{Chemotherapy}}{\underbrace{\mu b\frac{K}{K+sc}(1-e^{\zeta W}){\displaystyle \sum _{k=1}^{N_{chemo}}{\rho }_g(i_k){\psi }_{i_k}^{chemo}(t)c}}}.\end{array}& (10)\end{array}$

2.1.1.3 Two–Phase ECM macro–scale dynamics: F(x,t) and l(x,t)

The micro–scale mass distribution of fibre ECM phase determines a spatial orientation of ECM fibres at micro–scale level which represents their naturally emerging spatial bias for withstanding incoming cell forces (6). With this orientation, while deferring more consistent details for a later subsection, the ECM fibre phase is therefore represented as a macroscopic vector field θ_f(x,t) whose Euclidean norm stands for the amount of fibres at a given macro–scale point (x,t), and so F(x,t): = ∥θ_f(x,t)∥₂ (6, 7). Further, to incorporate the impact of treatment on the each of the two ECM phases, we build on the dynamics of the fibre and non–fibre ECM components introduced in Suveges et al. (7, 33) by considering also the decay effects that the chemo and radio therapies bring about, namely:

$\begin{array}{ll}\frac{\partial F}{\partial t}=-Fc({\beta }_F+{\beta }_{FChemo}+{\beta }_{FRadio}),& (11)\end{array}$

$\begin{array}{ll}\frac{\partial l}{\partial t}=-lc({\beta }_l+{\beta }_{lChemo}+{\beta }_{lRadio}),& (12)\end{array}$

where β_FChemo, β_FRadio and β_lChemo, β_lRadio are the corresponding constant decay rates due to the chemo and radio therapies on the ECM fibres and non-fibres phases, respectively.

2.1.1.4 Summary of the full macro–scale model

In summary, the full model for the macro–scale dynamics is:

$\begin{array}{ll}\{\begin{array}{l}\begin{array}{l}\frac{\partial c}{\partial t}=\underset{\text{Diffusion}}{\underbrace{\nabla \nabla :[{\mathbb{D}}_T(x)c]}}-\underset{\text{Adhesion}\hspace{0.17em}\text{interactions}}{\underbrace{\nabla [c\mathcal{A}(x,t,\mathbf{\text{u}},{\theta }_f)]}}+P(\mathbf{\text{u}})-Q(\mathbf{\text{u}})\\ -\underset{\text{Radiotherapy}}{\underbrace{{\displaystyle \sum _{m=1}^{N_{radio}}\overline{r}}(j_m){\psi }_{j_m}^{radio}(t)c}}-\underset{\text{Chemotherapy}}{\underbrace{\mu b\frac{K}{K+sc}(1-e^{\zeta W}){\displaystyle \sum _{k=1}^{N_{chemo}}{\rho }_g(i_k){\psi }_{i_k}^{chemo}(t)c}}},\end{array}\hfill \\ \frac{\partial F}{\partial t}=-Fc({\beta }_F+{\beta }_{FChemo}+{\beta }_{FRadio}),\hfill \\ \frac{\partial l}{\partial t}=-lc({\beta }_l+{\beta }_{lChemo}+{\beta }_{lRadio}),\hfill \\ 0=\nabla \cdot ({\mathcal{D}}_{\sigma }\nabla )\sigma -d_{\sigma }c\sigma ,\hfill \end{array}& (13)\end{array}$

in the presence of zero–flux boundary conditions for the cancer, fibre and non–fibre ECM phases, as well as, Dirichlet boundary condition for the nutrients.

2.1.2 Micro–scale dynamics within the bulk and at the tumour boundary

In this section, we focus on the micro–scale processes that contribute to cancer invasion. We first discuss the rearrangement of ECM fibres by cancer cells. ECM fibres are important for providing structural support to tissues. Cancer cells can rearrange ECM fibres using matrix–degrading enzymes (MDEs), such as matrix–metalloproteinases, which allows them to create new pathways for invasion. We then discuss the cell–scale proteolytic process at the edge of the tumour, whereby cancer cells secrete MDEs that degrade the ECM, allowing for further tumour invasion. Finally, we discuss the naturally arising double feedback loop that connects the micro–scale and macro–scale. In this loop, the micro–scale interactions between cancer cells and the ECM influence the macro–scale growth and spread of the tumour. The macro–scale growth and spread of the tumour, in turn, influences the micro–scale interactions between cancer cells and the ECM (6, 7, 33).

2.1.2.1 Micro–scale dynamics of ECM fibres and their macro–scale implications

As described in Shuttleworth and Trucu (6); Suveges et al. (7, 33), the macroscopic ECM fibres alongside their ability to withstand incoming forces are represented through the vector field θ_f (x,t) that at each spatio–temporal node (x,t) is non–locally induced from their micro–scale configuration, which is given with full details in Appendix 4 in Supplementary Material. This way, the global macro–scale oriented ECM fibre θ_f(x,t) characteristics (including its Euclidean magnitude which represent the amount of fibres at (x,t), namely F(x,t):= ∥θ_f (x,t)∥₂, arise and are fully determined from the micro–scale distribution of ECM fibres, providing this way a fibres bottom–up micro–to–macro scales link.

However, there exists also a macro–to–micro scales fibres top–bottom link, which is triggered by the movement of cancer cells through the ECM fibre distribution that cause the rearrangement of the ECM micro–fibres on each micro–domain δY(x). Specifically, the fibre rearrangement process is triggered by the macro–scale cancer cell spatial flux

$\begin{array}{ll}ℱ(x,t):={\mathbb{D}}_T(x)\nabla c+c\nabla \cdot {\mathbb{D}}_T(x)-c\mathcal{A}(x,t,\mathbf{\text{u}},{\theta }_f),& (14)\end{array}$

which is balanced by the oriented macro–scale ECM fibre θ_f (x,t), resulting in a rearrangement flux

$\begin{array}{ll}r(\delta Y(x),t):=w(x,t)ℱ(x,t)+(1-w(x,t)){\theta }_f(x,t).& (15)\end{array}$

with w(x,t):= c(x,t)/(c(x,t) + F(x,t)) being an appropriate mediating weight taking into account the amount of cells transported at (x,t) relative to the overall amount of cells and fibres at (x,t). This acts uniformly on the mass distribution of micro–fibre on each micro–domain δY(x), and induces a reallocation of the mass distribution of micro–fibres within both δY(x) and its adjacent neighbouring micro–domains, as described in Shuttleworth and Trucu (6); Suveges et al. (7, 33).

2.1.2.2 MDEs boundary micro–dynamics and its links to the macro–dynamics

Besides the bulk micro–dynamics that involve the ECM fibres, another key micro–dynamics for tumour invasion is the one involving the proteolytic activity that occurs on the invasive edge of the tumour, enabled by the MDEs (secreted by the cancer cells close to the tumour interface) and transported within the surrounding cell–scale peritumoural tissue neighbourhood. Consequently, this MDE micro–dynamics cause degradation of the peritumoural ECM, thereby inducing alterations in the morphological contours of the tumour boundary (7, 33).

This boundary micro–scale MDEs proteolytic activity is explored via the approach initially introduced in Trucu et al. (3), whereby the emergent spatio–temporal dynamics of MDEs on a micro–scale neighbouring envelope B(∂Ω(t),ϵ/2) of cell–scale thickness ϵ > 0, enabled by a bundle 𝒫(t) of overlapping cubic micro–domains $ϵY(z):=B_{\Vert \cdot {\Vert }_{\infty }}(z,ϵ/2)$, $\forall z\in \partial \Omega (t)$, i.e.,

$$𝒫(t):={\left\{ϵY(z)\right\}}_{z\in \partial \Omega (t)}$$

and

$$B(\partial \Omega (t),ϵ/2):={\displaystyle \underset{ϵY\in 𝒫(t)}{\cup }ϵY}$$

with $B_{\Vert \cdot {\Vert }_{\infty }}(z,ϵ/2)$ representing the ∥·∥_∞−ball of radius ϵ/2. This facilitates the decomposition of the overarching MDE micro–process occurring on ${\displaystyle {\cup }_{ϵY\in 𝒫(t)}ϵY}$ into an assembly of proteolytic micro–dynamics occurring on each distinct ϵY. Consequently, at any macroscopic time t₀ ∈ [0,T] during the tumour progression, this decomposing bundle 𝒫(t₀) enable us to explore the MDEs micro–dynamics on each individual micro–domain ϵY ∈𝒫(t₀), where a source of MDEs emerges naturally at micro–scale on the inner cancer side ϵY ∩ Ω(t) as result of collective contributions of the macroscopic distribution of cancer cells that arrives during the macro–dynamics within a close proximity, i.e., within distance γ_h > 0 from ∂Ω(t), which secretes the MDEs. Therefore, mathematically, on a small micro–scale time–length Δt > 0 and at each micro–scale spatio–temporal node (y,τ) ∈ ϵY × [0,Δt], this source of MDEs induced at the micro–scale by the macro–dynamics is expressed through the non–local term:

$\begin{array}{ll}h(y,\tau )=\{\begin{array}{ll}\frac{\mathbf{\text{B}}{\Vert \cdot \Vert }_{\infty } {\stackrel{{\displaystyle \int }}{(y,{\gamma }_h)}\text{∩Ω}(t_0)}^{c(x,t_0+\tau )dy}}{\lambda (\mathbf{\text{B}}(y,{\gamma }_h)\text{∩Ω}(t_0))},\hspace{1em}& y\in \text{ϵ}Y\cap \text{Ω}(t_0),\\ 0,& y\notin ϵY\backslash (\text{Ω}(t_0)+\left\{z\in Y{|\Vert z\Vert }_2<\rho \right\})\end{array}& (16)\end{array}$

where 0 < ρ < γ_h is a small mollification range, B(y,γ_h) represents the ∥·∥_∞ ball of radius γ_h which is centred at a micro–node y ∈ ϵY. Furthermore, in the presence of this source of MDEs on each of the micro–domains ϵY$\in \mathcal{P}$(t₀), the MDEs molecular mass–transport across the tumour interface takes place on each ϵY. Thus, denoting the MDEs density with m(y,τ), ∀(y,τ) ∈ ϵY × [0, Δt], this MDEs transport is assumed here to have a diffusive character and is expressed mathematically as

$\begin{array}{ll}\frac{\partial m}{\partial \tau }=D_m\text{Δ}m+h(y,\tau ),\text{\hspace{1em}\hspace{1em}on }ϵY\times [0,\text{Δ}t],& (17)\end{array}$

with D_m > 0 being a constant diffusion coefficient of the MDEs, while this diffusion process is assumed to take place with: (1) null initial conditions, as this is considered to occur with no molecular memory; and (2) with flux zero boundary conditions as we assume no MDEs molecular transport across the boundary of ∂ϵY.

Finally, as this source is induced and determined directly by the macro–scale cancer cell population c(·,·), this gives rise to a top–down link from the macro–scale to the MDE micro–scale dynamics. On the other hand, as detailed in Trucu et al. (3), the pattern of peritumoural ECM degradation that the MDEs micro–dynamics cause at micro–scale on each boundary micro–domain ϵY∈$\mathcal{P}$(t₀) determines the direction of tumour boundary relocation and enables to characterise this macro–scale movement of the cancer interface through rigorously derived movement laws that specifies precisely at each x ∈ ∂Ω(t₀) the associated relocation direction and magnitude. This ultimately results in a new evolved tumour macro–domain Ω(t₀ + Δt), and this way a bottom–up link is established between the boundary MDEs micro–dynamics and the macro–dynamics.

2.2 Reconstruction of the cancer–cell distribution within the oedema

It has been demonstrated that GBM cells invade the surrounding tissue via the peritumoural oedema which is populated by phenotypically distinct cancer cells that persist in the area following surgical intervention (9). While these cells typically remain untreated or survive the chemoradiotherapy treatment, these are not detectable on MRI scans, and contribute to tumour recurrence. Thus, to gain a better understanding of the tumour relapse process after surgery, it is of interest to explore whether there is any correlation between the shape of the distribution of GBM cells that remain within oedema right after surgery and the extent of the subsequent tumour relapse. Several numerical experiments that we carried out (as those shown in Figures 1, 2) suggest the following hypothesis, namely:

Figure 1

Figure 1. Comparative 3D simulations featuring the mollifier distribution: (A) k_R = 5 with no treatment; (B) k_R = 5 with treatment; (C) k_R = 20 with no treatment; and (D) k_R = 20 with treatment. All simulations captured at the final macro–micro stage 45.

Figure 2

Figure 2. Comparative 3D simulations featuring the Gaussian distribution: (A) $k_{\tilde{\sigma }}=10$ with no treatment; (B) $k_{\tilde{\sigma }}=10$ with treatment; (C) $k_{\tilde{\sigma }}=100$ with no treatment; and (D) $k_{\tilde{\sigma }}=100$ with treatment. All simulations captured at the final macro–micro stage 45.

H: a distribution of GBM cells within the oedema that has most cells mass concentrated within the immediate proximity of the cavity edge leads to a more limited spread and a slower progression of the tumour relapse.

This hypothesis also aligns with clinical findings suggesting that surgical resection removes a substantial number of cancer cells, leaving the remaining cells more dispersed throughout the oedema (41, 42).

In the following, hypothesis H will be examined on two relevant oedema cancer cell distribution types. Furthermore, in both cases, we propose a clinical data assimilation approach, by which we aim to reconstruct the particular shape of the cancer cell distribution that enables the predictive computational modelling solution for the post–surgery GBM relapse to match the available MRI imaging data.

Two possible post–surgery oedema cancer cell distribution scenarios: In the following, we explore hypothesis H by considering two possible scenarios for the post–surgery oedema cancer cell spatial distribution, namely one that is compactly supported strictly within Ω(0) and one that carries non–zero cell mass density distributed at any point in Ω(0). Specifically, denoting by n(x) the usual outward unit normal vector to the surgical cavity edge Γ, ∀x ∈ Γ, we assume that:

on the positive side of the normal direction associated to any x ∈ Γ, represented here parametrically by

$$d_x:=x+\upsilon n(x), \upsilon \ge 0,$$

the shape of immediate post–surgery cancer cell distribution remaining within the oedema along d_x, denoted here by $c_{oedema}^{d_x}$, is of either of the following two types:

case 1: a smooth compact support mollifier–type distribution of support radius $$R(n(x),k_R)$$ centred at x, which is given by

$$c_{oedema}^{d_x}(\upsilon ):=R{(n(x),k_R)}^{-1}{\psi }_1\left(\frac{\upsilon }{R(n(x),k_R)}\right),\upsilon \in [0,q(n(x),0)],$$

where ψ₁(·) is the 1D standard symmetric mollifier given in Appendix 5 in Supplementary Material, while, for any t ≥ 0, q(n(x),t) denotes the distance along line d_x between Γ and ∂Ω(t), with $R(n(x),k_R):=\frac{q(n(x),0)}{k_R}$ while $k_R>1$ represents an uniform scaling constant applied at each x ∈ Γ controls the cancer cells distribution spread in the normal direction described by n(x), see Figure 3A);

Figure 3

Figure 3. Schematic showing from top right to bottom left: a GBM tumour, radio and chemotherapy being applied to the surgically removed area, the surviving cells inside the oedema and finally the two scenarios of cancer cell distributions, (A) the mollifier, and (B) the Gaussian distributions, used in the simulations, where Γ represents the edge of the surgical cavity.

case 2: a Gaussian distribution centred at x and of standard deviation $\tilde{\sigma }(n(x),k_{\tilde{\sigma }}$), which is given by

$$c_{oedema}^{d_x}(\upsilon )\propto {\mathcal{N}}_{d_x}(0,\tilde{\sigma }(n(x),k_{\tilde{\sigma }})),\upsilon \in [0,q(n(x),0)],$$

where by ${\mathcal{N}}_{d_x}(0,\tilde{\sigma }(n(x),k_{\tilde{\sigma }}))$ we denote here the family of normal distributions along d_x, with $\tilde{\sigma }(n(x),k_{\tilde{\sigma }}):=\frac{q(n(x),0)}{k_{\tilde{\sigma }}}$, while $k_{\tilde{\sigma }}>1$ represents an uniform scaling constant applied at each x ∈ Γ controls the standard deviation, see Figure 3B).

For each of the two cases, we explore the correlation between the extent of significant tumour spread within oedema (characterised in case 1 by $R(n(x),k_R)$ and in case 2 by $\tilde{\sigma }(n(x),k_{\tilde{\sigma }})$) and the extent of tumour invasion post–surgery. A smaller $R(n(x),k_R)$ and $\tilde{\sigma }(n(x),k_{\tilde{\sigma }})$ corresponds to a higher concentration of cells near the cavity’s edge, with density decreasing as we move further away from it, as evident in Figure 3 and the upper–right image of Figure 4. Finally, we take advantage of available MRI scans to identify suitable values for $R(n(x),k_R)$ and $\tilde{\sigma }(n(x),k_{\tilde{\sigma }})$ that enable the closest possible match between the computed solutions and the imaging data.

Figure 4

Figure 4. Sagittal views of possible examples of initial conditions when applying them to the pre–surgical oedema: (A) the mollifier and (B) the Gaussian distribution.

Reformulation as least square minimisation problem: In order to assimilate available MRI data to identify appropriate values of parameters controlling the degree of spread of the residual cancer cells within oedema (namely, $R(n(x),k_R)$ and $\tilde{\sigma }(n(x),k_{\tilde{\sigma }})$ for case 1 and case 2, respectively), we proceed by conceptualising this as a minimisation problem. Indeed, to achieve this, to address simultaneously both cases, we consider the mapping Z(n(x),·): (1,∞) → (0,q(n(x),0)) that is defined at each ξ ∈ (1,∞) by

$\begin{array}{ll}Z(n(x),\xi ):=\{\begin{array}{ll}R(n(x),\xi ),\hfill & \text{\hspace{1em}for case }1\hfill \\ \tilde{\sigma }(n(x),\xi ),\hfill & \text{\hspace{1em}for case }2\hfill \end{array}& (18)\end{array}$

with $R(n(x),\xi ):=\frac{q(n(x),0)}{\xi }$ in case 1, and $\tilde{\sigma }(n(x),\xi ):=\frac{q(n(x),0)}{\xi }$ in case 2. In this context we aim to identify the point of minimum ${\xi }_{min}$ (representing the optimal controller parameters ${\overline{k}}_R$ and ${\overline{k}}_{\tilde{\sigma }}$) in case 1 and case 2, respectively) that minimises the following $\xi -$ dependant distance

$\begin{array}{ll}dist(c_{Z(n(x),\xi )},MRI):=\underset{i=1,\dots ,N_{data}}{\max }\Vert c_{Z(n(x),\xi )}(\cdot ,t_i)-MR{I}_i{\Vert }_2,& (19)\end{array}$

where ${\left\{t_i\right\}}_{i=1,\dots ,N_{data}}$ are the macroscopic times at which the corresponding MRI scans ${\left\{MR{I}_i\right\}}_{i=1,\dots ,N_{data}}$ will have been recorded. Moreover, MRI_i is the previously obtained “volume of interest” (VOI). Here, $c_{Z(n(x),\xi )}(\cdot ,t_i)$ represents the spatial density of the computed solution evaluated at t_i that is obtained for a guessed initial condition $c_0^{guess}(\xi ;d_x,\upsilon )$ that corresponds to $\xi \in (1,\infty )$. Finally, for each $\xi \in (1,\infty )$, the guessed initial condition $c_0^{guess}(\xi ;d_x,\upsilon )$ is defined in each of the two cases as:

case 1: $c_0^{guess}(\xi ;d_x,\upsilon ):=R{(n(x),\xi )}^{-1}{\psi }_1(\frac{\upsilon }{R(n(x),\xi )})$, $\upsilon \in [0,q(n(x),0)]$,

case 2: $c_0^{guess}(\xi ;d_x,\upsilon )\propto {\mathcal{N}}_{d_x}(0,\tilde{\sigma }(n(x),\xi ))$, $\upsilon \in [0,q(n(x),0)]$.

2.3 Clinical data assimilation

2.3.1 Acquisition of clinical data

The clinical data used for this study was acquired from one single GBM patient who received different treatments at Ninewells Hospital between 2017 and 2021, chosen due to their prolonged survival, giving us access to multiple MRI scans which can be used to improve our mathematical model. Ethical approval was obtained from the local Caldicott Guardian, Integrated Research Application System (IRAS)(project ID: 309957), Tayside Research and Development Committee (project ID: 2022NH01) and Research Ethics Committee (REC) (Ref: 22/NS/0021). To be included in the study, patients had to be over 16 years old but no older than 85 years old, with histologically confirmed GBM, and have undergone multiple pre–operative and post–operative MRI scans and received standard NHS chemotherapy and radiotherapy treatments. Patients with a limited number of MRI scans were excluded.

2.3.2 Brain imaging, preprocessing and segmentation

The MRI scans were conducted using NHS GE 1.5 Tesla scanners and included multiple pre–operative and post–operative scans for the selected patient. The scans consisted of T1–weighted (T1), T2–weighted (T2), contrast–enhanced T1–weighted (with Gadolinium) (T1+C), diffusion–weighted imaging (DWI) (for specific dates), and T2–FLAIR sequences.

A single typical patient from the series was used for the calculations described here. The patient received initial surgery, followed by chemoradiotherapy with Temozolomide (TMZ) at 130 mg per day concurrently with radiotherapy at a total of 60 Gy distributed equally in 30 total fractions, following the Stupp protocol, and in addition, adjuvant TMZ at a dose of 265–325 mg (6 cycles) after the initial radiotherapy and chemotherapy treatments. Due to recurrence, visible seven months after the completion of concurrent chemoradiotherapy and adjuvant TMZ, the patient also received Lomustine at 160 mg, Procarbazine at 150–200 mg, and Vincristine at 1–2 mg (6 cycles) (this treatment is named PCV). The delivery of the radiotherapy and chemotherapy can be seen in Figure 5, where we also considered all chemotherapeutic drugs functionally equivalent (TMZ, Lomustine, Procarbazine and Vincristine), adjusting only the dosage based on the treatment plan. With the purpose to simulate this treatment delivery, we need to modify the model such that every computational macro–micro stage corresponds to a certain amount of real time. This is later explained in Section 3.1, but clearly showcased in the bottom schematic of Figure 5.

Figure 5

Figure 5. Visualising treatment dynamics: Top image: These graphs depict the radiotherapy and chemotherapy delivery for this patient. The horizontal axis represents the computational stages of the treatment simulation, with every five stages corresponding to one actual day. The vertical axis represents the intensity of the radiation and chemotherapy doses delivered at each stage. Bottom image: Schematic showcasing the full timeline between the first and second surgery from the case considered of the treatments administered (radiotherapy, TMZ, adjuvant TMZ and PCV), the MRI scans dates and the corresponding dates of the comparison between our simulations and the MRI scans. Moreover, we have aligned the temporal progression (weeks and months) of the treatment phases with the corresponding stages of our computational model, providing a clear overview of the timeline and key milestones.

The patient underwent two more surgeries, and MRI scans were taken before and after each of the surgeries, as well as after the completion of the different treatments. When the patient was not undergoing any treatments, MRI scans were conducted every three to four months.

The MRI scans were first pre–processed using Statistical Parametric Mapping (SPM–12, http://www.fil.ion.ucl.ac.uk/spm/). This pre–processing involved reslicing, normalising and finally segmentation of the T1 scan to obtain the white and grey matter densities. As Diffusion Tensor Imaging (DTI) scans were not obtained for this patient, we modified a standard DTI scan from a healthy volunteer from the IXI Dataset (http://brain-development.org/ixi-dataset/), which was warped to match the anatomy of the T1 scan of the GBM patient, hence we were able to infer brain fibre tract directions for the GBM brain.

In Figure 6, on the top right, a T1–weighted scan (T1) is shown and on the top left a T1 scan with gadolinium contrast (T1+C), which outlines the tumour as gadolinium is taken up by the invasive edge of the tumour. The GBM proliferating edge is observed as enhancing in the T1+C and hypo- to iso- intense in T1, as seen in Figure 6. The region which is hyperintense in the T2 scan (and T2–FLAIR) and it is non–enhancing in T1+gadolinium represents the oedema, as seen in the bottom row of Figure 6.

Figure 6

Figure 6. Figure showing GBM and oedema in different axial MRI scans before the first surgery: (A) GBM tumour in a T1+C, (B) GBM tumour in a T1 without contrast, (C) oedema in a T1+C and (D) oedema in a T2 weighted scan.

After pre–processing has been completed, tumour segmentation was done using MRIcroGL, version v1.2.20220720 (www.nitrc.org). The segmentation was performed manually, under the supervision of NHS Consultant neurosurgeons KHI and MO, who specialise in the treatment of GBM. The scans were processed on a axial (transverse) slice–by–slice basis, as seen in Figure 7, for the post–contrast T1 and T2 sequences. These enabled the exploration of important characteristics, referred to as “volumes of interest” (or simply “VOI”), one for the pre–surgical tumour and another one for the oedema before surgery, which were given in the form of binary masks (i.e., individual indicator matrices of zeros and ones that give the footprints of the tumour and oedema) and that were later used in our mathematical model.

Figure 7

Figure 7. Figure showing the pre–surgical axial MRI scans and corresponding volumes of interest for the patient in question: (A) T1+C scan, (B) Pre–surgical GBM VOI, (C) T2 scan and (D) oedema VOI.

3 Results

The numerical approach employed in this work to tackle both the macro–scale and micro–scale dynamics, as well as the top–down and bottom up links between the scales, builds on a sequence of multiscale modelling and computational works introduced in Trucu et al. (3); Shuttleworth and Trucu (6, 32); Suveges et al. (7, 33), and extends these through the introduction of a new governing equation for capturing nutrients dynamics. Moreover, to identify the shape of the remaining post–surgery oedema cancer cell population distribution that lead to GBM relapse, the 3D computational modelling platform developed here is coupled with a least–square–type clinical data assimilation approach using post–surgical MRI scans.

Similar to the methodology outlined in Suveges et al. (33), we utilise the successive over–relaxation method for solving the nutrients Equation 1. For the rest of the macro–scale dynamics in (13), we follow similar steps as in Suveges et al. (7, 33) and employ the method of lines with the following details. Specifically, the spatial operators (i.e., the diffusion and adhesion operators) are addressed as follows: (a) for diffusion we implement a symmetric finite difference scheme based on convolution, as detailed in Suveges et al. (33); and (b) for adhesion we utilise a convolution–driven approach employing a fifth–order weighted essentially non–oscillatory (WENO5) finite difference scheme (43–46), also elaborated upon in Suveges et al. (33). Finally, the time marching is ensured through a predictor corrector scheme introduced in Shuttleworth and Trucu (6) and further detailed in Suveges et al. (7, 33).

3.1 Treatment scheduling

One of the primary aims of our work is to accurately replicate the treatment regimen and dosages administered to a specific patient, which in this case, revolves around the time span bridging the first and second surgery, during which various treatment modalities were employed throughout this entire duration.

The comprehensive timeline for this patient extends to 900 days, encompassing the period between the first and second surgery. During this span, chemotherapy and radiotherapy were administered, and MRI scans were conducted on specific dates, as shown in the bottom schematic of Figure 5. In order to forecast the possibility of relapse and tumour spread based on this patient’s treatment timeline, we need to simulate the treatment process over the course of these 900 days. To achieve this goal, it is important to demarcate the computational macro–micro stages and steps meticulously.

To precisely capture the daily dynamics of the patient’s treatment, the model was modified with a temporal discretisation scheme. Here, every five computational stages represent one actual day, resulting in 4500 stages. This discretisation comes from splitting the macro–scale time interval (i.e., the treatment duration) [0,T_f] into smaller intervals {[kΔt, (k+1)Δt]}_k=0,k_max. Each such increment, which encompasses both the macro–dynamics that takes place on Ω(kΔt) over the time period [kΔt, (k+1)Δt] and the micro–dynamics at its boundary (influenced by the “top–down” links (explained previously in greater detail in section 2.1.2) on each of the boundary micro–domains ϵY ∈ B(∂Ω(kΔt),ϵ/2)) constitutes a “stage k”. As described in Trucu et al. (3)Alzahrani et al. (47), these micro–dynamics at the boundary dictate the precise direction and displacement magnitude for the relocation of each of the points on ∂Ω(kΔt)), progressing this way the stage k tumour domain Ω(kΔt) into the newly obtained domain Ω((k + 1)Δt). With this method, we can match the exact treatment for each day and compare our simulations with the MRI scans taken on those specific dates.

3.2 Initial conditions

The initial micro–fibre distribution within a micro–domain δY (x) is considered here to be the one introduced in Suveges et al. (7), which in brief can be summarised as follows. On one hand, if x ∈ Y is located in the grey matter zone, random straight narrow 3D–stripes (i.e., narrow equal–square cross–section parallelepipedic bars that fit within δY (x)) are distributed until the ratio of the cumulative stripe volume occupied 35% out of the entire δY (x). On the other hand, if x is located in the white matter, a predefined set of aligned straight narrow 3D–stripes is distributed within δY (x) until the volume is filled up to the same percentage, i.e., up to 35%. We also incorporated information about the white and grey matter tracts from the T1+C scan into the micro–scale fibre distribution (7). For the non–fibre ECM phase, we have the following initial condition:

$\begin{array}{ll}l(x,0)=\text{min }\{h(x_1,x_2,x_3),1-c(x,0)\},& (20)\end{array}$

where for any $x:=(x_1,x_2,x_3)\in Y$ we have:

$$h(x)=\frac{1}{2}+\frac{1}{4}\text{sin }{(7\pi y_1(x)y_2(x)y_3(x))}^3\cdot \text{sin }(7\pi y_1(x)/y_2(x)/y_3(x)),$$

with:

$$\begin{array}{lll}{y}_1(x)\hfill & :=\hfill & \frac{1}{3}(x_1+1.5),\hfill \\ y_2(x)\hfill & :=\hfill & \frac{1}{3}(x_2+1.5),\hfill \\ y_3(x)\hfill & :=\hfill & \frac{1}{3}(x_3+1.5).\hfill \end{array}$$

Lastly, the initial condition for the nutrients is set to: σ(x,0) = 0.4.

3.3 Numerical simulations

This section presents the results of 3D numerical simulations of the multiscale model of GBM tumour growth. The parameter values used in the simulations are taken from Table S1 in Supplementary Material. Any modifications made to the values are stated in the text.

To display the evolution of the tumours at time 45Δt, we show four panels for each simulation. The first three panels show the tumour in the coronal, axial, and sagittal planes, respectively. The final panel shows a 3D image of the brain with the embedded tumour alongside the 3D tumour in isolation.

The figures below show the evolution of GBM tumours with different cancer cell distributions in the oedema, under the application (or not) of radiotherapy and chemotherapy. The densities of the main tumour and the ECM are shown in the top–right and bottom–right corners of each of the three classical–views panels, respectively.

Now, to initialise our simulations, we use the segmented masks for both the pre–surgical oedema and tumour, which are subtracted in order to create a surgical cavity, as depicted in Figure 3. Next, we apply either a cancer cell distribution within the modified oedema mask of the shape of a mollifier–type distribution or a Gaussian–type distribution. Moreover, the treatment used on this specific patient is also being applied at the simulation, as shown in Figure 5.

The figures below show the results of applying the mollifier distribution with different values for $k_R$, in Figure 1 and the Gaussian distribution with different values for $k_{\tilde{\sigma }}$, in Figure 2. Finally, we compare the results which showed a reduction in tumour size, as shown in Figure 8.

Figure 8

Figure 8. Simulations that showed the least tumour progression with the mollifier (top) and Gaussian distribution (bottom), respectively.

The results of the simulations are consistent with clinical data, which have shown that the highest concentration of cancer cells in recurrent GBM patients are located at the resection margin (41, 42), hence using the oedema mask, and applying either a mollifier or Gaussian distribution of cancer cells within it, can lead to clinically relevant results with the appropriate values for $k_R$ or $k_{\tilde{\sigma }}$, respectively.

Figure 1 illustrates the results from two experiments. In the first experiment, rows A) and B) used the parameter value of k_R = 5. Row A) depicts the results obtained without applying any treatment, whilst row B) shows the simulation when the treatment from Figure 5 was applied throughout the macro–micro stages.

The second experiment, showcased in rows C) and D), used k_R = 20. Similarly to the first experiment, row C) presents the results without any treatment, while row D) showcases the simulation with the treatment applied. Observe that applying the treatment, in Figure 1 rows B) and D), it highly reduces the densities and spread of the tumour, but there are still residual cancer cells left, mostly around the surgical cavity. Finally, observe that increasing the value of $k_R$ leads to less spread, when comparing the top two rows (A and B) with the bottom two (C and D).

Similarly to the previous case, Figure 2 displays the simulations using the Gaussian distribution with no treatment being applied in rows A) and C), whilst rows B) and D) are the simulations with chemoradiotherapy. Moreover, a value of $k_{\tilde{\sigma }}=10$ is used for both rows A) and B) and $k_{\tilde{\sigma }}$ = 100 for both rows C) and D). When using the Gaussian distribution for the residual cancer cells within oedema after surgery, we observe a similar morphology to the previous case. As with the mollifier distribution experiment, increasing $k_{\tilde{\sigma }}$ and applying the treatment also leads to less tumour growth and spread. Nonetheless, this still leads to a bigger tumour, and with much more spreading potential than in the mollifier case.

Furthermore, we performed experiments with different parameter values and found that the most compact and least invasive tumour spread was obtained when applying the chemoradiotherapy treatment to an initial maximum cancer cells density of 0.1, followed by applying the mollifier distribution to it within the oedema, with k_R = 30. As shown in Figure 8 top row this approach leads to barely any growth, and the tumour remains stable throughout the stages. Moreover, within the same scenario but considering the Gaussian distribution of cells within oedema with $k_{\tilde{\sigma }}$ =100, showcased in Figure 8 bottom row, this also leads to less spread than in the previous experiment from Figure 2, but as showcased in the 3D panel of Figure 8 (bottom row) for the Gaussian distribution simulation, the tumour is larger and spreads more than in the mollifier case from Figure 8 (top row).

Finally, during the course of various experiments, we observed an intriguing outcome. When applying the mollifier distribution to a specific set of values, the resulting outcome closely resembled the extent of the tumour from an MRI scan taken 881 days into the patient’s treatment, as evidenced by a visual comparison between the top–right image of our simulation and an actual MRI scan of the patient, as shown in Figure 9. This discovery guided us toward the subsequent phase of our goal: the comparative analysis of our simulations with MRI scans from this particular patient, enabled by the modification of $k_R$ or $k_{\tilde{\sigma }}$ according to the optimisation procedure from Section 2.2, so that our simulations can closely match the tumour extent observed in the given MRI scans from the considered patient.

Figure 9

Figure 9. Visual comparison between one of our simulations and the MRI scan of the patient, taken 881 days after the first surgery.

3.4 Comparison between our simulations and the MRI data

A key objective of this work is to predict the growth dynamics of GBM tumours whilst incorporating a range of pre–operative and post–operative MRI scans from the specific patient into our analysis. This pursuit, crucial in the field of neuro–oncology, demands a thorough examination of the treatments received by these patients plus the analysis of the MRI scans. This examination is carried out through rigorous comparisons between our computational simulations and the existing MRI scans.

Starting from the baseline parameters set from Table S1 in Supplementary Material, we use the optimisation procedure described in Section 2.2 to approximate the optimal k_R or $k_{\tilde{\sigma }}$ (depending on the case considered), which ensure a good match between the MRI scans and the simulations. Then to properly compare our simulations to MRI scans, we start by aligning the simulation data with the corresponding MRI scan taken at a specific time in our timeline. Comparing the outlined tumour volume from the MRI, outlined under the supervision of KHI and MO, with our simulated cancer density, we calculate the absolute difference following the methods described in Section 2.2, such that Equation 19 is satisfied. If, at any time point, the cancer growth exceeds a set threshold and the disparities between the actual and predicted data are significant, the simulation is halted. Subsequently, k_R or $k_{\tilde{\sigma }}$ are adjusted in a dyadic fashion until the simulation closely matches the real data, meeting our predefined threshold, which we set to be 5000 voxels, i.e., if the difference between our simulation and the VOIs (from MRI scan data) is larger than 5000, then the simulation will stop.

This iterative refinement process ensures that our simulations accurately represent tumour dynamics observed in MRI scans, thereby advancing towards enhancing the reliability and applicability of our computational models.

3.4.1 Utilising post–surgical MRI scans for more realistic tumour simulations

In earlier stages of our research, we focused solely on the initial oedema volume before surgery and the main tumour size before any operation took place, as shown in Figure 3. However, while this approach was methodologically sound, it falls short when attempting to replicate the evolving changes observed in later MRI scans of the patient. These changes occur as the patient’s anatomy undergoes significant transformations due to surgery, as visually depicted in Figure 10. The patient in question experienced a noticeable reduction in the size of the original tumour site after surgical intervention. This reduction was followed by a recurrence of a smaller tumour, as shown in Figure 11A. Consequently, starting our simulations solely based on the initial tumour outline inevitably leads to a tumour size pattern that exceeds our expectations. To address this methodological challenge, we introduced a novel element to our initial conditions: the post–surgical cavity MRI scan VOI. This post–surgical MRI scan provides a clear view of the changes in brain anatomy following surgery, as depicted in Figure 10. By incorporating this post–operative anatomical data into our computational framework, we bridge the crucial gap between the pre–operative and immediate post–operative states. This enables a more precise and anatomically and physiologically realistic simulation of GBM tumour growth dynamics in the context of surgical interventions.

Figure 10

Figure 10. Image showing: (A) axial view of the T1+C pre–operative MRI scan, and (B) axial view of the post–operative MRI scan of the same patient, showcasing the anatomical differences of the brain structure due to the surgical intervention.

Figure 11

Figure 11. (A) Superimposition illustrating the spatial alignment of the pre–surgical original tumour (depicted in blue), the post–surgery surgical cavity (highlighted in green), and the recurrent tumour preceding the second surgical procedure (presented in red). (B) Schematic illustrating the dynamics of the three masks: the pre–surgical oedema mask (in dark yellow), the initial tumour mask (in blue), and the surgical cavity mask (in green). The brighter green region that is not overlapping with the oedema mask, is designated as zero, i.e., no cancer will be located in this area.

The mathematical implementation of this innovative volume addition requires a robust framework. We introduce three distinct masks: the pre–surgical oedema mask, the initial tumour mask, and the surgical cavity mask, as illustrated in Figure 11B. We introduce two new constants which indicate the presence or absence of cancer cells within these masks, represented as α_it and α_sc. Following the notation described in Section 2.2, we employ the mollifier distribution within the different masks to articulate this operation as follows:

● For the initial tumour mask: $c_{it}^{d_x}(\upsilon ):=R{(n(x),k_R)}^{-1}{\psi }_1(\frac{\upsilon }{R(n(x),k_R)}){\alpha }_{it},\upsilon \in [0,q(n(x),0)],$

● For the surgical cavity mask: $c_{sc}^{d_x}(\upsilon ):=R{(n(x),k_R)}^{-1}{\psi }_1(\frac{\upsilon }{R(n(x),k_R)}){\alpha }_{sc},\upsilon \in [0,q(n(x),0)].$

On the other hand, when using a Gaussian distribution, we have the following equations:

● For the initial tumour mask: $c_{it}^{d_x}(\upsilon )\propto {\mathcal{N}}_{d_x}(0,\tilde{\sigma }(n(x),k_{\tilde{\sigma }})){\alpha }_{it},\upsilon \in [0,q(n(x),0)],$

● For the surgical cavity mask: $c_{sc}^{d_x}(\upsilon )\propto {\mathcal{N}}_{d_x}(0,\tilde{\sigma }(n(x),k_{\tilde{\sigma }})){\alpha }_{sc},\upsilon \in [0,q(n(x),0)].$

As shown in the schematic diagram in Figure 11B, it is clear that the surgical cavity is slightly more elongated than the oedema mask. Therefore, regions where these volumes do not overlap are defined by setting their values to zero.

In essence, this mathematical framework equips our computational model with the ability to smoothly incorporate the complex interactions between the oedema, initial tumour, and surgical cavity VOIs. This leads to a more physiologically accurate simulation of GBM tumour growth dynamics, especially in the context of surgical procedures and chemoradiotherapy treatments.

The outcome of the least square minimisation outlined in Section 2.2 led to the identification of the appropriate k_R, that is then used in the mollifier–type distribution of cancer cells in the oedema for the considered patient. With the identified parameter k_R = 20, the numerical simulation is shown in Figure 12.

Figure 12

Figure 12. Simulation with the mollifier distribution using the three introduced masks (the pre–surgical tumour VOI, the pre–surgical oedema VOI, and the surgical cavity VOI), with k_R = 20.

Notably, the tumour in these results is noticeably smaller and has a more compact spatial distribution. Importantly, this tumour is completely surrounded and confined within the boundaries outlined by the surgical cavity mask, which closely matches the patient’s MRI scans, as illustrated in Figure 13A, which represents an overlapping of our simulation (simulation of Figure 12 at stage 44, in green) and the recurred GBM cancer (in red), overlaid over the corresponding MRI scan slice taken 881 days into the treatment (or over 28 months, as depicted in the bottom schematic of Figure 5). The simulation of our model closely aligns with real–world clinical observations for this particular MRI slice. This strong agreement demonstrates the model’s effectiveness and its potential for predicting relevant outcomes. While Figure 13A showcases a high degree of accuracy, it is important to acknowledge that not all MRI slices achieve this level of precision, as shown in Figures 13B, C.

Figure 13

Figure 13. Top row: (A) Overlay depicting the spatial compatibility between our computational simulation (highlighted in green) and the MRI scan demonstrating GBM that recurred prior to the second surgery of the patient (emphasised in red), aligning with the data from the final MRI examination. (B, C) These overlays also involve our computational simulation (in green) and the recurred GBM pre–second–surgery MRI scan (in red). However, in this case, the alignment lacks a high degree of accuracy. Bottom two rows: Time evolution of the tumour recurrence (in red) and the simulation (in green), from month 6 (Comparison 1) into the treatment until month 26 (Comparison 8) from Figure 5. For each of these scans, a comparison was made between the real tumour growth and our simulations.

This achievement marks a significant milestone in our effort to accurately replicate the complexities of GBM tumour growth in the presence of surgical interventions, treatment administrations, and cancer cell distributions within the oedema. This mathematical modelling contributes to our understanding of the clinical management of this very challenging medical condition.

4 Discussion

GBM, an extremely aggressive brain tumour with a low 5–year survival rate of only 7.2%, poses significant challenges in terms of treatment (1, 2). In the search for better therapies, mathematical modelling has emerged as a valuable approach. Despite established treatments like the Stupp protocol, GBM almost always recurs, driven by its invasive nature and peritumoural oedema infiltration, in some cases (9–12). Mathematical models provide a promising way to understand the complexities of GBM.

Our study investigates the connection between the swelling around the tumour (peritumoural oedema) and the distribution of GBM cells within the oedema, whilst using MRI data. Building upon the 3D multiscale moving–boundary framework we introduced earlier, we have incorporated the treatment history of a specific patient from Ninewells Hospital. By simulating how tumours typically grow, our research sets the stage for future experiments using MRI data and treatment histories collected from GBM patients. Ultimately, we aim to develop a mathematical model that incorporates the effects of chemoradiotherapy and investigates the distribution of GBM cells within the oedema with greater accuracy, whilst also taking into account the anatomical changes of the brain due to surgery.

In each simulation, we initiate the process by segmenting the oedema and pre–surgical tumour masks obtained from the MRI scans of the specific GBM patient. Crucially, we meticulously replicate and take account of the exact treatment protocol administered to this patient in our simulations, as showcased in Figure 5. Furthermore, we investigate two scenarios for how cancer cells are distributed within the oedema: the mollifier and Gaussian distributions. The resulting figures show various outcomes based on different values for k_R and $k_{\tilde{\sigma }}$. These simulations closely resemble what doctors see in real clinical cases, where recurrent GBM often has the highest concentration of cancer cells at the edge of the surgically removed area (41, 42, 48, 49). As illustrated in Figures 1, 2, decreasing k_R and $k_{\tilde{\sigma }}$ respectively, corresponds to increased tumour aggressiveness. Notably here, each of the two distributions considered for the cancer cell population within the oedema decays towards the outer boundary of the oedema - an aspect that is consistent with evidence derived from clinical patient histopathology samples (48–50). Moreover, the mollifier distribution does go to zero as we approach the outer oedema boundary, however, the Gaussian does not. This is due to the fact that oedema represents a topologically compact region (subset) in ℝ³ where only the mollifier can die out to zero, while the Gaussian distribution remains always strictly positive (above a minimal threshold level).

Our experiments, involving a range of parameter combinations and the application of the chemoradiotherapy treatment, have shown that the most controlled and least invasive tumour growth occurs when we start with a maximum cancer cell density of 0.1 and use the mollifier distribution, with k_R = 30, to arrange the cancer cells within the oedema, as observed in Figure 8.

What is particularly noteworthy is that our simulations closely match MRI scans taken years into the treatment, as shown in Figure 9, suggesting good agreement between our model and real–world data. This promising finding motivated us to improve our methods to predict GBM growth dynamics by incorporating various pre–and post–operative MRI data along with treatment effects. Initially, we only focused on the initial pre–operative tumour and oedema regions. However, this approach fell short when trying to capture the dynamic changes that occur after surgery, as observed in Figure 10. To address this limitation, we introduced the most immediate and highest quality post–operative MRI data, such that the surgical cavity can be outlined without any complications, into our framework, bridging the gap between the pre–operative and post–operative states. This integration involved three masks (the pre–surgical oedema, the initial tumour, and a future representation of the surgical cavity), as shown in Figure 11B, regulated by constants and applying the mollifier or Gaussian distributions, according to the methods described in Section 2.2. This framework allowed for more precise simulations, as evidenced by the reduction in tumour size and spatial distribution, which closely matched the patient’s MRI scan taken 881 days after the first surgery, as seen in Figure 13A.

In conclusion, our model represents a significant advancement in our ability to predict how GBM tumours behave following surgery, treatment administration and the distributions of cancer cells within the oedema. By incorporating pre–operative and post–operative MRI scans and carefully considering patient treatment histories, we have developed a robust framework that accurately replicates the complex dynamics of GBM progression. This achievement not only enhances our understanding of this challenging disease but also opens up significant possibilities in the field of clinical management.

In terms of modelling challenges and limitations, our model has many unknown variables, which are either estimated from the relevant literature and/or are patient specific. Gathering more data from experimental studies that accurately measure the different parameters involved in GBM growth and invasion can help the model become more accurate. Furthermore, as evidenced in Figure 13, bottom images, the accuracy of the evolution of our simulations, which is being overlaid over the tumour recurrence over time, still needs to be further improved. While this work is meant to be a proof of concept, this model needs validation for its predictive potentials. This is something that we are currently working on, aiming to apply the model to other individuals from a larger patient cohort to explore its level of transferability.

This research reflects current advancements in GBM research by providing valuable insights into mathematical modelling and its potential to predict this aggressive disease. By translating these insights into improved treatments, we hope this work will lead to a significantly improved outlook for GBM patients.

Data availability statement

Permission for access to the data is available via request to the NHS Tayside Caldicott Guardian, and NHS Tayside Ethics and R&D Department. Requests to access the datasets should be directed to NHS Tayside Caldicott Guardian, and NHS Tayside Ethics and R&D Department.

Ethics statement

The studies involving humans were approved by Caldicott Guardian, Integrated Research Application System (IRAS)(project ID: 309957), Tayside Research and Development Committee (project ID: 2022NH01) and Research Ethics Committee (REC) (Ref: 22/NS/0021). The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and institutional requirements.

Author contributions

AM: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. SS: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. MO: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. KHI: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. JSD: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. DT: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. AM, DT, JDS, and KHI would like to acknowledge the generous funding received from the Ninewells Cancer Campaign (NCC) - Fraser Fellowships Doctoral Training Programme (DTP) in Precision Cancer Medicine, which fully funds this project.

Acknowledgments

We would like to acknowledge the generous help provided by Dr. Jennifer MacFarlane in gathering the GBM data.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2024.1447010/full#supplementary-material

References

1. Burri SH, Gondi V, Brown PD, Mehta MP. The evolving role of tumor treating fields in managing glioblastoma: Guide for oncologists. Am J Clin Oncol. (2018) 41:191–6. doi: 10.1097/COC.0000000000000395

PubMed Abstract | Crossref Full Text | Google Scholar

2. Wu W, Klockow JL, Zhang M, Lafortune F, Chang E, Jin L, et al. Glioblastoma multiforme (gbm): An overview of current therapies and mechanisms of resistance. Pharmacol Res. (2021) 171:105780. doi: 10.1016/j.phrs.2021.105780

PubMed Abstract | Crossref Full Text | Google Scholar

3. Trucu D, Lin P, Chaplain MAJ, Wang Y. A multiscale moving boundary model arising in cancer invasion. Multiscale Model Simulation. (2013) 11:309–35. doi: 10.1137/110839011

Crossref Full Text | Google Scholar

4. Michor F, Beal K. Improving cancer treatment via mathematical modeling: Surmounting the challenges is worth the effort. Cell. (2015) 163:1059–63. doi: 10.1016/j.cell.2015.11.002

PubMed Abstract | Crossref Full Text | Google Scholar

5. Malinzi J, Eladdadi A, Sibanda P. Modelling the spatiotemporal dynamics of chemovirotherapy cancer treatment. J Biol Dynamics. (2017) 11:244–74. doi: 10.1080/17513758.2017.1328079

PubMed Abstract | Crossref Full Text | Google Scholar

6. Shuttleworth R, Trucu D. Multiscale modelling of fibres dynamics and cell adhesion within moving boundary cancer invasion. Bull Math Biol. (2019) 81:2176–219. doi: 10.1007/s11538-019-00598-w

PubMed Abstract | Crossref Full Text | Google Scholar

7. Suveges S, Hossain-Ibrahim K, Steele JD, Eftimie R, Trucu D. Mathematical modelling of glioblastomas invasion within the brain: A 3d multi-scale moving-boundary approach. Mathematics. (2021) 9:2214. doi: 10.3390/math9182214

Crossref Full Text | Google Scholar

8. Stupp R, Mason WP, van den Bent MJ, Weller M, Fisher B, Taphoorn MJB, et al. Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. N Engl J Med. (2005) 352:987–96. doi: 10.1056/NEJMoa043330

PubMed Abstract | Crossref Full Text | Google Scholar

9. Lemée J-Me.a. Intratumoral heterogeneity in glioblastoma: don’t forget the peritumoral brain zone. Neuro-oncology. (2015) 17:1322–32. doi: 10.1093/neuonc/nov119

PubMed Abstract | Crossref Full Text | Google Scholar

10. Ringel F, Pape H, Sabel M, Krex D, Bock HC, Misch M, et al. Clinical benefit from resection of recurrent glioblastomas: results of a multicenter study including 503 patients with recurrent glioblastomas undergoing surgical resection. Neuro-Oncology. (2016) 18:96–104. doi: 10.1093/neuonc/nov145

PubMed Abstract | Crossref Full Text | Google Scholar

11. Chen W, Wang Y, Zhao B, Liu P, Liu L, Wang Y, et al. Optimal therapies for recurrent glioblastoma: A bayesian network meta-analysis. Front Oncol. (2021) 11:641878. doi: 10.3389/fonc.2021.641878

PubMed Abstract | Crossref Full Text | Google Scholar

12. Mizuhata M, Takamatsu S, Shibata S, Sakurai T, Minamikawa R, Yamazaki M, et al. Patterns of failure in glioblastoma multiforme following standard (60 gy) or short course (40 gy) radiation and concurrent temozolomide. Japanese J Radiol. (2023) 41:660–8. doi: 10.1007/s11604-023-01386-2

PubMed Abstract | Crossref Full Text | Google Scholar

13. Qin X, Liu R, Akter F, Qin F, Xie Q, Li Y, et al. Peri-tumoral brain edema associated with glioblastoma correlates with tumor recurrence. J Cancer. (2021) 12:2073–82. doi: 10.7150/jca.53198

PubMed Abstract | Crossref Full Text | Google Scholar

14. Niyazi M, Andratschke N, Bendszus M, Chalmers AJ, Erridge SC, Galldiks N, et al. Estro-eano guideline on target delineation and radiotherapy details for glioblastoma. Radiother Oncol. (2023) 184:109663. doi: 10.1016/j.radonc.2023.109663

PubMed Abstract | Crossref Full Text | Google Scholar

15. Yalamarty SSK, Filipczak N, Li X, Subhan MA, Parveen F, Ataide JA, et al. Mechanisms of resistance and current treatment options for glioblastoma multiforme (gbm). Cancers. (2023) 15:2116. doi: 10.3390/cancers15072116

PubMed Abstract | Crossref Full Text | Google Scholar

16. Hatzikirou H, Deutsch A, Schaller C, Simon M, Swanson K. Mathematical modelling of glioblastoma tumour development: a review. Math Models Methods Appl Sci. (2005) 15:1779–94. doi: 10.1142/s0218202505000960

Crossref Full Text | Google Scholar

17. Swanson K. A mathematical modelling tool for predicting survival of individual patients following resection of glioblastoma: a proof of principle. Br J Cancer. (2008) 98:113–9. doi: 10.1038/sj.bjc.6604125

PubMed Abstract | Crossref Full Text | Google Scholar

18. Rockne R, Rockhill JK, Mrugala M, Spence AM, Kalet I, Hendrickson K, et al. Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach. Phys Med Biol. (2010) 55:3271–85. doi: 10.1088/0031-9155/55/12/001

PubMed Abstract | Crossref Full Text | Google Scholar

19. Lipkova J, Angelikopoulos P, Wu S, Alberts E, Wiestler B, Diehl C, et al. Personalized radiotherapy design for glioblastoma: Integrating mathematical tumor models, multimodal scans, and bayesian inference. IEEE Trans Med Imaging. (2019) 38:1875–84. doi: 10.1109/TMI.2019.2902044

PubMed Abstract | Crossref Full Text | Google Scholar

20. Plaszczynski S, Grammaticos B, Pallud J, Campagne J-E, Badoual M. Predicting regrowth of low-grade gliomas after radiotherapy. PloS Comput Biol. (2023) 19:1–16. doi: 10.1371/journal.pcbi.1011002

PubMed Abstract | Crossref Full Text | Google Scholar

21. He L, Zhang H, Li T, Yang J, Zhou Y, Wang J, et al. Distinguishing tumor cell infiltration and vasogenic edema in the peritumoral region of glioblastoma at the voxel level via conventional mri sequences Acad Radiol. (2024) 31:1082–90. doi: 10.1016/j.acra.2023.08.008

PubMed Abstract | Crossref Full Text | Google Scholar

22. Lê M, Delingette H, Kalpathy-Cramer J, Gerstner ER, Batchelor T, Unkelbach J, et al. Personalized radiotherapy planning based on a computational tumor growth model. IEEE Trans Med Imaging. (2017) 36:815–25. doi: 10.1109/TMI.2016.2626443

PubMed Abstract | Crossref Full Text | Google Scholar

23. Yin A, Moes DJAR, van Hasselt JGC, Swen JJ, Guchelaar H-J. A review of mathematical models for tumor dynamics and treatment resistance evolution of solid tumors. CPT Pharmacometrics Syst Pharmacol. (2019) 8:720–37. doi: 10.1002/psp4.12450

PubMed Abstract | Crossref Full Text | Google Scholar

24. Clatz O, Sermesant M, Bondiau P-Y, Delingette H, Warfield S, Malandain G, et al. Realistic simulation of the 3-d growth of brain tumors in mr images coupling diffusion with biomechanical deformation. IEEE Trans Med Imaging. (2005) 24:1334–46. doi: 10.1109/TMI.2005.857217

PubMed Abstract | Crossref Full Text | Google Scholar

25. Suarez C, Maglietti F, Colonna M, Breitburd K, Marshall G. Mathematical modeling of human glioma growth based on brain topological structures: Study of two clinical cases. PLoS ONE. (2012) 7:e39616. doi: 10.1371/journal.pone.0039616

PubMed Abstract | Crossref Full Text | Google Scholar

26. Subramanian S, Gholami A, Biros G. Simulation of glioblastoma growth using a 3d multispecies tumor model with mass effect. J Math Biol. (2019) 79:941–67. doi: 10.1007/s00285-019-01383-y

PubMed Abstract | Crossref Full Text | Google Scholar

27. Hormuth DA, Al Feghali KA, Elliott AM, Yankeelov TE, Chung C. Image-based personalization of computational models for predicting response of high-grade glioma to chemoradiation. Sci Rep. (2021) 11:5820. doi: 10.1038/s41598-021-87887-4

PubMed Abstract | Crossref Full Text | Google Scholar

28. Subramanian S, Ghafouri A, Scheufele KM, Himthani N, Davatzikos C, Biros G. Ensemble inversion for brain tumor growth models with mass effect. IEEE Trans Med Imaging. (2023) 42:982–95. doi: 10.1109/TMI.2022.3221913

PubMed Abstract | Crossref Full Text | Google Scholar

29. Hunt A, Surulescu C. A multiscale modeling approach to glioma invasion with therapy. Vietnam J Mathematics. (2016) 45:221–40. doi: 10.1007/s10013-016-0223-x

Crossref Full Text | Google Scholar

30. Conte M, Dzierma Y, Knobe S, Surulescu C. Mathematical modeling of glioma invasion and therapy approaches via kinetic theory of active particles. Math Models Methods Applied Sci. (2023) 33:1009–51. doi: 10.1142/S0218202523500227

Crossref Full Text | Google Scholar

31. Urcun S, Baroli D, Rohan P-Y, Skalli W, Lubrano V, Bordas SP, et al. Non-operable glioblastoma: Proposition of patient-specific forecasting by image-informed poromechanical model. Brain Multiphysics. (2023) 4:100067. doi: 10.1016/j.brain.2023.100067

Crossref Full Text | Google Scholar

32. Shuttleworth R, Trucu D. Cell-scale degradation of peritumoural extracellular matrix fibre network and its role within tissue-scale cancer invasion. Bull Math Biol. (2020) 82:65. doi: 10.1007/s11538-020-00732-z

PubMed Abstract | Crossref Full Text | Google Scholar

33. Suveges S, Eftimie R, Trucu D. Re-polarisation of macrophages within collective tumour cell migration: A multiscale moving boundary approach. Front Appl Math Stat. (2022) 7:799650. doi: 10.3389/fams.2021.799650

Crossref Full Text | Google Scholar

34. Dataset IXI Dataset. IXI brain imaging dataset. IXI Repository. (2024).

Google Scholar

35. Painter KJ, Hillen T. Mathematical modelling of glioma growth: The use of diffusion tensor imaging (dti) data to predict the anisotropic pathways of cancer invasion. J Theoretical Biol. (2013) 323:25–39. doi: 10.1016/j.jtbi.2013.01.014

PubMed Abstract | Crossref Full Text | Google Scholar

36. Hillen T, Painter J, K. C, Swan A, Murtha AD. Moments of von mises and fisher distributions and applications. Math Biosci and Eng. (2017) 14:673–94. doi: 10.3934/mbe.2017038

PubMed Abstract | Crossref Full Text | Google Scholar

37. Mardia KV, Jupp PE. Directional Statistics. Hoboken, NJ, USA: Wiley (1999). doi: 10.1002/9780470316979

Crossref Full Text | Google Scholar

38. Hanahan D. Hallmarks of cancer: new dimensions. Cancer Discovery. (2022) 12:31–46. doi: 10.1158/2159-8290.CD-21-1059

PubMed Abstract | Crossref Full Text | Google Scholar

39. Engwer C, Hillen T, Knappitsch M, Surulescu C. Glioma follow white matter tracts: a multiscale dti-based model. J Math Biol. (2014) 71:551–82. doi: 10.1007/s00285-014-0822-7

PubMed Abstract | Crossref Full Text | Google Scholar

40. Bashkirtseva I, Ryashko L, López Á.G, Seoane JM, Sanjuán MAF. The effect of time ordering and concurrency in a mathematical model of chemoradiotherapy. Commun Nonlinear Sci Numerical Simulation. (2021) 96:105693. doi: 10.1016/j.cnsns.2021.105693

Crossref Full Text | Google Scholar

41. Oh J, Sahgal A, Sanghera P, Tsao MN, Davey P, Lam K, et al. Glioblastoma: patterns of recurrence and efficacy of salvage treatments. Can J Neurological Sci. (2011) 38:621–5. doi: 10.1017/s0317167100012166

PubMed Abstract | Crossref Full Text | Google Scholar

42. Petrecca K, Guiot M-C, Panet-Raymond V, Souhami L. Failure pattern following complete resection plus radiotherapy and temozolomide is at the resection margin in patients with glioblastoma. J Neuro-Oncol. (2013) 111:19–23. doi: 10.1007/s11060-012-0983-4

PubMed Abstract | Crossref Full Text | Google Scholar

43. Liu X-D, Osher S, Chan T. Weighted essentially non-oscillatory schemes. J Comput Phys. (1994) 115:200–12. doi: 10.1006/jcph.1994.1187

Crossref Full Text | Google Scholar

44. Jiang G-S, Shu C-W. Efficient implementation of weighted eno schemes. J Comput Phys. (1996) 126:202–28. doi: 10.1006/jcph.1996.0130

Crossref Full Text | Google Scholar

45. Zhang S, Shu C-W. A new smoothness indicator for the WENO schemes and its effect on the convergence to steady state solutions. J Sci Comput. (2006) 31:273–305. doi: 10.1007/s10915-006-9111-y

Crossref Full Text | Google Scholar

46. Kim D, Kwon JH. A high-order accurate hybrid scheme using a central flux scheme and a WENO scheme for compressible flowfield analysis. J Comput Phys. (2005) 210:554–83. doi: 10.1016/j.jcp.2005.04.023

Crossref Full Text | Google Scholar

47. Alzahrani T, Eftimie R, Trucu D. Multiscale modelling of cancer response to oncolytic viral therapy. Math Biosci. (2019) 310:76–95. doi: 10.1016/j.mbs.2018.12.018

PubMed Abstract | Crossref Full Text | Google Scholar

48. Wang X, Liu X, Chen Y, Lin G, Mei W, Chen J, et al. Histopathological findings in the peritumoral edema area of human glioma. Histology Histopathology. (2015) 30:1101–9. doi: 10.14670/HH-11-607

PubMed Abstract | Crossref Full Text | Google Scholar

49. Menna G, Marinno S, Valeri F, Mahadevan S, Mattogno PP, Gaudino S, et al. Diffusion tensor imaging in detecting gliomas sub-regions of infiltration, local and remote recurrences: a systematic review. Neurol Rev. (2024) 47:301. doi: 10.1007/s10143-024-02529-3

PubMed Abstract | Crossref Full Text | Google Scholar

50. Yamahara T, Numa Y, Oishi T, Kawaguchi T, Seno T, Asai A, et al. Morphological and flow cytometric analysis of cell infiltration in glioblastoma: a comparison of autopsy brain and neuroimaging. Brain Tumor Pathology. (2010) 27:81–7. doi: 10.1007/s10014-010-0275-7

PubMed Abstract | Crossref Full Text | Google Scholar