A Glioblastoma PDE-ODE model including chemotaxis and vasculature

In this work we analyse a PDE-ODE problem modelling the evolution of a Glioblastoma, which includes chemotaxis term directed to vasculature. First, we obtain some a priori estimates for the (possible) solutions of the model. In particular, under some conditions on the parameters, we obtain that the system does not develop blow-up at finite time. In addition, we design a fully discrete finite element scheme for the model which preserves some pointwise estimates of the continuous problem. Later, we make an adimensional study in order to reduce the number of parameters. Finally, we detect the main parameters determining different width of the ring formed by proliferative and necrotic cells and different regular/irregular behaviour of the tumor surface.


Introduction
Among the group of brain tumors, the Glioblastoma (GBM) is the most aggressive form with a survival of a little more than one year [26]. Moreover, GBM differs from many solid tumors in the sense that they grow infiltratively into the brain tissue, there exists an important presence of necrosis and they produce a high proliferation tumor cells. For all these reasons, GBM is one of the cancer types with more interest in the mathematical oncology community (see [1,4,32] and references therein).
Some studies about the morphology of GBM are based in the magnetic resonance images (MRI) in order to obtain results related to prognosis and survival (see [23,28,29,30]). Specifically, Molab 1 group classifies the GBM depending on the width of the tumor ring and/or the tumor surface regularity (see [28,30] respectively). The study of [28] concludes that tumors with slim ring have better prognostic, specifically 7 months of more survival than tumors with thick ring. In [30], the survival of patients in relation to the surface growth, regular or irregular, of the GBM, show that tumors with a regular surface have better prognostic, more than 5 moths of survival, than tumor with irregular surface.
In [39], the authors use the Fisher-Kolmogorov equation to reproduce the infiltrative characteristic of the GBM. However, more complex mathematical models are also built to simulate phenomena such that the tumor ring and the regularity surface of the GBM. One model appears in [27] where the tumor ring is studied by a PDE-ODE system of two equations (proliferative tumor and necrosis). In [11,12], the authors present a PDE-ODE system with three equations (proliferative tumor, necrosis and vasculature) which is able to capture different behaviours of tumor ring and regularity surface of the GBM via a nonlinear diffusion tumor increasing with vasculature.
In this paper, we present a PDE-ODE system, also with three equations (tumor, necrosis and vasculature) and we study the biological behaviours of the GBM such as the tumor ring volume, studied in [27,28], and the regularity surface considered in [30]. Unlike the system considered in [11,12], we have included a chemotaxis term. This term has been already introduced to model the movement of some populations towards a higher concentration of the chemical substance or another living organism, see for instance the reviews given in [5,10,15,16,31] and the references herein. Specifically, in this paper, we have included the chemotaxis term modelling the movement of tumor to vasculature.
Some previous chemotactic PDE-ODE models have been extensively studied in the literature, see for instance [6,33,34,35] where the authors model the cells movement with a parabolic-ODE system. Specifically, in [35] a system of PDEs is considered using a probabilistic framework of reinforced random walks. The authors analyse various combinations of taxis and local dynamics giving examples of aggregation, blow-up and collapse. Later, in [33], some analytical and numerical results which support the numerical observations of [35] are presented using a similar model than in [35]. Moreover, in [3,6] a model of tumor inducing angiogenesis is proposed consisting of a equation with chemotaxis and haptotaxis term, and two nonlinear ODEs. Finally, in [34] a stochastic system related to bacteria and particles of chemical substances is discussed where the position of each particle is described by a equation of a chemotaxis system.
Several works such as [19,36,37,38] have shown existence results for systems of three differential equations modelling cancer invasion. In [36] the global existence and boundedness of solution for a parabolic-parabolic-ODE system with nonlinear density-dependent chemotaxis and haptotaxis and logistic source is deduced. Furthermore, in [37], the authors have proved global existence of solutions for a parabolic-elliptic-ODE system with chemotaxis, haptotaxis and logistic growth. The study of existence of solutions for the chemotaxis and haptotaxis model with nonlinear diffusion is presented in [38]. The global existence of solution and its asymptotic behaviour are studied in [19] for a parabolicparabolic-ODE system modelling the cells invasion process.
Recently, a PDE-ODE model with chemotaxis is studied in [18] obtaining asymptotic stability results using a proper transformation and energy estimates. Another PDE-ODE with chemotaxis problem is considered in [24], see also [25], modelling the evolution of biological species and they obtain analytical results concerning the bifurcation of constant steady states and global existence of solutions for a range of initial data. In [14] a parabolic-ODE problem is analysed, and it is shown that, under several conditions, any stationary solution is locally stable.
In this paper, we investigate the following parabolic PDE-ODE system in (0, T f ) × Ω (Ω ⊆ R 3 is a bounded and regular domain and T f > 0 corresponds to the final time) endowed with non-flux boundary condition on the boundary ∂Ω where n is the outward unit normal vector to ∂Ω and initial conditions at time t = 0: Here, T (t, x), N (t, x) and Φ(t, x) represent the tumor and necrotic densities and the vasculature concentration at the point x ∈ Ω and time t > 0, respectively.
The nonlinear reactions functions f i : R 3 → R for i = 1, 2, 3 have the following form  The functions P (Φ, T ), S (Φ, T ), R (Φ, T ) and Q (Φ, T ) appearing in (1.4) are adimensional factors with the following biological meaning: Thus, these factor functions P (Φ, T ), S (Φ, T ), R (Φ, T ) and Q (Φ, T ) must satisfy the following modelling conditions: (1.5) and, We assume along the paper the following assumptions on the initial data (1.10) In order to obtain some estimates of the solutions of (1.1)-(1.3) (see (2.1)), we define the following truncated system of (1.1):  b) Assuming that there exists a constant C 1 > 0 such that c) Assuming additionally that there exist constants C i > 0 for i = 2, 3, 4 such that for all 0 ≤ Φ ≤ K and T ≥ 0, and ∇T is bounded in L 2 0, T f ; L 2 (Ω) .
By Theorem 1.1 a), for any (T, N, Φ) solution of (1.11), we deduce that T + = T , N + = N and Φ K + = Φ and then, . Hence, we obtain the following crucial corollary: The existence of solutions of problem (1.11) is out of the scope of this paper. It is an interesting open problem that could be treated in a forthcoming paper.
2. In Section 3, we design a Finite Element numerical scheme, computing Φ h k , T h k , N h k as an approximation of (Φ(t k , ·), T (t k , ·), N (t k , ·)) where t k is a partition of the time interval (0, T f ) and h is the mesh size. To build the scheme, we will use the change of variable in the PDE equation with chemotaxis, T = e κ ν Φ u, similar to the used in [8,9,20], in order to obtain an equivalent system with diffusion for the new variable u.  3. A parametric study through numerical simulations is made in order to detect different behaviours for the ring width and the regularity of the surface of the tumor.
The outline of the paper is as follows. In Section 2, we prove Theorem 1.1. In Section 3 we build a numerical scheme which preserves the a priori estimates of the continuous model given in Theorem 1.1 a). Later, in Section 4, we show a possible example of the dimensionless reaction functions of the system satisfying the hypotheses given in (1.6)-(1.9) and (1.12)-(1.16) and we make an adimensionalization of the model. Section 5 is dedicated to show, by means of some numerical simulations, the different behaviour of the ring width-volume and the regularity surface with respect to the dimensionless parameters. Finally, the more technical part of the proof of Theorem 1.1 b), obtained via an Alikakos' argument, is given in an Appendix.
Proof. Let (T, N, Φ) be a solution of (1.11). Since one can rewrite f 1 (T + , N + , Φ K + ) = T + f 1 (T + , N + , Φ K + ), multiplying the first equation of (1.11) by T − = min {T, 0} and integrating in Ω, we get We repeat the same argument for the other two equations of (1.11) using now that To obtain the upper bound Φ ≤ K, we multiply the third equation of (1.11) by Since Lemma 2.2. Any solution of (T, N, Φ) satisfies the estimates: Proof. Let (T, N, Φ) be a solution of (1.11). Integrating in Ω the first equation of (1.11) and using that P (Φ, T ) , S (Φ, T ) ≥ 0, we obtain that Thus, Rewriting P (Φ, T ) T = P (Φ, T ) P (Φ, T ) T and applying Young's inequality for the right side, we get, Hence, using that P (Φ, T ) ≤ 1, we conclude that Integrating in (0, t) for 0 < t ≤ T f , we obtain that To prove (2.3), we integrate the second equation of (1.11) in Ω × (0, t), with 0 < t ≤ T f , where we have used (1.5). Thus, using that Φ ≤ K and the bound obtained for T in (2.2), we get (2.3).

Proof of Theorem 1.1 b)
In order to obtain the L ∞ estimate for T , firstly we make a change of variable such that we rewrite the diffusion term and chemotaxis term as an unique diffusion term depending on the new variable. In fact, we consider: with u = e w and χ = κ ν .
Thus, the first equation of (1.1) changes to and the boundary condition (1.2) to ∇u · n = 0. Proof. To obtain the L ∞ estimates for T and N , taking into account the L ∞ estimates for Φ, it suffices that u be L ∞ . The proof of u is based in L p estimates with an Alikakos' argument. Let (T, N, Φ) be a solution of (1.11). We multiply (2.5) by u p−1 (for any p ≥ 2), and analyse term by term: • Time derivative term: and the second term of the right side of (2.7) can be expressed as Hence, from (2.7) and (2.8), • Nonlinear diffusion term: (2.10) • Reaction term: Rewriting in (2.9) the function Φ t as f 3 e χ Φ u, N, Φ and adding (2.9), (2.10) and (2.11), we get: (2.12) Due to hypothesis (1.13) and (1.12), it is easy to see in (2.12) that, Using now that 0 ≤ Φ ≤ K, (1.5) and (1.13) we obtain that with C > 0 independent of p (along the proof, we will denote by C different constants independent of p).
Using the auxiliary variable w = u p/2 , we can rewrite (2.14) as follows Thus, applying Gronwall's lemma, we deduce for p = 2 that Now, using the following equivalent norms with constants independent of p We are going to apply the following Gagliardo-Nirenberg interpolation inequality ([13, Theorem 10.1]) with ε > 0 and n the dimension of Ω (in this case n = 3). Applying (2.18) for ε = ν C p in the right hand side of (2.17), we deduce that Finally, due to (2.16), we can deduce that Hence, we obtain that (2.22) Following a similar argument to used by Alikakos in [2] (see Appendix), from (2.22) we can obtain that As consequence, T is bounded in L ∞ (0, T f ; L ∞ (Ω)).

Proof of Theorem 1.1 c)
Let (T, N, Φ) be a solution of (1.11). Taking gradient in the second and third equation of (1.11), (2.24) Using the change of variable T = e χ Φ u as in Lemma 2.3, we deduce that and we know from Lemma 2.3 that ∇ u is bounded in L 2 0, T f ; L 2 (Ω) . Taking into account that T and Φ are bounded in L ∞ (0, T f ; L ∞ (Ω)), it holds that Thus, rewriting (2.23) and (2.24) in terms of ∇ u, multiplying (2.23) and (2.24) by ∇ Φ and ∇ N respectively and integrating in Ω, we deduce and Using now Cauchy-Schwarz and Young's inequalities in (2.26) and (2.27) and adding them, it holds that with C i > 0 for i = 1, 2. Since ∇u is bounded in L 2 0, T f ; L 2 (Ω) , applying Gronwall's Lemma, it holds that ∇ N and ∇ Φ are bounded in L ∞ 0, T f ; L 2 (Ω) .

A FE numerical scheme
In this Section, we are going to design an uncoupled and linear fully discrete scheme to approach (1.1)-(1.3) by means of an Implicit-Explicit (IMEX) Finite Difference in time and P 1 continuous finite element with "mass-lumping" in space discretization. This scheme will preserve the pointwise estimates that appear in Lemma 2.1 considering acute triangulations.
Now we introduce the hypotheses required along this section. a) Let 0 < T f < +∞. We consider the uniform time partition with polygonal or polyhedral lipschitz-continuous boundary.
b) Let {T h } h>0 be a family of shape-regular, quasi-uniform triangulations of Ω formed by acute Nsimplexes (triangles in 2D and tetrahedral in 3D with all angles lowers than π/2), such that c) Conforming piecewise linear, finite element spaces associated to T h are assumed for approximating and its Lagrange basis is denoted by {ϕ a } a∈N h .
Let I h : C 0 Ω → N h be the nodal interpolation operator and consider the discrete inner product Before building the numerical scheme, we will transform the first equation of (1.1) into a non-linear diffusion equation throughout the change of variable T = u e χ Φ as in Lemma 2.3. Therefore, the first equation of (1.1) changes to: Thus, we consider the following linear uncoupled numerical scheme for (3.1) jointly with (1. h ∈ N h in a decoupled way (first Φ, then u and finally N ) satisfying We have denoted and similarly for δ t N k+1 h and δ t Φ k+1 h . The approximation of the initial conditions are taken as where we consider for simplicity that T 0 , N 0 , Φ 0 ∈ C 0 Ω with u 0 = e −χ Φ 0 T 0 .
Finally, the functions f i k h for i = 1, 2, 3 in (3.3), (3.4) and (3.5), have the following definitions: The functions P k h , S k h , R k h and Q k h in (3.7)-(3.9), are the corresponding dimensionless factors

Proof of Theorem 1.2
In this part, we are going to get a priori energy estimates for the fully discrete solution u k+1 h , N k+1 h and Φ k+1 h (and hence, for T k+1 h ) of (3.3), (3.4) and (3.5) which are independent of (h, k).
The following result is based on the hypothesis of acute triangulations to get a discrete maximum principle, see [7]. In fact, we arrive at discrete version of Lemma 2.1.
Multiplying (3.5) by (Φ k+1 h (a)) − and using that Φ k h (a) ≥ 0, it holds that: (3.10) Indeed, using the form of f 3 k h given in (3.9), the following estimates hold Adding the last two inequalities, one has On the other hand, since in every node a ∈ N , due to the form of f 3 k h given in (3.9) the following estimates hold Thus, adding the last two inequalities, we obtain that  Choosing where we have used in the left hand side that and that in every node a ∈ N h , On the other hand, we can make the following decomposition in the diffusion term Hence, using that u k+1 is a positive function and that ∇ϕ a · ∇ϕ a ≤ 0 ∀a = a ∈ N h (owing to the hypothesis of acute triangulation), we deduce, .

(3.15)
Adding (3.15) in (3.14), it holds that On the other hand, by using that in every node a ∈ N h , due to the form of f 1 k h given in (3.7), the following estimates hold Then, adding the last two inequalities, we obtain that Finally, for (3.4) it is easy to obtain that In addition, f 2 k h (a) ≥ 0 in every node a ∈ N h due to the form of f 2 k h given in (3.8). Hence,

Adimensionalization
Here, we simplify the number of the parameters of (1.1) and present the simulations according to the dimensionless parameters. For that, we consider as one possible example of the dimensionless factors P (Φ, T ), S (Φ, T ), Q (Φ, T ) and R (Φ, T ) appearing in (1.4) satisfying the hypotheses (1.12)-(1.16), the following ones: and These factors P , S, Q and R satisfy the conditions (1.5)-(1.9). To verify (1.12), observe that Moreover, in (4.2) and (4.3) we consider a regularization in the denominator since without this regularization, the partial derivatives of (4.2) and (4.3) degenerate in (0, 0).
For the adimensionalization, we start studying the carrying capacity parameter K > 0. We consider the change of variables T = T K , N = N K and Φ = Φ K passing the normalized capacity equal to 1.
Now, we consider the diffusion parameter ν and the tumor proliferation parameter ρ. We know that ρ is related to the time variable while ν is related to the spatial variable. Thus, we can make the following change of the independent variables: (4.5) Applying these changes in (1.1), it holds that (4.7) Hence, we obtain the following dimensionless parameters: Thus, we reduced three parameter from the original model (1.1): ρ, ν and K.
Finally, the adimensionalizated system is the following: (4.8)

Numerical Simulations
In this section, we will show some numerical simulations in order to detect which parameters of (4.8) are more important in the behaviour of the ring width between necrosis and tumor and the regular or irregular growth of the surface of a GBM.
For the numerical simulations we will use the uncoupled and linear fully discrete scheme defined in (3.3)-(3.5) by means of an Implicit-Explicit (IMEX) Finite Difference in time approximation and P 1 continuous finite element with "mass-lumping" in space.
We will use the computational domain, Ω = (−9, 9) × (−9, 9), the final time, T f = 500, the structured triangulation, {T h } h>0 of Ω such that Ω = We consider along the work necrosis zero initially and initial tumor given by: For the vasculature, we will take different initial conditions depending on the kind of tumor growth considered.

Ring width
Here, we present some numerical simulations according to the tumor-ring. Based on the study [28], we know that tumors with a thick tumor ring have the worst prognosis as we can see in the  In order to measure different rings, we will compare the density of tumor with respect to the density of tumor and necrosis. In every simulation, we will change the value of one parameter and testing how the tumor growth changes.
Since the subjects of study are tumor and necrosis, we change the parameters of the tumor and necrosis equations, these are, κ and α. Then, in all the simulations the value of γ and δ are fixed (see Table 3).

Tumor Ring quotient
We will start studying the ratio between proliferative tumor density, Ω T dx and total tumor density, Ω (T + N ) dx and we consider the different values of κ and α given in Table 4. In fact, we compute the following "ring quotient" (RQ) coefficient: Thus, if RQ is near to zero, there exists a high density of necrosis (which implies slim tumor ring) whereas if RQ is close to one, there is not enough necrosis in comparison with proliferative tumor density (which means thick tumor ring).  Hence, the best configurations to obtain a slim (resp. thick) ring would be choose a big (resp. small) α.

Density tumor growth
In Figure 4, we compute the total tumor Ω (T + N ) dx for the values of κ and α given in Table 4.  We can see in Figure 4 how the variation in the parameters κ and α produces changes in the total tumor density. In fact, the total tumor decreases with respect to κ and α.
Therefore, we conclude that α is the most important parameter in order to change the tumor ring and both κ and α have relevance to change the total density in the tumor growth.

Regularity surface
In this case, we will test if our model (4.8) can develop different regularities for the tumor surfaces. Now, we will base our results on the study published in [30] where appears the following survival curve: From Figure 5, the authors conclude that tumors with a regular surface have better prognosis than tumors with irregular surface.
Along this Section, we simulate the tumor growth with the initial tumor defined in Figure 1, necrosis zero and the vasculature distributed in different zones as in Figure 6: Thus, the question is if the chemotaxis term (of tumor going to the vasculature) implies tumor growth with regular or irregular surface. We remember that the chemotaxis term in (4.8) is defined by κ ∇ · (T ∇Φ) with κ > 0. Now, we want to detect which parameter is more relevant changing the regularity of the tumor surface, showing some simulations in which we move the value of one of them and observe how the tumor changes. For this, it is important the interaction between tumor and vasculature. Then, we will move the parameters which appear in tumor and vasculature equations, κ, α, γ and δ.
For these parameters we take the values of

Regularity Surface quotient
The pictures of Figure 7 show the quotient between the area occupied by the total tumor (tumor and necrosis) and the area of ratio the smallest circle containing the tumor. Thus, we present these computations for the different values of κ, α, γ and δ chosen in Table 5. In fact, we compute the following "surface quotient" (SQ) coefficient: where (T + N ) min and R max are defined as follows: Thus, we will deduce that if SQ is near to zero, tumor will have an irregular surface whereas if SQ is close to one, tumor will have a regular surface.   Remark 5.1. By the size of mesh considered, at the beginning of the pictures given in Figure 7, the value of SQ is larger than 1 and it is observed oscillations in the graphs of SQ. Indeed, if we consider a mesh size smaller, these initial values of SQ and the oscillations can be reduced. In order to check this, we show an example of SQ versus time for different κ considering a mesh size smaller: However, we think that it is not necessary the use of a mesh size smaller since we obtain the same behaviour (in average) for the mesh considered initially and with this mesh, we reduce the computational time.
We see in Figs 7a − 7d how our model differentiates two kinds of tumor growth changing the parameter κ, see Figure 7a, and with lower variation for α, see Figure 7b. On the other hand, we do not appreciate changes in the variation of parameters γ and δ for the irregularity of tumor growth as we see in Figs 7c − 7d.

Area
Once we have identified that the more important parameters for the regularity surface are firstly κ and later α, we measure the area of total tumor for these parameters as in Table 5:  We see in Figure 9 how the largest area corresponds to the smallest α = 10 and the smallest area holds for the highest α = 100. In the case of variation of κ, Figure 9a, a similar influence in the total tumor area for κ = 1 and κ = 10 is observed.
Thus, we have obtained a higher variation of total area for the different values of α than for κ, see Figure 9. Nevertheless, in the simulation of the "surface quotient" (SQ), we obtained more variation between the different values of κ that for the different values of α, see Figure 7. Hence, the factor which modifies this change is R max , defined by 5.4. In fact, R max will change more with the variation of κ than for the variation of α.

Tumor growth
Here, we examine the tumor growth for κ = 10 in five times step in order to see the variation in space of tumor. For this growth, the rest of parameters take the fixed values showed in Table 5. We observe an irregular tumor growth for κ = 10 when time increases. These results are in concordance with Figure 7a, where we observed a great irregularity for κ = 10 and with Figure 9a, where the area of the tumor for κ = 10 is increasing.
Finally, we conclude that κ is the more relevant parameters in the irregular surface of tumor and α is the most important parameter for total area in the tumor growth.

Discussion
Summarizing the results obtained with respect to the ring width and the regularity surface for the chemotactic and dimensionless system (4.8) related to GBM growth model, we deduce that this model can capture these two properties varying some parameters. Moreover, we have proved that the parameters more relevant according to the tumor growth are κ and α.
For the tumor ring, where the vasculature is uniformly distributed, the results show that the hypoxia parameter α is the most relevant coefficient as we can observe in Figs 3a − 3b.
In the case of regularity surface, where the vasculature is non-uniformly distributed, the parameter which produces more irregularity in the tumor surface is the chemotaxis parameter κ, see Figs 7a − 7d.