# 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.

There are no comments yet.

## Authors

• 2 publications
• 3 publications
• 2 publications
12/23/2020

### Theoretical and numerical analysis for a hybrid tumor model with diffusion depending on vasculature

In this work we analyse a PDE-ODE problem modelling the evolution of a G...
07/19/2020

### Convergence and positivity of finite element methods for a haptotaxis model of tumoral invasion

In this paper, we consider a mathematical model for the invasion of host...
07/27/2020

### Fractional semilinear optimal control: optimality conditions, convergence, and error analysis

We analyze an optimal control problem for a fractional semilinear PDE; c...
08/26/2019

### Well-posedness and discrete analysis for advection-diffusion-reaction in poroelastic media

We analyse a PDE system modelling poromechanical processes (formulated i...
05/22/2020

### A short note on plain convergence of adaptive least-squares finite element methods

We show that adaptive least-squares finite element methods driven by the...
12/09/2019

### Convective transport in nanofluids: regularity of solutions and error estimates for finite element approximations

We study the stationary version of a thermodynamically consistent varian...
12/17/2015

### Kauffman's adjacent possible in word order evolution

Word order evolution has been hypothesized to be constrained by a word o...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## Abstract

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.

Mathematics Subject Classification.
Keywords: Glioblastoma, Chemotaxis, PDE-ODE system, Numerical scheme.

The authors were supported by PGC2018-098308-B-I00 (MCI/AEI/FEDER, UE).

## 1 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 [19]. 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, 24] 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 [17, 21, 22, 23]). Specifically, Molab

group classifies the GBM depending on the width of the tumor ring and/or the tumor surface regularity (see

[21, 23] respectively). The study of [21] concludes that tumors with slim ring have better prognostic, specifically months of more survival than tumors with thick ring. In [23], 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 moths of survival, than tumor with irregular surface.

In [31], 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 [20] where the tumor ring is studied by a PDE-ODE system of two equations (proliferative tumor and necrosis). In [9, 10], 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) but, with chemotaxis term (tumor directed to vasculature) and we study the biological behaviours of the GBM such as the tumor ring volume, studied in [20, 21], and the regularity surface considered in [23], and compare the results with those in [9, 10].

Some previous chemotactic PDE-ODE models have been extensively studied in the literature, see for instance [5, 25, 26, 27] where the authors model the cells movement with a parabolic-ODE system. Specifically, in [27] 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 [25], some analytical and numerical results which support the numerical observations of [27] are presented using a similar model than in [27]. Moreover, in [3, 5] a model of tumor inducing angiogenesis is proposed consisting of a equation with chemotaxis and haptotaxis term, and two nonlinear ODEs. Finally, in [26] 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 [28, 29, 30] have shown existence results for systems of three differential equations modelling cancer invasion. In [28] 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 [29], 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 [30].

Recently, a PDE-ODE model with chemotaxis is studied in [13] obtaining asymptotic stability results using a proper transformation and energy estimates. Another PDE-ODE with chemotaxis problem is considered in [18] 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 this paper, we investigate the following parabolic PDE-ODE system in ( is a bounded and regular domain and corresponds to the final time)

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂T∂t−νΔTDiffusion+κ∇⋅(T∇Φ)Chemotaxis=f1(T,N,Φ)\par∂N∂t=f2(T,Φ)\par\par∂Φ∂t=f3(T,N,Φ)\par (1.1)

endowed with non-flux boundary condition on the boundary

 (−ν∇T+κT∇Φ)⋅n=0 (1.2)

where

is the outward unit normal vector to

and initial conditions at time :

 T(0,⋅)=T0,N(0,⋅)=N0,Φ(0,⋅)=Φ0inΩ. (1.3)

Here, and represent the tumor and necrotic densities and the vasculature concentration at the point and time , respectively.

The nonlinear reactions functions for have the following form

 (1.4)

The parameters in have the following description [12, 15, 16]:

The functions , , and appearing in are adimensional factors with the following biological meaning:

1. The tumor growth cells need space and a well amount of nutrients to grow. If this amount of nutrients per cell is suitable, the proliferation of tumor cells will occur. Hence, we introduce the tumor proliferation factor in as a volume fraction of the vasculature.

2. We consider the hypoxia as a decreasing term due to lack of vasculature. Hence, low vasculature produces more tumor destruction. Therefore, the factor must be a volume fraction of the lack of vasculature.

3. The vasculature growth factor will depend on the amount of tumor and the vasculature does not grow without tumor. Thus, will be a volume fraction of tumor.

4. The destruction of vasculature will increase with tumor and there will not be vascular destruction without tumor. In consequence, will be a volume fraction of tumor.

Thus, these factor functions , , and must satisfy the following modelling conditions:

 0≤P(Φ,T),S(Φ,T),Q(Φ,T),R(Φ,T)≤1∀(T,Φ)∈R2, (1.5)

and,

 P(Φ,T)=0 for Φ=0 and P(Φ,T)% increases if Φ increases, (1.6)
 S(Φ,T) increases if Φ decreases, (1.7)
 R(Φ,T)=0 for T=0 and R(Φ,T) % increases if T increases (at least for T≤K), (1.8)
 Q(Φ,T)=0 for T=0 and Q(Φ,T) % increases if T increases. (1.9)

We assume along the paper the following assumptions on the initial data

 0≤T0(x),N0(x),Φ0(x)≤K, a.e.x∈Ω. (1.10)

In order to obtain some estimates of the solutions of - (see ), we define the following truncated system of :

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂T∂t−νΔT+κ∇⋅(T+∇Φ)=f1(T+,N+,ΦK+)∂N∂t=f2(T+,Φ+)∂Φ∂t=f3(T+,N+,ΦK+) (1.11)

subject to and . We have denoted and and the same for and .

The main contributions of this work are the following:

1. ###### Theorem 1.1 (A priori estimates).
1. Any regular enough solution of the truncated problem - satisfies:

 0≤Φ≤K,T≥0andN≥0,a.e. in(0,Tf)×Ω

and

 T,N are bounded in L∞(0,Tf;L1(Ω)).
2. Assuming that there exists a constant such that

 C1P(Φ,T)≥R(Φ,T)Φ,∀0≤Φ≤K, and T≥0 (1.12)

and

 ρ≥κνγC1, (1.13)

then

 T,N are bounded in L∞(0,Tf;L∞(Ω)).
3. Assuming additionally that there exist constants for such that for all and ,

 ∣∣∂(R(Φ,T)Φ)∂Φ∣∣,∣∣∂(R(Φ,T)Φ)∂T∣∣≤C2, (1.14)
 ∣∣∂(Q(Φ,T)Φ)∂Φ∣∣,∣∣∂(Q(Φ,T)Φ)∂T∣∣≤C3 (1.15)

and

 ∣∣∂(S(Φ,T)T)∂Φ∣∣,∣∣∂(S(Φ,T)T)∂T∣∣≤C4, (1.16)

then

and

 ∇T is bounded in L2(0,Tf;L2(Ω)).

By Theorem 1.1 a), for any solution of , we deduce that , and and then, for and . Hence, we obtain the following crucial corollary:

###### Corollary 1.1.

If is a solution of the truncated problem , then is also a solution of - and satisfies the estimates of Theorem 1.1.

The existence of solutions of problem 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 as an approximation of where is a partition of the time interval and is the mesh size. To build the scheme, we will use the change of variable in the PDE equation with chemotaxis, , similar to the used in [7, 8, 14], in order to obtain an equivalent system with diffusion for the new variable .

###### Theorem 1.2 (Discrete version of Theorem 1.1 a)).

Scheme - has a unique solution satisfying the first pointwise estimates of Theorem 1.1 a), these are:

 0≤Φkh≤K,Tkh≥0andNkh≥0,inΩ. (1.17)

The design of a numerical scheme preserving the whole estimates of Theorem 1.1, and not only the estimates , remains as an open problem.

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 - and - 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.

## 2 A priori estimates of the solutions of (???)$−$(???)

### 2.1 Proof of Theorem 1.1 a)

###### Lemma 2.1.

Any solution of the truncated problem satisfy the following pointwise estimates:

 0≤Φ≤K,T≥0andN≥0,a.e. in(0,Tf)×Ω. (2.1)
###### Proof.

Let be a solution of . Since one can rewrite , multiplying the first equation of by and integrating in , we get

 12ddt∫Ω(T−)2dx+ν∫Ω∣∇T−∣2=∫ΩT−T+˜f1(T+,N+,ΦK+)dx=0, a.e. in(0,Tf).

Hence, since , then a.e. . We repeat the same argument for the other two equations of using now that

 Φ−f3(T+,N+,ΦK+)=0andN−f2(T+,Φ+)≤0.

To obtain the upper bound , we multiply the third equation of by and integrate in ,

Since , then . As , then a.e. .

###### Lemma 2.2.

Any solution of satisfies the estimates:

 (2.2)

 ∥N∥L∞(0,Tf;L1(Ω))≤C(ρ,α,δ,K,|Ω|,Tf). (2.3)
###### Proof.

Let be a solution of . Integrating in the first equation of and using that , we obtain that

 ddt∫ΩTdx=∫ΩρP(Φ,T)Tdx−∫ΩρP(Φ,T)T2Kdx−∫ΩρP(Φ,T)TN+ΦK≥0dx−
 −∫ΩαS(Φ,T)T≥0dx≤∫ΩρP(Φ,T)Tdx−1K∫ΩρP(Φ,T)T2dx.

Thus,

 ddt∫ΩTdx+1K∫ΩρP(Φ,T)T2dx≤∫ΩρP(Φ,T)Tdx.

Rewriting and applying Young’s inequality for the right side, we get,

 ddt∫ΩTdx+1K∫ΩρP(Φ,T)T2dx≤ρ(12K∫ΩP(Φ,T)T2dx+K2∫ΩP(Φ,T)).

Hence, using that , we conclude that

 ddt∫ΩTdx+ρ2K∫ΩP(Φ,T)T2dx≤ρK2|Ω|.

Integrating in for , we obtain that

 ∥T(t,⋅)∥L1(Ω)+ρ2K∫t0∫ΩP(Φ,T)T2dxdt≤TfρK2|Ω|,∀t∈(0,Tf)

whence we deduce .

To prove , we integrate the second equation of in , with ,

 ∥N(t,⋅)∥L1(Ω)≤α∫t0∫ΩTdxdt+δ∫t0∫ΩΦdxdt

where we have used . Thus, using that and the bound obtained for in , we get .

### 2.2 Proof of Theorem 1.1 b)

In order to obtain the estimate for , 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:

 w=log(T)−χΦ⇔T=eweχΦ=eχΦu (2.4)

with and .

Thus, the first equation of changes to

 (eχΦu)t−ν∇⋅(eχΦ∇u)=f1(eχΦu,N,Φ) (2.5)

and the boundary condition to

 ∇u⋅n=0. (2.6)
###### Lemma 2.3 (Proof of Theorem 1.1 b)).

Assume and . Then, given any solution of , it holds that is bounded in and is bounded in . Moreover, and are bounded in .

###### Proof.

To obtain the estimates for and , taking into account the estimates for , it suffices that be . The proof of is based in estimates with an Alikakos’ argument. Let be a solution of . We multiply by (for any ), and analyse term by term:

• Time derivative term:

 (eχΦu)tup−1=χΦteχΦup+1peχΦ(up)t (2.7)

and the second term of the right side of can be expressed as

 1peχΦ(up)t=1p(eχΦup)t−χpeχΦupΦt. (2.8)

Hence, from and ,

 (eχΦu)tup−1=1p(eχΦup)t+p−1pχΦteχΦup. (2.9)
• Nonlinear diffusion term:

 −ν∇⋅(eχΦ∇u)up−1=−ν∇⋅(eχΦ(∇u)up−1)+νeχΦ(p−1)up−2∣∇u∣2=−ν∇⋅(eχΦ(∇u)up−1)+νeχΦ(p−1)4p2∣∇(up/2)∣2. (2.10)
• Reaction term:

 f1(eχΦu,N,Φ)up−1=ρP(Φ,T)eχΦup(1−eχΦu+N+ΦK)−αS(Φ,T)eχΦup. (2.11)

Rewriting in the function as and adding , and , we get:

 1p(eχΦup)t−ν∇⋅(eχΦ(∇u)up−1)+νeχΦ(p−1)4p2∣∇(up/2)∣2+αS(Φ,T)eχΦup+(ρP(Φ,T)−(p−1p)χγR(Φ,T)Φ)eχΦup(eχΦu+N+ΦK)=(ρP(Φ,T)−(p−1p)χγR(Φ,T)Φ)eχΦup+χpeχΦupδQ(Φ,T)Φ. (2.12)

Due to hypothesis and , it is easy to see in that,

 ρP(Φ,T)−(p−1p)χγR(Φ,T)Φ≥0.

Using now that , and we obtain that

 1p(eχΦup)t−ν∇⋅(eχΦ(∇u)up−1)+νeχΦ(p−1)4p2∣∇(up/2)∣2+αS(Φ,T)eχΦup≤CeχΦup (2.13)

with . Integrating in , it holds that

 1pddt∫ΩeχΦupdx+ν(p−1)4p2∫ΩeχΦ∣∇(up/2)∣2dx+α∫ΩS(Φ,T)eχΦupdx≤C∫ΩeχΦupdx (2.14)

with independent of (along the proof, we will denote by different constants independent of ).

Using the auxiliary variable , we can rewrite as follows

 1pddt∥eχΦ2w∥2L2(Ω)+4ν(p−1)p2∥eχΦ2∇w∥2L2(Ω)≤C∥eχΦ2w∥2L2(Ω). (2.15)

Thus, applying Gronwall’s lemma, we deduce for that

Now, using the following equivalent norms with constants independent of

 ∥z∥2L2(Ω)≤∥eχΦ2z∥2L2(Ω)≤eχK∥z∥2L2(Ω), (2.16)

multiplying by and using that for any , we obtain that

 ddt∥eχΦ2w∥2L2(Ω)+2ν∥∇w∥2L2(Ω)≤Cp∥w∥2L2(Ω). (2.17)

We are going to apply the following Gagliardo-Nirenberg interpolation inequality (

[11, Theorem 10.1])

 (2.18)

with and the dimension of (in this case ). Applying for in the right hand side of , we deduce that

 ddt∥eχΦ2w∥2L2(Ω)+ν∥∇w∥2L2(Ω)≤Cp2∥w∥2L1(Ω). (2.19)

Using in but now for , it holds that

 (2.20)

Finally, due to , we can deduce that

 ddt∥eχΦ2w∥2L2(Ω)+C1∥eχΦ2w∥2L2(Ω)≤C(p2+1)∥w∥2L1(Ω) (2.21)

where .

Hence, we obtain that

 (2.22)

Following a similar argument to used by Alikakos in [2] (see Appendix), from we can obtain that

 u is bounded in L∞(0,Tf;L∞(Ω)).

As consequence, is bounded in .

Since and and are bounded in we obtain that is bounded in . ∎

### 2.3 Proof of Theorem 1.1 c)

Let be a solution of . Taking gradient in the second and third equation of ,

 (∇Φ)t=γ[(∂(R(Φ,T)Φ)∂Φ∇Φ+∂(R(Φ,T)Φ)∂T∇T)(1−T+N+ΦK)−−R(Φ,T)ΦK(∇T+∇N+∇Φ)]−δ(∂(Q(Φ,T)Φ)∂Φ∇Φ++∂(Q(Φ,T)Φ)∂T∇T), (2.23)
 (∇N)t=α(∂(S(Φ,T)T)∂Φ∇Φ+∂(S(Φ,T)T)∂T∇T)++δ(∂(Q(Φ,T)Φ)∂Φ∇Φ+∂(Q(Φ,T)Φ)∂T∇T). (2.24)

Using the change of variable as in Lemma 2.3, we deduce that

 ∇T=χeχΦu∇Φ+eχΦ∇u=χT∇Φ+eχΦ∇u (2.25)

and we know from Lemma 2.3 that is bounded in . Taking into account that and are bounded in , it holds that

 |∇T|≤C(|∇Φ|+|∇u|).

Thus, rewriting and in terms of , multiplying and by and respectively and integrating in , we deduce

 12ddt∥∇Φ∥2L2(Ω)≤C1∥∇Φ∥2L2(Ω)+C2∫Ω|∇u||∇Φ|dx+C3∫Ω|∇N||∇Φ|dx, (2.26)

and

 12ddt∥∇N∥2L2(Ω)≤C4∫Ω|∇Φ||∇N|dx+C5∫Ω|∇u||∇N|dx, (2.27)

with for . In and we have applied the inequality

 ∫Ωv|∇u||∇Φ|dx≤∥v∥L∞(Ω)∫Ω|∇u||∇Φ|dx

with since , and are bounded in .

Using now Cauchy-Schwarz and Young’s inequalities in and and adding them, it holds that

 12ddt(∥∇Φ∥2L2(Ω)+∥∇N∥2L2(Ω))≤ˆC1(∥∇Φ∥2L2(Ω)+∥∇N∥2L2(Ω))+ˆC2∥∇u∥2L2(Ω), (2.28)

with for . Since is bounded in , applying Gronwall’s Lemma, it holds that

 ∇Nand∇Φare bounded in L∞(0,Tf;L2(Ω)).

Finally, using in and , we obtain that

 (∇N)tand(∇Φ)t% are bounded in L2(0,Tf;L2(Ω)).

is bonded in .

## 3 A FE numerical scheme

In this Section, we are going to design an uncoupled and linear fully discrete scheme to approach - by means of an Implicit-Explicit (IMEX) Finite Difference in time and 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.

1. Let . We consider the uniform time partition

 (0,Tf]=Kf−1⋃k=0(tk,tk+1],

with where and is the time step. Let or a bounded domain with polygonal or polyhedral lipschitz-continuous boundary.

2. Let be a family of shape-regular, quasi-uniform triangulations of formed by acute N-simplexes (triangles in D and tetrahedral in D with all angles lowers than ), such that

where , with being the diameter of . We denote the set of all the nodes of .

3. Conforming piecewise linear, finite element spaces associated to are assumed for approximating :

 Nh={nh∈C0(¯¯¯¯Ω):nh|K∈P1(K),∀K∈Th}

and its Lagrange basis is denoted by .

Let be the nodal interpolation operator and consider the discrete inner product

 (nh,¯¯¯nh)h=∫ΩIh(nh⋅¯¯¯