Simulating the behaviour of glioblastoma multiforme based on patient MRI during treatments

Glioblastoma multiforme is a brain cancer that still shows poor prognosis for patients despite the active research for new treatments. In this work, the goal is to model and simulate the evolution of tumour associated angiogenesis and the therapeutic response to glioblastoma multiforme. Multiple phenomena are modelled in order to fit different biological pathways, such as the cellular cycle, apoptosis, hypoxia or angiogenesis. This leads to a nonlinear system with 4 equations and 4 unknowns: the density of tumour cells, the O2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {O}_{2}$$\end{document} concentration, the density of endothelial cells and the vascular endothelial growth factor concentration. This system is solved numerically on a mesh fitting the geometry of the brain and the tumour of a patient based on a 2D slice of MRI. We show that our numerical scheme is positive, and we give the energy estimates on the discrete solution to ensure its existence. The numerical scheme uses nonlinear control volume finite elements in space and is implicit in time. Numerical simulations have been done using the different standard treatments: surgery, chemotherapy and radiotherapy, in order to conform to the behaviour of a tumour in response to treatments according to empirical clinical knowledge. We find that our theoretical model exhibits realistic behaviours.


Introduction
Glioblastoma multiforme (GBM) is the deadliest and most frequent brain tumour. Despite the research of new treatments, patients still show poor prognosis in the long run: only 5% of patients survive 5-year post-prognosis.
Usually, patients undergo emergency surgery (if the tumour is reachable), then the treatment consists in radiotherapy plus concomitant and adjuvant Temozolomide (TMZ) therapy (Stupp et al. 2005). More efficient therapies remain a major preoccupation to cure GBM. Among them, immunotherapy is more and more a subject of research for gliomas (Lim et al. 2018;Kamran et al. 2018) and could improve the current prognosis for GBM patients.
Mathematics has been used for developing models matching the behaviour of glioma tumour cells in recent years. Some models use a spherical tumour growth approach using Partial Differential Equations (PDEs) (Papadogiorgaki et al. 2013;Stein et al. 2007;). Other models approach it using an elasticity (Subramanian et al. 2019) or using evolutionary game theoretical model (Basanta et al. 2011).
When a patient gets diagnosed with GBM, tumour cells have already achieved enough tumour promotion mechanisms in order to evade the immune system and to proliferate in the brain. In that sense, we chose to model the GBM growth based on the process of tumour associated angiogenesis.
Angiogenesis is the ensemble of phenomena that allow the formation of new blood vessels from pre-existing blood vessels. Those physiological processes happen not only for cancer patients, but tumours have the ability to use angiogenesis in their favour as a tumour promoter (Kim and Lee 2009). A simplification of the processes used by tumour cells to induce angiogenesis is proposed in Fig. 1. For their growth, tumour cells rely on nutrients and O 2 provided by blood vessels. During tumour growth, the tumour core lacks O 2 inducing hypoxia in the tumour core. Hypoxia prevents most tumour cellular activities, acting like a tumour suppressor process. To fight hypoxia, hypoxic tumour cells produce proangiogenic factors such as Vascular Endothelial Growth Factors (VEGF) that are the main factors produced in GBM. Proangiogenic factors promote angiogenesis, meaning that more blood vessels are produced, and so more nutrients and O 2 are provided to the tumour cells. Angiogenesis mathematical models have already been developed: using PDEs (Vilanova et al. 2017;Mantzaris et al. 2004;Schugart et al. 2008;Colin et al. 2015), some adding stochastic parts in the modelling (Travasso et al. 2011), or working at a mesoscopic scale (Spill et al. 2015).
However in this work, we consider an enhanced tumour associated angiogenesis model by working on Magnetic Resonance Imaging (MRIs) data based on a real patient and by modelling the behaviour of GBM growth through the treatments usually administered to patients. Indeed, MRIs are required to certify the diagnosis of GBM (Villanueva-Meyer et al. 2017), and it is easier nowadays to get information from MRI as some deep learning techniques can be used to extract medical data (Lundervold and Lundervold 2019). With tools like CaPTK (Bakas et al. 2017;Pati et al. 2020), it is possible to perform segmentation of GBM tumours based on MRIs. Recent studies also show that information on the tumour cells behaviour can be acquired with immunohistochemistry data, for example by identifying GBM subtypes (Orzan et al. Fig. 1 Endothelial cells carry blood vessels providing nutrients and O 2 in the brain. Due to tumour growth, hypoxic tumour cells are induced by a lack of O 2 . Hypoxic cells produce proangiogenic factors, mainly Vascular Endothelial Growth Factors, that enhance the creation of new blood vessels 2020) but we will not consider those different subtypes in this work. Working on MRIs is numerically challenging because on MRIs we cannot have constrained mesh to solve our equations on. Finite volume scheme based on TPFA (Two Point Flux Approximation) cannot ensure the positivity of numerical solutions. It is then needed to use more sophisticated numerical schemes in order to ensure the positivity of the solutions. Our approach is based on a CVFE (Control Volume Finite Element) scheme in which nonlinear numerical Gudonov fluxes are used to ensure the positivity.
To best describe the tumour behaviour from a patient diagnosis, we have to model the different treatments given. Currently, patients with GBM are treated using surgery and radiotherapy with adjuvant TMZ (chemotherapy), we will only consider those treatments in our model. Chemotherapy and surgery were first used in PDEs models around gliomas in Tracqui et al. (1995), Woodward et al. (1996) but more robust models have been developed: on chemotherapy with hypoxic cells (Hinow et al. 2009), on chemotherapy with heterogeneous drug delivery (Swanson et al. 2003), on surgery and radiotherapy with an haptotaxis model (Enderling et al. 2010), on chemotherapy and radiotherapy on low-grade glioma (Ribba et al. 2012), on using anti-angiogenesis drug with grow-or-go model (Saut et al. 2014) and metastatic spreading (Benzekry et al. 2012), on anti-invasive agents on solid cancers (Ribba et al. 2006) and even on immunotherapy in gliomas (Banerjee et al. 2015). Choosing to model those treatments will allow us to compare their theoretical impacts on the GBM growth with empirical clinical knowledge.
From now, the objectives of this work are -building a model based on angiogenesis and on the treatment schedule of a patient with GBM (Sect. 2), -converting the model into a numerical scheme (Sect. 3) and proving its adequacy with biological knowledge, e.g., positivity of solutions, upper-boundedness of cell population and existence of numerical solution (Sect. 4), -performing numerical simulations based on a patient's diagnosis. This includes the building of a mesh fitting specifically the patient's brain and its tumour, choosing the biological parameter in the model, and observing the impact of treatments on the tumour behaviour (Sect. 5), -discussing the reliability of the theoretical model and its usefulness (Sect. 6).

The anisotropic degenerate nonlinear angiogenesis model
We propose in this work a new model on angiogenesis inspired by the works of Enderling et al. (2010) and Hinow et al. (2009). Those works involve reactionadvection-diffusion equations around tumour cells and nutrients. In order to exhibit angiogenesis, two quantities are added into our model: endothelial cells that release O 2 in the brain and VEGF (Vascular Endothelial Growth Factor) that are produced by hypoxic tumour cells to promote the formation of new endothelial cells. In this section, we will describe the angiogenesis model and highlight all the biological processes.
Let Ω be an open bounded polygonal and connected subset of R 2 and T f > 0 a fixed finite time. We denote represents the area in the brain where the tumour is developing, here it is obtained from a slice of a pre-surgery axial MRI of a patient. ∂Ω is then the border of the brain around the skull and the ventricles if they are on the MRI (it depends on the location of the tumour). The behaviour of tumour cells during their spread and treatments is chosen to be described by the set of equations We associated with (1a)-(1d) homogeneous zero-flux boundary conditions These conditions model the no-exchange between the brain and the rest of the body. For each quantity, we associated an initial condition on Ω given by , w = u, c, u e , v. (3) In the model (1a)-(1d), u is the ratio between the number of tumour cells per cm 2 and the maximum tissue capacity u max , c is the concentration in O 2 in µ mol cm −2 , u e is the ratio between the number of endothelial cells per cm 2 and the maximum tissue capacity u max and v is the concentration in VEGF in µ mol cm −2 . The sum of the two cellular populations is u T = u + u e . This modelling enables the cell populations to be normalized between 0 and 1, e.g., 0 ≤ u, u e ≤ 1 (this will be proven in Sect. 4). Biologically, we would also expect to have 0 ≤ u + u e ≤ 1 because the maximum tissue capacity is for the total cell population and not only for each cell type independently. Even though our modelling choices do not allow us to prove this property, it is obtained in our simulations (see Sect. 5). The functions a(·), χ · (·) and f u T (·) are the cell-density dependent functions for diffusion, chemotaxis and growth rate respectively. More specifically, the functions a(·), χ · (·) and f u T (·) are used to degenerate the model. This is crucial in order to prove the positivity and upperboudedness of solutions (proof in Sect. 4). g(·) is the O 2 -dependent VEGF production by tumour cells function, e.g., tumour cells produce VEGF if and only if they are not in a hypoxic area. h(·) is an O 2 -dependent threshold allowing mitosis for tumour cells under normoxic conditions, e.g., tumour cells cannot undergo mitosis if they lack O 2 . Λ 1 (x) and Λ 3 (x) are the medium-dependent diffusion tensors for tumour cells and endothelial cells respectively. The diffusion of cells depends on the white matter, the grey matter and the post-surgical area. D 2 and D 4 are the isotropic diffusion tensors associated with O 2 and VEGF respectively. ρ 1 is the growth rate of tumour cells, α 2 is the production rate of O 2 by endothelial cells, ρ 3 is the growth rate of endothelial cells and α 4 is the production rate of VEGF by tumour cells. β 1 and β 3 are the apoptosis rates of tumour cells and endothelial cells respectively, β 2 and β 4 are the degradation rates of O 2 and VEGF respectively. γ 2 is the consumption rate of O 2 by tumour cells and γ 4 is the consumption rate of VEGF by endothelial cells. The map T treat (·, ·) represents the loss of tumour cells due to treatments.
There are a lot of parameters previously described, that appears in our modelling. Although some values can be extracted from the literature, most parameters have to be seen as patient dependent. Indeed, GBM regroups a lot a different tumour types, and so, different behaviours occur depending on the tumour type. As a consequence, most parameters will need to be tuned per patient. We are not focusing on this issue in this paper, and we will assume that the parameters used fit the patient.
Surgery, chemotherapy and radiotherapy are the treatments commonly used against GBM, all of them can be modelled in (1a)-(1d). Surgery is performed in emergency as soon as GBM is diagnosed. However, some patients cannot undergo surgery due to the non-accessibility of the tumour, in that case only a biopsy is done. If surgery is performed, the tumour core is removed and resection goes the largest possible without damaging healthy tissues.
In our model, we suppose that surgery is performed at a time t surg that can be either the initial time or a random time in ]0, T f [. Surgery is performed numerically by resetting to 0 the value of u, c, u e and v in the surgical operation area right after t surg .
Chemotherapy is the use of a drug designed against a tumour cell population. The molecule commonly used in GBM is TMZ at a daily dose of 75 mg/m 2 (Stupp et al. 2005). The chemotherapy part of treatments is modelled by where k c (t) = 1 if and only if the chemotherapy is effective at the time t and D che is the death rate of tumour cells due to the chemotherapy. The problem of chemotherapy is the possible existence of tumour-resistant cells that are not affected by the drug. Those cells keep proliferating, usually causing relapses for patients. For this work, we do not consider subtypes of cancer cells and suppose that all tumour cells are altered by the drug. This assumption is strong because patients treated with TMZ usually do not respond to TMZ (Lee 2016). However, TMZ is not used for newly diagnosed GBM for its alkylating properties, but only as an adjuvant of radiotherapy. Indeed, we know that the use of chemotherapy enhances the performances of radiotherapy because this drug makes the tumour cells more sensitive to radiation (Stupp et al. 2005).
Radiotherapy occurs for almost all cancer treatments because the use of radiotherapy depends mainly on the location of the tumour and its spatial spread and not on the cancer type. Indeed, radiotherapy works by sending a dose of radiation at a local position (the tumour location), that radiation cause micro-breaks into the DNA of irradiated cells. Normal cells can repair those DNA breaks, but tumour cells often cannot and die. Radiotherapy can cause side effects if normal cells are altered by the radiation, explaining why each cancer type has a specific guideline to dose the quantity of radiation allowed.
With GBM, radiotherapy is done by administering the dose of radiation in small fractions. The number of fractions and the cumulative dose depend on the patient clinical status, for example its age or its WHO performance status (ANOCEF 2018). The efficiency of radiotherapy depends on a lot of parameters: the number of fractions per day N f rac , the dose administered D rad , the duration of one irradiation τ , the time between irradiation Δτ , the DNA damaged rate μ and some sensitivity parameters α, β. In France, it is recommended to perform radiotherapy 5 days a week for 6 weeks with 30 fractions of 2G γ per day(ANOCEF 2018).
Because multiple small fractions are delivered with GBM, we model radiotherapy based on (Nilsson et al. 1990) where k r (t) = 1 if and only if the radiotherapy is effective at the time t and We use the term k chemo (t) to model the enhanced efficiency of radiotherapy if chemotherapy is applied concomitant with radiotherapy. So if chemotherapy is not administered at the time t, we set k chemo (t) to 1 and if radiotherapy and chemotherapy are concomitant then k chemo (t) = Υ > 1, Υ is the efficiency rate improvement of using radiotherapy with chemotherapy.
(A7) g(·) is a piecewise function that allows the production of VEGF by tumour cells only if the tumour cells are in a hypoxic environment and h(·) is the Heaviside step function around the hypoxia threshold, e.g. c necro is the threshold under which cells start to necrose and c hypo is the threshold under which cells lack of O 2 to be able to function normally.
(A8) The treatment map T treat : R + × R → R + is positive, piecewise in time and in space. In this work the available treatments are surgery, chemotherapy and radiotherapy. The map T treat (·, ·) models chemotherapy and radiotherapy and can be reconstructed using Eqs. (4) and (5) by T treat (t, y) = T chemo (t, y) + T radio (t, y). Moreover, ∀t ∈ R + , ∀y ∈ R − : T treat (t, y) = 0.
In order to ensure positivity of the solution in its discrete form, we use the following set of functions defined on R for the same ideas as in Cancès and Guichard (2016), Cancès et al. (2017), Foucher et al. (2018. We define the weak formulation associated to (1a)-(3). Definition (Weak Solution): Under assumptions (A1)-(A8), we say that the set of measurable functions (u,c,u • − • − Existence of the weak solution is not proved in this work. However, many works exist on Keller-Segel systems. Since our model consists of two Keller-Segel systems (coupled with a reaction term), we refer to (Bendahmane and Karlsen 2007) that proved the existence and Hölder regularity of a weak solution in a Keller-Segel system.

The nonlinear CVFE scheme for system
In this section, we discretize the model (1a)-(3) and introduce the numerical approximations necessary to perform simulations on a brain mesh, e.g., an unconstrained mesh fitting the shape of a brain.
The discretization of (1a)-(1d) is chosen following the work of Foucher et al. (2018), it uses two types of approximations: a conforming Finite Elements (FE) approximation for the diffusion terms, namely second terms in system (1a)-(1d), and a decentred Finite Volume (FV) scheme for the convection terms, namely the third terms in Eqs. (1a) and (1c). The FE approximations are done over a primal triangular mesh and the FV approximations are done over a dual barycentric mesh.
Let T be a conforming triangulation of the domain Ω, we denote by ϑ the set of vertices and E the set of edges in ρ J is the regularity of the mesh, where ρ J diameter of the incircle of the triangle J . For every vertex K ∈ ϑ, we denote by x K its coordinates, E K the set of edges having the vertex K as an extremity and T K the set of triangles that have K as a vertex. If two vertices K and L are joined by an edge, then we denote this edge by σ K L .
For every vertex K ∈ ϑ, we associate its dual element ω K constructed by connecting the barycentres of the triangles in T K with the barycentres of the edges in E K , the 2dimensional Lebesgue measure of ω K is m K . We denote by M the dual-mesh and H T the P 1 -finite element space on Ω defined by We associate H T with its canonical basis (Φ K ) K ∈ϑ . Furthermore, we consider the discrete control volume space χ M on Ω defined by In this paper, we choose a uniform time discretization with a time step δt = where N is a non-negative integer, and we set t n = nδt, for all n ∈ 0, N + 1 .
For a given (w n K ) K ∈ϑ,n∈ 0,N +1 , there exists a unique finite element reconstruction w T ,δt and a unique constant piecewise reconstruction w M ,δt such that The nonlinear CVFE scheme for the discretization of system (1a)-(3) is given by the following set of equations and In the above system, we have used a Finite Element approximation for the diffusion fluxes where the stiffness coefficients are given by We define the intervals to build a Godunov approximation for a n+1 Those terms are useful to ensure the positivity of the quantities. The functions μ 1 and μ 3 are approximated using an upwind scheme where the functions μ (1) ↓ and μ (3) ↑ are given by The description of all variables, coefficients, functions, spaces and mesh components are summed up in the supplementary Tables 5-9.

Discrete properties
In this section, we are interested in proving the good properties of the numerical resolution of (14)-(17). First, we show the positivity of all quantities and the upper-boundedness of cell populations. Then we introduce the iterative algorithm used to solve (14)- (17) in practice. Finally, we prove the existence of the numerical solution with the discrete energy estimates associated.

Positivity and upper-boundedness of quantities
We want to prove that the numerical scheme (14)- (17) is biologically trustable. Here, we will prove that the density of tumour cells and endothelial cells are minimized by 0 and maximized by 1 whereas the concentration in O 2 and VEGF are only positive. Equation (14) and (16) have similar terms, we prove the positivity of their respective quantity in Proposition 1. (16)) is positive. Proof We will show the result only for u n but the same steps are followed for showing the positivity of u n e . This proof works by induction on n, let's suppose that for a n ∈ 0, N we have u n

Proposition 1 (Positivity of tumour cells concentration)
L . Then by multiplying the equation of (14) associated to K by −u − K , we get According to assumptions (A5) and (A8), the functions f u T (·) and T treat (t, ·) are extended by zero outside of [0, 1], which implies that (u K − u L ) ≥ 0 due to the definition of u K and Λ (1) + K L a n+1 K L ≥ 0, so we have the positivity of the third term in (23), e.g.
For the fourth term in (23) we have μ (1) According to the definitions a n+1 K L in (19) and μ n+1 K L in (21), we have Λ , which gives the positivity of the fourth term in (23), e.g.
With inequalities (24) and (25), we can then conclude that the left-hand side of (23) is positive. However So all the terms on the left side of the above inequality are null because they were non-negative, implying that m K δt u − 2 K = 0, giving the result of the proposition. Now, we aim to prove that the density of tumour cells and endothelial cells are upperbounded by 1. This maximization is meaningful because tissues have a maximum capacity available. As for Proposition 1, this upper-boundedness is feasible because our model is degenerated. Proposition 2 summarizes this property. (16)) is upper-bounded by 1.

Proposition 2 (Boundedness of tumour cells concentration)
Proof We will show the result only for u n , but the same steps are followed for showing that u n e is upper-bounded. This proof works by induction on n. The idea is to show that the maximum value of {u n+1 K , K ∈ ϑ} = u n+1 K is upper-bounded by 1. This statement Proof We will show the result only for c n but the same steps are followed for showing the positivity of v n . This proof is very similar to the one of Proposition 1, but we need to use the transformations defined in Eq. (7) to ensure the monotony of the diffusion term. The detailed proof of this proposition is available in appendix B.
Something essential for the rest of the paper is to notice that Propositions 1, 2 and 3 remain true even if the system (14)-(17) becomes semi-implicit in time rather than fully implicit. In other words, we could consider a semi-implicit scheme by replacing c n+1 by c n in (14), u n+1 by u n in (15) and (17), v n+1 by v n in (16), u n+1 e by u n e in (15) and (17) and Proposition 1,2 and 3 would still hold. Knowing that, we are going to approximate the solution of the implicit scheme (14)-(17) with an iterative algorithm based on the semi-implicit version of this system. That way, biological properties are still ensured during the numerical resolution.

The algorithm to get the discrete solution (14)-(17)
In this work, we propose an iterative algorithm in order to get the discrete solution of the implicit numerical scheme (14)- (17). The main idea of this algorithm is to solve the nonlinear system (14)-(17) with an iterative method. The usual resolution of the full non-linear system (14)-(17) relies on a classical Newton method involving simultaneously the four unknowns at each iteration. We would rather solve four non-linear systems corresponding to each equation (14) to (17), consequently, diminishing the memory footprint. Moreover, it is more efficient to solve the four equations, corresponding to Eqs. (14)-(17), rather than the coupled system. That's why in this work we develop an iterative method that decouples the system (14)-(17), simplify the numerical resolution and converges to the solution of the implicit numerical scheme. Let's suppose that ∇c 0 and ∇v 0 are in (L 2 (Ω)) 2 and consider a given solution Π n T = (u n K , c n K , u n e,K , v n K ) K ∈ϑ for a given n. We look for the solution Π n+1 T at t n+1 as the limit of the iterative process described in Algorithm 1.
/* m is the number of iteration of the process */ while The iterative process has not converged do Observe that at each iteration m, the system (26)-(29), in Algorithm 1, is noncoupled.
Let τ (n) :X = (ũ,c,ũ e ,ṽ) → X = (u, c, u e , v) be the application that for the vectorX = (ũ,c,ũ e ,ṽ) ∈ (R #ϑ ) 4 associates the solution X = (u, c, u e , v) of (26)-(29). Replacing (u m , c m , u m e , v m ) byX and (u m+1 , c m+1 , u m+1 e , v m+1 ) by X , the iterative method is equivalent to Note that the numerical scheme (26)-(29) is consistent, e.g., if the sequence (X m ) m converges, then X m → m→∞ Π n+1 T . We consider that X m+1 = (u n+1 , c n+1 , u n+1 e , v n+1 ) when the relative error follows with I = 10 −3 being a numerical threshold. To find X m+1 from X m we solve (26)-(29) using Newton's method with N = 10 −11 as the error of convergence. Observe that choosing to solve the semi-implicit scheme (26)-(29) instead of the coupled system (14)- (17) allows solving four smaller systems rather than one system four times wider. Finally, we use a conjugate gradient method ending with a numerical threshold G of 10 −13 to solve each step of the Newton method.
Let's remember that 0 ≤ u n K , u n e,K ≤ 1 and c n K , v n K ≥ 0 for all K ∈ ϑ and n ≥ 0, so the solutions of the iterative method follow 0 ≤ u m , u m e ≤ 1 and c m , v m ≥ 0 for all m ≥ 0 due to Propositions 1-3.
To prove the existence of the solution of (26)-(29) and that the application τ (n) is well-defined, we will use the following energy estimates on the iterative method.
Lemma 2 There exists C 5 (D 2 , θ T ), C 6 (D 4 , θ T ) such that Proof The proof relies on the fact that and that the quantity is non-negative because p(·) is non-decreasing, we can then use the same arguments as Lemmas 3.1-3.3 in Cancès and Guichard (2016) to show that The same arguments are used for the inequality in v.
Now we want to prove the existence of the solution of (26)-(29). The proof is heavy and relies on the previous energy estimates and the Brouwer fixed-point theorem. The existence of the solution of (26)-(29) and the properties ensured are given by Proposition 4. Let's consider the application: ) is a solution of (26)-(29).
The aim in the first step is to prove that ∃k ≥ 0, ∀ y 2 > k : (W (y), y) ≥ 0. To do so, we develop (W (y), y) in four terms, corresponding to each equation of (26)-(29), and use the energy estimates given prior to this proposition. We show that (W (y), y) ≥ C 1 y 2 2 −C 2 y 2 − C 3 . Where C 1 , C 2 and C 3 are strictly non-negative constants.
Finally, we apply the Brouwer fixed-point theorem to the application S : y ∈ B(0, k) → −k W (y) W (y) ∈ B(0, k) to show that the function W has a zero point. The detailed proof of Proposition 4 is available in appendix C.
To sum up, we have shown that the implicit numerical scheme (14)-(17) was the limit of the iterative semi-implicit numerical scheme (26)-(29). Moreover, we have proved that this system ensured the biological properties expected at each iteration of Algorithm 1. And that the existence of the solution of (26)-(29) was ensured during the iterative process.

Numerical simulations
Our numerical simulations are based on the MRIs from the patient C3L_16 in the CPTAC-GBM database (TCIA 2018). This patient is a 60-year-old male (BMI of 28.81 and BSA of 2.07 m 2 ) diagnosed with a Glioblastoma of 4.3 cm in his parietal lobe who died 77 days after diagnosis.
There are different types of MRI from this patient that can be used to extract information as shown in Fig. 2: a T1 highlights the white matter in white while grey matter is not highlighted, a T1-Gado separates well the tumour core and edema from the brain, T2 and FLAIR show the enhancing tumour area.    Fig. 2. b Segmentation of the brain on the primal mesh T according to Fig. 2: the grey matter in dark blue, the white matter in light blue, the tumour core in green, the edema in yellow and the enhancing tumour in red. P1, P2 and P3 are three points used for investigating the local tumour cell density in the areas of the tumour (colour figure online) region of interests, e.g., on the outside of the brain, of the tumour core, of the edema and of the enhancing tumour. Once those vertices have been defined, we have used the software Triangle (version 1.6) (Shewchuk 1996(Shewchuk , 2002 to get the entire triangulation T . The primal and dual mesh obtained are represented in Fig. 3a. Information on the primal mesh T are given in Table 1. With the information readable in the MRIs of Fig. 2, we can extract data manually in the brain: the white and grey matter locations and the tumour segmentation. The segmentation of the brain is displayed in Fig. 3b in which we have, in different shades of blue, the locations of grey matter (in dark blue) and the white matter (in light blue) and the whole tumour which is divided into three parts: the tumour core with hypoxic cells in green, the edema in yellow and the enhancing tumour in red. For simplification, we suppose before surgery that the diffusion in the whole tumour behaves like the white matter.
Even though more detailed anisotropic expressions have been developed for the diffusion in the white matter, for example in (Painter and Hillen 2013), we choose here an easier expression not to rely on Diffusion Tensor Imaging (DTI) data. The diffusion in the white matter is then set to be 5 times faster than in the grey matter. In the post-surgical area, we set the diffusion being 50 times faster than in the grey matter. We have The different coefficients in (1a)-(1d) following the assumption (A4) are given in Table 2.
Following the timeline of a patient with GBM (ANOCEF 2018), we set the treatments as follows: -Surgery is performed on day 1, e.g., 24 h after the start of the simulation. Surgery removes every cell and every protein in the necrotic, enhancing and whole tumour area as segmented in Fig. 3b. -Radiotherapy and chemotherapy start both on the day 14 and last for a period of 6 weeks. Experiments in Stupp et al. (2005) found a 33% increase in median survival comparing radiotherapy only and radiotherapy with TMZ, but there is no study using chemotherapy only on patients. Υ is chosen to be 4 3 even though the value is overestimated.
The value of all parameters related to treatments are summarized in Table 3. Most of them come from the work of (Powathil et al. 2007) with the current state of treatments in (ANOCEF 2018).
Six simulations have been done to solve (1a)-(3), depending on the treatments administered to the patient: with or without surgery, chemotherapy and radiotherapy only (experiments 1 to 4) and then the concomitant use of chemotherapy and radiotherapy with or without surgery (experiments 5 and 6). The simulations' settings are given in Table 4 using the coefficients from Tables 2 and 3. For the initial conditions, we set the values of u 0 and u e 0 from the MRIs of the patient given in Fig. 2. According to the location on the MRIs of the tumour core and edema, we set u 0 to 12×10 4 cells cm −2 in the tumour core and to 5.5×10 7 cells cm −2  Fig. 4 a The initial number of tumour cells per cm 2 (u 0 ). b The initial number of endothelial cells per cm 2 (u e 0 ) on the right in the edema. In the enhancing tumour area, we set the tumour cells population with u 0 (x) = 1.4 × 10 8 exp(− 4 45 |x − x c |) cells cm −2 , where x c represents the position of the tumour's centre, and outside those areas we set u 0 to 0. According to the location of the grey and white matter, we set u e 0 to be 7.17 × 10 7 cells cm −2 in the grey matter and to 2.39 × 10 5 cells cm −2 in the white matter. The initial spatial distributions of tumour and endothelial cells are shown in Fig. 4. O 2 and VEGF concentrations cannot be seen on MRIs required for GBM diagnosis, such as those presented in Fig. 2. We have attributed their initial concentrations by using the previous expressions of u 0 and u e 0 , and then, calculating the stationary solutions of Eqs. (15) and (17), e.g., without the partial temporal term. The computed initial concentrations in O 2 and VEGF are shown in Fig. 5. We would expect their spatial distributions to be smoother around the tumour, but during the first times of the simulation, those distributions smooth welly.
Indeed, during the simulations, the O 2 concentration stabilizes between 74 and 104 µ mol cm −2 and VEGF are produced punctually by the hypoxic tumour cells. Without any treatment the tumour keeps growing inside the brain, inducing the growth of the hypoxic tumour core, the edema and the enhancing tumour area. Moreover, without In order to follow up the impact of treatments on the global tumour behaviour, we have displayed in Fig. 6 the total number of tumour cells and in Fig. 7 the total amount of VEGF through time.
We also display in Figs. 8, 9 and 10 the number of tumour cells per cm 2 at three points in the brain, namely at P1, P2 and P3. P1 is located in the tumour core, P2 in the edema and P3 in the enhancing tumour area as shown in Fig. 3b. With those figures, we analyse the local impact of treatments on tumour cells.
You can find a video of the 6 simulations related in Table 4 following this link https://www.youtube.com/watch?v=vJkMJ5bNoWA.
The tumour growth behaviour depends highly on the treatments used, e.g., surgery, chemotherapy and radiotherapy. We will analyse the impact of each of them in the following paragraphs.
Using surgery on day 1 allows decreasing intensely the number of tumour cells in the brain while stopping tumour induced angiogenesis because there is no more hypoxia in the brain. It remains only tumour cells in areas where no evidence of existence were detectable on the different MRIs. However those cells keep growing, because the brain does not lack of O 2 anymore, and may induce relapse. This is why in Fig. 6 the trajectory of the surgery only curve starts with a big drop on day 1 but comes back to the no treatment curve later on. We can even see in Fig. 7 that tumour cells start induced angiogenesis again around day 70.
In conclusion, surgery is not enough to stop the spread of the tumour growth. However, it stops locally the spread as there is no tumour cells remaining at P1 and P2 after the surgery (which is why in Figs. 8 and 9 the curves where surgery was used are not displayed). P3 is in the surgical area too, but also on the boundary of the Fig. 6 Number of the total tumour cells in the brain through time for different treatment schedules. If performed, surgery is done on the first day, chemotherapy (TMZ) is administered from day 14 to day 56 and radiotherapy is administered 5 days out of 7 from day 14 to day 54. The curves "No treatment" and "TMZ" are really close and difficult to distinguish at this scale Fig. 7 Total amount of VEGF, in µ mol, in the brain through time for different treatment schedules. If performed, surgery is done on the first day, chemotherapy (TMZ) is administered from day 14 to day 56 and radiotherapy is administered 5 days out of 7 from day 14 to day 54. The curves "Radiotherapy" and "TMZ+Radiotherapy" are really close and difficult to distinguish at this scale. Spikes observable before day 5 are due to poor initial condition knowledge enhancing tumour area. So some tumour cells near P3 remain after the surgery, which explains why in Fig. 10, there are still tumour cells located at P3 after the surgery.
The use of TMZ only does not affect the global behaviour of the tumour growth. Indeed, in Fig. 6 the trajectories of the "TMZ" curve and the "no treatment" curve are very similar. However, TMZ does impact the local behaviour of tumour growth at two  Fig. 3b. If performed, surgery is done on the first day, chemotherapy (TMZ) is administered from day 14 to day 56 and radiotherapy is administered 5 days out of 7 from day 14 to day 54. Treatments using surgery are not displayed because no tumour cell remains after surgery at P1 Fig. 9 Number of the tumour cells per cm 2 at P2 through time for different treatment schedules. The location of P2 (the edema) is shown in Fig. 3b. If performed, surgery is done on the first day, chemotherapy (TMZ) is administered from day 14 to day 56 and radiotherapy is administered 5 days out of 7 from day 14 to day 54. Treatments using surgery are not displayed because no tumour cell remains after surgery at P2  Fig. 3b. If performed, surgery is done on the first day, chemotherapy (TMZ) is administered from day 14 to day 56 and radiotherapy is administered 5 days out of 7 from day 14 to day 54 levels. On one hand, TMZ is an alkylating agent and induces DNA breaks into the tumour cells. Those DNA breaks do kill tumour cells, but the death rate is really low. On the other hand, this decrease in tumour cells density implies that the hypoxic area is shrinking. We can see this shrink in Fig. 7 as less VEGF are produced during the chemotherapy. Due to this shrinking, more tumour cells can reproduce. In conclusion, there is a balance between the death rate in tumour cells and the fact that more tumour cells can reproduce, which does not change the global behaviour. Locally, tumour cells can either be in a slightly higher number than without any treatment as at P2 or in a slightly fewer number as at P1 and P3. It must be remembered that TMZ is used concomitant with radiotherapy in practice and that no patient can be prescribed TMZ only to cure their newly diagnosed GBM.
Lastly, radiotherapy affects a lot the behaviour of the tumour growth in our simulations. Its use decreases drastically the number of tumour cells in the irradiated area, as shown in Figs. 6,7,8,9 and 10. The death rate of tumour cells due to radiotherapy is also enhanced when combined with TMZ as it was exhibited in Stupp et al. (2005). Figures 8, 9 and 10 describe the local density of tumour cells depending on the treatment. In that regard, the behaviour of the tumour can drastically differ between treatments and can be counter-intuitive knowing the global behaviour.
For example, in Figs. 8 and 9, the tumour cell density decreases in the case without treatment. This has to do with the angiogenesis process and more specifically on whether the tumour cells are in a hypoxic area or not. Inspecting Eq. (14), the reaction term has two behaviours:

Fig. 11
Infinity-norm of the normalized total cell population u T between 0 and 1 through time. The function t → u T (t, .) L ∞ (Ω) is displayed for each simulation. A zoomed portion is given to see every trajectory -locally we are in a hypoxic or necrotic area (P1 and P2 with no treatment). This implies that h(c) f u T (u) = 0 and so that the reaction term is strictly negative, e.g., the density of tumour cells is decreasing through time. -locally we are in a normoxic area (P3 with no treatment). This implies that h(c) f u T (u) ≥ 0 and so that the reaction term is usually positive, e.g., the density of tumour cells is increasing in time.
This reasoning also explains how the density of tumour cells can be larger at P1 and P2 in Figs. 8 and 9 with radiotherapy than without. In our simulations, we can see in Fig. 6 that the use of radiotherapy and chemotherapy gives the minimum number of tumour cells around the 50th day rather than with the use of surgery, radiotherapy and chemotherapy. This observation can be explained due to the fact that the use of surgery enhances the proliferation of the remaining tumour cells, because there are no more hypoxic tumour cells. Those cells migrate further into the brain, escaping the irradiated area. Yet, in the long run there are fewer tumour cells in the brain when using surgery, TMZ and radiotherapy than TMZ and radiotherapy which implies a longer survival time when using all treatments.
As introduced in Sect. 2, there is a maximum tissue capacity that would imply that 0 ≤ u + u e ≤ 1. Even though our modelling choices do not allow us to prove that property, we do have it in practice in our simulations. Please check the following link https://youtu.be/dw4K8X5nLcg, as it leads to a video displaying the cell population density (u+u e ) in the brain during all simulations. The infinity-norm of the normalized total cell density over time is also displayed in Fig. 11.

Discussion
It must be remembered that all results presented in this work rely on the behaviour of our theoretical model with its theoretical treatments.
This model is strongly suited for patients with GBM. Indeed the modelling uses the usual treatments administered and numerical simulations are done on a mesh with a geometry fitting the brain and the tumour. Moreover, we have proved in this work the existence of the numerical solution, and so, the numerical resolution is ensured. Finally, the global behaviour of the tumour cells fits empirical knowledge. In all the simulations, after the delivery of all the treatments, the tumour starts proliferating again until being slowed by the hypoxic tumour core.
Even with those strengths, the model has some weaknesses. Specifically, with the relationship between treatments and patients' response. Indeed, in our model, if we wanted to enhance the death rate of tumour cells due to radiotherapy, we would increase the irradiation dose. However, a higher dosage would impact the healthy cells in the brain that would not be able to repair their DNA breaks as before. Also, surgery in our model seems to amplify the spread of the remaining tumour cells. It would then be better to only use radiotherapy with concomitant TMZ as evidenced in Fig. 6 but if surgery is not performed then there is a necrotic tumour core in the brain that could have negative impacts on the surrounding healthy tissues. For example, the O 2 delivery would decrease for the healthy tissues in favour of supplying the hypoxic tumour cells.
Also, we decided in this work to consider only one population of tumour cells that react homogeneously to treatments. However, some tumour heterogeneity could be introduced in our model as we are working on Glioblastoma Multiforme, e.g., a cancer with intra-tumour and inter-tumour heterogeneity. For example, tumour heterogeneity could be on the efficiency of each treatment, on the diffusion rate or migration rate of each cell phenotype. Ultimately, it would be necessary to consider subpopulations of tumour cells that would not follow the exact same equation as (1a) depending on the origin of the heterogeneity.
In this work, modifying the value of the different coefficients in (1a)-(1d) is the only way in our model to exhibit different growth behaviours, like having a higher tumour proliferation, having a faster VEGF production or having a higher tumour diffusion rate. This approach works when the coefficients in (1a)-(1d) are known and well identified on a patient. Nowadays, a lot of information can be retrieved from the patient diagnosis: using immunohistochemistry-based algorithm (Orzan et al. 2020), analysing the extracellular vesicles situated in the GBM micro-environment (Simon et al. 2020) or by determining the GBM subtypes based on the OMS description (Louis et al. 2016). However, this information is not linked explicitly to the different coefficients, and so, the coefficients have to be extrapolated from the information and then adapted to fit the growth of the patient tumour.
A solution to find the coefficients of a patient could be to set some coefficients as unknowns and write the system (1a)-(3) as an optimization problem based on the knowledge of the solution at different time steps. However, this method shows weaknesses as there are more unknowns than equations which implies the need for more data that are not available from the diagnosis, so we cannot do simulations or predictions after the diagnosis.
To solve this problem, we could set the coefficients, not as unknowns, but as temporal functions and use Kalman and Luenberger filters as in Collin et al. (2020) to fit the model to the patient through time. This method allows having a unique model that will adapt to each patient during the simulations, matching simulations and data simultaneously.
Funding F.A is a recipient of a Ph.D. fellowship from the Ministère de la recherche, A.A.S is supported by chaire INSERM-Ecole Centrale de Nantes, we are also supported by IEA-CNRS and HIPHOP project.
Code availability Custom code.

Conflict of interests
The authors declare that they have no conflict of interests.

A Proof of Proposition 2
We will show the result only for u n but the same steps are followed for showing that u n e is upper-bounded. This proof works by induction on n, let's suppose that for a n ∈ 0, N we have u n K ≤ 1, ∀K ∈ ϑ. Let's u K = u n+1 H and ∀L ∈ ϑ, u L = u n+1 L . Then by multiplying the equation of (14) associated with K by (1 − u K ) − , we get Using proposition 1, we know that u K ≥ u L ≥ 0, ∀L ∈ ϑ. If Λ K L < 0, we have a n+1 K L (1−u K ) − = 0 due to the fact that a(·) is set to zero outside of [0, 1]. So knowing that a n+1 K L ≥ 0, we have then observing that (u K − u L ) ≥ 0, we conclude that the third term in (39) is positive, e.g.
The function f u T (·) is extended by zero outside of [0,1], implying that f u n+1 So whenever the sign of Λ (1) and, having a n+1 K L ≥ 0, we conclude that the fourth term in (39) is positive, e.g.

B Proof of Proposition 3
We will show the result only for c n but the same steps are followed for showing the positivity of v n . This proof works by induction on n, let's suppose that for a n ∈ 0, N we have c n and ∀L ∈ ϑ, c L = c n+1 L . Then by multiplying the equation of (15) associated with K by −c − K , we get Using the results of Proposition 1 and 2, we have u n+1 e , u n+1 If D K L < 0 then η n+1 K L c − K = 0. Observing that μ n+1 K L ≥ 0 and that the function p(·) is non-decreasing on R, we have ( p(c K ) − p(c L )) ≤ 0, meaning that Then The terms on the left side of (42) being non-negative, we conclude that they are all null. In particular c − 2 K = 0, which completes the proof of the proposition.

C Proof of Proposition 4
This works by induction on m.
The aim in the first step is to prove that ∃k ≥ 0, ∀ y 2 > k : (W (y), y) ≥ 0. So we develop the expression of (W (y), y) = (W u (y u ), y u ) + (W c (y c ), y c ) + (W u e (y u e ), y u e ) + (W v (y v ), y v ), y = (y u , y c , y u e , y v ) to build inequalities. We have With Lemma 1, we know that σ K L ∈E Λ (1) K L a m+1 K L (y K −y L ) 2 ≥ 0, σ K L ∈E Λ (3) K Lã m+1 K L (y K −y L ) 2 ≥ 0 and with Lemma 2 we have σ K L ∈E D K L η m+1 K L ( p(y K )− p(y L ))(y K − y L ) ≥ 0, σ K L ∈E D (4) K Lη m+1 K L ( p(y K ) − p(y L ))(y K − y L ) ≥ 0 because p(·) is non decreasing on R.