A flux-limited model for glioma patterning with hypoxia-induced angiogenesis
AA flux-limited model for glioma patterning withhypoxia-induced angiogenesis
Pawan Kumar and Christina SurulescuTU Kaiserslautern, Felix-Klein-Zentrum f ¨ur Mathematik,Paul-Ehrlich-Str. 31, 67663 Kaiserslautern, Germany (kumar, [email protected])
September 22, 2020
Abstract
We propose a model for glioma patterns in a microlocal tumor environment under theinfluence of acidity, angiogenesis, and tissue anisotropy. The bottom-up model deductioneventually leads to a system of reaction-diffusion-taxis equations for glioma andendothelial cell population densities, of which the former infers flux limitation both in theself-diffusion and taxis terms. The model extends a recently introduced [34] descriptionof glioma pseudopalisade formation, with the aim of studying the effect ofhypoxia-induced tumor vascularization on the establishment and maintenance of thesehistological patterns which are typical for high grade brain cancer. Numerical simulationsof the population level dynamics are performed to investigate several model scenarioscontaining this and further effects.
Keywords:
Glioblastoma, pseudopalisades, hypoxia-induced angiogenesis, flux-limited diffusionand taxis, kinetic transport equations, equilibrium-based moment closure
Classification of glioma, the most common type of primary brain tumors, typically comprises fourgrades, according to the degree of malignancy [32, 39]. The highest grade IV (including the mostaggressive type glioblastoma multiforme, shortly GBM) corresponds to characteristic histologicalpatterns called pseudopalisades; they exhibit garland-like hypercellular structures surrounding necroticregions usually centered around one or several sites of vasoocclusion [10, 11, 43, 52], and are associatedwith poor patient survival prognosis [32].More specifically, vascular occlusion and thrombosis [10, 43] result into hypoxia, which promotesmigration of glioma cells away from the highly acidic area, as well as tissue necrotization therein.Hypoxia around the pseudopalisading cells triggers the process of capillary formation, which involves acascade of coordinated steps [37, 47]. It typically starts with tumor cells overexpressinghypoxia-inducible regulators of angiogenesis, including vascular endothelial growth factor a r X i v : . [ q - b i o . T O ] S e p VEGF) [11]. The latter acts as a chemoattractant and proliferation promotor for endothelial cells(ECs) lining the surrounding blood vessels. The chemotactic migration of ECs towards VEGFgradients is accompanied by formation of new capillaries off existing ones, thus leading tomicrovascular hyperplasia which is common in GBM [9, 11] and which facilitates the growth of theneoplasm [11]. Understanding the mechanisms of tumor malignancy can contribute towards developingand tuning therapeutic strategies, e.g. by targeting angiogenesis [4] or by tumor alkalization [53], usedin a (neo)adjuvant way to other approaches (such as surgery and radiotherapy).Among the few mathematical models explicitly addressing glioma pseudopalisade formation (we referto [12, 13] for agent-based approaches and to [3, 34, 36] for continuous descriptions), the latterproposed a multiscale approach to deducing a reaction-(myopic)diffusion equation with repellentpH-taxis for glioma cell density coupled with a reaction-diffusion PDE for the acidity (expressed asproton concentration) produced by tumor cells. The derivation of this system started with describingthe evolution of the glioma distribution function on the mesoscopic scale, in the KTAP (kinetic theoryof active particles) framework [5]; thereby, the activity variable represented the amount of gliomatransmembrane units occupied with protons. A parabolic scaling led to the announcedreaction-diffusion-taxis (RDT) equation on the so-called macroscopic level of the glioma population ,including in its terms influences of acidity and tissue coming from the lower modeling scales. The twoequations for glioma and acidity were capable to qualitatively reproduce pseudopalisade-like patterns.In this work we aim to use a related approach and to extend the model in [34] such as to account forhypoxia-triggered angiogenesis and its effects on the pseudopalisade patterns. We thereby continue inthe line of previous works [14, 16–18, 21, 23, 24, 29] dealing with cancer cell migration inheterogeneous tissues.The rest of the paper is structured as follows: Section 2 is dedicated to setting up the micro-mesoformulation in the KTAP framework and to deduce a macroscopic system of differential equationscharacterizing the dynamics of glioma density (RDT) with flux-limited diffusion and repellentpH-taxis, EC density (RDT) with linear diffusion and chemotaxis towards VEGF gradients,concentrations of acidity and VEGF (both with linear diffusion and tumor- and EC-dependent sourceterms). In Section 3 we perform numerical simulations for various scenarios of the macroscopic modelobtained in Section 2 and interpret the results. Finally, Section 4 provides a discussion and an outlookon further problems of interest in this context. We begin by introducing some notations for the various variables and functions involved in our modelingprocess: • variables: time t ≥
0, position x ∈ R N , velocity of glioma cells v ∈ V = s S N − ⊂ R N , velocity ofECs ϑ ∈ Θ = σ S N − ⊂ R N ; we assume the cell speeds s , σ > S N − denotes theunit sphere in R N ; • ˆ v = v | v | , ˆ ϑ = ϑ | ϑ | are unit vectors denoting the directions of vectors v ∈ V and ϑ ∈ Θ , respectively; • p ( t , x , v ) : (mesoscopic) density function of glioma cells and M ( t , x ) = (cid:82) V p ( t , x , v ) d v :macroscopic density of glioma; We will use in the sequel the term ’macroscopic’ to designate quantities only depending on time and space,not to be confounded with the true biological size scale on which pseudopalisades are observed and which is -inthat acception- microscopic. Thus, in the following ’microscopic’ means single-cell level, while ’mesoscopic’refers to distributions of cells depending on time, space, and further kinetic and/or activity variables. w ( t , x , ϑ ) : (mesoscopic) density function of ECs and W ( t , x ) = (cid:82) Θ w ( t , x , ϑ ) d ϑ : macroscopicdensity of ECs. • q ( x , θ ) : (known) directional distribution function of tissue fibers with orientation θ ∈ S N − . Itholds that (cid:82) S N − q ( x , θ ) d θ =
1; for ω = (cid:82) V q ( x , ˆ v ) d v = s N − it holds then that (cid:82) V q ( x , ˆ v ) ω d v = • E q ( x ) : = (cid:82) S N − θ q ( x , θ ) d θ : mean fiber orientation. We also denote ˜ E q ( x ) : = (cid:82) V v q ( x , ˆ v ) ω d v and callit mean fiber direction; • V q ( x ) : = (cid:82) S N − ( θ − E q ( x )) ⊗ ( θ − E q ( x )) q ( x , θ ) d θ : variance-covariance matrix for orientationdistribution of tissue fibers, and, correspondingly, ˜ V q ( x ) = (cid:82) V ( v − ˜ E q ( x )) ⊗ ( v − ˜ E q ( x )) q ( x , ˆ v ) ω d v ; • h ( t , x ) : concentration of protons (acidity), a macroscopic quantity; • g ( t , x ) : concentration of VEGF, also a macroscopic quantity; • λ >
0: constant glioma turning rate; • η ( ϑ , g ( t , x )) >
0: EC turning rate.
In the spirit of [17] we consider single cell (microscale) dynamics in the form d x dt = v , d v dt = − α S ( v , h , M ) , (2.1)with S ( v , h , M ) = B ( v ) γ ∇ h (cid:113) K h / X + | ∇ h | + γ ∇ M (cid:113) K M / X + | ∇ M | = : B ( v ) φ ( h , M ) , (2.2)where γ , γ ≥ K h , K M , X > ), and the tensor B ( v ) : = | v | I N − v ⊗ v models biomechanical cellstress; thereby v ⊗ v represents active cell stress, while −| v | I N is the isotropic part. Hence, similarlyto [14, 17, 21], the velocity (and in particular the direction) of glioma cells is also influenced by aweighted sum of gradients (each inferring some limitation): increasing gradients of protonconcentration have a repelling effect, and the cells also try to avoid too crowded regions. The constant α > h and M ; for further details we referto [17]. Notice that (2.1) actually implies d | v | dt =
0, which is in line with our assumption of s = | v | being constant.Correspondingly, on the mesoscale we have the following kinetic transport equation (KTE) for theevolution of glioma cell distribution function p ( t , x , v ) : ∇ ( t , x , v ) · (( , v , − α S ( v , h , M )) p ) = p t + ∇ x · ( v p ) − α ∇ v · ( S ( v , h , M ) p ) = L [ λ ] p + P ( M , h , W ) , (2.3)where L [ λ ] p ( t , x , v ) = − λ p ( t , x , v ) + λ (cid:82) V K ( v , v (cid:48) ) p ( t , x , v (cid:48) ) d v (cid:48) is a turning operator with constantturning rate λ > K . As in [28] and many subsequent works [16–18, 23, 24, 34] wechoose K ( v , v (cid:48) ) : = q ( x , θ ) ω , meaning that the reorientation of cells is due to contact guidance by tissue, with X to be chosen in Paragraph 3.1 ith its orientation distribution q . Hence, L [ λ ] p ( t , x , v ) = λ ( M q ( x , ˆ v ) ω − p ( t , x , v )) . Moreover, as e.g.,in [17] we assume p to be compactly supported in the velocity space.The term P ( M , h , W ) models glioma cell proliferation, which is supposed to depend on the amount ofcells already present (irrespective of their orientation) and on favorable (nutrients provided by ECs ofdensity W ) and unfavorable (acidity h ) factors in the environment. More detailed descriptions accountfor proliferative interactions mediated by kernels depending on kinetic or activity variables, as in [17,24]or even differentiate between moving and proliferating cells [25, 30]. Here we would rather focus on themacrolevel influences and choose a uniform velocity kernel to characterize such interactions, i.e. P ( M , h , W ) = (cid:18) − MK M (cid:19) (cid:90) V (cid:18) µ (cid:18) − hK h (cid:19) + µ W / K W + W / K W (cid:19) | V | p ( t , x , v (cid:48) ) d v (cid:48) = µ | V | M (cid:18) − MK M (cid:19) (cid:18) − hK h (cid:19) + µ | V | M W / K W + W / K W , (2.4)where µ , µ > K M and K W represent carrying capacities of glioma and endothelial cells, respectively, while K h > M from the KTEformulation, in [34] was used a parabolic upscaling ; the microscale dynamics in that setting focusedon the evolution of an activity variable quantifying the amount of transmembrane units occupied byprotons. In [17], where velocity dynamics akin to (2.1) was considered (however without limitations asin the denominators of S ), it was yet a parabolic scaling leading to the reaction-diffusion-taxis PDE forglioma evolution of the macroscale. Instead, we will consider here a moment closure approach similarto that in [14, 29], ensuring the closure by way of an equilibrium distribution.For that purpose we introduce the following quantities: ( M U )( t , x ) : = (cid:90) V v p ( t , x , v ) d v , (2.5) P ( t , x ) : = (cid:90) V ( v − U ( t , x )) ⊗ ( v − U ( t , x )) p ( t , x , v ) d v , (2.6)commonly designated as momentum and pressure tensor, respectively. Thereby, U represents theso-called ensemble (macroscopic) cell velocity.Integrating (2.3) w.r.t. v we get ∂ t M + ∇ x · ( M U ) = µ M (cid:18) − MK M (cid:19) (cid:18) − hK h (cid:19) + µ M W / K W + W / K W (2.7)Multiplying (2.3) by v and integrating again w.r.t. v we get ∂ t ( M U ) + ∇ x · (cid:90) V v ⊗ v pd v + α (cid:90) V S ( v , h , M ) pd v = λ M ( ˜ E q − U ) (2.8)From (2.6), P = (cid:90) V v ⊗ v pd v − M U ⊗ U , with the observation that a hyperbolic one leads to a PDE which is too much drift-dominated and not able toreproduce pseudopalisade patterns ence due to (2.2) we can rewrite (2.8) as ∂ t ( M U ) + ∇ x · ( P + M U ⊗ U ) − α ( P + M U ⊗ U − s M I N ) φ ( h , M ) = λ M ( ˜ E q − U ) . (2.9)System (2.8), (2.9) is not closed, since the pressure tensor P involves a second order moment of p andis therefore unknown. In order to close the system we proceed as e.g., in [29] and assume that it is closeto equilibrium, which also dominates the second order moments. If we also assume that at equilibriumwe can neglect cell proliferation (as it indeed happens much slower than the motility-related processes),then in virtue of (2.3) we can write p equil = q ω M , (2.10)which also leads to U equil = M (cid:90) V v p equil d v = ˜ E q (2.11)and P equil = (cid:90) V ( v − ˜ E q ) ⊗ ( v − ˜ E q ) q ω Md v = M ˜ V q , (2.12)hence at equilibrium the momentum of glioma population is driven by the cell ensemble aligning theirmotion to the average fiber orientation ˜ E q , and the ’population pressure’ is generated by thevariance-covariance matrix of the orientation distribution q / ω of tissue surrounding the cells.Thus, we get from (2.9) ∂ t ( M U ) + ∇ x · ( M U ⊗ U ) = − ∇ x · P equil + α (cid:0) P equil + M U ⊗ U − s M I N (cid:1) φ ( h , M ) + λ ( M ˜ E q − M U ) . (2.13)The middle term on the right hand side of (2.13) stems from the description of forces acting on the cells;we suppose that the tensor M U ⊗ U therein is itself relating to the equilibrium situation, thus replace itby M ˜ E q ⊗ ˜ E q . On the other hand, the left side in (2.13) can be interpreted as directional derivative ofthe momentum M U in the direction U of the ’cell flow’. As in [29] we consider that the cell flux relaxesquickly to its equilibrium, so that M U = M ˜ E q − λ ∇ x · P equil + αλ (cid:0) P equil + M ˜ E q ⊗ ˜ E q − s M I N (cid:1) φ ( h , M ) . (2.14)If, furthermore, we assume that the tissue is undirected , this translates into ˜ E q = , which simplifies theabove expression of the momentum: M U = − ∇ x · ( D T ( x ) M ) + α ( D T ( x ) − s λ I N ) M φ ( h , M ) , (2.15)with D T ( x ) : = λ (cid:82) V v ⊗ v q ( x , ˆ v ) ω d v denoting as in previous works (see e.g. [23, 34]) the so-called tumordiffusion tensor.Plugging (2.15) into (2.7) then leads to the macroscopic PDE for glioma evolution ∂ t M = ∇∇ : ( D T M ) + α ∇ · (cid:18)(cid:18) s λ I N − D T (cid:19) φ ( h , M ) M (cid:19) + µ M (cid:18) − MK M (cid:19) (cid:18) − hK h (cid:19) + µ MW / K W + W / K W , (2.16) as suggested by the simulation results in [34] hich is of course merely an approximation, in view of the many assumptions made above.The myopic diffusion term in (2.16) is made of a drift and a ’classical’ diffusion term: ∇∇ : ( D T M ) = ∇ · ( ∇ · D T M ) + ∇ · ( D T ∇ M ) , hence using the expression of φ ( h , M ) we can rewrite (2.16) in the form of a PDE in which both self-diffusion (first term in (2.17)) and repellent pH-taxis (third term in (2.17)) are (at least partially) flux-limited: ∂ t M = ∇ · D T + αγ (cid:18) s λ I N − D T (cid:19) M (cid:113) K M / X + | ∇ M | ∇ M + ∇ · ( ∇ · D T M )+ ∇ · αγ (cid:18) s λ I N − D T (cid:19) M ∇ h (cid:113) K h / X + | ∇ h | + µ M (cid:18) − MK M (cid:19) (cid:18) − hK h (cid:19) + µ MW / K W + W / K W . (2.17)Notice that in order to have a genuine repellent pH-taxis the tensor s λ I N − D T should be positive definite,which also ensures a nondegenerating, proper diffusion . Under the hypothesis of undirected tissue, i.e.for E q = this amounts to I N − V q = I N − (cid:82) S N − θ ⊗ θ q ( θ ) d θ being positive definite. This is, forinstance, the case if q is the so-called peanut distribution [38]: q ( x , θ ) = N | S N − | trace D water ( x ) θ T D water ( x ) θ , where D water ( x ) denotes the (positive definite and usually diagonalized) water diffusion tensor. Sameapplies to our choice of q in Subsection 3.2.A much simplified situation is encountered when the stress tensor B ( v ) has only an isotropic part, whilethe active cell stress is neglected. The above deduction then directly yields instead of (2.17): ∂ t M = ∇ · D T + αγ s λ I N M (cid:113) K M / X + | ∇ M | ∇ M + ∇ · ( ∇ · D T M )+ ∇ · αγ s λ M ∇ h (cid:113) K h / X + | ∇ h | + µ M (cid:18) − MK M (cid:19) (cid:18) − hK h (cid:19) + µ MW / K W + W / K W . For yet another choice of the right hand side in (2.1) leading in a different, but related way to a flux-limited reaction-diffusion-taxis PDE for glioma evolving in a tissue (also under the influence of acidityand vascularization) we refer to [21].
We turn now to obtaining a macroscopic description for the evolution of ECs. For this purpose we canreproduce the derivation made in [16], the only difference being that ECs are supposed to respond hereto (gradients of) VEGF concentration rather than proliferating glioma cells. Therefore, we will onlysketch the main steps and refer for details to [16]. without being, however, a necessary condition therefor Indeed, the latter assumption would be inappropriate in the present context of glioma patterns, where most ofthe cells are migrating to avoid hypoxic regions, only starting to proliferate when there is enough vascularization. e start with the KTE formulation for mesoscopic dynamics of EC distribution function w ( t , x , ϑ ) : w t + ∇ x · ( ϑ w ) = L w [ η ] w + µ w ( W , g ) w , (2.18)with turning operator L w [ η ] = − η ( ϑ , g ) w ( t , x , ϑ ) + (cid:90) Θ | Θ | η ( ϑ (cid:48) , g ) w ( t , x , ϑ (cid:48) ) d ϑ (cid:48) , (2.19)involving a uniform turning kernel and the turning rate η ( ϑ , g ) = η e − a ( g ) D t g (2.20)depending upon the pathwise gradient D t g = g t + ϑ · ∇ g of VEGF concentration. The coefficientfunction a ( g ) takes into account the way in which ECs perceive VEGF and respond to it,correspondingly adapting their turning. Following [16] we choose a ( g ) = χ a K R ( K R + g / K g ) , with χ a , K R > y ∗ ( g ) = g / K g g / K g + K R representing theequilibrium of EC receptor binding dynamics to g ; the constant χ a scales changes in the turning rateper unit of change in dy ∗ / dt .The coefficient µ w ( W , g ) in the proliferation term is chosen as µ w ( W , g ) = µ W (cid:18) − WK W (cid:19) gK g , (2.21)with µ W , K W , K g > t (cid:32) ε t , x (cid:32) ε x performed in the usual way in combination with a Hilbertexpansion w = w + ε w + ε w + ... , a linearization of η , and the assumption that EC proliferation ismuch slower than migration then leads to the reaction-diffusion-chemotaxis PDE for the leading termin the (macroscopic) Hilbert expansion W = W + ε W + ε W + ... (for details compare [16]): ∂ t W − ∇ · ( D EC ∇ W ) + ∇ · ( χ a ( g ) W ∇ g ) = µ w ( W , g ) W , (2.22)where D EC ( x ) : = σ N η I N and χ a ( g ) : = σ a ( g ) N I N = η a ( g ) D EC ( x ) .In the sequel we will (formally) approximate the evolution of W by that of its leading term W and willuse the notation W instead of W . Equations (2.17) and (2.22) are coupled with the following equations set directly on the macroscopicscale and describing the dynamics of VEGF and acidity: g t = D g ∆ g + γ g M hK h M + K M − ζ g g WK W , (2.23) h t = D h ∆ h + γ h MK M + M − ζ h h WK W , (2.24)where, D i , γ i and ζ i , for i = h , g , are positive constants representing diffusion coefficients, production,and uptake rates for the involved quantities.The above PDE system needs to be supplied with initial conditions; they will be addressed in Section3. Furthermore, since we are interested in the solution behavior inside a bounded region (representinge.g., a slice of a histologic sample), we restrict our equations (hitherto written for x ∈ R N ) to a boundeddomain Ω ⊂ R N and consider no-flux boundary conditions for M , W , g , h . They are naturally assumedfor g , h and can be obtained (at least for W ) during the upscaling process, as detailed in [17]. Numerical simulations
This section is dedicated to numerical simulations of the macroscopic system obtained in Section 2 forthe quantities M , W , g , h . We consider the following rescalings of the hitherto dimensional quantities:˜ M : = MK M , ˜ h : = hK h , ˜ W : = WK W , ˜ g : = gK g , ˜ t : = µ t , ˜ x : = x (cid:114) µ D h , ˜ D T : = D T D h , ˜ D EC : = D EC D h , ˜ D g : = D g D h , ˜ D s = s λ D h , ˜ µ c : = µ µ , ˜ µ e : = µ W µ , ˜ χ W : = χ a η , ˜ γ g : = γ g µ K g , ˜ ζ g : = ζ g µ , ˜ γ h : = γ h µ K h , ˜ ζ h : = ζ h µ , ˜ φ ( ˜ h , ˜ M ) : = γ ∇ ˜ x ˜ h (cid:113) + | ∇ ˜ x ˜ h | + γ ∇ ˜ x ˜ M (cid:112) + | ∇ ˜ x ˜ M | , X = (cid:115) D h µ . After dropping the tildes for simplicity of notation we therewith get the following nondimensionalizedsystem: M t = ∇∇ : ( D T M ) + α ∇ · (( D s I N − D T ) φ ( h , M ) M ) + M ( − M ) ( − h ) + µ c MW + W (3.1a) W t = ∇ · ( D EC ∇ W ) − ∇ · (cid:18) χ W K R ( K R + g ) D EC W ∇ g (cid:19) + µ e gW ( − W ) (3.1b) h t = ∇ h + γ h M + M − ζ h hW (3.1c) g t = D g ∇ g + γ g Mh + M − ζ g gW , (3.1d)whose parameters are given by the above scalings and the values provided in Table 1. Typically the brain structure of each patient is assessed via diffusion tensor imaging (DTI) and thetherewith provided water diffusion tensor D water already occurring in Subsection 2.1 above in thepeanut distribution. The space scale on which glioma patterns are observed is, however, smaller thanthe best resolution of DTI, which does not go below the size of a voxel (ca. 1 mm ): pseudopalisadesare structures with a medium width of 200 − µ m ; for more details see the introduction of [34] andreferences therein. Therefore, as in that reference, we use an artificially created tissue with acorresponding water diffusion tensor as introduced in [38], in order to analyze the effect of tissueanisotropy. We consider the water diffusion tensor as D water ( x , y ) = (cid:18) . − d ( x , y )
00 0 . + d ( x , y ) , (cid:19) (3.2)where d ( x , y ) = . e − . ( x − ) − . e − . ( y − ) . A combination of uniform and von Mises-Fisher distributions is considered for the orientation distribution of tissue fibers: q ( x , θ ) = δ π + ( − δ ) (cid:18) π I ( k ( x )) (cid:19) e k ( x ) ϕ ( x ) · θ + e − k ( x ) ϕ ( x ) · θ ith δ ∈ [ , ] being a weighting coefficient, ϕ denotes the eigenvector corresponding to the leadingeigenvalue of D water ( x ) and I is the modified Bessel function of first kind of order 0. Furthermore, θ = ( cos ξ , sin ξ ) for ξ ∈ [ , π ] , and k ( x ) = κ FA ( x ) , where κ ≥ FA ( x ) = | λ − λ | (cid:113) λ + λ denotes the fractional anisotropy in 2D [38] with λ i ( i = ,
2) denoting the eigenvalues of D water ( x ) . Parameter Meaning Value Reference K M glioma carrying capacity 0 . − . µ m this work, [1] K h acidity threshold for cancer celldeath 10 − . mol/L [48] s speed of glioma cells 15 − µ m/h [20, 42] λ turning frequency coefficient 360 1/h [23, 45] γ h proton production rate 3 . · − mol /(L · h) this work, [35] ζ h proton removal rate 3 . · − − . · − /h this work D h acidity diffusion coefficient 1 . · − . · µ m /h this work, [34] µ glioma growth/decay rate, influ-enced by acidity ( . − . ) · − /h [22, 46] µ glioma growth rate, influenced byECs ( − . ) · − /h this work α advection constant 10 [17] γ weight coefficient related to acidicstress 3 · − − · − this work γ weight coefficient related to cellpopulation pressure 10 − − − this work σ speed of ECs 15 − µ m/h [16, 19] η turning rate of ECs 36 /h this work µ w EC proliferation rate ( . − . ) · − /h this work, [16] χ a coefficient of chemotactic sensi-tivity of ECs 0 . − . K W carrying capacity of ECs 0 . − . µ m this work, [2] D g diffusion coefficient of VEGF 3 . · − . · µ m /h [41], [44] K g VEGF threshold value 8 · − − − mol/L this work, [26] γ g VEGF production rate 3 . · − mol/(L · h) [41] ζ g VEGF uptake rate 3 . · − − . · − /h this work Table 1:
Parameters (dimensional quantities) .3 Initial conditions We consider the following initial densities and concentrations for the quantities involved: M ( , x ) = . · K M · (cid:18) e − ( x − ) − ( y − ) ( ) + e − ( x − ) − ( y − ) ( ) + e − ( x − ) − ( y − ) ( ) (cid:19) (3.4a) h ( , x ) = − e − ( x − ) − ( y − ) ( ) + − e − ( x − ) − ( y − ) ( ) + − . e − ( x − ) − ( y − ) ( . ) (3.4b) g ( , x ) = − e − ( x − ) − ( y − ) ( ) + − e − ( x − ) − ( y − ) ( ) + − e − ( x − ) − ( y − ) ( . ) (3.4c) W ( , x ) = . · K W · sin ( . π y ) (cid:18) e − ( x − ) . + e − ( x − ) . (cid:19) . (3.4d) (a) Tumor (b)
Acidity (c)
ECs (d)
VEGF
Figure 1:
Initial conditions
We solve system (3.1) with no-flux boundary conditions and with the initial data given above on asquare domain [ , ] × [ , ] (in µ m ) by using a finite difference method. The diffusion parts ofacidity , ECs and VEGF are discretized using a standard central difference scheme, while anonnegativity-preserving discretization scheme, as proposed in [49], is used for the anisotropicdiffusion of glioma cells. We apply a first order upwind scheme for the taxis terms in glioma and ECequations. For time discretization, an IMEX scheme is used for the acidity, EC, and VEFG equations,thereby treating the diffusion parts implicitly and discretizing the taxis and source terms with anexplicit Euler method. We use an explicit Euler scheme for each term in the glioma equation. We addressed in [34] the case with isotropic tissue ( δ =
1) as well as that with anisotropy ( δ = . k = Experiment 1
To start with, we solve the full system (3.1) with no-flux boundary conditions andwith initial data as in Subsection 3.3. The results are shown in Figure 2. The garland-shaped gliomapattern forms (and is well visible after 60 simulated days) due to the high proton concentration in itscentral region, in spite of ECs being attracted by the glioma producing VEGF. In fact, ECs spreadover almost the entire domain of simulation, coming from the left and right margins; the EC plot at here and subsequently we will mean by ’acidity’ the concentration of protons, keeping in mind its relationshipto the corresponding pH value; the latter is usually refer to ’acidity’ in medical practice =
60 days exhibits large densities (except in the narrow central region, where the cells did not yetarrived), a phenomenon known as hyperplasia. Subsequently, ECs have accumulated enough near thetumor to uptake acidity and release nutrients, which leads among others to the tumor cells being ableto move away from the acidic area and toward more favorable regions (due to γ (cid:29) γ the influenceof acidity dominates the correction of self-diffusion), where they start producing acidity and VEGF.The pseudopalisades develop larger diameters and involve glioma garlands with lower densities. Thecancer cells are followed by ECs, which again reduce the proton concentration, thus changing the gliomabehavior. The pseudopalisades are disrupted and at later times both ECs and glioma cells can repopulatethe previously too acidic regions. Experiment 2
Same as
Experiment 1 , but with γ = γ , i.e. the flux-limited self-diffusion andpH-taxis being equally weighted. The results are illustrated in Figure 3. The behavior of solutioncomponents is qualitatively similar to the previous scenario, except the pseudopalisade-like patternspersisting for a longer time, being succeeded by a more uniform tumor- and EC spread occupying thewhole domain (also migrating into the formerly most hypoxic areas) and exhibiting higher cell densities. Experiment 3
Here we want to assess the effects of extending the model in [34]: M t = ∇∇ : ( D T M ) + ∇ · ( G ( h ) M D T ∇ h ) + M ( − M ) ( − h ) (3.5a) h t = ∇ h + γ h M + M − ζ h h , (3.5b)where G ( h ) = λ k D ( h + k D ) ( h + k D + λ ) , λ , λ , k D > M t = ∇∇ : ( D T M ) + α ∇ · (( D s I N − D T ) φ ( h , M ) M ) + M ( − M ) ( − h ) (3.6a) h t = ∇ h + γ h M + M − ζ h h . (3.6b)for the right part of Figure 4.The plots on the left column exhibit enhanced tumor growth and spread in the presence of ECs andVEGF. Earlier stages feature larger glioma densities (and therewith associated proton concentrations)distributed in a ring-shaped manner, which suggests more pronounced pseudopalisades. At later times,in the presence of vascularization, the tumor cells are not only able to proliferate, but also to reoccupythe formerly too acidic areas near the center of the simulation domain, while still responding to thespatial distribution of protons. This leads (as observed also in the previous numerical experiments) topatterns no longer having the pseudopalisade-like appearance. Similar effects can be observed in theplots on the right column, although less pronounced, due to the smaller structural difference betweenthe two models compared therein. In fact, at later times there are only minimal differences between theplots on the left and on the right columns, which suggests that on the long run the effects of fluxlimitation in pH-taxis and additional self-diffusion from (3.6a) are increasingly diminished incomparison to the simple pH-taxis from (3.5), most probably because of the diffusion dominance in allthese models.Figure 7 shows differences between 1D tumor and acidity patterns for the model comparison donein the left column of Figure 4 and for two different choices of the diffusion coefficient D T (which isin this situation a scalar function of position): constant (upper row) and involving strong (over twointervals) or milder (at a finite number of positions) degenerations. The plots suggest on the one hand hat the availability of vasculature permits an enhanced tumor spread in regions with higher diffusivity,but also enables invasion in areas where the cells modeled by (3.5) infer strongly degenerate diffusionand therewith associated annihilation of pH-taxis. Moreover, the model version (3.5) predicts high cellaccumulation near the interface with highest diffusivity difference, whereas model (3.1) eludes suchtendency potentially leading to singularities. This is presumably due to the flux-limitations in the taxisand additional self-diffusion terms (for more comments on this issue we refer to Section 4). Experiment 4
Effects of flux limitations. We compare the full model (3.1) with its counterpartinvolving in the motility terms only myopic diffusion and pH-taxis without flux saturation, i.e.characterizing glioma density dynamics by M t = ∇∇ : ( D T M ) + α ∇ · ( γ ( D s I N − D T ) M ∇ h ) + M ( − M ) ( − h ) + µ c MW + W , (3.7)the rest of the equations in (3.1) remaining unchanged. The plots in Figure 5 illustrate the differencebetween solution components obtained with these two settings, in the above order. At a relatively earlystage (30 days), with (3.1a) there are more glioma cells at the tumor edge, with an obvious preferencefor moving away from the more acidic, central area. The proton concentration is higher towards thedomain margins, and with time passing (plots at 60 days) it becomes consistently larger all over thedomain in the model version using (3.7), in which diffusion is dominating over taxis. In that setting thetumor cell density gets larger, too, which also applies to EC density; although (3.1) features (slightly)more glioma diffusion, however the flux-limited pH-taxis prevents a faster moving away from acidicregions. Instead, (3.7) ensures a wider spread of glioma, favored by the larger magnitude of thepH-tactic cell fluxes, which lets the tumor cells move faster and therewith produce the attracting VEGFat more locations. After a while (about further 30 days) the previous trend begins to get reversed, as faras glioma density is concerned: the flux limitation keeps the tumor more confined, but the acidity is stilldiffusing fast; the expression of VEGF is enhanced, which subsequently leads to hyperplasia, tumorenhancement, and revival of tiny glioma populations at more distant sites (in right half of thesimulation domain). At the respective ’cancer spots’ the cells stay more confined in the flux-limitedcase, while the model variant with (3.7) allows the cells to migrate faster (green-yellow margins intumor difference plot at 90 days). The tumor cells accordingly produce acid and VEGF, attracting evenmore ECs at their main sites. Much later (300 days) (3.1) predicts larger densities of glioma cells in thecenter of the domain and near the boundary (the former because of the protons being increasinglyremoved by arriving ECs, and the latter due to the mutual proliferation support of the two cell types).The setting with (3.7) leads to pseudopalisades with a larger diameter, being preserved for a longerwhile. Still later (390 days and more) the differences become smaller and smaller in magnitude, whilethe areas with more or less tumor cells become mixed up, with the pseudopalisades disappearing andgiving way to more or less heterogeneous patterns where the glioma cells and ECs are spread over thewhole domain. The simulations (observe in particular the plots at 390 and 420 days) also suggest thatthe motion of ECs is dominated by diffusion (and its limitation near tumor cells) rather than chemotaxistowards VEGF.Figure 6 illustrates 1D tumor patterns obtained with model (3.1) and with (3.7) replacing (3.1a). Twodifferent choices of the glioma diffusion coefficient D T ( x ) (the same as in Figure 7) are considered.Model (3.1) enables glioma diffusion even in the region where D T becomes zero, due to theflux-limited self-diffusion, hence the tumor cells can invade a wider area. In both model versions, as aconsequence of VEGF availability, more ECs gather in the areas with larger glioma density, which inturn favors proton removal and tumor proliferation, followed again by glioma diffusion. For thestrongly degenerating D T the setting without flux limitation leads, as in the case investigated in Figure7, to higher cell accumulations near the interface with highest diffusivity difference, while (as long as M is bounded) the cell flux involving φ ( h , M ) in (3.1) is bounded even if the gradients of h and M are a y s a y s a y s a y s a y s a y s a y s (a) Tumor (b)
Acidity (c)
ECs (d)
VEGF
Figure 2:
Evolution of the solution components at several times in the framework of
Experiment 1 (full model (3.1)). a y s a y s a y s a y s a y s a y s a y s (a) Tumor (b)
Acidity (c)
ECs (d)
VEGF
Figure 3:
Evolution of the solution components at several times in the framework of
Experiment 2 (equally weighted influence of flux-limitedrepellent pH-taxis and self-diffusion). a y s a y s a y s a y s a y s a y s a y s (a) Tumor difference (b)
Acidity difference (c)
Tumor difference (d)
Acidity difference
Figure 4: Experiment 3 : effect of angiogenesis. Left columns: Difference between tumor density and acidity concentration, as obtained withfull system (3.1), and those computed from the model variant (3.5) with proton uptake rate ζ h = . · − /h. Right columns: Same differences,but between computations done with (3.1) and with (3.6). a y s a y s a y s a y s a y s a y s a y s (a) Tumor difference (b)
Acidity difference (c)
EC difference (d)
VEGF difference
Figure 5: Experiment 4 : effect of flux-limited motility terms. Differences between tumor density and proton concentration computed with(3.1) (in the parameter setting of
Experiment 2 ) and those obtained for a model with glioma motility terms only involving myopic diffusionand pH-taxis without flux saturation (i.e. with glioma dynamics as in (3.7)). ot; thus this model version should be less prone to singular structure formation. It also allows forsharper fronts of the cell mass, as can be seen in the right plot of the second row of Figure 6. Figure 6:
1D patterns for two different choices of the glioma diffusion coefficient D T . 2nd row: tumor patterns for model (3.1). 3rd row:tumor patterns when (3.1a) is replaced by (3.7). 4th row: differences between respective patterns in rows 2 and 3. 5th and 6th rows: differencesbetween acidity and EC densities computed with model (3.1) and with (3.7) replacing (3.1a). igure 7: Differences between tumor and acidity patterns (1D) for the model comparison done in the left column of Figure 4 (framework of
Experiment 3 ). Different choices are made for the glioma diffusion coefficient D T ( x ) , as shown in the rightmost column. The vast majority of continuous models for cell migration involve reaction-diffusion(-taxis) PDEs onthe macroscopic scale, irrespective to whether those equations were obtained by a direct formulation ora more or less formal deduction from lower scales. Thereby, the diffusion is in some (most oftenconstant) proportion to the density gradient, which is associated to infinite speed of propagation - anunrealistic feature, accompanied by smearing the sharp profile of the cell mass advancing into itssurroundings; possible overcrowding effects cannot be properly captured either. These and other issuesmotivated the introduction of models with flux saturation also in this context. We refer e.g., to [6, 40]for formal and respectively rigorous derivations of chemotaxis models from KTEs, and to [15, 21, 31]for settings specifically relating to glioma invasion. Among the latter, [21] contains a deduction ofmacroscopic equations with flux-saturated diffusion, haptotaxis, and chemotaxis; the passage from themesoscopic KTAP framework to the RDT system is different from the one performed here, yet stillformal. Very recently, [54] started from a KTE formulation similar to that in [21], hence also related tothe one in the present work, and obtained by yet another upscaling method (with rigorousconvergences) an RDT with myopic diffusion for the zeroth order approximation of the solution. ThatPDE did not involve flux-limitation, but the equation deduced for the first order correction did.The model proposed here is still a bottom-up approach connecting subcellular and mesoscopicdynamics with the macroscopic evolution of cell populations in response to cues from the tumormicroenvironment, leading to patterns which are qualitatively similar to those observed in histologicalimaging of glioblastoma. The setting is an extension of the previous one introduced in [34], since itaccommodates the description of tumor vascularization in response to hypoxia. The main novelty hereis the way of obtaining flux-limited pH-taxis and (supplementary) self-diffusion, namely uponconsidering on the single cell scale the joint effect of cell stress, chemical gradients, and populationpressure to describe forces acting on the cells, cf. (2.1), (2.2). The approach is akin to that employedin [14, 17, 21], however it involves several modifications, as [14, 17] did not contain any fluxlimitations, while [21] does, but models differently the single cell dynamics (allowing, among others,for non-constant speeds) and does not take into account cell reorientations (by way of a turningoperator) in response to local tissue anisotropy and nonlocal orientation distribution of fibers and cells. ur method also differs from [6, 40], where the flux limitation stems from the choice of signal responsefunction involved in the turning operator and depending on the directional derivative of the tacticsignals (chemoattractants). We included in (2.2) only effects of pH-taxis and population pressure, butin modeling glioma migration through a highly anisotropic tissue one could also involve the influenceof its gradients (via macroscopic tissue density, as in [21] or mesoscopic spatial and orientationaldistribution of tissue fibers, as in [17]). There are two ’sources’ of diffusion in our model: the myopicdiffusion comes from the description of cell reorientations via the turning operator, while theself-diffusion term with flux limitation originates from Newton’s second law in (2.1), (2.2). The latterdiffusion part has the effect of eluding diffusion degeneration when the tensor D T nullifies as aconsequence of the mesoscopic tissue distribution q ( x , ˆ v ) vanishing locally. (Strongly) degeneratingmyopic diffusion, alone or in combination with taxis have been related to high cell aggregates andpossible singularity formation even in lower dimensions [27, 50, 51], while flux limitation couldalleviate such behavior, at least under certain conditions, see e.g., [7, 8]. This tendency was alsoobserved in Figure 7, in the case with strongly degenerating D T . Figure 6 showed a sharper tumor cellprofile at the interface with large diffusivity drop, while the model without flux limitation generated amore uniform pattern with a tendency of smearing the interface and faster filling the areas of lower celldensity. We have to stress, nevertheless, that our deduction of macroscopic PDEs is merely formal; asmentioned above, a rigorous one was done in [54], but for a simplified model and involving onlylimitation(s) of signal gradient(s) of the form ∇ S + | ∇ S | instead of our choice encoded in φ ( h , M ) .The obtained macroscopic model is able to reproduce, at least qualitatively, the histological patternsobserved in patient biopsies. Indeed, the typical pseudopalisades due to severe hypoxia are formedeven when (incipient) angiogenesis is present, followed by a disruption of the ’garlands’ asconsequence of acidity removal by vasculature (motile ECs) within the pseudopalisades and therewithassociated repopulation with tumor cells of the formally acidic sites, which was not or very weaklynoticed in simulations of the previous model [34]. The flux limitation does not seem to be necessary forthe latter to happen, however it does influence the duration of pattern formation and maintenance, and(to a lesser extent) their appearance. We also noticed (cf. Figures 2 and 3) that the weights assigned to(flux-limited) pH-taxis and self-diffusion affect the patterns in terms of shape, duration, andpersistence.The macroscopic system (3.1) features two types of taxis: a repellent, flux-limited pH-taxis of gliomacells moving in the opposite direction of the proton concentration gradient ∇ h , and chemotaxis ofendothelial cells following the gradient of VEGF. As such, the setting does not fall into any class ofmodels with multiple taxis as reviewed in [33], since each cell population is performing only one typeof taxis, but it raises mathematical questions which are at least partially related to those models. Theproton production by tumor cells gives the pH-taxis a direct character, while VEGF, the chemotacticcue of ECs, is produced by glioma cells in presence of acidity, while glioma proliferation is favorablyinfluenced by ECs. The latter ensure uptake of protons as well as VEGF, hence exert an indirect andalso a direct influence on their tactic signal. The mathematical analysis of system (3.1) with appropriateinitial and boundary conditions is far from standard; this is not only due to the myopic diffusion and thenonlinearities arising from the source terms and from couplings, but especially to the flux-limitedmotility terms in the glioma reaction-diffusion-taxis PDE. Results about well-posedness and long timebehavior of solutions to such systems are unknown. We merely provided a couple of examples ofsimulation-based, 1D pattern assessments in Figures 7 and 6 - with no pretension of completeness orrigor, but only for the purpose of giving a flavor of the patterns arising from such model. Beyond that, atraveling wave analysis, particularly of glioma and EC dynamics, would be of interest, boththeoretically and from the application viewpoint. cknowledgement PK acknowledges funding by DAAD in form of a PhD scholarship. CS was funded by BMBF in theproject
GlioMaTh
References
Scientific reports , 6:37283, 2016.[4] T. T. Batchelor, D. A. Reardon, J. F. de Groot, W. Wick, and M. Weller. Antiangiogenic therapy forglioblastoma: Current status and future prospects.
Clinical Cancer Research , 20(22):5612–5619,November 2014.[5] N. Bellomo, A. Bellouquid, L. Gibelli, and N. Outada.
A Quest Towards a Mathematical Theoryof Living Systems . Birkh¨auser, 2018.[6] N. Bellomo, A. Bellouquid, J. Nieto, and J. Soler. Multiscale biological tissue models andflux-limited chemotaxis for multicellular growing systems.
Math. Models Methods Appl. Sci. ,20(7):1179–1207, 2010.[7] N. Bellomo and M. Winkler. A degenerate chemotaxis system with flux limitation: maximallyextended solutions and absence of gradient blow-up.
Comm. Partial Differential Equations ,42(3):436–473, 2017.[8] N. Bellomo and M. Winkler. Finite-time blow-up in a degenerate chemotaxis system with fluxlimitation.
Trans. Amer. Math. Soc. Ser. B , 4:31–67, 2017.[9] D. J. Brat, A. Castellano-Sanchez, B. Kaur, and E. G. Van Meir. Genetic and biologic progressionin astrocytomas and their relation to angiogenic dysregulation.
Advances in anatomic pathology ,9(1):24–36, 2002.[10] D. J. Brat, A. A. Castellano-Sanchez, S. B. Hunter, M. Pecot, C. Cohen, E. H. Hammond,S. N. Devi, B. Kaur, and E. G. Van Meir. Pseudopalisades in glioblastoma are hypoxic, expressextracellular matrix proteases, and are formed by an actively migrating cell population.
Cancerresearch , 64(3):920–927, 2004.[11] D. J. Brat and E. G. Van Meir. Vaso-occlusive and prothrombotic mechanisms associated withtumor hypoxia, necrosis, and accelerated growth in glioblastoma.
Laboratory investigation ,84(4):397–405, 2004.[12] Y. Cai, J. Wu, Z. Li, and Q. Long. Mathematical modelling of a brain tumour initiation and earlydevelopment: a coupled model of glioblastoma growth, pre-existing vessel co-option, angiogenesisand blood perfusion.
PloS one , 11(3):e0150296, 2016.[13] A. Caiazzo and I. Ramis-Conde. Multiscale modelling of palisade formation in gliobastomamultiforme.
Journal of theoretical biology , 383:145–156, 2015.
14] A. Chauvi`ere, T. Hillen, and L. Preziosi. Modeling cell movement in anisotropic and heterogeneousnetwork tissues.
Networks & Heterogeneous Media , 2(2):333, 2007.[15] M. Conte, S. Casas-Tint`o, and J. Soler. Modeling invasion patterns in the glioblastoma battlefield.bioRxiv:2020.06.17.156497, 2020.[16] M. Conte and C. Surulescu. Mathematical modeling of glioma invasion: acid-and vasculaturemediated go-or-grow dichotomy and the influence of tissue anisotropy. arXiv preprintarXiv:2007.12204 , 2020.[17] G. Corbin, C. Engwer, A. Klar, J. Nieto, J. Soler, C. Surulescu, and M. Wenske. Modeling gliomainvasion with anisotropy-and hypoxia-triggered motility enhancement: from subcellular dynamicsto macroscopic pdes with multiple taxis. arXiv preprint arXiv:2006.12322 , 2020.[18] G. Corbin, A. Hunt, A. Klar, F. Schneider, and C. Surulescu. Higher-order models for gliomainvasion: From a two-scale description to effective equations for mass density and momentum.
Mathematical Models and Methods in Applied Sciences , 28(09):1771–1800, 2018.[19] A. Czirok. Endothelial cell motility, coordination and pattern formation during vasculogenesis.
Wiley Interdisciplinary Reviews: Systems Biology and Medicine , 5(5):587–602, 2013.[20] W. Diao, X. Tong, C. Yang, F. Zhang, C. Bao, H. Chen, L. Liu, M. Li, F. Ye, Q. Fan, et al. Behaviorsof glioblastoma cells in in vitro microenvironments.
Scientific reports , 9(1):1–9, 2019.[21] A. Dietrich, N. Kolbe, N. Sfakianakis, and C. Surulescu. Multiscale modeling of glioma invasion:from receptor binding to flux-limited macroscopic PDEs. in preparation.[22] S. Eikenberry, T. Sankar, M. Preul, E. Kostelich, C. Thalhauser, and Y. Kuang. Virtualglioblastoma: growth, migration and treatment in a three-dimensional mathematical model.
CellProliferation , 42(4):511–528, 2009.[23] C. Engwer, T. Hillen, M. Knappitsch, and C. Surulescu. Glioma follow white matter tracts: amultiscale dti-based model.
Journal of mathematical biology , 71(3):551–582, 2015.[24] C. Engwer, A. Hunt, and C. Surulescu. Effective equations for anisotropic glioma spread withproliferation: a multiscale approach and comparisons with previous settings.
Mathematicalmedicine and biology: a journal of the IMA , 33(4):435–459, 2016.[25] C. Engwer, M. Knappitsch, and C. Surulescu. A multiscale model for glioma spread includingcell-tissue interactions and proliferation.
Mathematical Biosciences & Engineering , 13(2):443–460, 2016.[26] J. L. Gevertz and S. Torquato. Modeling the effects of vasculature evolution on early brain tumorgrowth.
Journal of Theoretical Biology , 243(4):517–531, 2006.[27] T. Hillen, K. Painter, and M. Winkler. Anisotropic diffusion in oriented environments can leadto singularity formation.
European Journal of Applied Mathematics , 24(3):371–413, December2013.[28] T. Hillen. M 5 mesoscopic and macroscopic models for mesenchymal motion.
Journal ofmathematical biology , 53(4):585–616, 2006.[29] T. Hillen and K. J. Painter. Transport and anisotropic diffusion models for movement in orientedhabitats. In
Dispersal, individual movement and spatial ecology , pages 177–222. Springer, 2013.
30] A. Hunt and C. Surulescu. A multiscale modeling approach to glioma invasion with therapy,.
Vietnam Journal of Mathematics , 45:221–240, 2017.[31] Y. Kim, S. Lawler, M. Nowicki, E. Chiocca, and A. Friedman. A mathematical model forpattern formation of glioma cells outside the tumor spheroid core.
Journal of Theoretical Biology ,260(3):359–371, 2009.[32] P. Kleihues, F. Soylemezoglu, B. Sch¨auble, B. Scheithauer, and P. Burger. Histopathology,classification and grading of gliomas.
Glia , 5:211–221, 1995.[33] N. Kolbe, N. Sfakianakis, C. Stinner, C. Surulescu, and J. Lenz. Modeling multiple taxis: tumorinvasion with phenotypic heterogeneity, haptotaxis, and unilateral interspecies repellence.
Discreteand Continuous Dynamical Systems B , in print, arXiv:2005.01444.[34] P. Kumar, J. Li, and C. Surulescu. Multiscale modeling of glioma pseudopalisades: contributionsfrom the tumor microenvironment. arXiv preprint arXiv:2007.05297 , 2020.[35] G. Martin and R. Jain. Noninvasive measurement of interstitial pH profiles in normal and neoplastictissue using fluorescence ratio imaging microscopy.
Cancer Research , 54(21):5670–5674, 1994.[36] A. Mart´ınez-Gonz´alez, G. Calvo, L. P´erez Romasanta, and V. P´erez-Garc´ıa. Hypoxic cell wavesaround necrotic cores in glioblastoma: a biomathematical model and its therapeutic implications.
Bulletin of Mathematical Biology , 74(12):2875–2896, 2012.[37] M. Onishi, T. Ichikawa, K. Kurozumi, and I. Date. Angiogenesis and invasion in glioma.
Braintumor pathology , 28(1):13–24, 2011.[38] K. Painter and T. Hillen. Mathematical modelling of glioma growth: The use of diffusion tensorimaging (DTI) data to predict the anisotropic pathways of cancer invasion.
Journal of TheoreticalBiology , 323:25–39, 2013.[39] A. Perry and P. Wesseling. Histologic classification of gliomas. In
Handbook of ClinicalNeurology , pages 71–95. Elsevier, 2016.[40] B. Perthame, N. Vauchelet, and Z. Wang. The flux limited Keller-Segel system; properties andderivation from kinetic equations. arXiv:1801.07062, 2018.[41] M. Plank, B. Sleeman, and P. Jones. A mathematical model of tumour angiogenesis, regulatedby vascular endothelial growth factor and the angiopoietins.
Journal of theoretical biology ,229(4):435–454, 2004.[42] S. Prag, E. Lepekhin, K. Kolkova, R. Hartmann-Petersen, A. Kawa, P. Walmod, V. Belman, H.Gallagher, V. Berezin, E. Bock, and N. Pedersen. Ncam regulates cell motility.
Journal of CellScience , 115(2):283–292, 2002.[43] Y. Rong, D. L. Durden, E. G. Van Meir, and D. J. Brat. ‘pseudopalisading’necrosis in glioblastoma:a familiar morphologic feature that links vascular pathology, hypoxia, and angiogenesis.
Journalof Neuropathology & Experimental Neurology , 65(6):529–539, 2006.[44] M. Shamsi, M. Saghafian, M. Dejam, and A. Sanati-Nezhad. Mathematical modeling of thefunction of warburg effect in tumor microenvironment.
Scientific reports , 8(1):1–13, 2018.[45] M. Sidani, D. Wessels, G. Mouneimne, M. Ghosh, S. Goswami, C. Sarmiento, W. Wang, S. Kuhl,M. El-Sibai, J. Backer, et al. Cofilin determines the migration behavior and turning frequency ofmetastatic cancer cells.
The Journal of Cell Biology , 179(4):777–791, 2007.
46] A. Stein, T. Demuth, D. Mobley, M. Berens, and L. Sander. A mathematical model ofglioblastoma tumor spheroid invasion in a three-dimensional in vitro experiment.
BiophysicalJournal , 92(1):356–365, 2007.[47] M. C. Tate and M. K. Aghi. Biology of angiogenesis and invasion in glioma.
Neurotherapeutics ,6(3):447–457, 2009.[48] B. Webb, M. Chimenti, M. Jacobson, and D. Barber. Dysregulated pH: a perfect storm for cancerprogression.
Nature Reviews Cancer , 11(9):671–677, 2011.[49] J. Weickert.
Anisotropic diffusion in image processing , volume 1. Teubner Stuttgart, 1998.[50] M. Winkler. Singular structure formation in a degenerate haptotaxis model involving myopicdiffusion.
Journal de Math´ematiques Pures et Appliqu´ees , 112:118–169, 2018.[51] M. Winkler and C. Surulescu. Global weak solutions to a strongly degenerate haptotaxis model.
Comm. Math. Sci. , 15:1581–1616, 2017.[52] F. Wippold, M. L¨ammle, F. Anatelli, J. Lennerz, and A. Perry. Neuropathology forthe neuroradiologist: palisades and pseudopalisades.
American journal of neuroradiology ,27(10):2037–2041, 2006.[53] X. Zhang, K. Ding, J. Wang, X. Li, and P. Zhao. Chemoresistance caused by the microenvironmentof glioblastoma and the corresponding solutions.
Biomedicine & Pharmacotherapy , 109:39–46,2019.[54] A. Zhigun and C. Surulescu. A novel derivation of rigorous macroscopic limits from a micro-mesodescription of signal-triggered cell migration in fibrous environments. in preparation., 109:39–46,2019.[54] A. Zhigun and C. Surulescu. A novel derivation of rigorous macroscopic limits from a micro-mesodescription of signal-triggered cell migration in fibrous environments. in preparation.