MModeling glioma invasion with anisotropy- and hypoxia-triggeredmotility enhancement: from subcellular dynamics to macroscopicPDEs with multiple taxis
G. Corbin, C. Engwer, A. Klar, J. Nieto, J. Soler, C. Surulescu, and M. WenskeJune 23, 2020
Abstract
We deduce a model for glioma invasion making use of DTI data and accounting for the dy-namics of brain tissue being actively degraded by tumor cells via excessive acidity production,but also according to the local orientation of tissue fibers. Our approach has a multiscale charac-ter: we start with a microscopic description of single cell dynamics including biochemical and/orbiophysical effects of the tumor microenvironment, translated on the one hand into cell stressand corresponding forces and on the other hand into receptor binding dynamics; these lead onthe mesoscopic level to kinetic equations involving transport terms w.r.t. all kinetic variablesand eventually, by appropriate upscaling, to a macroscopic reaction-diffusion equation for gliomadensity with multiple taxis, coupled to (integro-)differential equations characterizing the evolu-tion of acidity and macro- and mesoscopic tissue. Our approach also allows for a switch betweenfast and slower moving regimes, according to the local tissue anisotropy. We perform numericalsimulations to investigate the behavior of solutions w.r.t. various scenarios of tissue dynamics andthe dominance of each of the tactic terms, also suggesting how the model can be used to performa numerical necrosis-based tumor grading or support radiotherapy planning by dose painting. Wealso provide a discussion about alternative ways of including cell level environmental influencesin such multiscale modeling approach, ultimately leading in the macroscopic limit to (multiple)taxis.
Glioma is a common type of primary, fast growing brain tumor with a poor prognosis, themedian survival time with the most frequent (and most aggressive) form called glioblastomemultiforme amounting at 60 weeks, in spite of modern treatment involving resection, radio-, andchemotherapy (see e.g., [91] and references therein). An exhaustive removal of the tumor is ingeneral impossible, due to the rapid advancement of glioma cells through their cycle and to thediffuse tumor infiltration. This leads to serious clinical challenges and to rather modest treatmentoutcomes. The correct assessment of tumor margins and of the related GTV, CTV, and PTV is therefore of utmost importance. Noninvasive medical imaging techniques like MRI and CTare valuable diagnostic tools, however they only provide a macroscopic classification of neoplasticregions and the surounding oedema, without being able to evaluate the actual tumor extent whichis heavily influenced by the biological, mechanical, and chemical processes on single cell andsubcellular levels. Moreover, every patient has a different brain structure and there is abundantevidence that glioma follow the white matter tracts made up of myelinated axon bundles [32, 33].Therefore, the personalized prediction of the tumor burden is necessary for an enhanced therapyplanning and mathematical models are called upon to provide such information via numericalsimulations, thereby relying on medical data, of which diffusion tensor imaging (DTI) is quitecommon [80]. It measures the spatial diffusion of water molecules by MRI per volume element(voxel). The information obtained in this way can be expressed in the form of a diffusion tensorcontaining the full (apparent) diffusion information along six directions. One way of visualizing GTV=gross tumor volume; CTV= clinical target volume; PTV=planning target volume a r X i v : . [ q - b i o . T O ] J un iffusion tensors is by the so-called fractional anisotropy index FA which is a measure for thelocal tissue alignment and can directly be obtained from DTI data; it is a scalar value betweenzero and one, calculated by using the eigenvalues of the apparent diffusion tensor. A value closeto one means high anisotropy, i.e., a strong preference for a specific direction, whereas a verysmall value corresponds to the nearly isotropic case. For more details we refer e.g., to [10], seealso [25]. For more about DTI and its visualization approaches along with advantages, drawbacksand extensions we refer e.g., to [2, 24, 44].Several approaches to modeling glioma invasion have been considered so far; they range fromdiscrete formulations [36,45,46] over hybrid settings [31,82,93] combining decriptions of individualcells moving on a lattice with PDEs describing the evolution of some stimulus (e.g., a chemoat-tractant) and up to pure continuum models coupling several types of differential equations. Thelatter category comprises in turn several classes of models, according to the scales taken intoaccount and to the type of PDEs employed. Pure macroscopic models using reaction-diffusionequations (possibly with space- and time-dependent diffusion coefficients, or accounting for theinteraction with the surrounding tissue by letting the diffusion coefficient be proportional to thewater diffusion tensor assessed by DTI) have been introduced e.g., in [13, 43, 53, 85] to charac-terize the evolution of glioma density; in [21] some biomechanical effects have been included aswell. While these models impose the precise form of the equations and their coefficients, otherapproaches [39, 72] deduce them by using a scaling argument from a mesoscopic setting relyingon kinetic transport equations for the glioma density function depending -supplementary to timeand position- also on the cell velocities. These models have been further extended in [25–27,41] totake into account subcellular level effects on the mesoscopic and ultimately macroscopic behaviorof the tumor. Concretely, these multiscale settings included on the microscopic level receptorbinding to tissue fibers, which was related to the mesoscale by way of an additional transportterm and by the turning rate. Appropriate scalings then led to the macroscale dynamics, whosededuced coefficients are still carrying subcellular level information. In [68] different types of scal-ing have been considered for this modeling approach, and a well-posedness result was providedfor a model also involving the time-space evolution of tissue. A more general analytical result inless regular function spaces for a complex micro-meso-macroscale model involving both chemo-and haptotaxis introduced on the mesolevel has been proved in [61].Here we propose a multiscale model starting from a kinetic transport equation formulation andaccounting, too, for receptor binding dynamics. Unlike in [25–27, 41] where the space-dependenttissue was given by DTI data at a fixed time point and the glioma cells were migrating on thisbrain structure, in the current model the tissue is allowed to evolve in time as well and the gliomacells degrade it according to their respective orientation to the tissue fibers, similarly to the wayintroduced in [38] and readdressed in [61,68]. Moreover, in this work we also include a character-ization of velocity changes via a transport term involving derivatives with respect to the velocityvariable, hence allowing for external forces to act on the cells. Specifically, these external actionsare decisively influenced by the spatial gradients of tissue fibers in the tumor surroundings andby the bias induced through the acidity gradient. Reference [19] considered a similar approachaccounting for a given chemotactic force acting on the cells and being proportional to the concen-tration gradient of some non-evolving chemoattractant. Here we pay more detailed attention tothe influence of time- and space-dependent acidity and tissue on the cell migration, proliferation,and depletion, thereby also looking into the effect of tissue fiber orientation, hence dealing withboth macroscopic and mesoscopic tissue evolution. By formal upscaling we then obtain a systemof ordinary and partial (integro-)differential equations with (myopic) diffusion, multiple taxis,and nonlinear source terms, which couples macroscopic dynamics of tumor cells, necrotic matter,acidity, and normal tissue with that of mesoscopic tissue. The rest of the paper is organized asfollows: in Section 2 we set up the micro- and mesoscopic descriptions of (sub)cellular gliomadynamics in the kinetic theory of active particles (KTAP) framework developed e.g., in [6, 7]and also provide descriptions of the other involved model variables: tissue, necrotic matter, andacidity. Section 3 is dedicated to deducing effective equations for the macroscopic glioma densityinfluenced by tissue and acidity dynamics. Formal parabolic and hyporbolic upscalings are per-formed and a low-order moment approximation is deduced for the characterization of mesoscopictissue, in order to reduce the complexity of the system. Numerical simulations are conductedin Section 4 for three different scenarios of increasing complexity, and the obtained results arecompared and commented. Finally, Section 5 provides a discussion of this work’s outcome, modelextensions, and perspectives, along with a review of various ways of multiscale modeling of tacticcell behavior in the KTAP framework. Set up of the micro-meso model
Our modeling approach aims at obtaining a macroscopic description of the tumor cell evolutionin interaction with the underlying tissue. The resulting PDE should carry in its coefficients infor-mation about relevant processes occuring on the lower levels (subcellular/microscopic, individualcells/mesoscopic) and should be deduced from mathematical descriptions of the dynamics on thelatter. On the macroscopic scale the two main components of cell motion in a tissue should beretrieved, namely diffusion and haptotaxis.We start by introducing some notations for the various variables and functions involved inthis work.
Notation 1 • Mechanical variables: time t ∈ R + , space x ∈ R N , velocity v ∈ V := s S N − , where S N − is the unit sphere in R N , and, as in [25–27, 38, 42], we assume the cell speed s > • ˆ v := v | v | a unitary vector denoting the direction of a vector v ∈ V . • Q ( t, x ): macroscopic tissue density. • Biological (activity) variable y : volume fraction of bound receptors on the cell membrane.Its dynamics is given by an ODE obtained for y by mass action kinetics of receptor bindingsand detachments (with rates k + and k − , respectively):˙ y = k + ( R − y ) Q ( t, x ) − k − y, y ∈ (0 , R ) . (2.1)This represents in our model the subcellular dynamics. Thereby, as in previous works[49, 68], we assume for simplicity that the total amount of cell membrane receptors R isconstant and we rescale y → y/R , thus obtaining ˙ y = G ( y, Q ) := k + (1 − y ) Q ( t, x ) − k − y, y ∈ Y := (0 ,
1) (2.2) • p ( t, x, v, y ): density function of tumor cells (mesoscopic tumor cell density). We will alsouse ¯ p ( t, x, y ) := (cid:82) V p ( t, x, v, y ) dv . • ρ ( t, x ) := (cid:82) Y ¯ p ( t, x, y ) dy = (cid:82) Y (cid:82) V p ( t, x, v, y ) dvdy : macroscopic density of tumor cells. • q ( t, x, θ ): directional distribution function of tissue fibers with orientation θ ∈ S N − , where S N − denotes the unit sphere in R N . Hence (cid:82) S N − q ( t, x, θ ) dθ = 1. • E q ( t, x ) := (cid:82) S N − θq ( t, x, θ ) dθ : mean fiber orientation. • V q ( t, x ) := (cid:82) S N − ( θ − E q ( t, x )) ⊗ ( θ − E q ( t, x )) q ( t, x, θ ) dθ : variance-covariance matrix fororientation distribution of tissue fibers. • h ( t, x ): concentration of acidity (protons), a genuinely macroscopic quantity. • n ( t, x ): density of necrotic tissue (also a macroscopic quantity). • λ ( x, y ): cell turning rate.Next we characterize innovations of cell density function due to changes in velocity and ac-tivity variable under the influence of the tissue structure. Depending on their nature, these areincorporated into the model in different ways:(i) As transport parts of the mesoscopic equation for the cell density function p : a div y termaccounts for the binding of cell receptors to tissue fibers, while a div v term describes biome-chanical effects, such as cell stress and directionality dictated by the spatial variations in(mesoscopic) tissue density and macroscopic acidity concentration.(ii) As a jump process, being modeled as an integral operator L on the right-hand side of theequation. The kernel of that operator will depend as in previous works [25–27,38,42] on thetissue structure, by way of the fiber directional distribution q ( t, x, ˆ v ), as will be explainedbelow. Observe that the activity variable y belongs to an open set; indeed, y = 0 would mean that no receptors are bound(which is not possible for cells living in tissue), while y = 1 would correspond to all receptors being occupied (whichactually never happens, as it would impede migration, proliferation, and even cell survival, due to anoikis). For anactual mathematical justification we refer to [49, 61]. iii) As a source term describing proliferation/decay and being modeled as another integral oper-ator P . This can be done for instance similarly to the approach in [26], where the interactionof cells with the surrounding tissue was considered to be at the onset of cell survival andmitosis. Here we modify that approach in order to account for the unfavorable effect ofacidity on these processes. The concrete form of the corresponding term P ( p, Q, h, ρ ) willbe specified later.Thus, the equation for p takes the following form: ∂ t p + v ·∇ x p + ∂ y ( G ( y, Q ) p ) + α div v (cid:16) B ( v )( ∇ x q ( t, x, ˆ v ) − ∇ x h ( t, x )) p (cid:17) = λ ( x, y ) L ( p ) + P ( p, Q, h, ρ ) , (2.3)where we used the notation B ( v ) := | v | I N − v ⊗ v , and where the transport coefficients areobtained from the subcellular (microlevel) dynamics : dxdt = v ; dvdt = α B ( v )( ∇ x q ( t, x, ˆ v ) − ∇ x h ( t, x )); dydt = G ( y, Q ) , (2.4)the latter equation in (2.4) being just (2.1) and the middle equation characterizing the dynamicsof v influenced by the space gradients of the directional distribution of tissue fibers and of theacidity. Thereby, the tensor v ⊗ v represents the active cell stress, while −| v | I N is the isotropicpart. Indeed, concerning cell velocities we are foremost paying attention to the direction of thevector v and consider that it is mainly influenced by the fiber distribution (and in particular byits spatial variations), as well as by the acidity gradient. The latter is known to drive the tumorcells and seems to have a bidirectional effect: Low extracellular pH favorizes both their migrationand proliferation [29,79,86], but acidosis can also inhibit cell proliferation, induce stress response,and apoptosis, see e.g. [69]. Specifically in glioblastoma multiforme (GBM) it can trigger motionin the opposite direction of the acidity gradient. Indeed, GBM cells form typical, garland-likestructures called pseudopalisades, which are centered around the highly hypoxic occlusion site ofa capillary [14, 89]. The presence of such histological patterns is an indication for poor prognosisof patient survival [51]. In order not to complicate the exposition too much we will only considerhere the latter situation, with the repellent effect of acidity, hence the motion of glioma cellsbeing biased towards −∇ x h .The dimensionless parameter α will play an important role, since it incorporates informationfrom the different scales on which the tumor density, fiber directional distribution or acidity arecorrelated. While the tumor has volume dimensions, the path of travel made up of targetingfibers has a thinner structure close to the flat one. We can assume that the relationship withacidity is similar in scale. Therefore, the gradients of q and h w.r.t. x contribute with this newscale parameter α among the agents of this process, which must be, as a consequence, a largequantity.Observe that the dynamics of v in (2.4) ensures that d | v | dt , hence the cell speed | v | = s isconstant (in line with the above assumption V = s S N − ), hence only the direction matters.For the evolution of acidity concentration h we also need to provide an equation; it will takethe form of a reaction-diffusion PDE describing the facts that protons diffuse through the entireavailable space with diffusion parameter D H , acidity is produced by cancer cells, and may inferdecay, e.g. through uptake by vasculature or normal cells. We denote by b > a >
0. In (2.5) the constant h denotes an acidity thresholdwhich is critical for non-cancerous matter; the proton buffering is enhanced when that thresholdis exceeded, and it is limited by the availability of normal tissue and cells, which decreases underhypoxia (as explained later on): h t = D H ∆ h + aρ ρ − bQ ( h − h ) . (2.5)The turning operator L on the right hand side in (2.3) is taken as a relaxation-type operatordescribing the velocity-jump discontinuities in the cell density function p due to contact guidance, Strictly speaking we understand the ODEs for v and y as being written for tissue q , Q and acidity h , scaled withtheir reference quantities (e.g., tissue carrying capacity K Q and acidity threshold h introduced below), but for thesake of simplifying the notation we do not explicitly write those denominators. Same applies to the complementarysolution components involved in the source terms of the forthcoming equations for cell density, acidity, and necroticmatter. hat is, these velocity alterations are assumed to be caused by the cells interacting with tissuefibers and adapting their motion to biochemical (biomechanical, etc.) cues available at sites intheir direct proximity. More precisely, it is given by L ( p ) := ¯ p ( t, x, y ) M ( t, x, v ) − p ( t, x, v, y ) , (2.6)where M ( t, x, v ) depends on the fiber directional distribution q : M ( t, x, v ) := q ( t, x, ˆ v ) ω , ω := (cid:90) V q ( t, x, ˆ v ) d v, thus (cid:82) V M ( t, x, v ) d v = 1. Notice that this corresponds to considering, as in [25–27, 38, 72], aturning operator of type L ( p ) := − p ( t, x, v, y ) + (cid:90) V ˜ K ( v, v (cid:48) ) p ( t, x, v (cid:48) , y ) dv (cid:48) with the turning kernel defined by ˜ K ( v, v (cid:48) ) = q ( t,x, ˆ v ) ω , that is, the random migration from anyvelocity v (cid:48) to a new one v is given by the (normalized) orientational distribution of tissue fibers onthis last velocity v , so through contact guidance. We recall that this turning operator is multipliedby the cell turning rate λ ( x, y ). It will play an essential role in our attempt to characterize thealternation between slow and enhanced migration of tumor cells, which is mainly due to the tissueanisotropy. To capture this dependence we use the fractional anisotropy index FA mentioned inSection 1. Henceforth we deal with a turning rate of the form λ ( x, y ) = κyF A ( x ) + y , which is in agreement with the assumption of reduced turning in highly aligned regions and withits dependence on the receptor binding state y , the latter with a certain saturation modeled bythe Monod type factor. Thereby, κ is a constant that refers to the maximum turning frequency.Finally, we introduce the proliferation/decay operator. Similarly to the approach in [26], bothproliferation and decay of tumor cells originate (at least partially) in the interaction of cells withthe surrounding tissue, but here we modify that approach in order to account for the unfavorableeffects of acidity on these processes. Concretely, P ( p, Q, h, ρ ) is given by P ( p, Q, h, ρ ) := µ ( ρ, h ) (cid:90) Z χ ( x, z, z (cid:48) ) p ( t, x, v, z (cid:48) ) Q ( t, x ) dz (cid:48) , (2.7)where χ ( x, z, z (cid:48) ) is a kernel with respect to z and represents the transition from the state z (cid:48) tothe state z during a proliferation-favorable glioma-tissue interaction.For the tissue dynamics we adopt the approach in [38] which was also employed in [49,61,68],with slight modifications. We look for the evolution of q ( t, x, θ ), hence want an equation of theform ∂ t q ( t, x, θ ) = τ ( θ, p, q ) , with τ satisfying the natural conditions: τ ( · , · ,
0) = 0 , (cid:90) S N − τ ( θ, · , · ) dθ = 0 , ∀ θ ∈ S N − . We thus compare -as in [38]- the fiber orientation θ with the cell direction ˆ v , and introducethe operator:Π[ p ]( t, x, θ ) := (cid:40) ρ ( t,x ) (cid:82) Y (cid:82) V | θ · ˆ v | p ( t, x, v, y ) d v d y, for undirected tissue ρ ( t,x ) (cid:82) Y (cid:82) V θ · ˆ v p ( t, x, v, y ) d v d y, for directed tissue, (2.8)which averages the projection of cell movement direction on the fiber orientation. It modelsthe fact that cells preferentially degrade fibers which are nearly orthogonal to their movementdirection (for more details see [38]). It holds that 0 ≤ Π[ p ] ≤ − ≤ Π[ p ] ≤ undirected means that the fibers making up thetissue are symmetrical all along their axes, i.e. there is no ’up’ and ’down’ on such fibers, whichtranslates into symmetry of the orientational distribution: q ( t, x, − θ ) = q ( t, x, θ ) , ∀ θ ∈ S N − . f this is not the case, then the tissue is said to be directed . These tissue properties reflecton the cell motility: While undirected fibers do not impose any supplementary bias on the(re)orientation of a cell, in a directed tissue the cells will have a preferred direction of advancementalong the fibers. As DTI data do not provide any information about such directionality, we willconsider in the following both types of tissue, which will lead to two different ways of performingthe transition from the lower scales (micro and meso) to the macroscopic description of tumordynamics in interaction with their fibrous and biochemical surroundings.If we denote by g ( t, x, θ ) the mesoscopic density of tissue fibers with orientation θ we obtainfor the evolution of this quantity the equation ∂ t g ( t, x, θ ) = r D ( h )(Π[ p ]( t, x, θ ) − ρ ( t, x ) g ( t, x, θ ) , (2.9)where r D is a nonnegative quantity characterizing the efficiency of fiber degradation (e.g., itcan be considered to be proportional to the amount of matrix degrading enzymes expressed bythe cells, thus to the acidity h , since MDEs are actually known to be enhanced by an acidicextracellular pH, see [16, 47]). We then consider r D to be a monotonically increasing function of h with r D ( h ) = r ( h − h ) + .It would be important to have information about the directional distribution function q ( t, x, θ ),hence we express the latter with respect to the mesoscopic tissue density g ( t, x, θ ). Following [38],this relationship takes the form q ( t, x, θ ) = g ( t, x, θ ) (cid:82) S N − g ( t, x, θ ) dθ , g (cid:54) = 0 . The macroscopic tissue density Q ( t, x ) represents the volume fraction of tissue fibers, irre-spective of their orientation, thus we also have the following relationship: g ( t, x, θ ) = q ( t, x, θ ) Q ( t, x ) , t > , x ∈ R N , θ ∈ S N − . (2.10)Integrating (2.9) w.r.t. θ ∈ S N − we get ∂ t Q = r D ( h ) ρQ (cid:18)(cid:90) S N − Π[ p ]( t, x, θ ) q ( t, x, θ ) dθ − (cid:19) = r D ( h ) ρ Q (cid:90) S N − (cid:0) Π[ p ] − (cid:1) q dθ. (2.11)In the sequel the average of Π[ p ] w.r.t. the distribution of fiber orientations will be denotedas in [38] by A [ p ]( t, x ) := (cid:90) S N − Π[ p ]( t, x, θ ) q ( t, x, θ ) dθ, (2.12)hence (2.11) becomes ∂ t Q = r D ( h ) ρQ ( A [ p ] − . (2.13)Likewise, we can deduce an equation for the fiber directional distribution function: ∂ t q ( t, x, θ ) = r D ( h ) ρ ( t, x ) q ( t, x, θ ) (Π[ p ]( t, x, θ ) − A [ p ]( t, x )) . (2.14)Eventually, the decay of both tissue and glioma (primarily due to acidity) generates necrotictissue, large amounts of which are an indicative of poor survival prognosis [35, 62]. Therefore,the detection and assessment of necrosis is an important issue in therapy. Here we describe thedynamics of necrotic tissue density by ∂ t n = r D ( h ) ρQ (cid:90) S N − (1 − Π[ p ]( t, x, θ )) q ( t, x, θ ) dθ + F ( h, ρ, Q ) (2.15)where F ( h, ρ, Q ) describes glioma death due to hypoxia and is to be specified later.Thus, the full model for cell-tissue interactions and response to acidity is given by (2.3), (2.5),(2.13), (2.14), (2.15), supplemented with appropriate initial and boundary conditions.As in previous works [25–27, 49] we assume p to be compactly supported in the ( v, y )-space.The initial conditions and the boundary conditions w.r.t. space will be addressed in the followingsections. The wellposedness of this problem (especially in less regular function spaces) is nottrivial and also not the aim of this paper. When the space domain is the whole R N the approachin [48, 49, 61] could be used as a starting point in order to address this issue. shortly MDEs Towards population behavior: Upscaling
Our objective is to assess the macroscopic evolution of the tumor cell population interactingwith the brain tissue. As the latter exhibits anisotropy variability, with highly aligned regionsalternating with areas of isotropic fiber distribution, the invasive behavior of glioma cells will becorrespondingly -and locally- dominated by diffusion with or without drift. Therefore, it is desir-able that our model includes a switch between these two kinds of motion, in the sense that theinfluence of drift together with diffusion can be potentiated against pure diffusion. These effectsare built in via the turning rate λ ( x, y ) given in Section 2 above. Indeed, its dependence on the(local) fractional anisotropy F A ( x ) and the amount of receptors bound to their ligands on thetissue fibers will be able to capture the alternation between the epochs of higher- or less-alignedtissue.1. Mass conservation: (cid:82) V L ( p )( t, x, v, y ) d v = 0.2. Self–adjointness: L is self–adjoint in L ( dvM ( v ) ).3. Kernel of L : L ( p ) = 0 ⇔ p ∈ (cid:104) M ( v ) (cid:105) , thus Ker ( L ) = (cid:104) M ( v ) (cid:105) , the space gener-ated by M ( v ). In particular, using the self–adjointness, we know that M ⊥ ∈ (cid:104) M ( v ) (cid:105) ⊥ iff (cid:82) V M ⊥ ( v ) dv = 0.Before proceeding with the scaling we re-express (similarly to e.g. [25]) the subcellular dy-namics in a more convenient way: Let us consider y ∗ = f ( Q ) := k + Qk − + k + Q , the steady-state of(2.1), and denote by z := y − y ∗ the deviation of the current activity variable y from the steady state. Then z ∈ Z ⊆ [ − y ∗ , − y ∗ ]and the microscopic dynamics (2.4) turns into dxdt = v ; dvdt = α B ( v ) (cid:0) ∇ x q ( t, x, ˆ v ) − ∇ x h ( t, x ) (cid:1) ; dzdt = − ( k + Q + k − ) z − f (cid:48) ( Q )( Q t + v ·∇ x Q ) . Then (2.3) becomes ∂ t p + v ·∇ x p − ∂ z (cid:16) (( k + Q + k − ) z + f (cid:48) ( Q )( Q t + v · ∇ Q )) p (cid:17) + α div v (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) p (cid:17) = λ ( x, z ) L ( p ) + P ( p, Q, h, ρ ) , (3.16)where the turning rate in terms of z is given by λ ( x, z ) = κ ( z + f ( Q ( t,x ))) F A + z + f ( Q ( t,x )) . To simplify thesubsequent computations we linearize it as follows: λ ( t, x, z ) (cid:39) λ ( t, x ) + λ ( t, x ) z, where λ ( t, x ) = κf ( Q ) F A + f ( Q ) and λ ( t, x ) = κF A ( F A + f ( Q )) . Remark 3.1
Notice that the fractional anisotropy FA changes dynamically, i.e. depend bothon t and x , since tissue evolution (here this means degradation) leads to local modifications ofthe water diffusion tensor and, correspondingly, of its eigenvalues. Indeed, the apparent diffusiontensor can even vanish locally, in which situation the cell diffusivity degenerates. In hithertoworks [22,23,25–27,41–43,72] FA has been assumed to be time-independent, motivating that theresolution of DTI data does rarely go below voxels with sizes of 1 mm , which is rather roughcompared with the size of glioma cell bodies (15-60 µm , [56]). This simplifies the numericalhandling and also avoids problems related to the possibly singular behavior of solutions to themacroscopic PDE, as systems with degenerating myopic diffusion and haptotaxis can lead toblow-up even in 1D [87, 88]. In Section 4, however, we will consider several simulation scenarios,including, in turn, evolving and time-stationary tissue, and compare their outcome. The parabolic scaling corresponds formally to the change of variables t → ε t, x → εx. Weperform it and thereby rescale as in [26] the source term P ( p, Q, h, ρ ) with ε , in order to let itact on the correct time scale. Correspondingly, equation (3.16) becomes: ε ∂ t p + ε ∇ x · ( vp ) − ∂ z (cid:16) (( k + Q + k − ) z + f (cid:48) ( Q )( ε Q t + εv · ∇ Q )) p (cid:17) + αε div v (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) p (cid:17) = λ ( t, x, z ) L ( p ) + ε µ ( ρ, h ) Q ( t, x ) (cid:90) Z χ ( x, z, z (cid:48) ) p ( t, x, v, z (cid:48) ) dz (cid:48) . (3.17) ext we consider the moments of p with respect to v and especially with z : m ( t, x, v ) := (cid:90) Z p ( t, x, v, z ) dz, m z ( t, x, v ) := (cid:90) Z zp ( t, x, v, z ) dz, ρ z ( t, x ) = (cid:90) V m z ( t, x, v ) dv ¯ p ( t, x, z ) := (cid:90) V p ( t, x, v, z ) dv, ρ ( t, x ) := (cid:90) V m ( t, x, v ) dv = (cid:90) Z ¯ p ( t, x, z ) dz and neglect the higher order moments with respect to z by assuming very small deviations ofthe receptor binding dynamics from the steady-state, i.e. by assuming z to be very small. Asthe subcellular dynamics is very fast in comparison to cell motion and proliferation, this is areasonable assumption.Then from (3.17) we obtain the moment equations: ε ∂ t m + ε ∇ x · ( vm ) + αε div v (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) m (cid:17) = λ ( t, x )( M ( v ) ρ − m )+ λ ( t, x )( M ( v ) ρ z − m z ) + ε µ ( ρ, h ) Q ( t, x ) (cid:90) Z (cid:90) Z χ ( x, z, z (cid:48) ) p ( t, x, v, z (cid:48) ) dz (cid:48) dz (3.18a) ε ∂ t m z + ε ∇ x · ( vm z ) + αε div v (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) m z (cid:17) + ( k + Q + k − ) m z + εf (cid:48) ( Q ) v · ∇ Qm + ε f (cid:48) ( Q ) mQ t = λ ( t, x )( M ( v ) ρ z − m z ) + ε µ ( ρ, h ) Q ( t, x ) (cid:90) Z (cid:90) Z zχ ( x, z, z (cid:48) ) p ( t, x, v, z (cid:48) ) dz (cid:48) dz (3.18b)Performing the usual Hilbert expansion p = (cid:80) k ε k p k and consequently the expansion of themoments: m = (cid:88) k ε k m k , m z = (cid:88) k ε k m zk , ρ = (cid:88) k ε k ρ k , ρ z = (cid:88) k ε k ρ zk (3.19)and identifying the powers of ε we obtain from (3.18a) and (3.18b) the following relationships: (cid:15) terms:0 = λ ( M ( v ) ρ − m ) + λ ( M ( v ) ρ z − m z )( k + Q + k − ) m z = λ ( M ( v ) ρ z − m z ) ⇒ ρ z = m z = 0 and M ( v ) ρ = m .(cid:15) terms: ∇ x · ( vm ) + α div v (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) m (cid:17) = λ ( M ( v ) ρ − m ) + λ ( M ( v ) ρ z − m z ) ∇ x · ( vm z ) + α div v (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) m z (cid:17) + ( k + Q + k − ) m z + f (cid:48) ( Q ) v · ∇ Qm = λ ( M ( v ) ρ z − m z ) . Integrating with respect to v we obtain (for undirected tissue) ρ z = 0 , m z = − M ( v ) ρ k + Q + k − + λ f (cid:48) ( Q ) v · ∇ Q, (3.20a) m = 1 λ (cid:104) − ∇ x · ( vM ( v ) ρ ) − α div v (cid:16) B ( v )( ∇ x q ( t, x, ˆ v ) − ∇ x h ) M ( v ) ρ (cid:17) − λ m z (cid:105) (3.20b)+ M ( v ) ρ . (3.20c) (cid:15) terms from (3.18a), after expanding µ about ρ and integrating the equation with respect to v : ∂ t ρ + ∇ x · (cid:90) V vm dv = µ ( ρ , h ) Qρ . From (3.20) we can compute (cid:90) V vm dv = 1 λ (cid:104) − ∇ x · ( (cid:90) V v ⊗ vM ( v ) dvρ ) − α (cid:90) V v div v (cid:16) B ( v )( ∇ x q ( t, x, ˆ v ) − ∇ x h ( t, x )) M ( v ) ρ (cid:17) dv + λ f (cid:48) ( Q ) k + Q + k − + λ (cid:90) V v ⊗ vM ( v ) dv ∇ Qρ (cid:105) . e denote by ˜ E q ( t, x ) := (cid:90) V v q ( t, x, ˆ v ) ω d v = s (cid:90) S N − θq ( t, x, θ ) dθ = s E q ( t, x ) (3.21)and with D T ( t, x ) := (cid:90) V ( v − ˜ E q ) ⊗ ( v − ˜ E q ) M ( t, x, v ) d v = (cid:90) V ( v − ˜ E q ) ⊗ ( v − ˜ E q ) q ( t, x, ˆ v ) ω d v (3.22)the so-called tumor diffusion tensor . In the situation with undirected tissue we have E q ( t, x ) = 0 , thus in this case D T ( t, x ) = s V q ( t, x ), where V q denotes the variance-covariance matrix w.r.t.the fiber orientation distribution q . With the notation g ( Q, λ )( t, x ) := f (cid:48) ( Q ) k + Q + k − + λ , we have (cid:90) V vm dv = − λ ∇ x · ( D T ρ ) + λ λ g ( Q, λ ) D T ∇ Qρ + αλ Σ( t, x ) ρ , (3.23)with Σ( t, x ) := S ( t, x ; q ) − S ( t, x ; h ), where S ( t, x ; q ) := (cid:90) V B ( v ) ∇ x q ( t, x, ˆ v ) M ( t, x, v ) dv = 1 ω (cid:90) V B ( v ) ∇ x q ( t, x, ˆ v ) q ( t, x, ˆ v ) dv (3.24a) S ( t, x ; h ) := (cid:90) V B ( v ) M ( t, x, v ) dv ∇ x h ( t, x ) = s ( I N − V q ( t, x )) ∇ x h ( t, x ) . (3.24b)Then, the limiting macroscopic equation for the tumor cell density takes the form ∂ t ρ = ∇ · (cid:104) λ ( x ) (cid:16) ∇ · ( D T ρ ) − λ ( x ) g ( Q, λ ) D T ∇ Qρ − αS ( t, x ; q ) ρ + αS ( t, x ; h ) ρ (cid:17)(cid:105) + µ ( ρ , h ) Qρ . (3.25)The obtained equation is of drift-diffusion type; the first term on the right hand side representsa non-Fickian, so-called myopic diffusion: the cells spread out according to information availablein their immediate surroundings. The next two terms feature a sign opposite to the myopicdiffusion and represent drift corrections of the diffusive part in the direction of the macroscopicand mesoscopic tissue gradients ∇ Q and ∇ q , respectively. They could therefore be interpretedas haptotaxis-like terms. The first of them carries through the function g the influence of thesubcellular receptor binding dynamics to the surrounding tissue, while the remaining drift termaccounts for the (mesoscopic) stress exerted by the cells on the tissue, the latter being describedby the directional distribution of fibers q . The drift term involving S in (3.25) has the samesign with the diffusion, and contains a bias towards the opposite direction of acidity gradient ∇ h , thus describes a repellent pH-taxis. The information about the brain tissue structure, whichis decisive for personalized predictions of the tumor space-time evolution, is contained in (3.25)via the tumor ’diffusion’ tensor D T and the vector S ( t, x ; q ), which in view of (2.4) can be seenas a biomechanical interaction force between tumor cells and tissue. Thereby, the stress tensor B ( v ) contributes (together with the kernel M ( t, x, v ) = q ( t,x, ˆ v ) ω ) to the (mesoscopic) haptotacticsensitivity of the cells, which is expressed here in a velocity-averaged way. The last term in(3.25) encodes proliferation and decay of glioma embedded in tissue. The function µ ( ρ , h ) canbe correspondingly chosen, e.g. in the form µ ( ρ , h ) = η (1 − ρ − n )( h T − h ) + − γh, (3.26)which would correspond to a logistic type of growth and a hypoxia-triggered necrosis. Thereby η, γ > h T denotes an acidity threshold which is critical for the cancer cells: This designation is somewhat abusive, since here -unlike previous works [25–27,72]- the turning rate is not involvedin its expression, due to its dependence on the position. When writing out the myopic diffusion it becomes evidentthat the ’actual’ diffusion tensor is D T /λ . hen the proton concentration exceeds it, they become hypoxic and cease proliferation. Thechoice of µ establishes the form of the last source term in (2.15). Thus, with (3.26) we have inthe equation for necrotic tissue F ( h, ρ, Q ) = γhρQ. (3.27)The role of the turning rate λ is twofold: it connects the cell reorientations to the recep-tor binding kinetics and it also captures the effect of tissue anisotropy on the cells migratingthrough the tissue. While λ influences all motility terms, λ is specific for the haptotaxis com-ponent including subcellular level effects. The factor λ is independent on y and increasing withthe space-varying fractional anisotropy F A , the effect of which is particularly accentuated inthe term multiplied with λ . This means that λ with the inherent y -dependence provides ananisotropy-triggered switch between myopic diffusion with ’mesolevel’ haptotaxis and migratorymode with enhanced haptotaxis, the latter supplementary providing bias in the direction of thetissue gradient. Observe that an almost isotropic tissue lets λ be near constantly κ and nearlyturns off λ and therewith the influence of subcellular dynamics. In this case, without the effectof S ( t, x ; q ) there would only be myopic diffusion, as e.g., in [72]. So far we considered the space variable x ∈ R N , however the brain occupies a well delimitedregion inside the skull. Therefore we consider a bounded space domain Ω ⊂ R N and assume itto have a smooth enough boundary. Through the rescaling x → εx the domain on which (3.25)holds is ˜Ω = ε Ω, having outer unit normal vector ν ( x ) at x ∈ ∂ ˜Ω. In order to determine thecorresponding boundary conditions on ∂ ˜Ω we assume that there is no normal mass flux acrossthe boundary [59], which translates into the mesoscopic no-flux condition [74] (cid:90) V vp ( t, x, v, z ) · ν ( x ) dv = 0 , for all x ∈ ∂ ˜Ω , t > . (3.28)This condition actually means that the normal component of the macroscopic ensemble ve-locity U ( t, x ) = (cid:82) V vm ( t, x, v ) dv of the tumor across the space boundary vanishes (the tumorcannot leave the brain). Following [74] we write the boundary of the phase space as ∂ ˜Ω × V × Z = (Γ + ∪ Γ − ∪ Γ ) × Z, whereΓ ± := { ( x, v ) ∈ ∂ ˜Ω × V : ± v · ν ( x ) > } , Γ := { ( x, v ) ∈ ∂ ˜Ω × V : v · ν ( x ) = 0 } . We assume that Γ has zero measure w.r.t. the Lebesgue measure on ∂ ˜Ω × V and considerthe trace spaces L ± := L (Γ ± × Z ; | v · ν ( x ) | dσ ( x ) dvdz ) . Moreover, p is supposed to be regular enough so that we can define the traces p | Γ ± × Z ∈ L ± ,and that for a fixed t > p | ∂ ˜Ω × V × Z ( t, x, v, z ) = lim ˜ x ∈ ˜Ω˜ x → x p ( t, ˜ x, z ) , for each x ∈ ∂ ˜Ω . Assuming that a regular Hilbert expansion is valid in ˜Ω we can therefore compute the trace bysimply passing to the corresponding limit in the Hilbert expansions for p ( t, x, v, z ) and accordinglyalso for the moments (3.19). Thus, the no-flux condition (3.28) becomes (cid:90) V v ( p ( t, x, v, z ) + εp ( t, x, v, z )) dv · ν ( x ) + O ( ε ) = 0 , x ∈ ∂ ˜Ω , z ∈ Z, t > . This condition should not depend on ε , therefore we should have (cid:90) V vp j ( t, x, v ) · ν ( x ) dv = 0 for all j ≥ , x ∈ ∂ ˜Ω , z ∈ Z, t > . ndeed, for j = 0 we already have this condition satisfied for undirected tissue: (cid:90) Z (cid:90) V vp ( t, x, v, z ) dvdz · ν ( x ) = (cid:90) V vm ( t, x, v ) dv · ν ( x )= (cid:90) V v q ( t, x, ˆ v ) ω ρ ( t, x ) dv · ν ( x )= ρ ˜ E q ( t, x ) · ν ( x ) = 0 . Next we look at the condition for j = 1, more precisely at its integration w.r.t. z : (cid:90) Z (cid:90) V vp ( t, x, v, z ) dvdz · ν ( x ) = (cid:90) V vm ( t, x, v ) dv · ν ( x ) (3.23) = (cid:16) − λ ∇ x · ( D T ρ ) + λ λ g ( Q, λ ) D T ∇ Qρ + αλ Σ( t, x ) ρ (cid:17) · ν ( x ) , which leads to a typical no-flux boundary condition (cid:16) − λ ∇ x · ( D T ρ ) + λ λ g ( Q, λ ) D T ∇ Qρ + αλ Σ( t, x ) ρ (cid:17) · ν ( x ) = 0 on ∂ ˜Ω , t > h by wayof S ( t, x ; h ) in Σ( t, x ). For this ab initio macroscopic equation we can directly impose a no-fluxcondition: D H ∇ x h = 0 on ∂ ˜Ω , t > . (3.30)In previous works [25–27, 38–40] the quantities relating to the tissue, i.e. q and (where ap-plicable) Q were assumed to be time-invariant, whereby Q itself was estimated from the DTIdata, as in [25–27, 42]. No tissue degradation was accounted for, and one could also include aproliferation term, e.g. like in [26]. In this much simplified situation only coupling (3.25) with(2.5) and using the no-flux boundary conditions given above, the system takes the form of aKeller-Segel problem with some supplementary drift terms, all of which are linear in ρ .When the full evolution of the tissue becomes relevant, then (3.25) and (2.5) have to besupplementally coupled with the dynamics of q , as given by (2.14), and with that of Q , given by(2.11). The resulting (reaction-)diffusion-taxis system then characterizes cell, acidity, and tissuedynamics evolving on two scales (macroscopic and mesoscopic, respectively). Thereby, we dealwith the macroscopic cell density in the first PDE, while (2.14) and (2.11) involve the mesoscopicquantity p ( t, x, v, y ), which is inconvenient both for the analysis and the numerics.Remark that, thanks to property 3 of L , the Hilbert expansion is equivalent to splitting p as p ( t, x, v, y ) = ¯ p ( t, x, y ) M ( t, x, v ) + εM ⊥ ( t, x, v ), with (cid:82) V M ⊥ ( t, x, v ) dv = 0, i.e. the Chapman-Enskog expansion used in [38]. The leading order of the operator in (2.8) leads toΠ a [ q ]( t, x, θ ) (cid:39) (cid:26) (cid:82) S N − | θ · θ (cid:48) | q ( t, x, θ (cid:48) ) dθ (cid:48) for undirected tissue (cid:82) S N − θ · θ (cid:48) q ( t, x, θ (cid:48) ) dθ (cid:48) for directed tissue, (3.31)along with correspondingly rewriting (2.12) as A [ q ]( t, x ) = (cid:90) S N − Π a [ q ]( t, x, θ ) q ( t, x, θ ) dθ. (3.32)Therewith we obtain from (2.14) the equation ∂ t q ( t, x, θ ) = r D ( h ) ρ ( t, x ) q ( t, x, θ ) (cid:16) Π a [ q ]( t, x, θ ) − A [ q ]( t, x ) (cid:17) (3.33)and from (2.13) ∂ t Q = r D ( h ) ρQ ( A [ q ]( t, x ) − . (3.34)Thus, in order to determine the dynamics of tumor cells in interaction with the tissue they egrade and with the acidity they produce, we have to solve the macro-meso system ∂ t ρ = ∇ x · (cid:104) λ ( t, x ) (cid:16) ∇ x · ( D T ρ ) − λ ( t, x ) g ( Q, λ ) D T ∇ Qρ − αS ( t, x ; q ) ρ + αS ( t, x ; h ) ρ (cid:17)(cid:105) + µ ( ρ, h ) Qρ (3.35a) ∂ t q ( t, x, θ ) = r D ( h ) ρ ( t, x ) q ( t, x, θ ) (cid:16) Π a [ q ]( t, x, θ ) − A [ q ]( t, x ) (cid:17) (3.35b) ∂ t Q = r D ( h ) ρQ ( A [ q ] −
1) (3.35c) h t = D H ∆ h + aρ ρ − bQ ( h − h ) (3.35d) ∂ t n = r D ( h ) ρQ (1 − A [ q ]) + F ( h, ρ, Q ) , (3.35e)where D T is as in (3.22) with E q = 0, with S and S as in (3.24), Π a [ q ] as in (3.31), A [ q ]as in (3.32), µ as in (3.26), F as in (3.27), with boundary conditions (3.29), (3.30), and withgiven initial conditions ρ (0 , x ), q (0 , x, θ ), Q (0 , x ), h (0 , x ), and n (0 , x ). These can be the tumorcell distribution (or an approximation of it) observed at diagnosis, the directional distribution offibers obtained via DTI, some estimate of the macroscopic volume fraction of the tissue (e.g., mostsimply F A , as in [23, 25]), some (estimated) acidity distribution at diagnosis, and the necrotictissue distribution, respectively. A tumor segmentation of the diagnosis image could be useful inassessing the latter.The mathematical handling of system (3.35) is nontrivial, both with respect to well posednessand numerics. The equations connect two modeling levels (macroscopic and mesoscopic) and thecouplings via q involved in the coefficients of all terms on the right hand side of (3.35a) render theproblem highly nonlinear. Moreover, the equation for ρ features (along with the myopic diffusion)three types of taxis: • macroscopic haptotaxis towards ∇ Q , • a new kind of mesoscopic haptotaxis (term with S ), where the bias is actually given by ∇ x q , and • pH-taxis, describing chemorepellence due to acidity (term with S ). The system (3.35) is macroscopic with respect to the cell dynamics: by an asymptotic expansionof p ( t, x, v, y ) around the local Maxwellian q ( t, x, θ ), we derived a model for the local cell density ρ ( t, x ). However, the tissue is modeled by the mesoscopic quantity q ( t, x, θ ). To reduce the levelof detail in the tissue dynamics to match the rest of the model, we derive a low-order momentapproximation to (3.35b). For simplicity, we only consider the two-dimensional case.Recall the peanut distribution [38, 72]: q ( θ ) = 1 (cid:104) θ (cid:62) D W θ (cid:105) θ (cid:62) D W θ = 1 π tr D W θ (cid:62) D W θ. We interpret this as a P -approximation [15] [75, III.5] under some additional constraints. Asecond-order basis of monomials is given by aaa = (cid:0) θ x θ y θ x θ x θ y θ y (cid:1) . We drop the constant function from the basis, because it would introduce a linear dependence: θ x + θ y = 1. The P ansatz is defined by the linear combination of basis functions q ( θ ) = ααα · aaa, wherein ααα is a vector of multipliers such that the moment constraints uuu := (cid:104) aaaq (cid:105) = (cid:104) aaa q (cid:105) = (cid:10) aaaaaa (cid:62) (cid:11) k are fulfilled. Due to the symmetry constraints (cid:104) θq (cid:105) = 0, we drop the first-order monomials θ x , θ y from the basis. The remaining basis functions are then: bbb = (cid:0) θ x θ x θ y θ y (cid:1) . ith this choice it is easy to associate the multipliers βββ in the ansatz q = βββ · bbb ( θ )with the components of the normalized water diffusion tensor: βββ := (cid:0) β xx β xy β yy (cid:1) = 1 π tr D W (cid:0) D W xx D W xy D W yy (cid:1) . Moreover, the moments := (cid:104) bbbq (cid:105) are directly related to the components of the (normalized)tumor diffusion tensor P = (cid:10) θθ (cid:62) q (cid:11) : := (cid:0) w xx w xy w yy (cid:1) = (cid:0) P xx P xy P yy (cid:1) . We can translate multipliers to moments with the moment constraints = (cid:104) bbb q (cid:105) = (cid:10) bbbbbb (cid:62) (cid:11) βββ . Thetransfer matrix is given by H := (cid:10) bbbbbb (cid:62) (cid:11) = π . With these tools in hand, we derive a P -approximation of the tissue dynamics (3.35b) throughthe following steps: Insert the ansatz q into (3.35b), multiply by the basis bbb and integrate over S to obtain ∂ t = r D ( h ) ρ ( (cid:104) bbb q Π [ q ] (cid:105) − [ q ]) . (3.36)It remains to calculate the moments (cid:104) bbb q Π [ q ] (cid:105) . Inserting the definitions of q and Π, we obtain (cid:104) bbb q Π [ q ] (cid:105) = (cid:90) S (cid:90) S bbb ( θ ) βββ · bbb ( θ ) βββ · bbb ( θ (cid:48) ) | θ · θ (cid:48) | d θ (cid:48) d θ. With the tensor P k,ij = (cid:90) S (cid:90) S b k ( θ ) b i ( θ ) b j ( θ (cid:48) ) | θ · θ (cid:48) | d θ (cid:48) d θ we can write the k -th component (cid:104) bbb q Π [ q ] (cid:105) k as (cid:104) bbb q Π [ q ] (cid:105) k = βββ (cid:62) P k βββ = (cid:62) (cid:0) H − (cid:62) P k H − (cid:1) Finally, we obtain A [ q ] from the identity θ x + θ y = 1: A [ q ] = (cid:104) q Π [ q ] (cid:105) = (cid:10) ( θ x + θ y ) q Π [ q ] (cid:11) = (cid:104) bbb q Π [ q ] (cid:105) xx + (cid:104) bbb q Π [ q ] (cid:105) yy . The tensor P k,ij can be precomputed once and for all with a high-order quadrature. To evaluateeach component of the right-hand side of (3.36) we need to compute only a quadratic form at runtime. Note that the normalization (cid:104) q (cid:105) = 1 results in the loss of an additional degree of freedom.It holds tr P = 1, therefore we can reconstruct w yy = 1 − w xx and only need the evolve the twomoments w xx , w xy in (3.36). Remark:
In three dimensions, the previous considerations are completely analogous for thebasis bbb = (cid:0) θ x θ y θ z θ x θ y θ x θ z θ y θ z (cid:1) . Here we also aim to investigate the effect of reducing diffusivity that might not be experimentallyconsistent. In order to achieve this objective and for the sake of completeness, we propose in thissubsection to perform a hyperbolic limit that will provide us with a macroscopic vision in whichthe terms of transport and potentials win the battle over diffusion. We will take advantage ofmost of the calculations in Subsection 3.1, that we will omit in part so as not to be repetitive.Finally, in order to try to model a series of effects of a priori minor influence according to theexperiments, not contemplated in the model variables, we carry out a ”small” extension to the econd order of the hyperbolic expansion and compare this double development on the scale withthe parabolic approach.In the following we address a version of (3.16) where we neglect the dependence of the turningrate on the subcellular dynamics, i.e. we consider it to be of the form λ ( x ) = κF A ( x )+1 and areinterested in the situation of directed tissue.We thus consider the kinetic transport equation ∂ t p + v ·∇ x p − ∂ z (cid:16) (( k + Q + k − ) z + f (cid:48) ( Q )( Q t + v · ∇ Q )) p (cid:17) + α div v (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) p (cid:17) = λ ( x ) L ( p ) + P ( p, Q, h, ρ ) , (3.37)with L ( p ) = ¯ p ( t, x, z ) q ( x, ˆ v ) ω − p ( t, x, v, z ) having the same properties as in Subsection 3.1 andperform a hyperbolic scaling t → εt , x → εx , while the proliferation term is rescaled as previouslywith ε . This leads to ε∂ t p + εv ·∇ x p − ∂ z (cid:16) (( k + Q + k − ) z + εf (cid:48) ( Q )( Q t + v · ∇ Q )) p (cid:17) + αε div v (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) p (cid:17) = λ ( x ) L ( p ) + ε P ( p, Q, h, ρ ) . (3.38)We now consider a Chapman-Enskog expansion (equivalent here to the Hilbert one, as statedin Subsection 3.1), i.e., a decomposition of p into a Ker ( L )-component and a Ker ( L ) ⊥ -part,as follows: p ( t, x, v, z ) = ¯ p ( t, x, z ) q ( t, x, ˆ v ) ω + εp ⊥ ( t, x, v, z ) , (3.39)where p ⊥ ∈ (cid:10) qω (cid:11) ⊥ verifies (cid:82) V p ⊥ ( t, x, v, z ) dv = 0. This decomposition leads to the correspondingexpansion for the moments of p , in particular: m ( t, x, v ) = ρ ( t, x ) q ( t,x, ˆ v ) ω + εm ⊥ ( t, x, v, z ), andthen, integrating (3.38) w.r.t. z ∈ Z we obtain ε∂ t ρ qω + ερ∂ t ( qω ) + ε ∂ t m ⊥ + εv · ∇ x ( ρ qω ) + ε v · ∇ x m ⊥ + αε ∇ v · (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) m (cid:17) = ελ ( x ) L ( m ⊥ ) + ε µ ( ρ, h ) Qm. (3.40)Integrate w.r.t. v and divide by ε to get ρ t + ∇ x · ( ρ ˜ E q + ε (cid:90) V vm ⊥ dv ) = ερµ ( ρ, h ) Q, (3.41)where as before ˜ E q ( t, x ) = (cid:82) V v q ( t,x, ˆ v ) ω dv = s E q ( t, x ).Clearly (3.41) is drift-dominated; however, it is worth computing the O ( ε ) correction withrespect to the pure drift. From (3.40) and (3.41) follows λ ( x ) L ( m ⊥ ) = ε qω ρµ ( ρ, h ) Q − qω ∇ x · (cid:16) ρ ˜ E q + ε (cid:90) V vm ⊥ dv (cid:17) + ρ ∂ t qω + ε∂ t m ⊥ + v · ∇ x ( ρ qω + εm ⊥ )+ α ∇ v · (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) m (cid:17) − ερµ ( ρ, h ) Q. (3.42)Now observe that the integral w.r.t. v of the right hand side in (3.42) vanishes, so that wecan take the (pseudo)inverse of L giving (at leading order) m ⊥ (cid:39) − λ ( x ) (cid:104) qω ( v − ˜ E q ) · ∇ ρ + ρ (cid:16) v · ∇ x qω − qω ∇ x · ˜ E q (cid:17) + ρ ∂ t qω (3.43)+ α ∇ v · (cid:16) B ( v )( ∇ x q (ˆ v ) − ∇ x h ) m (cid:17)(cid:105) , (3.44)hence (cid:90) V vm ⊥ dv (cid:39) − λ ( x ) (cid:16) (cid:90) V v ⊗ ( v − ˜ E q ) qω dv ∇ ρ (3.45)+ ρ (cid:16) ∇ x · (cid:90) V v ⊗ v qω dv − ˜ E q ∇ x · ˜ E q (cid:17) + ρ∂ t ˜ E q − αρ Σ( t, x ) (cid:17) = − λ ( x ) (cid:16) ∇ x · ( D T ρ ) − ρ ˜ E q ∇ x · ˜ E q + ρ∂ t ˜ E q − αρ Σ( t, x ) (cid:17) , (3.46) here D T ( t, x ) and Σ( t, x ) := S ( t, x ; q ) − S ( t, x ; h ) are as introduced in Section 3.1 and weobserve that (cid:90) V v ⊗ ( v − ˜ E q ) qω dv = (cid:90) V ( v − ˜ E q ) ⊗ ( v − ˜ E q ) qω dv = D T . Together with (3.41) this leads to the macroscopic PDE ρ t + ∇· ( ρ ˜ E q ) = ε ∇· (cid:32) λ ( x ) (cid:16) ∇· ( ρ D T ( t, x ))+ ρ ( ∂ t ˜ E q − ˜ E q ∇· ˜ E q − α Σ( t, x )) (cid:17)(cid:33) + ερµ ( ρ, h ) Q, (3.47)which is drift-dominated. Notice that for E q = 0, i.e. if the tissue is undirected, the correctionterm in (3.47) has the same form as the right hand side of (3.25) - except for the middle termtherein, which got lost when integrating (3.38) w.r.t. z . In fact, it was the assumption of theturning rate not depending on z which effaced the whole influence of subcellular dynamics. Theeffect of this is not having the macroscopic haptotaxis term in the ε -correction on the righthand side. Nevertheless we still get the pH-taxis and (mesoscopic) haptotaxis correction termscontained in ε ∇ · (cid:16) αλ ( x ) ρ Σ( t, x ) (cid:17) = ε ∇ · (cid:16) αλ ( x ) ρ ( S ( t, x ; q ) − S ( t, x ; h )) (cid:17) . Coupling (3.47) with the equations (3.35b) and (3.35d) describing tissue and acidity evolution,respectively, leads to challenges similar to those in Section 3.1.We need to prescribe boundary conditions. To this aim we consider again a bounded, suffi-ciently smooth domain ˜Ω = ε Ω ⊂ R N and start from a mesoscopic no-flux condition of the form(3.28). Using the Chapman-Enskog expansion (3.39) of p ( t, x, v, z ) we rewrite this condition as (cid:16) ρ ( t, x )˜ E q ( x ) + ε (cid:90) Z (cid:90) V vp ⊥ ( t, x, v, z ) dvdz (cid:17) · ν ( x ) = 0 , for all x ∈ ∂ ˜Ω , t > , which in virtue of (3.45) leads to (cid:16) − ρ ˜ E q + ελ ( x ) (cid:16) ∇ · ( ρ D T ( t, x )) + ρ ( ∂ t ˜ E q − ˜ E q ∇ · ˜ E q − α Σ( t, x )) (cid:17)(cid:17) · ν ( x ) = 0 , x ∈ ∂ ˜Ω , t > , (3.48)meaning that the normal macroscopic flux and its ε -correction vanish on the spatial boundary.Thus, in the case with directed tissue and when no subcellular dynamics are taken intoaccount, the macro-meso system to be solved becomes ρ t + s ∇ · ( ρ E q ) = ε ∇ · (cid:32) s λ ( x ) (cid:16) ∇ · ( ρ V q ( t, x )) + ρ ( 1 s ∂ t E q − E q ∇ · E q − αs Σ( t, x )) (cid:17)(cid:33) + ερµ ( ρ, h ) Q (3.49a) h t = D H ∆ h + aρ ρ − bQ ( h − h ) (3.49b) ∂ t q ( t, x, θ ) = r D ( h ) ρ ( t, x ) q ( t, x, θ ) (Π[ p ]( t, x, θ ) − A [ p ]( t, x )) , (3.49c) ∂ t Q = r D ( h ) ρQ ( A [ p ] −
1) (3.49d) ∂ t n = r D ( h ) ρQ (1 − A [ p ]) + F ( h, ρ, Q ) , (3.49e)with Π[ p ] and A [ p ] as given in (2.8) and (2.12), respectively. Observe that upon assuming thetisue to be undirected (i.e. that E q = 0), (3.49a) takes the form of (3.35a), which was obtained byparabolic scaling, with the difference of the right hand side of the glioma dynamics being scaledhere by ε .As for the parabolic limit, in (3.49) the dependence of the featured operators on the mesoscopiccell distribution function p is inconvenient, so we use the same approach as in Subsection 3.1.1to approximate Π[ p ] and A [ p ] by (3.31) and (3.32), respectively. Notice that in the present case(directed tissue) we can rewrite Π a [ q ]( t, x, θ ) = θ · E q ( t, x ) A [ q ]( t, x ) = E q ( t, x ) · E q ( t, x ) onsequently, (3.49c) can be replaced by ∂ t q ( t, x, θ ) = r D ( h ) ρ ( t, x ) q ( t, x, θ )( θ · E q ( t, x ) − E q ( t, x )) , hence we can compute ∂ t E q = r D ( h ) ρ ( V q + E q ⊗ E q − E q I N ) · E q and plug it into (3.49a). Let us denote T q := 1 s r D ( h ) ρ ( V q + E q ⊗ E q − E q I N ) · E q − E q ∇ · E q . Then altogether the system becomes ρ t + ∇ · ( ρ E q ) = ε ∇ · (cid:32) s λ ( x ) (cid:16) ∇ · ( ρ V q ( t, x )) + ρ ( T q − αs Σ( t, x )) (cid:17)(cid:33) + ερµ ( ρ, h ) Q (3.50a) h t = D H ∆ h + aρ ρ − bQ ( h − h ) (3.50b) ∂ t q = r D ( h ) ρq ( θ · E q − E q ) (3.50c) ∂ t Q = r D ( h ) ρQ ( E q −
1) (3.50d) ∂ t n = r D ( h ) ρQ (cid:0) − E q (cid:1) + F ( h, ρ, Q ) , (3.50e)with θ ∈ S N − , t >
0, and x ∈ ˜Ω. The functions µ ( ρ, h ) and F ( h, ρ, Q ) are still those given in(3.26) and (3.27), respectively. The equations are completed by initial conditions as in Subsection3.1.1 and by the no-flux boundary conditions ∇ h · ν = 0 on ∂ ˜Ω and (cid:16) − sρ E q + s ελ ( x ) (cid:16) ∇ · ( ρ V q ( t, x )) + ρ ( T q ( t, x ) − αs Σ( t, x )) (cid:17)(cid:17) · ν ( x ) = 0 , x ∈ ∂ ˜Ω , t > . (3.51) In this section the equations obtained by parabolic upscaling are solved numerically. Simulationsare conducted for three scenarios of increasing complexity. In each scenario we describe a set ofassumptions simplifying the equations. This allows for directly investigating the relative impor-tance of the different effects included in our model. The results are presented in Subsection 4.3.For convenience we state in brevity the exact macroscopic systems of PDEs for the consideredscenarios:
Scenario 1:
We assume that neither q nor Q depend on time, thus can be assessed from the DTI data.We follow previous works, starting from a pre-assigned directional distribution q (peanut)and independently estimate Q using the approach in [26,41,42] relying on Brownian motiondescription of water molecule diffusion. The system then reduces to (3.35a), (3.35d), merelyfeaturing myopic diffusion, pH-taxis, and transport terms with drift velocities computedfrom the data, together with the necrosis dynamics (3.35e) with r D ( h ) = 0, i.e. ∂ t ρ = ∇ x · (cid:104) λ ( x ) (cid:16) ∇ x · ( D T ρ ) − λ ( x ) g ( Q, λ ) D T ∇ Qρ − αS ( t, x ; q ) ρ + αh S ( t, x ; h ) ρ (cid:17)(cid:105) + Qρη (1 − ρ − n ) ( h T − h ) + h − γ hh ρQ (4.52a) ∂ t h = D H ∆ h + aρ ρ − bQ ( h − h ) (4.52b) ∂ t n = F ( h, ρ, Q ) . (4.52c) Alternative, more precise ways to assess the (macroscopic and/or mesoscopic) tissue density from DTI data wouldbe to perform density estimation by using statistical or variational methods, see e.g. [24, 30, 50, 90]. Here we do notaddress this issue, as our focus is on providing a continuous approach to modeling glioma growth and spread which isamenable to efficient computations and at the same time able to capture lower scale effects. cenario 2: Q evolves in time, meaning that the tissue is degrading irrespective of the fibre orientations(hence the anisotropy does not play any role in the degradation). Here A [ q ] is zero. Thiscorresponds to an indirect depletion, i.e. not upon contact with glioma cells, but ratherthrough the acidity they produce. Consequently, the macroscopic system (3.35) simplifiesto ∂ t ρ = ∇ x · (cid:104) λ ( x ) (cid:16) ∇ x · ( D T ρ ) − λ ( x ) g ( Q, λ ) D T ∇ Qρ − αS ( t, x ; q ) ρ + αS ( t, x ; h ) ρ (cid:17)(cid:105) + Qρη (1 − ρ − n ) ( h T − h ) + h − γ hh ρQ (4.53a) ∂ t Q = − r D ( h ) ρQ (4.53b) ∂ t h = D H ∆ h + aρ ρ − bQ ( h − h ) (4.53c) ∂ t n = r D ( h ) ρQ + F ( h, ρ, Q ) . (4.53d) Scenario 3:
Both q and Q evolve and the full system (3.35) is solved. The equations were discretized by the method of lines approach in all three cases. For thespatial discretization we used a vertex centered finite volume method. The mesh is naturallystructured because of the underlying medical dataset. The mesh width h = 2 mm is imposed bythe particular DTI dataset, which is publicly available [17]. The reconstructed diffusion tensorswere averaged with a component-wise arithmetic mean at the control volume interfaces. Theadvective fluxes were approximated by a first order upwind discretization. The no-flux barrierbetween the brain and the background were implemented by suppressing flux assembly at therelevant cell interfaces . The integration of the reaction terms coupling the evolving fields wasperformed by a first order midpoint rule at the cell centers. As an initial condition we used asmall Gaußian peak as an approximation of the tumor and took for the spatial distribution of theacidity field h the uniform steady-state of (3.35d). The initial condition for necrotic matter wassimply n = 0. The time discretization was realized by an implicit Euler scheme. The resultingnonlinear system was solved via Newton iteration which internally used a BiCGSTAB solverpreconditioned via AMG with an SSOR smoother [83]. The implementation was performed within the DUNE software framework [3,4], specifically withthe dune-pdelab discretization module [5,12]. The BiCGSTAB solver and the AMG were providedby the iterative solver template library dune-istl [11].
In Table 1 we present the parameters used throughout all numerical simulations, their meaning,and their respective sources. We can also assess the parabolic scaling parameter ε = sκX (cid:39) · − , where we use X = 0 . m as a reference length, corresponding to the side length of he plots ine.g., Figure 2). This value indicates that the parabolic limit is an adequate approximation to thekinetic model for the chosen parameters.Let us also highlight here the role of the scaling parameter α . It is influenced by the positiongradients of q and h (involved in Σ) and hence should be a large quantity, say of the order γ ε − γ forsome γ >
0. This is similar to what happens when we perform hyperbolic scaling in a fairly similarkinetic framework (see [8, 77]), where hydrodynamic limits with dominant drift are obtained, butwhere the small diffusion correction of a non-diffusive limit can be made explicit by displayingthe dependency of ε on all scaling parameters. In the present context γ should be a parameterrelated to smaller scales, such as the microtubule (MT) extension zones that are responsible forthe biochemical and biomechanical exchange of cells with their environment [34, 54]. The length here we mean by ’cell’ the numerical discretization element and not the biological entity f microtubules is between 5 and 8 cell units [76], whereas a glioma cell has a diameter of ca.15 µm = 15 × − mm [1]. This means that the regions of MTs that would correspond to thetumor front have a length around 100 µm = 0 . mm . Then, the parameter γ will represent theratio of such MT area with respect to the whole tumor; here we take γ = 0 . ε , the parameter α is taken here of the order of 10 .The coefficient functions involved in the simulations are summarized below:For all scenarios: F ( h, ρ, Q ) = γ hh ρQ,λ = κf ( Q ) F A + f ( Q ) , λ = κF A ( F A + f ( Q )) , f ( Q ) = k + Qk + Q + k − f (cid:48) ( Q ) = k + k − ( k + Q + k − ) , g ( Q, λ ) = f (cid:48) ( Q ) k + Q + k − + λ , r D ( h ) = r (cid:18) hh − (cid:19) + . For Scenario 1: Q ( x ) = 1 − (cid:16) tr ( D W ( x )) λ max (cid:17) , D T ( x ) = s ( N + 2) (cid:16) I N + 2 D W ( x ) tr ( D W ( x )) (cid:17) , with λ max denoting the maximum eigenvalue of the water diffusion tensor D W . parameter value meaning source s
50 ( µ m/h) (= 1.389e-8 m/s) cell speed [67] α (1) weight for advective terms estimated κ h h T D H . e −
03 ( mm /s ) diffusion coefficient of protons [60] a . e −
17 (mol cm ( cells · sec ) − ) acid production rate [65] b η k + k − γ r Table 1: Full set of parameters for the numerical simulations
We focus the presentation of the numerical results on those aspects where the new model andits extension, in particular by including tissue dynamics, produce new qualitative effects. Com-putations were conducted with a simulated time of 52 weeks. We first present line plots fromthe coordinate origin through the tumor to the upper right of the domain in order to visualizeand discuss the evolution of the involved fields. This helps to establish the general dynamicsof the model in Subsection 4.3.1. We present the final states of the involved fields in 4.3.2 andinvestigate in Subsection 4.3.3 which advective effect is playing a dominant role with respect tothe tumor extent and appearance. We then directly compare the solutions of the three consid-ered scenarios in Subsection 4.3.4, and address in Subsection 4.3.5 the grading of a tumor uponrelying on the numerical results, particularly with respect to the amount of necrotic matter. Wealso suggest in Subsection 4.3.6 how the simulations of the model could be used to perform dosepainting for radiation treatment.
Figure 1 shows the amplitude of the four fields ρ, n, Q for densities of tumor cells, necroticmatter, and macroscopic tissue, respectively, and for the proton concentration field h convertedto its corresponding pH -value (right scale) at three distinct points in time. The uppermost plotindicates the initial condition we chose for the ρ - and h fields, progressing then downwards withthe simulation time. s the tumor density ρ increases due to the proliferation (thus also leading to an increase inacidity concentration) the cells start to accumulate at the interface with tissue due to pH-taxisand meso- and macroscopic haptotaxis. The sustained expression of protons (lower pH) leadsto necrotic core formation and growth. The total neoplasm is considered to be made up by thedensities of necrotic matter n and living tumor cells ρ . At the final stage of the simulation, theformer represents the largest part of the neoplasm ( n > . ρ < .
2) of living tumor cells. This is consistent with high grade tumors typicallyobserved in clinical practice [35, 78]. . . . . x/cm ρ n Q . . x/cm p H h . . . . x/cm
26 weeks . . x/cm p H . . . . x/cm
52 weeks . . x/cm p H Fig. 1: Profiles of tumor density ρ , acidity h (as pH), necrosis n , and tissue density Q over 52 weeks.Scenario 3: the tissue is dynamically degraded depending on its local orientation and the acidityconcentration; a necrotic core is formed and develops, becoming the main part of the neoplasm.Available in color online. In Figure 2 we show the final states of the four fields of Scenario 3, overlain with the main axisof the water diffusion tensor D w from the DTI data sets, as well as the direction of the combinedadvective fields. The dominant S term describing repellent pH-taxis results in active transportof tumor mass away from the center. Similarly to the line plots presented earlier in Figure 1, wecan identify a reduction of active tumor cells and tissue degradation within the acidic core areaof the neoplasm, where necrosis is growing instead. x/cm y / c m ρ . x/cm y / c m pH . . .
26 7 8 9 10 11 126789101112 x/cm y / c m n . x/cm y / c m Q . Fig. 2: Dynamic fields of solution components ρ, pH = − log ( h ) , n, Q after 52 weeks in a computa-tion for Scenario 3. Red arrows: main axes of D W . Green arrows: total advection vector. Top left:density ρ of living tumor cells. Top right: pH value calculated from the acidity field pH = − log ( h ).Bottom left: density n of necrosis. Bottom right: macroscopic density Q of normal tissue. Availablein color online. 20 .3.3 Scenario 3: study of dominant advective effects To evaluate the relative importance of the four advective terms within the model, we investigatedtheir dynamically changing magnitudes. We include the dynamic scaling term 1 /λ ( x ) in thesepresentations. Figure 3 shows the fractional anisotropy of the underlying data set, the magnitudesof S , S , and the magnitude of all combined advective terms at 52 weeks simulation time. Forthe literature-based parameter set in Table 1 the chemotactic term S is orders of magnitudesstronger than the other advective fields in this late stage of tumor progression. In the bottomright plot of Figure 3 it can be observed that the dominant advective effect at the proliferationfront is given by S , while the other effects may only contribute significantly in the absence ofacidity gradients, i.e. beyond the visible tumor margins.We also investigated in Figure 4 the relative importance of the different advective effects bysuccessively excluding them in the simulation runs. That way, the total time-integrated effect ofinclusion and exclusion of the terms could be investigated. The repellent pH-taxis described bythe S -term was found to have the most dominant effect on the tumor density, but the haptotaxisterms also contribute to establishing low-density areas at the tumor edges. Given the infiltrativespread of glioma cells and post-treatment tumor recurrence which is mainly due to such relativelysmall cell aggregates located further away from the main tumor mass, such margins might berelevant for treatment planning, especially as far as radiotherapy is concerned (see also Subsection4.3.6). In the following we compare the simulation results for Scenarios 1-3 with identical initial condition,in order to find out whether there are any potentially relevant differences. Figure 5 shows thedensity ρ + n of the neoplasm after 52 weeks. As expected, Scenario 3 involving three-fold taxisand myopic diffusion predicts the largest tumor spread, with the most substantial difference beingthat w.r.t. Scenario 1. It is a modeling problem in itself to predict how living cells and necroticmatter influence the gray scale images in different ways. A direct comparison of Scenarios 2 and3 with medical images is therefore not possible without knowing the underlying mapping of thetwo species to the MRI gray scale imaging. However, the distinction of ρ and n allows deeperinsights into the dynamics of the tumor expansion. The two plots in the second row of Figure 5illustrate the differences between the densities ρ + n after 52 simulated weeks, computed uponusing Scenarios 3-2, and 3-1. Thus, the latter comparison shows that taking into account theevolution of normal tissue is of high relevance, while considering the dynamics of mesoscopic tissuedistribution q (i.e the step from Scenario 2 to 3) is contributing less (please note the differentscales of the color maps therein). Hence, the evolution of macroscopic tissue density Q seems tocontribute most to the prediction of tumor spread.On the other hand, the dynamics of mesoscopic tissue q contributes to substantially modifyingthe fractional anisotropy FA, which is a key observable in the DTI imaging of brain tumors. Ournew model allows to directly model the depletion of brain structure and the therewith associatedchanges in FA, as indicated in Figure 6. The two plots of the first row show the initial fractionalanisotropy and how it is diminished by the tumor growth after 52 weeks, again overlain by thedirection of the total advective fields. The plot on the second row illustrates the difference in FAbetween the final and the initial simulation times.The results indicate that tissue dynamics can be seen as a downstream effect, after tumor densityand acidity have changed from baseline. The differences in F A are therefore most prominentcloser to the tumor center, and not so much at the invasion edge. Scenario 3 is more complexand requires more effort to treat, but beyond the differences put in evidence in Figures 5 and 6 italso provides a framework for a careful description of macroscopic tissue dynamics - along withits effects on the other solution components.Still with the aim of investigating the influence of evolving tissue, we also compare a reducedversion of Scenario 3 with a previous model obtained in [72] and accounting for the effect of tissueanisotropy on the migration of glioma cells, however with fixed tissue. Specifically, we comparethe PDE ρ t = ∇∇ : (cid:18) λ D T ρ (cid:19) + ηρ (1 − ρ ) (4.54)which is a slight modification of that in [72], by adding the proliferation term and letting λ be x/cm y / c m F A . x/cm y / c m | S | · − x/cm y / c m | S | · − x/cm y / c m combined advection · − Fig. 3: Top left: fractional anisotropy of the camino dataset [17]. Contours: neoplasm density ρ + n from Scenario 3 after 52 weeks. Red arrows: main axis of D W . Green arrows: total advectionvector. Top right: magnitude of the S term. Bottom left: magnitude of the S term. Bottom right:magnitude of all advective terms (myopic drift, macroscopic haptotaxis, mesoscopic haptotaxis S ,repellent pH-taxis S ). The scales of the color maps differ between the plots. All advective termscontain the scaling factor 1 /λ ( x ). Available in color online.22 x/cm y / c m scenario 3: neoplasm: ρ ( x ) + n ( x ), all adv. fields, 52 weeks . x/cm y / c m scenario 3: no haptotaxis − · − x/cm y / c m scenario 3 no S − . . x/cm y / c m scenario 3 no S − . . Fig. 4: The effect of advective fields on the model predictions. Red arrows: main axis of D W . Greenarrows: total advection vector. All results are obtained for Scenario 3. Top left: simulation resultfor the density ρ + n of the neoplasm after 52 weeks. The rest of the plots show absolute differencesin ρ + n between the full model and Scenario 3 without one of the adective terms. Top right: nomacroscopic haptotaxis. Bottom left: no S advection. Bottom right: no S advection. The scalesof the color maps differ between the plots. Available in color online. nonconstant, but depend on x like in the rest of this paper, with the system ρ t = ∇ · (cid:18) λ ∇ · ( D T ρ ) − λ g ( Q, λ ) D T ∇ Qρ − αS ( q ) ρ (cid:19) + ηρ (1 − ρ ) (4.55a) q t = r ρq (Π a [ q ] − A [ q ]) (4.55b) Q t = r ρQ ( A [ q ] − , (4.55c) x/cm y / c m Scenario 1: ρ + n . x/cm y / c m Scenario 3: ρ + n . x/cm y / c m Scenario 3 - Scenario 2 − · − x/cm y / c m Scenario 3 - Scenario 1 − . . Fig. 5: Comparison of the scenarios. Contours: density ρ + n of the neoplasm computed in theframework of Scenario 3. Red arrows: main axis of D W . Green arrows: total advection vector.Top row: density ρ + n at the final time (52 weeks) for Scenarios 1 and 3. Bottom row: pointwisedifference in ρ + n between Scenarios 3 and 2 and difference between Scenarios 3 and 1. Note: thescales of the color maps differ between the plots. Available in color online.24 oth with the same initial conditions and no-flux boundary conditions and the same function λ . Figure 4.3.4 shows the computation results. Note that none of the two model makes anydistinction between active and necrotic tumor. The complete tumor mass is encoded in ρ . Asbefore, the evolution of tissue leads to substantial changes in FA, although the dynamics of acidityand its influence are absent in these settings. The simulations predict a larger area of gliomaspread when the cells are allowed to degrade the tissue; although there are rather small amountsof tumor cells in the enlarged regions, these might be relevant for tumor recurrence. Model(4.54) leads to higher cell densities at the core of the neoplasm, due to the migrating cells onlyperforming myopic diffusion and not being supplementary driven by haptotaxis. Tumor grade assessment is decisive for the treatment planning of glioma and for patient survivalprognosis. Typically, glioma are classified according to cell activity and tumor aggressiveness,upon relying on a series of indicators, including histological patterns [14, 89] and tumor composi-tion. Among these, grading by the amount of necrosis relative to the whole tumor volume (bothvisible via biomedical imaging) has been addressed e.g., in [35, 55, 78]. In fact, [35] observed that’location and volume of tumors were not statistically significant predictors of survival’, which isin accordance to most conducted studies, as reviewed e.g. in [37]. We therefore address here onlynecrosis-based grading. In [35] the percentage of necrotic matter visible on MRI is related tothe total visible tumor volume. An attempt at tumor grading by way of connecting the resultsof numerical simulations for our model to this definition proves to be difficult without in-depthinformation on how the tissue types (e.g., obtained by segmentation) affect the MRI shading.The numerical results of our model do not contain any thresholding or shading effects, thereforewe define the time-dependent grade G ( t ) ∈ [0 ,
1] of the simulated tumors via: G ( t ) := V n ( t ) V n ( t ) + V ρ ( t ) (4.56)where V n ( t ) and V ρ ( t ) denote integrals over the whole space domain Ω of the necrosis and livingcell densities n ( t, x ) and ρ ( t, x ), respectively. We investigate then the evolution of this quantityover time, upon being guided by the percentage classification in [35]: 0 < G < ≤ G < G ≥
50: grade 3, with the highest grade corresponding to the mostaggressive tumor and the poorest prognosis.Figure 8 illustrates the evolution of the tumor grade for the three considered scenarios. Wefound only tiny differences between the grades for Scenarios 2 and 3, whereas the difference intumor grade progression between Scenarios 1 and 2 is substantial. The definition (4.56) of grade G effectively skips grade 0, as there is no threshold on the visibility of the necrotic tissue density n . The fact that Scenario 1 fails to predict tumor progression into grade 3 in due (biomedicallyrealistic) time might be an indicator of Scenario 1 being insufficient for modeling glioma dynamics,thus endorsing the role played by tissue evolution. Beyond tumor grading based on assessment of necrotic matter and heterogeneity of the neoplasm,the model can provide valuable information about the extent and shape of the tumor and hencehelp establishing the contours for radiation treatment. The ability to predict space-time densitiesof clonogens as well as necrosis opens the way for more accurate tumor segmentation and dosepainting. Figure 9 shows a gray scale image of the computed neoplasm ( ρ + n ) for a simulated timeof 52 weeks. Level sets for the density of living glioma cells are shown, thus providing informationabout the regions where higher radiation doses should be applied in order to eradicate the activecells, at the same time sparing areas where their density is lower (either as a consequence ofnecrosis - in the tumor core, or due to less cells having invaded the healthy tissue at the tumormargins). We proposed in this paper a novel model for glioma invasion, which takes into account theinteraction with the biochemical and biophysical tumor microenvironment, represented here by x/cm y / c m Initial
F A . x/cm y / c m F A after 52 weeks . x/cm y / c m F A difference − . . Fig. 6: Tissue dynamics in Scenario 3. Contours: neoplasm density ρ + n from Scenario 3. Redarrows: main axis of D W . Green arrows: total advection vector. Top left: gray scale plot of F A at the start of the simulation. Top right:
F A after 52 weeks, with visible reduction of anisotropycloser to the tumor centre. Bottom left: difference in
F A between final and initial simulation times.Available in color online. 26 x/cm y / c m Reduced scenario 3: ρ . x/cm y / c m PH model: ρ . x/cm y / c m ρ difference − . . x/cm y / c m F A difference − . . Fig. 7: Comparison between models (4.55) and (4.54). Contours: tumor density ρ from the com-putation of (4.55) (purple: ρ = 0 .
1, blue: ρ = 0 .
25, green: ρ = 0 . D W . Green arrows: total advection vector. Top left: tumor density ρ at the final time (52 weeks)for model (4.55). Top right: tumor density ρ for model (4.54). Bottom left: difference in ρ com-puted with model (4.55) and with model (4.54). Bottom right: Difference in F A between end andbeginning of simulation, computed with model (4.55). Available in color online.27ig. 8: Grade of tumor (4.56) for the three considered Scenarios. G is essentially the same forScenarios 2 and 3, whereas in Scenario 1 (with the same parameter choice) it exhibits significantlyslower progress. Full information is assumed, so that there is no visibility threshold. Scenario 1predicts the progression from grade 1 towards grade 2 after about 28 weeks, whereas in Scenarios2 and 3 it needs only about 5 weeks. The progression to grade 3 occurs after about 16 weeks forScenarios 2 and 3. For Scenario 1 it does not happen within the simulated time of 52 weeks (recallthat the median survival time of patients with glioma is approx. 60 weeks [91]).28 x/cm y / c m Level sets for neoplasm density (Scenario 3) . Fig. 9: Level sets for density ρ of living tumor cells could be used for dose painting in radiotherapy.Gray scale image: density of neoplasm ρ + n , computation from the reduced Scenario 3, simulationof 52 weeks. Level sets: ρ = 0 . ρ = 0 .
075 (green). acidity acting as a repellent and, respectively, oriented tissue fibers guiding the cell migration.Both environmental components, most prominently the acidity, also influence the growth ordepletion of tumor cells.The modeling approach has a multiscale character: it starts from ODE descriptions of singlecell (position, velocity, activity variable) dynamics, which are translated into a kinetic transportequation for the space-time distribution function of cells sharing the same regimes of velocity andactivity variables, and eventually arrives at a macroscopic PDE for the evolution of the wholetumor as a population of such cells performing myopic diffusion and multiple taxis. The cellbehavior is thereby dynamically coupled to that of the environment, the latter featuring a systemof macroscopic and mesoscopic (integro-)differential equations with dependencies on space, time,and orientation of tissue fibers. We refer to [52] for a review of existing models of cell migrationwith multiple taxis and therewith associated challenges. The diffusion and taxis coefficients ofthe PDE obtained here for the cell density and the source terms of the equations for meso- andmacrocopic tissue and necrosis densities involve nonlocalities w.r.t. the orientations of tissuefibers and cell velocities; see [20] for a review of models for cell migration involving this and othertypes of nonlocality.There are several alternative ways to include cell level environmental influences in a KTAPmodeling framework, ultimately leading in the macroscopic limit to taxis terms:(i) by using in the operator describing velocity reorientations turning rates which depend onthe pathwise gradient of one or several signals, as done e.g. in [63, 70];(ii) by describing the dynamics of activity variables for the cell state, which translates intocorresponding transport terms in the kinetic equations and turning rates depending on theactivity variable(s), see [23, 25–27, 41, 42] for glioma interacting with tissue and e.g., [28, 92]for swimming bacteria biased by some chemoattractant;(iii) by accounting for biochemical and/or biophysical effects translated into cell stress and forcesacting on the cells and leading to transport terms w.r.t. the velocity variable in the kineticPDE on the mesolevel, see [19] for the case of cancer cells responding to a ’chemotaxis force’proportional to the gradient of a given chemoattractant concentration;(iv) by considering turning operators where the cell reorientation depends on the interactionwith (some of) the environmental cues and using equilibrium distributions for the momentclosure, as e.g. in [19]. hile (i) and (iv) have a rather mesoscopic flavor, (ii) and (iii) refer closer to microscopic, singlecell descriptions, therefore we employed in this work a combination of the latter: (ii) led to themacroscopic haptotaxis term, while (iii) with the description of cell stress and associated, speed-preserving forces acting on glioma cells yielded the terms characterizing repellent pH-taxis andmesoscopic haptotaxis. Notice that the coupling with evolving acidity and macro- and mesoscopictissue led to genuine taxis, unlike previous models [23, 25–27, 41, 42] where the tactic cue (e.g.,macroscopic tissue) was only space-dependent. To our knowledge dynamics of such complexity hasnot been treated in previous works, neither numerically nor from an analytical viewpoint. In [71]were performed interesting simulations on the kinetic level for a simpler model of cells migratingthrough the extracellular matrix with oriented fibers, which they degrade according to certainrules. Our closed macro-meso description (3.35) has the advantage of a lower dimensional phasespace, thus allowing -among others- to accommodate even further possibly relevant biologicaleffects with a rather detailed description via (ii) and/or (iii), without facing problems due to aprohibitive increase of the number of kinetic variables and therewith associated computationalcosts.Concerning possible relationships between (i)-(iv), it has been proved in [73] that (ii) actuallyimplied (iii) under certain assumptions on the subcellular dynamics for a model of bacterialrun-and-tumble swimming. In [57] such connection was investigated in a less rigorous way for amodel describing glioma pseudopalisade patterning under the influence of repellent pH-taxis andanisotropic, but non-evolving tissue.In order to pass from the KTAP framework to effective equations for the tumor spread weperformed here a hyperbolic as well as a parabolic upscaling, according to the underlying tissuebeing directed or not. The latter is still not clarified from a biological viewpoint; the study(via numerical simulations) performed in [57] for the mentioned model of glioma pseudopalisadesrevealed that such patterns which are invariably observed on histological samples of high gradetumors [14, 89] formed in due time when using the system obtained in the parabolic limit, whilethe hyperbolic scaling led to a drift-dominated system predicting only a rapid shift of the initialaggregate of glioma cells towards the leading tissue orientation, without pattern development.This suggests that brain tissue might be undirected and induced us here to perform numericalsimulations only for the equations obtained by parabolic scaling.Our model makes use of DTI data and accounts for the evolution of brain tissue, transferringbiomedical information through the directional distribution of tissue fibers. As mentioned before,this information belongs to both micro and macro worlds and they coexist in the macroscopiclimit. In Subsection 3.1 we put in evidence the versatility of this model w.r.t. the influence of(spatially varying) tissue anisotropy in characterizing the tumor spread, thereby also consideringthe impact of hypoxia. In fact, our study of dominant advective effects suggests that pH-taxisis the dominating drift - at least as far as this model and its parameterization are concerned.In this work we used the approaches (ii) and (iii) to describe cell-tissue interactions and (iii)for the response of cells to acidity. This led to the anisotropy-triggered switch between myopicdiffusion with ’meso-haptotaxis’ and the migratory behavior with the supplementary ’macro-haptotaxis’. Instead, using (ii) for cell-acid interplay and correspondingly define the turning ratewould conduct to a hypoxia-triggered switch. Several other options of combining (ii) and/or (iii)can be conceived as well, raising further questions about the relationship of these alternativeapproaches. How appropriate a particular model is for a specific problem eventually depends ona multitude of factors and can properly be assessed when validated with adequate patient data.The current high costs of collecting such data still keeps a reliable validation out of reach, butthe unprecedented development of biomedical technology during the last decade gives hope of itbecoming available in the near future. Acknowledgement
GC, CE, AK, CS, and MW acknowledge funding by the Federal Ministry of Education andResearch BMBF in the project
GlioMaTh eferences [1] https://bionumbers.hms.harvard.edu/bionumber.aspx?s=n&v=0&id=108941.[2] P. Basser, Diffusion and diffusion tensor MR imaging:fundamentals. Magnetic ResonanceImaging of the Brain and Spine (S. Atlas ed.). Philadelphia: Lippincott Williams, 2008,pp. 1752–1767.[3] P. Bastian, M. Blatt, A. Dedner et al, A generic grid interface for parallel and adaptivescientific computing. Part I: abstract framework. Computing 82 (2008) 103–119.[4] P. Bastian, M. Blatt, A. Dedner et al, A generic grid interface for parallel and adaptivescientific computing. Part II: implementation and tests in DUNE. Computing 82 (2008)121–138.[5] P. Bastian, F. Heimann, S. Marnach, Generic implementation of finite element methodsin the Distributed and Unified Numerics Environment (DUNE). Kybernetika 46 (2010),294–315.[6] N. Bellomo. Modeling Complex Living Systems. Birk¨auser, Boston, 2008.[7] N. Bellomo, A. Bellouquid, J. Nieto, J. Soler, On the asymptotic theory from microscopicto macroscopic growing tissue models: an overview with perspectives. Math. Models Meth-ods Appl. Sci. 22 (2012), 1130001, 37 pages.[8] A. Bellouquid, J. Calvo, J. Nieto, J. Soler, Hyperbolic vs parabolic asymptotics in kinetictheory towards fluid dynamic models, SIAM J. Appl. Math. 73(4) (2013) 1327–1346.[9] S. Benedetto, R. Pulito, S. Geninatti Crich, G. Tarone, S. Aime, L. Silengo, J. Hamm,Quantification of the Expression Level of Integrin Receptor α v β in Cell Lines and MRImaging With Antibody-Coated Iron Oxide Particles. Magnetic Resonance in Medicine56 (2006) 711–716.[10] T. Beppu, T. Inoue, Y. Shibata, A. Kurose, H. Arai, K. Ogasawara, A. Ogawa, S. Naka-mura, H. Kabasawa, Measurement of fractional anisotropy using diffuson tensor MRI insupratentorial astrocytic tumors. J Neuro Oncol 63 (2003) 109–116.[11] M. Blatt, P. Bastian, The Iterative Solver Template Library. Proc. of the 8th intern. conf.on Applied Parallel Computing: State of the Art in Scientific Computing, 2007, 666–675.[12] M. Blatt, A. Burchardt, A. Dedner, Ch. Engwer, J. Fahlke, B. Flemisch, Ch. Gersbacher,C. Gr¨aser, F. Gruber, Ch. Gr¨uninger, D. Kempf, R. Kl¨ofkorn, T. Malkmus, S. M¨uthing,M. Nolte, M. Piatkowski, O. Sander, The Distributed and Unified Numerics Environment,Version 2.4. Archive of Numerical Software, 2016.[13] P. Bondiau, O. Clatz, M. Sermesant, P. Marcy, H. Delignette, M. Frenay, N. Ayache, Bio-computing: numerical simulation of glioblastoma growth using diffusion tensor imaging.Phys. Med. Biol., 53 (2008) 879–893.[14] D. J. Brat, A. A. Castellano-Sanchez, S. B. Hunter, M. Pecot, C. Cohen, E. H Hammond,S. N Devi, B. Kaur, E. G. Van Meir, Pseudopalisades in glioblastoma are hypoxic, expressextracellular matrix proteases, and are formed by an actively migrating cell population,Cancer Research 64 (2004) 920–927.[15] T.A. Brunner, J.P. Holloway, Two-dimensional time dependent Riemann solvers for neu-tron transport. J. comp. Phys. 210 (2005) 386-399.[16] L. Calorini, S. Peppicelli, and F. Bianchini, Extracellular acidity as favoring factor oftumor progression and metastatic dissemination. Exp. Oncol. 34 (2012) 79–84.[17] Camino: Open-Source Diffusion-MRI Reconstruction and Processing, 14th ScientificMeeting of the International Society for Magnetic Resonance in Medicine, 2006.[18] J. Chiche, M.C. Brahimi-Horn, J. Pouyss´egur, Tumour hypoxia induces a metabolic shiftcausing acidosis: A common feature in cancer. J. Cell. Mol. Med. 14 (2010) 771–794.[19] A. Chauvi`ere, T. Hillen, L. Preziosi, Modeling cell movement in anisotropic and hetero-geneous network tissues. Networks & Heterogeneous Media 2 (2007) 333–357.[20] L. Chen, K. Painter, C. Surulescu, A. Zhigun, Mathematical models for cell migration: anon-local perspective. Phil. Trans. R. Soc. B (in press) arXiv:1911.05200v1.
21] O. Clatz, M. Sermesant, P. Bondiau, H. Delignette, S. Warfield, G. Malandain, N. Ayache,Realistic simulation of the 3D growth of brain tumors in MRI images coupling diffusionwith biomechanical deformation. IEEE Trans. Med. Imaging, 24 (2005), 1334–1346.[22] M. Conte, L. Gerardo-Giorda, M. Groppi, Glioma invasion and its interplay with thenervous tissue: a multiscale model. J. Theor. Biol. 486 (2020): 110088[23] G. Corbin, A. Hunt, A. Klar, F. Schneider, C. Surulescu, Higher-order models for gliomainvasion: from a two-scale description to effective equations for mass density and momen-tum, Math. Mod. Meth. Appl. Sci. 28 (2018) 1771-1800.[24] M. Descoteaux, High angular resolution diffusion MRI: from local estimation to segmen-tation and tractography. Ph.D. Thesis, Universit´e Nice-Sophia Antipolis, 2008.[25] C. Engwer, T. Hillen, M. Knappitsch, C. Surulescu, Glioma Follow White Matter Tracts;a Multiscale DTI-based Model, J. of Math. Biol. 71 (2015) 551–582.[26] C. Engwer, A. Hunt, C. Surulescu, Effective equations for anisotropic glioma spread withproliferation: a multiscale approach, IMA J. Mathematical Medicine & Biology 33 (2016)435-459.[27] C. Engwer, M. Knappitsch, C. Surulescu, A multiscale model for glioma spread includingcell-tissue interactions and proliferation, Math. Biosc. Eng. 13 (2016) 443–460.[28] R. Erban, H.G. Othmer, From Individual to Collective Behavior in Bacterial Chemotaxis,SIAM J. Appl. Math. 65 (2004) 361–391.[29] V. Estrella, T. Chen, M. Lloyd, J. Wojtkowiak, H. H. Cornnell, A. Ibrahim-Hashim, K.Bailey, Y. Balagurunathan, J. M. Rothberg, B. F. Sloane, J. Johnson, R. A. Gatenby, R.J. Gillies, Acidity generated by the tumor microenvironment drives local invasion. CancerRes, 73 (2013) 1524–1535.[30] P.T. Fletcher, S. Joshi, Riemannian geometry for the statistical analysis of diffusion tensordata. Signal Processing 87 (2007) 250–262.[31] H.B. Frieboes, J.S. Lowengrub, S. Wise, X. Zheng, P. Macklin, E. Bearer, V. Cristini,Computer Simulation of Glioma Growth and Morphology, Neuroimage 37 (2007), S59S70.[32] A. Giese, L. Kluwe, H. Meissner, M.E. Berens, M. Westphal, Migration of human gliomacells on myelin, Neurosurgery 38 (1996), 755–764.[33] A. Giese, M. Westphal, Glioma invasion in the central nervous system, Neurosurgery, 39(1996) 235–252.[34] L. Gonz´alez-M´endez, A. C. Gradilla, I. Guerrero, The cytoneme connection: direct long-distance signal transfer during development. Development 146(9), dev174607 (2019).[35] M.A. Hammoud, R. Sawaya, W. Shi, P.F. Thall, N.E. Leeds: Prognostic significance ofpreoperative MRI scans in glioblastoma multiforme, J. Neurooncol. 27 (1996) 65–73.[36] H. Hatzikirou, A. Deutsch: Cellular automata as microscopic models of cell migration inheterogeneous environments, Curr. Top. Dev. Biol. 81 (2008) 401–434.[37] C. Henker, T. Kriesen, ¨A. Glass, B. Schneider, J. Piek, Volumetric quantification ofglioblastoma: experiences with different measurement techniques and impact on survival.J Neurooncol 135 (2017) 391–402.[38] T. Hillen, M mesoscopic and macroscopic models for mesenchymal motion, J. Math.Biol. 53 (2006) 585–616.[39] T. Hillen, K. Painter, Transport and anisotropic diffusion models for movement in orientedhabitats, Dispersal, Individual Movement and Spatial Ecology (M. Lewis, P. Maini & S.Petrovskii eds). Lecture Notes in Mathematics, vol. 2071. Berlin, Heidelberg, New York:Springer 2013, p. 46.[40] T. Hillen, A. Swan, The Diffusion Limit of Transport Equations in Biology. In: L. Preziosi,P. Ciarletta, T. Hillen, M. Chaplain, A. Pugliese, H. Othmer, D. Trucu (eds.), Mathemat-ical Models and Methods for Living Systems, LNM 2167, Springer, 2016.[41] A. Hunt, DTI-Based Multiscale Models for Glioma Invasion, PhD thesis, TU Kaiser-slautern, 2017.[42] A. Hunt, C. Surulescu, A multiscale modeling approach to glioma invasion with therapy,Vietnam J. Math. 45 (2017) 221–240.
43] A. Jbabdi, E. Mandonnet, H. Duffau, L. Capelle, K. Swanson, M. Pelegrini-Issac, R.Guillevin, H. Benali, Simulation of anisotropic growth of low-grade gliomas using diffusiontensor imaging. Mang. Res. Med., 54 (2005) 616–624.[44] B. Jellison, A. Field, J. Medow, M. Lazar, M. Salamat, A. Alexander, Diffusion tensorimaging of cerebral white matter: a pictorial review of physics, fiber tract anatomy, andtumor imaging patterns. Am. J. Neuroradiol., 25 (2004) 356–369.[45] A.R. Kansal, S. Torquato, G.R. Harsh IV, E.A. Chiocca, T.S. Deisbrock, Cellular au-tomaton of idealized brain tumor growth dynamics, BioSystems 55 (2000) 119–127.[46] A.R. Kansal, S. Torquato, G.R. Harsh IV, E.A. Chiocca, T.S. Deisbrock, Simulated BrainTumor Growth Dynamics Using a Three-Dimensional Cellular Automaton, J. Theor. Biol.203 (2000) 367–382.[47] Y. Kato, S. Ozawa, C. Miyamoto, Y. Maehata, A. Suzuki, T. Maeda, Y. Baba, Acidicextracellular microenvironment and cancer. Cancer Cell Int. 13 (2013) 8 pp.[48] J. Kelkel, C. Surulescu, On some models for cancer cell migration through tissue, Math.Biosci. Engrg. 8 (2011), 575–589.[49] J. Kelkel, C. Surulescu, A multiscale approach to cell migration in tissue networks, Math.Mod. Meth. Appl. Sci. 22 (2013), 1150017 (25 p.).[50] P.T. Kim, D. Richards, Deconvolution density estimation on the space of positive definitesymmetric matrices. In Nonparametric Statistics and Mixture Models: A Festschrift inHonor of Thomas P Hettmansperger (pp. 147-168). World Scientific Publishing Co., 2011.[51] P. Kleihues, F. Soylemezoglu, B. Sch¨auble, B.W. Scheithauer, P.C. Burger, Histopathol-ogy, classification and grading of gliomas. Glia 5 (1995) 211–221.[52] N. Kolbe, N. Sfakianakis, C. Stinner, C. Surulescu, J. Lenz, Modeling multiple taxis:tumor invasion with phenotypic heterogeneity, haptotaxis, and unilateral interspecies re-pellence, arXiv:2005.01444v1.[53] E. Konukoglu, O. Clatz, P. Bondiau, H. Delignette, N. Ayache, Extrapolation gliomainvasion margin in brain magnetic resonance images: suggesting new irradiation margins.Med. Image Anal., 14 (2010) 111–125.[54] T.B. Kornberg, L. Gilboa, Developmental biology: Nanotubes in the niche. Nature 523(2015) 292–293.[55] J.M. Kros, Grading of Gliomas: The Road From Eminence to Evidence, J. Neuropathol.Exp. Neurol. 70 (2011) 101–109.[56] L.Y. Kucheryavykh, Y.V. Kucheryavykh, K. Rolon-Reyes, S.N. Skatchkov, M.J. Eaton,L.A. Cubano, M. Inyushin, Visualization of implanted GL261 glioma cells in livingmouse brain slices using fluorescent 4-(4-(dimethylamino)-styryl)-N-methylpyridinium io-dide (ASP+). BioTechniques, 53 (2012) 305–309.[57] P. Kumar, C. Surulescu, Modeling glioma pseudopalisade patterning: a multiscale ap-proach, in preparation.[58] D.A. Lauffenburger, J.L. Lindermann, Receptors. Models for binding, trafficking and sig-naling. Oxford University Press, 1993.[59] M. Lemou and L. Mieussens, A new asymptotic preserving scheme based on micro-macroformulation for linear kinetic equations in the diffusion limit. SIAM J. Sci. Comput. 31(2008) 334–368.[60] D.R. Lide (ed). CRC handbook of chemistry and physics, vol. 85. CRC Press, 2004.[61] T. Lorenz, C. Surulescu, On a class of multiscale cancer cell migration models: well-posedness in less regular function spaces, Math. Mod. Meth. Appl. Sci. 24 (2014) 2383-2436.[62] D. Louis: Molecular Pathology of Malignant Gliomas, Annu. Rev. Pathol. 1 (2006) 97–117.[63] N. Loy, L. Preziosi, Kinetic models with non-local sensing determining cell polarizationand speed according to independent cues, J. Math. Biol. 80 (2020) 373–421.[64] N.K. Martin, E.A. Gaffney, R.A. Gatenby, R.J. Gillies, I.F. Robey, and P.K. Maini. Amathematical model of tumour and blood pHe regulation: The buffering system. Mathe-matical Biosciences, 230 (2011) 1–11.
65] G.R. Martin and R.K. Jain. Noninvasive measurement of interstitial ph profiles in normaland neoplastic tissue using fluorescence ratio imaging microscopy. Cancer Research, 54(1994) 5670–5674.[66] J. Mercapide, R. Lopez De Cicco, J.S. Castresana, and A.J.P. Klein-Szanto. Stromelysin-1/matrix metalloproteinase-3 (MMP-3) expression accounts for invasive properties of hu-man astrocytoma cell lines. International Journal of Cancer, 106 (2003) 676–682.[67] R. Milo and R. Phillips, Cell biology by the numbers, Garland Science, 2015.[68] J. Nieto and L. Urrutia: A multiscale model of cell mobility: From a kinetic to a hydro-dynamic description, J. Math. Anal. Appl. 433 (2016) 1055–1071.[69] T. Ohtsubo, X. Wang, A. Takahashi, K. Ohnishi, H. Saito, et al. p53-dependent inductionof WAF1 by a low-pH culture condition in human glioblastoma cells. Cancer Res. 57(1997) 3910–3913.[70] H.G. Othmer, T. Hillen, The diffusion limit of transport equations II: Chemotaxis equa-tions. SIAM Journal on Applied Mathematics, 62 (2003) 1220–1250.[71] K. Painter, Modelling cell migration strategies in the extracellularmatrix. J. Math. Biol.58 (2009) 511–543.[72] K. Painter, T. Hillen, Mathematical modelling of glioma growth: the use of diffusion tensorimaging (DTI) data to predict the anisotropic pathways of cancer invasion, J. Theor. Biol.323 (2013) 25–39.[73] B. Perthame, M. Tang, N. Vauchelet, Derivation of the bacterial run-and-tumble kineticequation from a model with biochemical pathway, J. Math. Biol. 73 (2016), pp.1161-1178.[74] R.G. Plaza, Derivation of a bacterial nutrient-taxis system with doubly degenerate cross-diffusion as the parabolic limit of a velocity-jump process. J. Math. Biol. 78 (2019) 1681–1711.[75] G.C. Pomraning, The Equations of Radiation Hydrodynamics. Dover books on physics(2005)[76] M. Portela, et al., Glioblastoma cells vampirize WNT from neurons and trigger aJNK/MMP signaling loop that enhances glioblastoma progression and neurodegenera-tion. PLoS Biol. 17(12) (2019).[77] D. Poyato, J. Soler, Euler-type equations and commutators in singular and hyperboliclimits of kinetic Cucker-Smale models, Math. Mod. Meth. Appl. Sci, 27(6) (2017), pp.1089–1152.[78] S.M. Raza, G.N. Fuller, C.H. Rhee, S. Huang, K. Hess, W. Zhang, R. Sawaya, Identifica-tion of Necrosis-Associated Genes in Glioblastoma by cDNA Microarray Analysis, Clin.Canc. Res. 10 (2004) 212–221.[79] C. Stock, A. Schwab, Protons make tumor cells move like clockwork. Pflugers Arch. 458(2009) 981–992.[80] P.C. Sundgren, Q. Dong, D. Gomez-Hassan, S.K. Mukherji, P. Maly, R. Welsh, Diffusiontensor imaging of the brain: Review of clinical applications, Neurocarciology, 46 (2004)339–350.[81] K.R. Swanson, R.C. Rockne, J. Claridge, M.A. Chaplain, E.C. Alvord, and A.R.A. An-derson. Quantifying the role of angiogenesis in malignant progression of gliomas: In silicomodeling integrates imaging and histology. Cancer Research, 71 (2011) 7366–7375.[82] M.L. Tanaka, W. Debinski, I.K. Puri, Hybrid mathematical model of glioma progression,Cell Prolif. 42 (2009) 637–646.[83] van der Vorst, Bi-cgstab: A fast and smoothly converging variant of bi-cg for the solutionof nonsymmetric linear systems. SIAM Journal on Scientific and Statistical Computing13 (1992) 631–644.[84] P. Vaupel, F. Kallinowski, P. Okunieff: Blood flow, oxygen and nutrient supply, andmetabolic microenvironment of human tumors: A review. Cancer Res. 49 (1989) 6449–6465.[85] C. Wang, J. Rockhill, M. Mrugala, D. Peacock, A. Lai, K. Jusenius, J. Wardlaw, T.Cloughesy, A. Spence, R. Rockne, E. Alvord Jr., K. Swanson, Prognostic significance ofgrowth kinetics in newly diagnosed glioblastomas revealed by combining serial imagingwith a novel biomathematical model. Cancer Res., 69 (2009) 9133–9140.
86] B.A. Webb, M. Chimenti, M.P. Jacobson, D.L. Barber, Dysregulated pH: A perfect stormfor cancer progression. Nature Rev. Cancer 11 (2011) 671–677.[87] M. Winkler, Singular structure formation in a degenerate haptotaxis model involvingmyopic diffusion. J. Math. Pures Appl. 112 (2018) 118–169.[88] M. Winkler, C. Surulescu, Global weak solutions to a strongly degenerate haptotaxismodel. Commun. Math. Sci. 15 (2017) 1581–1616.[89] F. J. Wippold, M. L¨ammle, F. Anatelli, J. Lennerz, A. Perry, Neuropathology for theneuroradiologist: palisades and pseudopalisades. American Journal of Neuroradiology 27(2006) 2037–2041.[90] R.K.W. Wong, T.C.M. Lee, D. Paul, J. Peng, Fiber direction estimation, smoothing, andtracking in diffusion MRI. Ann. appl. Stat. 10 (2016) 1137–1156.[91] M. Wrensch, Y. Minn, T. Chew, M. Bondy, M. S. Berger, Epidemiology of primary braintumors: Current concepts and review of the literature, Neuro-Oncology, 4 (2002) 278–299.[92] C. Xue, H.G. Othmer, Multiscale models of taxis-driven patterning in bacterial popula-tions, SIAM J. Appl. Math. 70 (2009) 133-169.[93] L. Zhang, B. Jiang, Y. Wu, C. Strouthos, P. Zhe Sun, J. Su, X. Zhou, Developing amultiscale, multi-resolution agent- based brain tumor model by graphics processing units,Theor. Biol. Medical Mod. 8:46 (2011).86] B.A. Webb, M. Chimenti, M.P. Jacobson, D.L. Barber, Dysregulated pH: A perfect stormfor cancer progression. Nature Rev. Cancer 11 (2011) 671–677.[87] M. Winkler, Singular structure formation in a degenerate haptotaxis model involvingmyopic diffusion. J. Math. Pures Appl. 112 (2018) 118–169.[88] M. Winkler, C. Surulescu, Global weak solutions to a strongly degenerate haptotaxismodel. Commun. Math. Sci. 15 (2017) 1581–1616.[89] F. J. Wippold, M. L¨ammle, F. Anatelli, J. Lennerz, A. Perry, Neuropathology for theneuroradiologist: palisades and pseudopalisades. American Journal of Neuroradiology 27(2006) 2037–2041.[90] R.K.W. Wong, T.C.M. Lee, D. Paul, J. Peng, Fiber direction estimation, smoothing, andtracking in diffusion MRI. Ann. appl. Stat. 10 (2016) 1137–1156.[91] M. Wrensch, Y. Minn, T. Chew, M. Bondy, M. S. Berger, Epidemiology of primary braintumors: Current concepts and review of the literature, Neuro-Oncology, 4 (2002) 278–299.[92] C. Xue, H.G. Othmer, Multiscale models of taxis-driven patterning in bacterial popula-tions, SIAM J. Appl. Math. 70 (2009) 133-169.[93] L. Zhang, B. Jiang, Y. Wu, C. Strouthos, P. Zhe Sun, J. Su, X. Zhou, Developing amultiscale, multi-resolution agent- based brain tumor model by graphics processing units,Theor. Biol. Medical Mod. 8:46 (2011).