Mathematical modeling of glioma invasion: acid- and vasculature mediated go-or-grow dichotomy and the influence of tissue anisotropy
MMathematical modeling of glioma invasion: acid- and vasculaturemediated go-or-grow dichotomy and the influence of tissue anisotropy
Martina Conte ∗ and Christina Surulescu BCAM - Basque Center for Applied MathematicsAlameda de Mazarredo, 14 - E-48009 Bilbao, Spain Technische Universit¨at Kaiserslautern, Felix-Klein-Zentrum f¨ur MathematikPaul-Ehrlich-Str. 31 - 67663 Kaiserslautern, Germany
July 27, 2020
Abstract
Starting from kinetic transport equations and subcellular dynamics we deduce a multiscale model forglioma invasion relying on the go-or-grow dichotomy and the influence of vasculature, acidity, and braintissue anisotropy. Numerical simulations are performed for this model with multiple taxis, in order toassess the solution behavior under several scenarios of taxis and growth for tumor and endothelial cells. Anextension of the model to incorporate the macroscopic evolution of normal tissue and necrotic matter allowsus to perform tumor grading.
Keyworks –
Multiscale modeling of glioma invasion, go-or-grow dichotomy and hypoxia-driven phenotypicswitch, multiple taxis, necrosis-based tumor grading.
Glioma is the most frequent type of primary brain tumor. It originates from glia cells and accounts for 78percent of malignant brain tumors, of which glioblastoma multiforme (GBM) is the most aggressive, beingcharacterized by a fast, infiltrative spread and high proliferation. These features make it very difficult to treat,and go along with a poor survival prognosis [34, 93].Like in most tumors, glioma development, growth, and invasion are influenced by a multitude of intrinsicand extrinsic factors. Among these, the phenotypic heterogeneity and the biological, physical, and chemicalcomposition of the tumor microenvironment play a decisive role. Experiments with cultures of glioma cellssuggest mutual exclusion of migratory and proliferative behavior, as reviewed e.g. in [7, 26, 94]; this is knownas go-or-grow dichotomy [27, 28]. Biological evidence indicates that migratory and proliferative processes sharecommon signaling pathways, suggesting a unique intracellular mechanism that regulates both phenotypes [26].Hypoxia is a prominent trait of tumor microenvironment and glioma neoplasms are no exception. It has beensuggested (see e.g. [44,94] and references therein) that it is putatively influencing the phenotypic switch betweenmigrating and proliferating behavior - along with other regulating factors, like angiogenesis, ECM productionand degradation, etc. Indeed, glioma cells have been observed to move away from highly hypoxic sites created,for instance, by the occlusion of a capillary and to form so-called pseudopalisade patterns, which are typical forglioblastoma [10,11,92]. They have garland-like shapes exhibiting central necrotic zones surrounded by stacks oftumor cells, most of which are actively migrating. As (tumor) cell proliferation is impaired at (too) low pH, thisseems to endorse the antagonistic relationship between (transiently) migratory and proliferative phenotypes [68].While cancer cells are able to survive in relatively acidic regions by anaerobic glycolysis, which confers them anadvantage against normal cells [35, 90], the large amounts of lactate and alanine they produce during this pro-cess can decrease the pH below critical levels. As a consequence, they initiate (re)vascularization by expressingpro-angiogenic factors, in order to provide adequate supply with blood-transported nutrients [32, 89]. They candeter proliferation for migration towards more favorable areas [65,90], and the above mentioned pseudopalisadeformation is just one aspect of this complex behavior.Previous models for glioma invasion have been proposed in (semi)discrete [9, 33, 53] or continuous frameworks.Most of the latter are purely macroscopic, describing the evolution of glioma cell density under the influenceof surrounding tissue, chemical signals, and/or vasculature, see e.g. [37, 46, 56, 64, 86] and the review [3]; many ∗ Corresponding author: [email protected] a r X i v : . [ q - b i o . T O ] J u l f them are versions or extensions of a model proposed by Murray [69]. More recent continuous models leavebehind the classical reaction-diffusion prototype and take into account advection bias of glioma cells in re-sponse to environmental cues. Some of these are directly set on the macroscopic scale and rely on balance ofmass/flux/momentum [13, 42, 52], others are obtained by more detailed descriptions from lower scale dynam-ics, as proposed in [71] and further developed in [14, 16, 19–21, 45, 57, 85]. In particular, some of the latterhave a multiscale character that is (partially) preserved during the upscaling process from kinetic transportequations (mesoscale) to reaction-diffusion-taxis PDEs (macroscale). As mentioned above, hypoxia is an es-sential factor in tumor evolution. Early mathematical models for cancer invasion and patterning under acidicconditions were introduced in [23, 66, 80]. Further PDE-based models, that characterize acid-mediated tumordevelopment, directly related or applicable to pH-influenced glioma spread have been proposed and investigatede.g. in [64] and [22, 67, 81], respectively. Instead, for stochastic multiscale settings see e.g. [38, 39, 41] and [40]for a review also addressing further related models. Tumor heterogeneity has been usually modeled in thecontinuous setting by describing the dynamics of the correspondingly defined subpopulations of cells, e.g. hy-poxic/normoxic/necrotic [37, 64, 86] or moving/proliferating [21, 25, 45, 67, 74, 82, 96]. An indirect accounting fora go-grow-recede heterogeneity under the influence of intra- and extracellular acidity was proposed and analyzedin [39]; its numerical simulations were able to explain a large variety of invasion patterns.The models in [21, 45] started from a mesoscopic description of glioma density functions for migrating, respec-tively proliferating cells in interaction with the anisotropic brain tissue, whereby the latter took into accountsubcellular dynamics of receptor binding to tissue fibers on the microscale. The cells were able to switch betweenthe two phenotypes, the corresponding rates being dependent on space and (in [45]) also on the doses of somechemotherapeutic agent administered in order to impair cell motility. In this note we extend those settings byincluding the dynamics of endothelial cells developing vasculature and of acidity, both being responsible for thephenotypic switch. The upscaling of the kinetic transport PDEs for glioma and endothelial cells leads in theparabolic limit to a system of reaction-advection-diffusion equations featuring several types of taxis, nonlinearmyopic diffusion of glioma, and highly complicated couplings between the variables of the model, which has amultiscale character due to (some of) the taxis coefficients encoding information from the lower modeling levels.The subsequent content is organized as follows. Section 2 is concerned with the setup of the model on thesubcellular and mesoscopic scales and with the deduction of the macroscopic PDEs. In Section 3 we concretizethe coefficient functions in preparation for the numerical simulations to be performed in Section 4. In order tofacilitate the evaluation of the tumor burden in relation to the necrotic and the normal tissue, Section 5 includesa model extension accounting for the dynamics of the latter. We conclude with a discussion of the results inSection 6. The Appendix contains the assessment of the model parameters and a nondimensionalization of themacroscopic PDEs. Relying on the go-or-grow dichotomy, we consider a tumor containing two mutually exclusive subpopulationsof glioma cells, which are either migrating or proliferating. The respective states are transient, the tumor canchange dynamically its composition, according to the signals received by the cells from their surroundings. Fromthe huge variety of chemical and physical cues present in the extracellular space and influencing the developmentand spread of cancer [31] we focus here on the effects of acidity and cell-tissue interactions. Since acidificationand angiogenesis are tightly interrelated and crucial for the tumor evolution, we also model vascularization, byway of endothelial cell dynamics. We develop a multiscale model upon starting from the subcellular level ofinteractions between cells, acidity, and tissue, setting up the corresponding kinetic transport equations (KTEs)for glioma cells of the two phenotypes and for endothelial cells, and performing a parabolic limit to deducethe macroscopic system of reaction-advection-diffusion PDEs for the involved quantities: total tumor burden(moving + proliferating cells), endothelial cells, and acidity (concentration of protons).
On the microscopic scale, we describe the interaction of glioma cells with the extracellular space, more preciselywith tissue fibers and protons. Cells exchange information with their environment through various transmem-brane entities, e.g. cell surface receptors and ion channels. We will use the former to account for cell-tissueinteractions and both for the cell-proton exchange. Indeed, additionally to ion channels and membrane trans-porters which have been extensively studied in the context of intra- and extracellular pH regulation, there alsoexist proton-sensing receptors [43], e.g. the G protein-coupled receptors (GPCRs) involved among others inregulating the migration and proliferation of cells in tumor development and wound healing [47, 91]. We ignorehere all intricate details about the intracellular machinery activated by receptor binding and channel openingand closing. Instead, we see the events of occupying such transmembrane units as triggering the cellular pro-cesses leading to migration, proliferation, and phenotypic switch.2onsidering such interactions between cells and soluble as well as unsoluble ligands follows the idea employedin [50, 51, 60] to build a micro-meso model for tumor invasion with chemo- and haptotaxis and revisitedin [14, 16, 19–21, 45] for cell-tissue interactions, in [57] for cell-proton interactions leading to pseudopalisadepatterns, or in [22] for a nonlocal micro-macro model with cell-tissue, cell-cell, and cell-protons interactions.The micro-macro framework proposed in [38, 39, 41] offers a related, yet different perspective on the interplaybetween intra- and extracellular acidity and tumor cells.We denote by y ( t ) the amount of receptors bound to tissue fibers and by y ( t ) that of transmembrane entities(ion channels, pH-sensing receptors) occupied by protons. The corresponding binding/occupying dynamics ischaracterized by simple mass action kinetics:¯ R − ( y + y ) + QQ ∗ k +1 (cid:10) k − y ¯ R − ( y + y ) + SS c, k +2 (cid:10) k − y . Here ¯ R is the total amount of receptors on a cell membrane, Q ( x ) denotes the macroscopic tissue density,depending on the position x ∈ R n , Q ∗ ( x ) is the reference tissue density (corresponding to its carrying capac-ity), while S ( t, x ) is the macroscopic concentration of protons, and S c, ( x ) the reference proton concentration.Accordingly, we get the ODE system˙ y = k +1 Q ∗ Q ( x ) ( ¯ R − ( y + y )) − k − y ˙ y = k +2 S c, S ( t, x ) ( ¯ R − ( y + y )) − k − y where k +1 and k − represent the attachment and detachment rates of cells to tissue, respectively, while k +2 and k − are the corresponding rates in the process of proton binding. As in [22], we define y := y + y to be thetotal amount of transmembrane entities occupied by tissue or protons, which allows us to lump together thetwo ODEs, into ˙ y = (cid:18) k +1 Q ∗ Q + k +2 S c, S (cid:19) ¯ R − y (cid:18) k +1 Q ∗ Q + k +2 S c, S (cid:19) − k − y − k − y . Assuming that k − = k − = k − , we get the microscopic equation for the subcellular dynamics˙ y = (cid:18) k +1 Q ∗ Q + k +2 S c, S (cid:19) ( ¯ R − y ) − k − y . (1)Since processes on this scale are much faster than those happening on the macroscopic time scale, they canbe assumed to equilibrate rapidly and moreover, that on this scale we can ignore the time dependence of S .Rescaling y/ ¯ R (cid:32) y will further simplify the notation. Then, the unique steady state of the above equation isgiven by: y ∗ = k +1 Q ∗ Q + k +2 S c, S k +1 Q ∗ Q + k +2 S c, S + k − =: ¯ f ( Q, S ) . The variable y can be seen as characterizing the ’internal’ cellular state. In the kinetic theory of active particles(see e.g. [5, 6]) it is called ’activity variable’. In the sequel, we will consider the mesoscopic densities p ( t, x, v, y )and r ( t, x, y ) of migrating and respectively non-moving (thus proliferating) glioma cells, hence both dependingon such activity variable y . Thereby, v ∈ V ⊂ R n is the cell velocity vector, with the space V to be closerexplained in Subsection 2.2 below.Further, we assume that the glioma cells follow the tissue gradient, but move away from highly acidic areas.Therefore, we look at the path of a single cell starting at position x and moving to position x with velocity v in the (locally) time-invariant density fields Q and S , so that Q ( x ) = Q ( x + vt ), while S ( x ) = S ( x − vt ).Denoting by z := y ∗ − y the deviation of y from its steady state, we have:˙ z = ∂ ¯ f∂Q v · ∇ x Q − ∂ ¯ f∂S v · ∇ x S − z (cid:18) k +1 Q ∗ Q + k +2 S c, S + k − (cid:19) , ∂ ¯ f∂Q = 1 Q ∗ k +1 k − (cid:16) k +1 Q ∗ Q + k +2 S c, S + k − (cid:17) ∂ ¯ f∂S = 1 S c, k +2 k − (cid:16) k +1 Q ∗ Q + k +2 S c, S + k − (cid:17) , thus the equation for z is given by:˙ z = − z (cid:18) k +1 Q ∗ Q + k +2 S c, S + k − (cid:19) + k − (cid:16) k +1 Q ∗ Q + k +2 S c, S + k − (cid:17) (cid:18) k +1 Q ∗ v · ∇ x Q − k +2 S c, v · ∇ x S (cid:19) =: G ( z, Q, S ) . (2)To simplify the writing, we denote B ( Q, S ) := (cid:18) k +1 Q ∗ Q + k +2 S c, S + k − (cid:19) . We model the mesoscale behavior of glioma and endothelial cells with the aid of kinetic transport equations(KTEs) describing velocity-jump processes and taking into account the subcellular dynamics. Concretely, weconsider the following cell density functions: • p ( t, x, v, y ) for moving glioma cells; • r ( t, x, y ) for non-moving, hence (in virtue of the go-or-grow dichotomy) proliferating glioma cells; • w ( t, x, ϑ ) for endothelial cells (ECs) forming capillaries,with the time and space variables t > x ∈ R n , velocities v ∈ V = s S n − and ϑ ∈ Θ := σ S n − , and activityvariable y ∈ Y = (0 , s > σ >
0, respectively, where S n − denotes the unit sphere in R n . As in [19–21,45], in the sequel we will work withthe deviation z = y ∗ − y ∈ Z ⊆ ( y ∗ − , y ∗ ) rather then with y .The corresponding macroscopic cell densities are denoted by M ( t, x ), R ( t, x ) and W ( t, x ), respectively, and weuse the notation N ( t, x ) := M ( t, x ) + R ( t, x ) for the space-time varying macroscopic total tumor burden.The kinetic transport equation for the motile glioma phenotype is given by: ∂ t p + ∇ x · ( vp ) + ∂ z ( G ( z, Q, S ) p ) = L p [ λ ( z )] p + β ( S ) qω r − α ( w, S ) p − l m ( N ) p , (3)where L p [ λ ( z )] p denotes the turning operator, describing velocity changes. In particular, such changes are dueto contact guidance: the cells have a tendency to align to the brain tissue anisotropy, mainly associated withwhite matter tracts. As in previous models of this type, L p [ λ ] p is a Boltzmann-like integral operator of theform L p [ λ ( z )] p = − λ ( z ) p + λ ( z ) (cid:90) V K ( x, v ) p ( v (cid:48) ) dv (cid:48) , (4)where λ ( z ) := λ − λ z ≥ z , while λ and λ are positive constants. The integral term describes the reorientation of cells from any previous velocity v (cid:48) to a new velocity v upon interacting with the tissue fibers. This is controlled by the turning kernel K ( x, v ),which we assume to be independent on the incoming velocity v (cid:48) . In particular, we consider that the dominatingdirectional cue is given by the orientation of tissue fibers and take as e.g. in [19, 71] K ( x, v ) := q ( x, ˆ v ) ω , whereˆ v = v | v | and q ( x, θ ) with θ ∈ S n − is the orientational distribution of the fibers, normalized by ω = s n − . Itencodes the individual brain structure obtained by diffusion tensor imaging (DTI); concrete choices will beprovided in Section 3. We assume the tissue to be undirected, hence q ( x, θ ) = q ( x, − θ ) for all x ∈ R n . For laterreference, we introduce the notations E q ( x ) := (cid:90) S n − θq ( x, θ ) dθ V q ( x ) := (cid:90) S n − ( θ − E q ) ⊗ ( θ − E q ) q ( x, θ ) dθ for the mean fiber orientation and the variance-covariance matrix for the orientation distribution of tissue fibers,respectively. Notice that the above symmetry of q implies E q = 0.4he term β ( S ) qω r in (3) describes the phenotypic switch r → p ; the rate β depends on S , since in a too acidicenvironment the cells are supposed to stop proliferation and migrate towards regions with higher pH (see Sec-tion 1). Therefore, we assume that β ( S ) is an increasing function of the proton concentration S and we require β ( S ) >
0, since there will always be some proliferating cells switching into a migratory regime .The terms α ( w, S ) p and l m ( N ) p model the phenotypic switch p → r due to environmental signals. Thus, theformer term depends on the proton concentration S and on the mesoscopic EC density w and describes theadoption of a proliferative phenotype when there are enough oxygen and nutrient supply, while the acidityremains below a certain threshold. The latter term models the switch to proliferation caused by the gliomacell population being too crowded to allow effective migration (but still allowing some limited proliferation). Inorder not to complicate too much the model we do not explicitly account for cell recession or quiescence.The evolution of proliferating tumor cells is characterized by the integro-differential equation ∂ t r = ( α ( w, S ) + l m ( N )) (cid:90) V p ( t, x, v, z ) dv + µ ( W, N, S ) (cid:90) Z χ ( x, z, z (cid:48) ) Q ( x ) Q ∗ r ( t, x, z (cid:48) ) dz (cid:48) − ( β ( S ) + γ ( S )) r , (5)where, additionally to the already described switch terms, we model intrinsic proliferation and death. µ ( W, N, S ) (cid:82) Z χ ( x, z, z (cid:48) ) Q ( x ) Q ∗ ( x ) r ( z (cid:48) ) dz (cid:48) describes as in [20, 45] proliferation triggered by cell receptor binding to tissue. Theproliferation rate µ ( W, N, S ) depends on the total macroscopic tumor density N , on the concentration of protons S , and on the vasculature. In the integral operator, the kernel χ ( x, z, z (cid:48) ) characterizes the transition from state z (cid:48) to state z during such proliferation-initiating interaction at position x . No further conditions are requiredon χ , we only assume that the nonlinear proliferative operator is uniformly bounded in the L -norm, which isreasonable, in view of space-limited cell division. The last term − γ ( S ) r describes acid-induced death of gliomacells when the pH value drops below a certain threshold.The KTE for ECs is given by: ∂ t w + ∇ x · ( ϑw ) = L w [ η ] w + µ W ( W, Q ) w , (6)where the turning operator L w [ η ] w describes changes in the orientation of ECs. It is well-established (seee.g. [31]) that tumor cells produce angiogenic signals acting as chemoattractants for ECs. Less known is whethersuch signals are expressed by proliferating rather than moving cells or by both phenotypes. There is, however,evidence that hypoxia induces production of VEGF and other angiogenic cytokines (see e.g. [12, 78, 95] andreferences therein). Since cancer cells are highly glycolytic (which leads to acidification), and increased glucosemetabolism is selected for in proliferating cells [61], we assume that pro-angiogenic signals are mainly producedby proliferating cells. With the aim of avoiding the introduction of a new variable for the concentration of suchchemoattractants, we let the ECs be attracted by the proliferating glioma cells as main sources therewith. Thistranslates into the following form of the turning operator acting on the right hand side in (6): L w [ η ] w = − η ( x, ϑ, R ) w ( t, x, ϑ ) + (cid:90) Θ | Θ | η ( x, ϑ (cid:48) , R ) w ( ϑ (cid:48) ) dϑ (cid:48) , where for the turning kernel modeling ECs reorientations we simply took a uniform density function over theunit sphere S n − and η ( x, ϑ, R ) = η ( x ) e − a ( R ) D t R (7)represents the turning rate of ECs. It depends on the macroscopic density of proliferating tumor cells R andon the pathwise gradient D t R = ∂ t R + ϑ · ∇ R . The coefficient function a ( R ) is related to the interactionsbetween ECs and proliferating glioma, more precisely ECs and pro-angiogenic signals produced by the latter.This can be described for instance via equilibrium of EC receptor binding. This way to include directional biasprovides an alternative to that using a transport term with respect to the activity variable z in (3) and it hasbeen introduced in [70] in the context of bacteria movement and recently used in [57] for a model for gliomapseudopalisade patterning. Under certain conditions, the relation between the two approaches was establishedrigorously for bacteria chemotaxis in [73] and investigated more formally for glioma repellent pH-taxis in [57].Finally, µ W ( W, Q ) w is the proliferation term for ECs. Besides the total population of ECs irrespective of theirorientation, µ W is supposed to depend on the available macroscopic tissue Q . More details about the concretechoices of the coefficient functions involved in (3), (5), and (6) will be provided in Section 3. The high dimensionality of the KTE system (3), (5), and (6) makes its numerical simulations too expensive.Moreover, clinicians are interested in the macroscopic evolution of the tumor along with its vascularization and otherwise the tumor would stay confined, which is not the case for glioblastoma N = R + M . The PDE for the proton concentration S is already macroscopic anddoes not need to be upscaled. We rescale the time and space variables as t → ε t and x → εx . Moreover,the proliferation terms in (5) and (6) and the death term in (5) will be scaled by ε in order to account forthe mitotic and apoptotic events taking place on a much larger time scale than migration and switching frommoving to non-moving regimes. Hence, (3), (5), and (6) become: ε ∂ t p + ε ∇ x · ( vp ) − ∂ z (cid:32)(cid:16) zB ( Q, S ) − εk − B ( Q, S ) (cid:18) k +1 Q ∗ v · ∇ Q − k +2 S c, v · ∇ S (cid:19) p (cid:33) = L p [ λ ( z )] p + β ( S ) qω r − α ( w, S ) p − l m ( N ) p (8) ε ∂ t r = ( α ( w, S ) + l m ( N )) (cid:90) V p ( v ) dv + ε µ ( W, N, S ) (cid:90) Z χ ( x, z, z (cid:48) ) Q ( x ) Q ∗ r ( t, x, z (cid:48) ) dz (cid:48) − ( β ( S ) + ε γ ( S )) r , (9) ε ∂ t w + ε ∇ x · ( ϑw ) = L w [ η ε ] w + ε µ W ( W, Q ) w , (10)with η ε ( x, ϑ, R ) = η ( x ) exp (cid:16) − a ( R )( ε ∂ t R + εϑ · ∇ R ) (cid:17) . (11)As in [19–21], we define the following moments: m ( t, x, v ) = (cid:90) Z p ( t, x, v, z ) dz, R ( t, x ) = (cid:90) Z r ( t, x, z ) dz, m z ( t, x, v ) = (cid:90) Z zp ( t, x, v, z ) dz,R z ( t, x ) = (cid:90) Z zr ( t, x, z ) dz, M ( t, x ) = (cid:90) V m ( t, x, v ) dv, W ( t, x ) = (cid:90) Θ w ( t, x, ϑ ) dϑ,M z ( t, x ) = (cid:90) V m z ( t, x, v ) dv and neglect the higher order moments w.r.t. the variable z , in virtue of the subcellular dynamics being muchfaster than the events on the higher scales, hence z (cid:28)
1. We assume the functions p and r to be compactlysupported in the phase space R n × V × Z and w to be compactly supported in R n × Θ.We first integrate equation (8) with respect to z , getting the following equation for m ( t, x, v ): ε ∂ t m + ε ∇ x · ( vm ) = − λ (cid:16) m − qω M (cid:17) + λ (cid:16) m z − qω M z (cid:17) + β ( S ) qω R − α ( w, S ) m − l m ( N ) m . (12)Integrating equation (9) with respect to z we get ε ∂ t R ( t, x ) = ( α ( w, S ) + l m ( N )) M ( t, x )+ ε (cid:90) Z µ ( W, N, S ) (cid:90) Z χ ( z, z (cid:48) , x ) r ( z (cid:48) ) Q ( x ) Q ∗ dz (cid:48) dz − (cid:0) β ( S ) + ε γ ( S ) (cid:1) R ( t, x ) . Using the fact that χ ( z, z (cid:48) , x ) is a probability kernel with respect to z , the previous equation for R ( t, x ) reducesto: ε ∂ t R = ( α ( w, S ) + l m ( N )) M + ε µ ( W, N, S ) Q ( x ) Q ∗ R − (cid:0) β ( S ) + ε γ ( S ) (cid:1) R . (13)Then, we multiply equation (8) by z and integrate it w.r.t. z , obtaining ε ∂ t m z + ε ∇ x · ( vm z ) − (cid:90) Z z∂ z (cid:20)(cid:18) zB ( Q, S ) − εk − B ( Q, S ) (cid:18) k +1 Q ∗ v · ∇ Q − k +2 S c, v · ∇ S (cid:19)(cid:19) p ( z ) (cid:21) dz = (cid:90) Z z L p [ λ ( z )] p ( z ) dz + β ( S ) qω R z − ( α ( w, S ) + l m ( N )) m z . The calculation of the integral term on the left hand side leads to the following equation for m z ( t, x, v ): ε ∂ t m z + ε ∇ x · ( vm z ) + B ( Q, S ) m z − εk − B ( Q, S ) (cid:18) k +1 Q ∗ v · ∇ Q − k +2 S c, v · ∇ S (cid:19) m = − λ m z + λ qω M z + β ( S ) qω R z − ( α ( w, S ) + l m ( N )) m z . (14)Applying the same procedure to equation (9) we obtain an equation for R z ( t, x ): ε ∂ t R z = ( α ( w, S ) + l m ( N )) M z + ε µ ( W, N, S ) Q ( x ) Q ∗ (cid:90) Z (cid:90) Z zχ ( x, z, z (cid:48) ) r ( z (cid:48) ) dz (cid:48) dz − (cid:0) β ( S ) + ε γ ( S ) (cid:1) R z . (15)6e consider Hilbert expansions for the previously introduced moments: m ( t, x, v ) = ∞ (cid:88) k =0 ε k m k , R ( t, x ) = ∞ (cid:88) k =0 ε k R k , m z ( t, x, v ) = ∞ (cid:88) k =0 ε k m zk R z ( t, x ) = ∞ (cid:88) k =0 ε k R zk ,M ( t, x ) = ∞ (cid:88) k =0 ε k M k , w ( x, t, ϑ ) = ∞ (cid:88) k =0 ε k w k , M z ( t, x ) = ∞ (cid:88) k =0 ε k M zk , W ( x, t ) = ∞ (cid:88) k =0 ε k W k . For the subsequent calculations it will be useful to Taylor-expand the coefficient functions involving any of w , W , R , M or N in the scaled equations (12)-(15) and (10): α ( w, S ) = α ( w , S ) + ∂ w α ( w , S ) ( w − w ) + 12 ∂ ww α ( w , S ) ( w − w ) + O ( | w − w | ) ,l m ( N ) = l m ( N ) + l (cid:48) m ( N ) ( N − N ) + 12 l (cid:48)(cid:48) m ( N ) ( N − N ) + O ( | N − N | ) ,µ ( W, N, S ) = µ ( W , N , S ) + ∂ W µ ( W , N , S )( W − W ) + ∂ N µ ( W , N , S )( N − N ) + O ( ε ) ,µ W ( W, Q ) = µ W ( W , Q ) + ∂ W µ W ( W , Q )( W − W ) + O ( | W − W | ) ,η ε ( x, ϑ, R ) = η ( x ) (cid:104) − εa ( R ) ϑ · ∇ R + ε (cid:16) − a ( R ) ∂ t R + 12 ( a ( R )) ( ϑ · ∇ R ) (cid:17) + O ( ε ) (cid:105) ,a ( R ) = a ( R ) + a (cid:48) ( R )( R − R ) + 12 a (cid:48)(cid:48) ( R )( R − R ) + O ( | R − R | ) . Then, equating the powers of ε in the scaled equations (12)-(15) and (10), we obtain: ε terms: 0 = − λ (cid:16) m − qω M (cid:17) + λ (cid:16) m z − qω M z (cid:17) + β ( S ) qω R − (cid:16) α ( w , S ) + l m ( N ) (cid:17) m , (16)0 = (cid:16) α ( w , S ) + l m ( N ) (cid:17) M − β ( S ) R , (17)0 = − (cid:16) B ( Q, S ) + λ + α ( w , S ) + l m ( N ) (cid:17) m z + λ qω M z + β ( S ) qω R z , (18)0 = (cid:16) α ( w , S ) + l m ( N ) (cid:17) M z − β ( S ) R z , (19)0 = η ( x )( S σn W − w ) , (20)where S σn := | Θ | = σ − n | S n − | . ε terms: ∇ x · ( vm ) = − λ (cid:16) m − qω M (cid:17) + λ (cid:16) m z − qω M z (cid:17) + β ( S ) qω R − (cid:16) α ( w , S ) + l m ( N ) (cid:17) m − (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) m , (21)0 = (cid:16) α ( w , S ) + l m ( N ) (cid:17) M + (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) M − β ( S ) R , (22) ∇ x · ( vm z ) = − B ( Q, S ) m z + k − B ( Q, S ) (cid:18) k +1 Q ∗ v · ∇ Q − k +2 S c, v · ∇ S (cid:19) m − λ (cid:16) m z − qω M z (cid:17) ++ β ( S ) qω R z − (cid:16) α ( w , S ) + l m ( N ) (cid:17) m z − (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) m z , (23)0 = (cid:16) α ( w , S ) + l m ( N ) (cid:17) M z + (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) M z − β ( S ) R z , (24) ∇ x · ( ϑw ) = η ( x ) (cid:16) S σn W − w (cid:17) + η ( x ) a ( R ) ϑ · ∇ R w − S σn η ( x ) a ( R ) (cid:90) Θ w ( ϑ (cid:48) ) ϑ (cid:48) dϑ (cid:48) · ∇ R . (25)7 terms: ∂ t m + ∇ x · ( vm ) = − λ (cid:16) m − qω M (cid:17) + λ (cid:16) m z − qω M z (cid:17) + β ( S ) qω R − (cid:16) α ( w , S ) + l m ( N ) (cid:17) m − (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) m − (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) m − (cid:16) ∂ ww α ( w , S ) w + l (cid:48)(cid:48) m ( N ) N (cid:17) m , (26) ∂ t R = (cid:16) α ( w , S ) + l m ( N ) (cid:17) M + (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) M + 12 (cid:16) ∂ ww α ( w , S ) w + l (cid:48)(cid:48) m ( N ) N (cid:17) M + (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) M + µ ( W , N , S ) Q ( x ) Q ∗ R − β ( S ) R − γ ( S ) R , (27) ∂ t w + ∇ x · ( ϑw ) = η ( x ) (cid:104)(cid:16) a ( R ) ∂ t R + a ( R ) ϑ · ∇ R + a (cid:48) ( R ) R ϑ · ∇ R −
12 ( a ( R ) ϑ · ∇ R ) (cid:17) w + a ( R ) ϑ · ∇ R w − w (cid:105) + S σn η ( x ) (cid:104) W − a ( R ) (cid:90) Θ ϑ (cid:48) w ( ϑ (cid:48) ) dϑ (cid:48) · ∇ R − a ( R ) ∂ t R W − (cid:90) Θ ϑ (cid:48) w ( ϑ (cid:48) ) dϑ (cid:48) · (cid:16) a ( R ) ∇ R + a (cid:48) ( R ) R ∇ R (cid:17) + 12 (cid:90) Θ ( a ( R ) ϑ (cid:48) · ∇ R ) w ( ϑ (cid:48) ) dϑ (cid:48) (cid:105) + µ W ( W , Q ) w = η ( x ) (cid:104)(cid:16) a ( R ) ∂ t R + a ( R ) ϑ · ∇ R + a (cid:48) ( R ) R ϑ · ∇ R −
12 ( a ( R ) ϑ · ∇ R ) (cid:17) w + a ( R ) ϑ · ∇ R w − w (cid:105) + S σn η ( x ) (cid:104) W − a ( R ) (cid:90) Θ ϑ (cid:48) w ( ϑ (cid:48) ) dϑ (cid:48) · ∇ R − a ( R ) ∂ t R W + 12 (cid:90) Θ ( a ( R ) ϑ (cid:48) · ∇ R ) w ( ϑ (cid:48) ) dϑ (cid:48) (cid:105) + µ W ( W , Q ) w , (28)due to (20).From (17) we get R ( t, x ) = α ( w , S ) + l m ( N ) β ( S ) M ( t, x ) (29)and from (19) we get R z ( t, x ) = α ( w , S ) + l m ( N ) β ( S ) M z ( t, x ) . (30)Integrating (18) with respect to v , we have0 = − B ( Q, S ) M z + β ( S ) R z − α ( w , S ) M z − l m ( N ) M z , from which by using (30) we obtain M z = 0 (31) R z = 0 . (32)These, together with (18), lead to m z = 0 . (33)8rom (16) and (18) we obtain:0 = − λ m + λ qω M + (cid:16) α ( w , S ) + l m ( N ) (cid:17) qω M − (cid:16) α ( w , S ) − l m ( N ) (cid:17) m ⇒ (cid:16) qω M − m (cid:17) (cid:16) λ + α ( w , S ) + l m ( N ) (cid:17) . Since ( λ + α ( w , S ) + l m ( N )) (cid:54) = 0 for any N , S, w , we obtain m = qω M . (34)From equation (20) we see that w = S σn W , (35)thus w depends on the (constant) speed σ , but not on the direction θ ∈ S n − .Now, turning to the equations stemming from the ε -terms, from (24) we get: R z ( x, t ) = α ( w , S ) + l m ( N ) β ( S ) M z ( x, t ) . (36)Integrating (23) with respect to v , we have upon using (34) and (33):0 = − B ( Q, S ) M z + k − B ( Q, S ) M s n E q · (cid:18) k +1 Q ∗ ∇ Q − k +2 S c, ∇ S (cid:19) + β ( S ) R z − (cid:16) α ( w , S ) + l m ( N ) (cid:17) M z . Using the fact that the tissue is undirected, i.e. E q = 0 and (36), the previous equation leads to: M z = 0 (37)and whence R z = 0 . (38)Now plugging these results into (23) we get0 = − B ( Q, S ) m z + k − B ( Q, S ) (cid:18) k +1 Q ∗ v · ∇ Q − k +2 S c, v · ∇ S (cid:19) m − λ m z − (cid:16) α ( w , S ) + l m ( N ) (cid:17) m z ⇒ m z [ B ( Q, S ) + λ + α ( w , S ) + l m ( N )] = k − B ( Q, S ) (cid:18) k +1 Q ∗ v · ∇ Q − k +2 S c, v · ∇ S (cid:19) m . We denote F ( Q, S ) := k − B ( S, Q ) [ B ( Q, S ) + λ + α ( w , S ) + l m ( N )] (39)and obtain therewith the following expression for m z : m z = F ( Q, S ) (cid:18) k +1 Q ∗ v · ∇ Q − k +2 S c, v · ∇ S (cid:19) m . (40)From equation (22) we have R = α ( w , S ) + l m ( N ) β ( S ) M + ∂ w α ( w , S ) w + l (cid:48) m ( N ) N β ( S ) M . (41)Using (21) we derive ∇ x · ( vm ) = − λ (cid:16) m − qω M (cid:17) + λ m z + (cid:16) α ( w , S ) + l m ( N ) (cid:17) qω M + (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) qω M − (cid:16) α ( w , S ) + l m ( N ) (cid:17) m − (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) m , L m [ λ + α ( w , S ) + l m ( N )] m := − (cid:16) λ + α ( w , S ) + l m ( N ) (cid:17) m + (cid:16) λ + α ( w , S ) + l m ( N ) (cid:17) qω M = ∇ x · ( vm ) − λ m z . (42)In order to get an explicit expression for m , we would like to invert the operator ¯ L m [ λ + α ( w , S )+ l m ( N )]. Ase.g. in [19, 71], we define it on the weighted L -space L q ( V ), in which the measure dv is weighted by q ( x, ˆ v ) /ω .In particular, L q ( V ) can be decomposed as L q ( V ) = < q/ω > ⊕ < q/ω > ⊥ . Due to the properties of the chosenturning kernel, ¯ L m [ λ + α ( w , S ) + l m ( N )] is a compact Hilbert-Schmidt operator with kernel < q/ω > . Wecan therefore calculate its pseudo-inverse on < q/ω > ⊥ .Thus, to determine m from (42) we need to check the solvability condition (cid:90) V [ ∇ x · ( vm ) − λ m z ] dv = 0 . This holds thanks to the above results and to the symmetry of q ( x, ˆ v ). Therefore, we obtain from (42) and (40) m = − λ + α ( w , S ) + l m ( N ) (cid:20) ∇ x · ( vm ) − λ F ( Q, S ) (cid:18) k +1 Q ∗ v · ∇ Q − k +2 S c, v · ∇ S (cid:19) m (cid:21) (43)and M = 0 . (44)On the other hand, (25) and (35) give: ∇ x · ( ϑw ) = − η ( x ) w + η ( x ) a ( R ) ϑ · ∇ R w + S σn η ( x ) W . (45)Likewise, we observe that the operator ¯ L w [ η ] w := − η ( x ) w + S σn η ( x ) W can be inverted, so that (45) leadsto w = − η ( x ) ∇ x · ( ϑw ) + w a ( R ) ϑ · ∇ R (46)and W = 0 . (47)From (27) we derive the following expression for β ( S ) ω R : β ( S ) ω R = 1 ω (cid:16) α ( w , S ) + l m ( N ) (cid:17) M + 12 ω (cid:16) ∂ ww α ( w , S ) w + l (cid:48)(cid:48) m ( N ) N (cid:17) M + 1 ω (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) M + µ ( W , N , S ) ω Q ( x ) Q ∗ R − γ ( S ) ω R − ω ∂ t R . (48)Plugging it into (26) we get: ∂ t m + ∇ x · ( vm ) = − λ (cid:16) m − qω M (cid:17) + λ (cid:16) m z − qω M z (cid:17) + qω (cid:16) α ( w , S ) + l m ( N ) (cid:17) M + q ω (cid:16) ∂ ww α ( w , S ) w + l (cid:48)(cid:48) m ( N ) N (cid:17) M + qω (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) M + qω µ ( W , N , S ) Q ( x ) Q ∗ R − qω γ ( S ) R − qω ∂ t R − (cid:16) α ( w , S ) + l m ( N ) (cid:17) m − (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) m − (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) m − (cid:16) ∂ ww α ( w , S ) w + l (cid:48)(cid:48) m ( N ) N (cid:17) m . (49)Integrating with respect to v we get ∂ t M + (cid:90) V ∇ x · ( vm ) dv = 12 (cid:16) ∂ ww α ( w , S ) w + l (cid:48)(cid:48) m ( N ) N (cid:17) M + (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) M + µ ( W , N , S ) Q ( x ) Q ∗ R − γ ( S ) R − ∂ t R − (cid:16) ∂ w α ( w , S ) w + l (cid:48) m ( N ) N (cid:17) M − (cid:16) ∂ ww α ( w , S ) w + l (cid:48)(cid:48) m ( N ) N (cid:17) M ∂ t M + (cid:90) V ∇ x · ( vm ) dv = µ ( W , N , S ) Q ( x ) Q ∗ R − γ ( S ) R − ∂ t R , (50)where (cid:90) V ∇ x · ( vm ) dv = (cid:90) V ∇ x · (cid:20) v (cid:18) − λ + α ( w , S ) + l m ( N ) (cid:18) ∇ x · ( vm ) − λ F ( Q, S ) (cid:18) k +1 Q ∗ v · ∇ Q − k +2 S c, v · ∇ S (cid:19) m (cid:19)(cid:19)(cid:21) dv = ∇ x · (cid:20)(cid:90) V − λ + α ( w , S ) + l m ( N ) v ⊗ v ∇ x (cid:16) qω M (cid:17)(cid:21) dv + ∇ x · (cid:34) λ F ( Q, S ) ω ( λ + α ( w , S ) + l m ( N )) (cid:90) V v ⊗ v q ( x, ˆ v ) dv (cid:18) k +1 Q ∗ ∇ Q − k +2 S c, ∇ S (cid:19) M (cid:35) . With the notation D T ( x ) := 1 ω (cid:90) V v ⊗ v q ( x, ˆ v ) dv = s (cid:90) S n − θ ⊗ θ q ( θ, x ) dθ = s V q ( x ) (51)and recalling that N ( t, x ) = M ( t, x ) + R ( t, x ), i.e., N = (cid:18) α ( w , S ) + l m ( N ) β ( S ) (cid:19) M , (52)we obtain from (50) the following macroscopic equation for N ( t, x ): ∂ t N − ∇ x · (cid:20) λ + α ( w , S ) + l m ( N ) ∇ x · (cid:16) β ( S ) β ( S ) + α ( w , S ) + l m ( N ) D T ( x ) N (cid:17)(cid:21) + ∇ x · λ F ( Q, S ) β ( S ) λ + α ( w , S ) + l m ( N ) D T ( x ) k +1 Q ∗ ∇ Q − k +2 S c, ∇ Sβ ( S ) + α ( w , S ) + l m ( N ) N = α ( w , S ) + l m ( N ) α ( w , S ) + l m ( N ) + β ( S ) N (cid:16) µ ( W , N , S ) Q ( x ) Q ∗ − γ ( S ) (cid:17) . (53)We denote ϕ ( w , N, S ) := β ( S ) β ( S ) + α ( w , S ) + l m ( N ) (54) (cid:37) ( w , N, S ) := ( λ + α ( w , S ) + l m ( N )) . (55)Observe that, due to (20), these are, in fact, purely macroscopic quantities.Integrating (28) with respect to ϑ ∈ Θ gives ∂ t W + ∇ x · (cid:90) Θ ϑw dϑ = µ W ( W , Q ) W , again due to (20). Recalling (46), this leads to the following macroscopic equation for the density W ofendothelial cells: ∂ t W − ∇ · (cid:16) D EC ∇ W (cid:17) + ∇ · (cid:16) χ a ( R ) W ∇ R (cid:17) = µ W ( W , Q ) W , (56)where D EC ( x ) := σ nη ( x ) I n and χ a ( R ) := σ a ( R ) n I n = η ( x ) a ( R ) D EC ( x ).The two macroscopic equations obtained in (53) and (56) for the evolution of the tumor and ECs, respectively,are coupled with a PDE for proton concentration dynamics: ∂ t S = D S ∆ S + g ( S, N , W , Q ) , (57)with the concrete form of g ( S, N , W , Q ) given in (65).11n view of (29), (41), (44), and (47), the ε -correction terms for N and W can be left out and ignoring the higherorder terms we get the following closed PDE system characterizing the macroscopic evolution of the tumorunder the influence of tissue, vasculature, and acidity: ∂ t N − ∇ · (cid:20) (cid:37) ( W, N, S ) ∇ · (cid:16) ϕ ( W, N, S ) D T ( x ) N (cid:17)(cid:21) + ∇ · (cid:20) λ F ( Q, S ) ϕ ( W, N, S ) (cid:37) ( W, N, S ) D T ( x ) (cid:18) k +1 Q ∗ ∇ Q − k +2 S c, ∇ S (cid:19) N (cid:21) = ϕ ( W, N, S ) α ( W, S ) + l m ( N ) β ( S ) N (cid:16) µ ( W, N, S ) QQ ∗ − γ ( S ) (cid:17) ,∂ t W = ∇ · (cid:16) D EC ∇ W (cid:17) − ∇ · (cid:16) η a ( R ) D EC W ∇ R (cid:17) + µ W ( W, Q ) W,∂ t S = D S ∆ S + g ( S, N, W, Q ) , (58)with F ( Q, S ) given in (39), the tumor diffusion tensor D T from (51), and with the coefficients ϕ and ρ from(54) and (55), respectively, whereby we used w = W | Θ | , in virtue of (20). This system has to be supplementedwith adequate initial conditions; its deduction has been carried out for x ∈ R n . For the numerical simulationsto be performed in Section 4, we will consider it to be set in a bounded, sufficiently regular domain Ω ⊂ R n and endow it with no-flux boundary conditions.The next section is dedicated to specifying the remaining coefficients of the above macroscopic system. To determine the tumor diffusion tensor D T ( x ) in (51) we need to provide a concrete form for the (mesoscopic)orientational distribution of tissue fibers q ( x, θ ). Several different choices are available in the literature, seee.g. [14, 71]. As in [45] we use the orientation distribution function (ODF): q ( x, θ ) = 14 π | D W ( x ) | ( θ T ( D W ( x )) − θ ) , where D W ( x ) is the water diffusion tensor obtained from processing the (patient-specific) DTI data.We choose for the macroscopic tissue density Q ( x ) the expression proposed in [20]: Q ( x ) = 1 − l c ( x ) h , (59)where h is the side length of one voxel of the DTI dataset and l c a characteristic length, estimated as l c ( x ) = (cid:115) tr ( D W ( x )) h l with l being the leading eigenvalue of the diffusion tensor D W .For the rate α ( w, S ) describing the cell phenotypic switch p → r from migration to proliferation, we choosea combination of an increasing function of w and a decreasing function of S . As explained in Subsection 2.2,this rate is influenced by the availability of nutrient (provided by vasculature) for sustaining the processesinvolved in the cell cycle and by the pH of the environment; the dynamical balance of these factors decides themigratory/mitotic fate. Recalling (35), we set α ( W, S ) = α S σn W/W c, S σn + S σn W/W c,
11 +
S/S c, , (60)where W c, and S c, are reference values for EC density and proton concentration, respectively.The second rate of phenotypic switch p → r l m ( N ), describing the influence of a crowed environment on cellphenotype changes, is chosen as l m ( N ) = l m, (1 + tanh( N/N c, − N ∗ /N c, )) , (61)where N ∗ represents a threshold value for glioma density: when it is exceeded, the cells are not able to moveanymore. The constant N c, denotes a reference value for the density of (moving and proliferating) glioma cells.12he switching rate r → p given by the function β ( S ) controls the acidity-triggered motility enhancement offormerly proliferating cells. We set β ( S ) = β [ ε + ( S/S c, − S T, /S c, ) + ] (62)with ε (cid:28) · ) + the positive part. S T, is the pH threshold which, when underrun, induces the cells toswitch from a proliferative to a migrative phenotype.Although tumor cells can live in an environment with substantially lower pH than that for normal tissue [23,90],when it drops below a certain threshold (which in terms of proton concentration we denote by S T, ), the cancercells become necrotic [24, 90]. This suggests the following cell death coefficient γ ( S ): γ ( S ) = γ ( S/S c, − S T, /S c, ) + . (63)The growth rate µ ( W, N, S ) can be also defined in different ways and it should be motivated by biologicalevidence. As, for instance, in [20], we choose a logistic-like function to describe the growth self-limitation anda growth enhancement factor depending on the vascularization W , along with a growth-limiting one due toacidification, like that employed in (60): µ ( W, N, S ) = µ N, (cid:18) − NK N − c e N e K N e (cid:19) WW c,
11 +
S/S c, . (64)Here, K N is the tumor carrying capacity and µ N, is a constant. The term − c e N e K Ne is related to the extensionof the model described in Section 5. In particular, c e = 0 when we consider the evolution of system (58), while c e = 1 when the dynamics of healthy tissue and necrotic tissues ( N e ) are included. K N e represents the carryingcapacity for the necrotic component.Similarly, for the term µ W ( W, Q ) describing proliferation of ECs we take µ W ( W, Q ) = µ W, (cid:18) − WK W (cid:19) QQ ∗ , with K W and Q ∗ the carrying capacities for ECs and healthy tissue, respectively.For the tactic sensitivity χ a ( R ) = η a ( R ) D EC of ECs towards (gradients of) proliferating glioma cells, we needto specify the function a ( R ) involved in the definition (7) of the EC turning rate. We choose a ( R ) = χ a K N ( K N + R ) , which corresponds to the rate of change of the expression ζ ( R ) = RR + K N representing the equilibrium of theinteractions between ECs and proliferating glioma cells R , scaled by a constant χ a that is used to account forchanges in the turning rate per unit of change in dζ/dt . Thereby, we assume that attachment and detachmentof ECs and R -cells happen with the same rates.The reaction term in the PDE for acidity dynamics is chosen as g ( S, N, W, Q ) = g s NN c, − g d (cid:18) WW c, + QQ ∗ (cid:19) S, (65)hence it has a source term for the production of protons by the tumor, and a decay term, as the protons arebuffered by healthy tissue and also uptaken by the capillary network.In Table 1 we report the range of the values for the constant parameters involved in system (58), as well as thereferences from which they were drawn. In order to simplify, we assume a constant coefficient η in the turningrate (7) for ECs. Further details on the estimation of the parameters and on the nondimensionalization of (58)are given in Appendix A and B. We perform 2D simulations of the resulting macroscopic system of coupled advection-diffusion-reaction equations(58). For the initial conditions we take a Gaussian-like aggregate of tumor cells centered at ( x ,N , y ,N ) =( − , R ( x, y ) = e − ( ( x − x ,N )2+( y − y ,N )2 ) λ turning frequency in L p [ λ ( z )] 0 .
001 (s − ) [79] λ turning frequency in L p [ λ ( z )] 0 .
001 (s − ) [19, 79] α phenotype switch rate p → r . s − ) [54, 72] β phenotype switch rate r → p . − ) [54, 72] l m, overcrowding switch rate p → r . − ) [54, 72] N ∗ optimal tumor density value for cell movement 0 . · K N [8, 74] s speed of tumor cells 0 . · − (mm · s − ) [18] k +1 attachment rate of tumor cells to tissue fibers 0 .
034 (s − ) [58] k +2 interaction rate between tumor cells and protons 0 .
01 (s − ) [58] k − detachment rate 0 .
01 (s − ) [58] µ N, tumor proliferation rate 9 . · − (s − ) [49] K N tumor carrying capacity ∼ (cells · mm − ) [1] γ acid-induced death rate for tumor cells 0 . · − (s − ) [86] S T, proton concentration threshold for r → p . · − (M) [88, 90] S T, proton concentration threshold for tumor cell death 3 . · − (M) [88, 90] η turning frequency of ECs in L w [ η ] 0 .
001 (s − ) [87] χ a duration between R -damped EC turnings 4 . σ speed of ECs 0 . · − (mm · s − ) [17] K W carrying capacity for ECs ∼ (cells · mm − ) [2] µ W, EC proliferation rate 0 . · − (s − ) [4, 30] D S diffusion coefficient of protons 0 . · − (mm · s − ) [59] g s proton production rate 2 . · − (M · mm · (cells · s ) − ) [62] g d proton removal rate 0 . · − (s − ) [63] Table 1: Model parameters and a ring-like profile for the migrating tumor cells centered at the same location: M ( x, y ) = 0 . e − (cid:16) √ ( x − x ,N ) +( y − y ,N ) − (cid:17) . The initial distribution of the total tumor population is given by N ( x, y ) = R ( x, y ) + M ( x, y ). For the ECswe consider x ,W = − W ( x, y ) = 0 . e − ( x − x ,W )20 . sin (cid:16) π y (cid:17) ∀ y ∈ [ − , x ,S , y ,S ) =( − , S ( x, y ) = 0 . e − ( ( x − x ,S )2+( y − y ,S )2 ) . The initial pH distribution is calculated considering that pH = − log ( S ). Figure 1 shows the plots for theinitial conditions on the entire 2D brain slice, zooming then in the region ¯Ω = [ − , × [ − , D T ( x ) is precalculated using DTI data and the ODF for the fiber distribution function (seeSection 3). The dataset was acquired at the Hospital Galdakao-Usansolo (Galdakao, Spain), and approved byits Ethics Committee: all the methods employed were in accordance to approved guidelines. A Galerkin finiteelement scheme for the spatial discretization of the equations for tumor cells, ECs, and proton concentration is14 igure 1: Initial conditions for Model (58).Figure 2: Healthy tissue density for Model (58). considered, together with an implicit Euler scheme for the time discretization. We present different simulationsaddressing several aspects.Firstly, we study the behavior of species involved in (58) for the parameters listed in Table 1. The correspondingsimulation results are shown in Figure 3, where the five columns report the evolution of the whole tumor mass( N ), the two subpopulations of proliferating ( R ) and migrating ( M ) tumor cells, the endothelial cells ( W ), andpH (computed from S ). The tumor spread, which seems to mainly depend on the choice of the parameters λ and s , as well as on the EC evolution, is rather slow, with a partial exchange between the two subpopula-tions of tumor cells in relation to pH at the core of the tumor mass. The tumor cells increasingly adopt theproliferating phenotype when they approach ECs, as they provide the necessary nutrient and oxygen supply tosustain glioma proliferation. On the other hand, ECs diffuse and grow, with a higher accumulation around thefirst of the three vessels situated in the upper part of the plots and where there is more healthy tissue availableto sustain their proliferation. They clearly exhibit tactic behavior toward the (pro-angiogenic growth factorsreleased by) proliferating tumor cells mainly located around the core of the neoplasm. In particular, the sub-plots for the evolution of ECs at later times (last two rows of Figure 3) show an increasing amount of high ECaggregates developing towards the tumor. This behavior can be associated with the phenomenon of microvas-cular hyperplasia and glomeruloid bodies. The latter are tumor-associated vascular structures that develop inthe presence of high levels of VEGF and are important histopathological features of glioblastoma multiforme [77].In Figure 4 we show comparisons between simulations done upon varying the parameters λ , referring to theturning rate of glioma cells, and the two speeds ( s, σ ) for tumor and endothelial cells, respectively, with the aimto illustrate how sensitive the model predictions are w.r.t. these parameters. The tumor and EC densities areplotted after 560 days of evolution for three different values of λ (expressed in s − ), i.e., 10 − , 10 − and 10 − ,and for four pairs ( s, σ ) of speed values (expressed in µ m · h − ), i.e., (15 , , , , λ accompanied by relatively large EC speed σ . A too small λ effects the (biologically rather unrealistic) shiftof tumor cells from their original location to the site of blood vessels, where they switch to the proliferative15 igure 3: Numerical simulation of Model (58) with parameters listed in Table 1. phenotype. The faster the glioma cells are, the more pronounced is this behavior, obviously dominated bymigration during its first stage and subsequent proliferation. Increasing λ by one or two orders of magnitudeleads to more realistic behaviors of tumor cells and ECs, with less sensitivity towards variations in λ . Naturally,wide-spread hyperplasia and pronounced tumor invasion occur for higher cell speeds.In order to test the effect of the go-or-grow dichotomy on the evolution of cell populations involved in system(58), we compare that model with a setting in which the tumor cells migrate and proliferate without deterringone of these phenotypes for the other. In particular, we do not differentiate between proliferating and migratingcells and accordingly let the ECs be biased by the density gradient of the whole tumor. Using a scaling argumentsimilar to the one described in Section 2.3, we get a system of three partial differential equations for tumorcells ( N ), ECs ( W ), and protons ( S ), which is analogous to (58) with α = 0 and l m, = 0. This choice ofparameters reduces the former coefficient functions to ϕ ( w , N, S ) = 1 and (cid:37) ( w , N, S ) = λ . We refer to thissetting as Model NGG . Figure 5 shows the solution behavior for this new setting. The initial conditions for thethree populations are the same as those shown in Figure 1. Comparing Figures 3 and 5 we observe that thego-or-growth model predicts -as expected- a slower tumor spread, with lower cell density, that, consequently,induces lower acidity concentrations in the environment, and the differences between the two settings becomingmore accentuated with increasing time. Moreover, the taxis driving ECs towards the tumor mass is strongerfor the case shown in Figure 5, and accumulations of ECs indicating microvascular hyperplasia are now earlierformed and become larger.To enable a direct assessment of the two settings we plot in Figure 6 the differences (at 400 and 560 days) betweenthe (overall) densities of tumor and endothelial cells for the model with go-or-grow and its
NGG counterpart(the quantities for the latter are marked by the index 1), as well as between the respective pH distributions, thelatter illustrated on a larger domain ˜Ω = [ − , × [ − , ", $) (30,25)(30,20)(20,15)(15,10) + , -. -/ -0 Figure 4: Comparison between the evolution of tumor (first row of each box) and endothelial cells (second row of eachbox) for three different values of λ (s − ): 10 − , 10 − and 10 − , and four pairs ( s, σ ) of speed values ( µ m · h − ): (15 , , , , We extend the above model by considering the evolution of macroscopic tissue density Q and that of necroticmatter N e , the latter including tissue as well as glioma cells degraded by hypoxia. The whole system consists17 igure 5: Numerical simulation of Model NGG , i.e., with simultaneously moving and proliferating cancer cells. Theemployed parameter values are listed in Table 1, except now α = 0 and l m, = 0. of equations (58) together with the following ODEs: ∂ t Q = − d Q ( S ) Q∂ t N e = d Q ( S ) Q + ϕ ( w, N, S ) α ( w, S ) + l m ( N ) β ( S ) γ ( S ) N, (66)where the coefficient d Q ( S ) models the pH-triggered tissue degradation occurring when a certain acidity thresh-old S T,Q is exceeded d Q ( S ) = d ,Q (cid:18) SS c, − S T,Q S c, (cid:19) + . Thereby d ,Q > · ) + means as usual the positive part. Details on theestimation of d Q, and S T,Q are given in Appendix A.Numerical simulations for the extended model (58), (66) are shown in Figure 8, corresponding to the initial con-ditions in Figure 7. Figure 9 illustrates the evolution of proliferating and migrating glioma cells. Although thequalitative behavior of glioma and endothelial cells is comparable with that shown in Figure 3, the degradationof tissue by environmental acidity affects both tumor and EC proliferation, as less healthy tissue is available tosustain growth. Therefore, the simulations show lower densities for both species and, particularly, tumor cellgrowth is affected by the reduction of vasculature and the depletion of healthy tissue, as clearly shown by theevolution of the proliferating tumor cells in Figure 9 (first row). Moreover, Figure 8 shows that lower tumor18 igure 6: Differences between the respective solution components of Model (58) and of
Model NGG at 400 (upper row)and 560 (lower row) days. Figure 7: Initial conditions for the extended Model (58), (66). density effects less proton extrusion, i.e. higher pH.The above model extension enables us to perform necrosis-based tumor grading, which is essential for assessingpatient survival and treatment planning. Other indicators of tumor aggressiveness are employed as well (e.g.histological patterns [92] or tumor size [75]), however, we focus here on grading by the amount of necrosisrelative to the whole tumor volume, in view of [29, 36], where the tumor volume by itself was found to have noinfluence on overall survival. Following [15], we define the time-dependent grade G ( t ) ∈ [0 ,
1] of the simulatedtumor via: G ( t ) := V Ne ( t ) V Ne ( t ) + V N ( t ) (67)where V Ne ( t ) and V N ( t ) denote, respectively, the fractions of necrosis and living cell densities in the visibletumor volume. They are defined as the integrals of the densities N e and N over the domain defined by the level19 igure 8: Numerical simulation of the extended Model (58), (66). Parameters are listed in Table 1 and the value of ECsspeed is here set to σ = 0 . · − mm · s − . sets of the tumor population for a detection threshold of 80% of the carrying capacity, which corresponds to thedetection threshold for T1-Gd images [86]. We represent in Figure 10 the time evolution of G , guided by thepercentage classification in [29], i.e., 0 < G < ≤ G < G ≥
50: grade 3.The highest grade corresponds to the most aggressive tumor and the poorest survival prognosis. In particular,we compare the effect of four different scenarios on necrosis-based tumor grading: the grey curves therein referto the extended model (58), (66)) involving vascularization, i.e. EC density W (solid line: go-or-grow, dottedline: Model NGG ), while the red curves illustrate the evolution of G for the corresponding variants of theextended model without EC dynamics (i.e., without W ).Figure 10 provides rich information. On the one hand, it shows that assuming the go-or-grow dichotomy leadsto slower progression of neoplasia, due to the cells deterring one phenotype for the other. Indeed, in the longrun, the full go-or-grow model with vascularization predicts a slower advancement towards high tumor grades.When EC dynamics are accounted for, the differences between the model with or without go-or-grow are rathersmall; when W is, instead, not included, then such differences increase. The vascularization seems to havea significant impact: if we focus e.g., on Model NGG (dotted curves) we see that during the first simulated14-15 weeks the vascularization ensures a higher percentage of necrosis, after which the situation reverses, withdifferences becoming larger while time is advancing. The early phase (which is supposed to correspond toa lower tumor grade) might seem somehow paradoxical when thinking about blood capillaries buffering theacidity and reducing necrosis. It is, however, well-known that a tumor usually develops angiogenesis when ithas reached a more advanced phase in its development and begun to get increasingly hypoxic, which leads toenhanced VEGF expression and capillary formation. The small amount of vessels prior to such stage is on theone side supporting the growth of tumor cells, while on the other side it is not able to buffer the ever increasingproton concentration triggered by the exuberant growth. Without a substantial enhancement of angiogenesis,the tumor will develop a larger necrotic component, thus receiving a higher grade. Therefore, omitting thedynamics of endothelial cells from the model might overestimate the tumor growth and spread, which for brain20 igure 9: Numerical simulation of the extended Model (58), (66). Here the evolution of the proliferating and of themigrating tumor cells is shown.Figure 10: Evolution of the grade function G ( t ) given in (67). Grey curves relate to the extended full model (i.e.,including dynamics of endothelial cells) in the case with go-or-grow (solid line, model (58), (66)) and without themigration-proliferation dichotomy (dotted line, Model NGG ). The red curves refer, respectively, to the same modelvariants, but without vascularization. In both cases we set µ N, = 5 . · − s − . tumors can have a significant therapeutic impact. To our knowledge this is the first continuous mathematical model with pH- and vasculature-induced phenotypicswitch between moving and proliferating cells where the two cell types are seen as distinct tumor subpopulationsevolving under mutual, direct and indirect interactions. As such, it can be seen as a further development of themodels in [21, 45]. This novel model not only includes on the macroscopic scale nonlinear, myopic self-diffusionand multiple taxis (haptotaxis, repellent pH-taxis) for one of the tactic populations (glioma cells), but alsofeatures taxis of the second population (ECs) towards only one part of the former: the ECs are following, in21act, the gradient of proliferating glioma. Thereby, the chemorepellent is produced by the population (gliomacells) which tries to avoid it, but at the same time is degraded by the other tactic actor (ECs), which is, inturn, directionally biased by (a part of) the first one. This renders the repellent acidity taxis of glioma cellsboth direct and indirect. The extended model involving dynamics of tissue and necrotic matter involves, too,an indirect kind of haptotaxis, with yet more complex couplings. These features raise several highly interestingquestions related to rigorous well-posedness and long term behavior of solutions. Among others, the simula-tions showed that for certain parameter combinations the solution can infer singularities. For a recent reviewof available models with multiple taxis we refer to [55].We used a multiscale approach starting from single cell dynamics for glioma interacting with (macroscopicallyrepresented) extracellular acidity and tissue, we wrote the corresponding kinetic transport equations for gliomaand EC density distribution functions on the mesoscopic scale, and employed a parabolic scaling to deduce thepopulation level behavior for the model variables: glioma density, acid concentration (or pH), EC density, and-for the extended model- densities of normal tissue and necrotic matter. Thereby, tumor cell migration is drivenas in previous works, e.g, [19, 21, 71], by a myopic diffusion term involving a tumor diffusion tensor which takesinto account the local tissue structure. Moreover, taxis terms carrying information from the microscopic scaledirect glioma migration towards increasing tissue gradients and away from highly acidic regions. We consid-ered a heterogeneous tumor and used the go-or-grow dichotomy asserting that glioma cells can either move orproliferate, the respective behavior being transient and switching according to nutrient availability (supplied byvasculature), pH (determined by proton production and buffering), as well as crowded environments.The deduced macroscopic motion of ECs is characterized, too, by a combination of diffusion and drift terms. Inorder not to increase too much the number of solution components, we described rather indirectly the responseof ECs to pro-angiogenic growth factors (e.g. VEGFs), upon biasing their tactic motion towards proliferatingtumor cells. This was realized in an alternative manner to that by which taxis terms have been obtained forglioma cells.The extension of the model by taking into account the macroscopic evolution of healthy tissue and necroticmatter opened the way for necrosis-based tumor grading, which is highly relevant for diagnosis and therapyplanning. A forthcoming work will be dedicated to several aspects of the latter issue.The validation of model predictions through comparison with adequate patient data would be a future taskto address. Although the available technology is potentially able to provide a large amount of the neededquantitative information, the standard clinical imaging and therewith associated treatment does not providesuch data sets, among others due to prohibitively high costs. With the advent of technological developmentand the advancement of personalized medicine, however, such data may become available in the near future.
Acknowledgements
The authors thank Prof. Juan Soler for the helpful discussions and his valuable advices on modelling issues. Theauthors also thank Dr. Marian Gomez Beldarrain from the Department of Neurology at the Galdakao-UsansoloHospital (Galdakao, Spain) for the DTI data. MC acknowledges funding by the Basque Government throughthe BERC 2018- 2021 program and by the Spanish State Research Agency through BCAM Severo Ochoaexcellence accreditation SEV-2017-0718. This project has received funding from the European Unions Horizon2020 research and innovation programme under the Marie Sk(cid:32)lodowska-Curie grant agreement No. 713673.The project that gave rise to these results received the support of a fellowship from la Caixa Foundation (ID100010434). The fellowship code is LCF/BQ/IN17/11620056. CS was supported by the Federal Ministry ofEducation and Research BMBF, project
GlioMaTh
A Parameter assessment
A.1 Turning rates and diffusion related parameters: λ , λ , s , σ , η , D S λ , λ : In [79] the authors presented experiments on the migratory behavior and turning frequency of metastaticcancer cells from rat mammary adenocarcinoma cell line, reporting values for λ in the range [0 . − .
1] s − .Considering the highly aligned brain structure which influences cell migration upon enhancing cell persistencein the favorable direction of motion, we assume a reduction of the turning likelihood and take the range[10 − − − ] s − for the parameter λ . The choice of λ is, unfortunately, rather imprecise, since we found nodata or references. In [19] variations of this parameter by ±
50% were tested in a similar context. As proposedthere, we consider the same order of magnitude for λ and λ . s : Different references are available for the tumor cell speed. In a recent work [18], four types of typical GBMcell lines were cultured in a microfabricated 3D model to study their in vitro behavior; according to that, we22onsider for glioma cell speed the range [10 . − µ m · h − . A further upper limit for this parameter can befound in [76], where a maximum speed of 54 − µ m · h − was reported for glioma cells. σ : For the average speed of ECs, in [17] the range σ ∈ [10 − µ m · h − was reported for individual cells inseveral culture conditions. In vivo, the registered mean speed for motile endocardial and endothelial cells is ofapproximately 20 µ m · h − . η : In [87], the authors analyzed the statistical properties of the random streaming behavior for endothelialcell cultures. In particular, they estimated a needed time T = 5m to alter cell polarity and influence the cellturning. With this value, they obtained a cell speed estimation in the range of 2040 µ m · h − within monolayers,in agreement with their empirical data about cell speed in [10 − µ m · h − . Following their results we consider η ∈ [0 . − . − . D S : From [59], we get an estimation of the diffusion coefficient for different ions; an average value between[0 . − · − mm · s − is reported. We test the model for different values of this parameter and we chose D S = 0 . · − mm · s − . D S was rescaled in the simulations in order to account for the fast dynamicscharacterizing proton evolution. A.2 Phenotypic switch related parameters: α , β , l m, , N ∗ , S T, , S T, , γ α , β , l m, : The estimation of reasonable ranges for the values of the phenotypic switch parameters α and β might be quite imprecise, since there are no specific data or references available. However, considering theexperiment done in [54, 72], it is possible to define a wide range for the duration of the glioma cell cycle thattranslates into α , β ∈ [0 . − · − s − . Analogous arguments apply for the rate l m, , for which there arenot estimations directly derived from biological data. For these reasons we test the model for several values of α , β , l m, with order of magnitude 10 − . N ∗ : Due to the lack of biological data, for the estimation of this parameter related to the total tumor densitylevel that still allows cells movement, we choose it to be proportional to the estimated tumor carrying capacity,referring to the range [0 . − . · K N proposed in [8, 74]. S T, , S T, : Following the pH range [6 . − .
3] proposed in [88, 90] for the brain tumor microenvironment, wechoose the thresholds values that determines the phenotypic switch from a proliferating to a migrating cells( S T, ) and the acid-mediated death of the resting cells ( S T, ). We set S T, = 1 . · − M (referring to apH= 6 . S T, = 3 . · − M (referring to a pH= 6 . γ : In [86], the authors assumed the cell necrosis rate to be proportional to their metabolic rate µ N, . Inparticular, this estimation appears reasonable, considering also our assumption in the parabolic scaling procedureabout having similar time scales for birth and death. Therefore, we set γ = µ N, / A.3 Adhesion related parameters: k +1 , k +2 , k − , χ a k +1 , k +2 : For the estimation of the attachment rates between tumor cell and ECM or protons, we referto [58]. In particular, for both the cell-ECM attachment rate and the cell-protons interaction rate, we set k +1 = k +2 = 10 (M s) − . Then, assuming that the main ECM component is collagen, with a molecular weight of ≈ Q ∗ (in Table 2), we deduce k +1 = 0 .
034 s − . Analogously,considering the reference value for the protons concentration S c, (in Table 2), we get k +2 = 0 .
01 s − . k − : For the estimation of the cells-ECM and cells-protons detachment rate, referring to [58], we set k − =0 .
01 s − . χ a : We estimate the parameter χ a by considering the values reported in [84] and [83] for the chemotacticsensitivity. In particular, in [84] the authors analyzed the chemotactic coefficient of migrating endothelial cellsin gradients of aFGF, measuring a maximum chemotactic response of 2600 cm · (M s) − at a concentration ofaFGF around 10 − M. Instead, in the further work [83], the authors analyzed the changes of this parameter inresponse to cell speed and persistence time. Taking into account the above described range for the EC speed,we get χ a ∈ [3 . − .
5] d.
A.4 Proliferation related parameters: µ N, , K N , µ W, , K W µ N, : For the estimation of glioma growth rate, we analyze the doubling times reported in [49] for severalglioma cell lines. There, the authors reported a range of variably between 21 .
1h and 46h, which translated into µ N, ∈ [0 . − . · − s − . K N , K W : Considering that the mean diameter of a glioma cell is around [12 − µ m [1], we estimate a valuefor the tumor carrying capacity of K N ∼ cells · mm − . In the same way, considering a mean diameter foran endothelial cell of [10 − µ m [2], we set K W ∼ cells · mm .23 W, : The doubling time of EC density has been estimated in several experiments to vary between the differentphases of an endothelial colony growth until the formation of a monolayer. In [4, 30], the authors gave a rangeof variability for the value of the EC doubling time of [3 −
13] d. This leads to µ W, ∈ [0 . − . · − s − . A.5 Production and consumption related parameters g ,N , g ,W , d ,Q , S T,Q g s : In [23], the authors estimated the rate of proton production due to tumor cell activity by fitting theirequation for proton dynamics (analogous to our PDE (57)) to a converted form of the data in [62]. In particu-lar, in [62] pH measurements were taken at a variety of points within both the tumor and surrounding healthytissue for four composite cases, giving a geometric mean for the production rate of 2 . · − M · mm · (s · cell) − . g d : Following [23, 63], for the rate of proton uptake by vasculature we consider the range of variability given by[0 . − . · − s − . As for D S , proton production and consumption rates were rescaled in the simulations inorder to account for the fast dynamics characterizing proton evolution. d ,Q : For the rate of tissue degradation due to the acidic environment, in [86] the authors proposed an esti-mation choosing the parameter such that 10% necrosis gives tissue a 50-day half-life. Starting from the valueproposed in [86], we test a wider range of possible estimations, that translates into d ,Q ∈ [0 . − .
07] d − . S T,Q : In [88], tissue pH values in normal brain and in brain tumors were reported. Specifically, consideringthat these values vary depending on the type of brain tissue (i.e., gray matter, white matter, cerebellum), theminimum pH requires for the normal cell activity is in the range of [6 . − . S T,Q = 1 . · − M (referring to a pH= 6 . B Nondimensionalization
To proceed with nondimensionalizing the system, we firstly observe that the variables N , W , and N e involved insystems (58) and (66) are expressed in cells/mm , Q in g/mm , while the concentration of protons S is given inmol/liter (=:M). The reference values we use for the nondimensionalization are listed in Table 2. In particular,we rescale the tumor, EC, and necrotic matter (dead cells and tissue) densities with respect to their carryingcapacities, i.e., N c, = K N , W c, = K W , and N e, = K N , assuming a similar carrying capacity for tumor cellsand necrosis. Parameter Description Value (units) SourceT time 1 (d)L length 0 .
875 (mm) N c, tumor cell density 10 (cell · mm − ) this work W c, EC density 10 (cell · mm − ) this work S c, Proton concentration 10 − (M) [88] Q ∗ healthy tissue density 10 − (mg · mm − ) [48] N e, density of necrotic matter 10 (cell · mm − ) this work Table 2: Reference variables for the nondimensionalization.
We nondimensionalize the partial differential equations introduced above as follows:˜ t = tT , ˜ x = xL , ˜ N = NK N , ˜ W = WK W , ˜ S = SS c, , ˜ Q = QQ ∗ , ˜ N e = N e K N . The proper scaling of the parameters involved in the macroscopic setting then reads24 α = α λ , ˜ l m, = l m, λ , ˜ N ∗ = N ∗ K N , ˜ β = β λ ˜ S T,j = S T,j S c, ( j = 1 , , ˜ k +1 = k +1 λ , ˜ k +2 = k +2 λ , ˜ k − = k − λ , ˜ µ N, = µ N, T, ˜ γ = γ T, ˜ λ = λ λ , ˜ D T = 1 λ TL D T , ˜ D EC = TL D EC , ˜ µ W, = µ W, T, ˜ χ a = χ a K N , ˜ D S = TL D S , ˜ g s = g s TS c, , ˜ g d = g d T, ˜ d ,Q = d ,Q T, ˜ S T,Q = S T,Q S c, Dropping the tilde (” ˜ ”) in the new variables and parameters, the differential equations in system (58) keepthe same form, with the following rescaled functions:˜ α ( ˜ W , ˜ S ) = ˜ α ˜ W W
11 + ˜
S , ˜ l m ( ˜ N ) = ˜ l m, (1 + tanh( ˜ N − ˜ N ∗ )) , ˜ ρ ( ˜ W , ˜ N , ˜ S ) = 1 + ˜ α ( ˜ W , ˜ S ) + ˜ l m ( ˜ N ) , ˜ β ( ˜ S ) = ˜ β ( ε + ( ˜ S − ˜ S T, ) + ) , ˜ ϕ ( ˜ W , ˜ N , ˜ S ) = ˜ β ( ˜ S )˜ β ( ˜ S ) + ˜ α ( ˜ W , ˜ S ) + ˜ l m ( ˜ N ) , ˜ a ( ˜ R ) = ˜ χ a (1 + ˜ R ) , ˜ µ ( ˜ W , ˜ N , ˜ S ) := ˜ µ N, (cid:16) − ˜ N − ˜ N e (cid:17) ˜ W
11 + ˜
S , ˜ µ W ( ˜ W , ˜ Q ) := ˜ µ W, (cid:16) − ˜ W (cid:17) ˜ Q, ˜ γ ( ˜ S ) = ˜ γ ( ˜ S − ˜ S T, ) + . ˜ g ( ˜ N , ˜ S, ˜ W , ˜ Q ) = ˜ g s ˜ N − ˜ g g ( ˜ W + ˜ Q ) ˜ S, ˜ B ( ˜ Q, ˜ S ) = (cid:16) ˜ k +1 ˜ Q + ˜ k +2 ˜ S + ˜ k − (cid:17) , ˜ F ( ˜ Q, ˜ S ) = ˜ k − ˜ B ( ˜ S, ˜ Q ) (cid:104) ˜ B ( ˜ Q, ˜ S ) + 1 + ˜ α ( ˜ W , ˜ S ) + ˜ l m ( ˜ N ) (cid:105) . References
Journal of TheRoyal Society Interface , 14(136):20170490, November 2017.[4] K. Alhazzani, A. Alaseem, M. Algahtani, S. Dhandayuthapani, T. Venkatesan, and A. Rathinavelu. An-giogenesis in cancer treatment: 60 years swing between promising trials and disappointing tribulations.
Anti-Angiogenesis Drug Discovery and Development: Volume 4 , 4:34, 2019.[5] N. Bellomo.
Modeling Complex Living Systems . Birkh¨auser Boston, 2008.[6] N. Bellomo, A. Bellouquid, J. Nieto, and J. Soler. Complexity and mathematical tools toward the modellingof multicellular growing systems.
Mathematical and Computer Modelling , 51(5-6):441–451, 2010.[7] M.E. Berens and A. Giese. “...those left behind.” Biology and oncology of invasive glioma cells.
Neoplasia ,1(3):208–219, 1999.[8] F. Billy, B. Ribba, O. Saut, H. Morre-Trouilhet, T. Colin, D. Bresch, J.-P. Boissel, E. Grenier, and J.-P. Flandrois. A pharmacologically based multiscale mathematical model of angiogenesis and its use ininvestigating the efficacy of a new cancer treatment strategy.
Journal of Theoretical Biology , 260(4):545–562, 2009.[9] K. B¨ottger, H. Hatzikirou, A. Chauviere, and A. Deutsch. Investigation of the migration/proliferationdichotomy and its impact on avascular glioma invasion.
Mathematical Modelling of Natural Phenomena ,7:105–135, 2012. 2510] D.J. Brat, A.A. Castellano-Sanchez, S.B. Hunter, M. Pecot, C. Cohen, E.H. Hammond, S.N. Devi, B. Kaur,and E.G. Van Meir. Pseudopalisades in glioblastoma are hypoxic, express extracellular matrix proteases,and are formed by an actively migrating cell population.
Cancer Research , 64(3):920–927, 2004.[11] D.J. Brat and E.G. Van Meir. Vaso-occlusive and prothrombotic mechanisms associated with tumorhypoxia, necrosis, and accelerated growth in glioblastoma.
Laboratory Investigation , 84(4):397–405, 2004.[12] S. Chouaib, Y. Messai, S. Couve, B. Escudier, M. Hasmim, and M.Z. Noman. Hypoxia promotes tumorgrowth in linking angiogenesis to immune escape.
Frontiers in Immunology , 3, 2012.[13] M.C. Colombo, C. Giverso, E. Faggiano, C. Boffano, F. Acerbi, and P. Ciarletta. Towards the personalizedtreatment of glioblastoma: Integrating patient-specific clinical data in a continuous mechanical model.
PLoS ONE , 10(7):e0132887, 2015.[14] M. Conte, L. Gerardo-Giorda, and M. Groppi. Glioma invasion and its interplay with nervous tissue andtherapy: A multiscale model.
Journal of Theoretical Biology , 486:110088, 2020.[15] G. Corbin, C. Engwer, A. Klar, J. Nieto, J. Soler, C. Surulescu, and M. Wenske. Modeling glioma invasionwith anisotropy- and hypoxia-triggered motility enhancement: from subcellular dynamics to macroscopicpdes with multiple taxis. arXiv:2006.12322.[16] G. Corbin, A. Hunt, A. Klar, F. Schneider, and C. Surulescu. Higher-order models for glioma invasion:From a two-scale description to effective equations for mass density and momentum.
Mathematical Modelsand Methods in Applied Sciences , 28(09):1771–1800, 2018.[17] A. Czirok. Endothelial cell motility, coordination and pattern formation during vasculogenesis.
WileyInterdisciplinary Reviews: Systems Biology and Medicine , 5(5):587–602, 2013.[18] W. Diao, X. Tong, C. Yang, F. Zhang, C. Bao, H. Chen, L. Liu, M. Li, F. Ye, Q. Fan, et al. Behaviors ofglioblastoma cells in in vitro microenvironments.
Scientific reports , 9(1):1–9, 2019.[19] C. Engwer, T. Hillen, M. Knappitsch, and C. Surulescu. Glioma follow white matter tracts: a multiscaleDTI-based model.
Journal of Mathematical Biology , 71(3):551–582, 2014.[20] C. Engwer, A. Hunt, and C. Surulescu. Effective equations for anisotropic glioma spread with proliferation:a multiscale approach.
IMA Journal of Mathematical Medicine and Biology , 33:435–459, 2016.[21] C. Engwer, M. Knappitsch, and C. Surulescu. A multiscale model for glioma spread including cell-tissueinteractions and proliferation.
Mathematical Biosciences and Engineering , 13(2):443–460, 2016.[22] C. Engwer, C. Stinner, and C. Surulescu. On a structured multiscale model for acid-mediated tumorinvasion: The effects of adhesion and proliferation.
Mathematical Models and Methods in Applied Sciences ,27:1355–1390, 2017.[23] R.A. Gatenby and E.T. Gawlinski. A reaction-diffusion model of cancer invasion.
Cancer Research ,56(24):5745–5753, 1996.[24] R.A. Gatenby, E.T. Gawlinski, A.F. Gmitro, B. Kaylor, and R.J. Gillies. Acid-mediated tumor invasion:a multidisciplinary study.
Cancer Research , 66(10):5216–5223, 2006.[25] P. Gerlee and S. Nelander. The impact of phenotypic switching on glioblastoma growth and invasion.
PLoSComputational Biology , 8(6):e1002556, 2012.[26] A. Giese, R. Bjerkvig, M.E. Berens, and M. Westphal. Cost of migration: Invasion of malignant gliomasand implications for treatment.
Journal of Clinical Oncology , 21(8):1624–1636, 2003.[27] A. Giese, L. Kluwe, Meissner H., Michael E., and M. Westphal. Migration of human glioma cells on myelin.
Neurosurgery , 38:755–764, 1996.[28] A. Giese, M.A. Loo, N. Tran, D. Haskett, S.W. Coons, and M.E. Berens. Dichotomy of astrocytomamigration and proliferation.
International Journal of Cancer , 67(2):275–282, 1996.[29] M.A. Hammoud, R. Sawaya, W. Shi, P.F. Thall, and N.E. Leeds. Prognostic significance of preoperativeMRI scans in glioblastoma multiforme.
Journal of Neuro-Oncology , 27(1):65–73, 1996.[30] D. Hanahan and J. Folkman. Patterns and emerging mechanisms of the angiogenic switch during tumori-genesis.
Cell , 86(3):353–364, 1996. 2631] D. Hanahan and R.A. Weinberg. Hallmarks of cancer: the next generation.
Cell , 144(5):646–674, 2011.[32] M.E. Hardee and D. Zagzag. Mechanisms of glioma-associated neovascularization.
The American Journalof Pathology , 181(4):1126–1141, 2012.[33] H. Hatzikirou, D. Basanta, M. Simon, K. Schaller, and A. Deutsch. Go or grow: the key to the emergenceof invasion in tumour progression?
Mathematical Medicine and Biology , 29(1):49–65, 2010.[34] M.A. Hayat. Introduction. In M.A. Hayat (ed.) , Tumors of the Central Nervous System, Volume 1 , pages3–8. Springer Netherlands, 2011.[35] M.G. Vander Heiden, L.C. Cantley, and C.B. Thompson. Understanding the Warburg effect: The metabolicrequirements of cell proliferation.
Science , 324(5930):1029–1033, 2009.[36] C. Henker, T. Kriesen, ¨A. Glass, B. Schneider, and J. Piek. Volumetric quantification of glioblastoma:experiences with different measurement techniques and impact on survival.
Journal of Neuro-Oncology ,135(2):391–402, 2017.[37] P. Hinow, P. Gerlee, L.J. McCawley, V. Quaranta, M. Ciobanu, S. Wang, J.M. Graham, B.P. Ayati,J. Claridge, K.R. Swanson, M. Loveless, and A.R.A. Anderson. A spatial model of tumor-host interaction:application of chemotherapy.
Mathematical Biosciences and Engineering , 6(3):521–546, 2009.[38] S.A. Hiremath and C. Surulescu. A stochastic multiscale model for acid mediated cancer invasion.
NonlinearAnalysis: Real World Applications , 22:176–205, 2015.[39] S.A. Hiremath and C. Surulescu. A stochastic model featuring acid-induced gaps during tumor progression.
Nonlinearity , 29(3):851–914, 2016.[40] S.A. Hiremath and C. Surulescu. Mathematical models for acid-mediated tumor invasion: from determin-istic to stochastic approaches. In
Multiscale models in mechano and tumor biology , volume 122 of
LectureNotes in Computational Science and Engineering , pages 45–71. Springer, Cham, 2017.[41] S.A. Hiremath, C. Surulescu, A. Zhigun, and S. Sonner. On a coupled SDE-PDE system modeling acid-mediated tumor invasion.
Discrete & Continuous Dynamical Systems - Series B , 23(6):2339–2369, 2018.[42] C. Hogea, C. Davatzikos, and G. Biros. An image-driven parameter estimation problem for a reac-tion–diffusion glioma growth model with mass effects.
Journal of Mathematical Biology , 56(6):793–825,2007.[43] P. Holzer. Acid-sensitive ion channels and receptors. In
Sensory Nerves , pages 283–332. Springer BerlinHeidelberg, 2009.[44] E. H¨oring, P.N. Harter, J. Seznec, J. Schittenhelm, H.-J. B¨uhring, S. Bhattacharyya, E. von Hattingen,C. Zachskorn, M. Mittelbronn, and U. Naumann. The “go or grow” potential of gliomas is linked to the neu-ropeptide processing enzyme carboxypeptidase e and mediated by metabolic stress.
Acta Neuropathologica ,124(1):83–97, 2012.[45] A. Hunt and C. Surulescu. A multiscale modeling approach to glioma invasion with therapy.
VietnamJournal of Mathematics , 45(1-2):221–240, 2016.[46] S. Jbabdi, E. Mandonnet, H. Duffau, L. Capelle, K.R. Swanson, M. P´el´egrini-Issac, R. Guillevin, andH. Benali. Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging.
MagneticResonance in Medicine , 54(3):616–624, 2005.[47] C.R. Justus, L. Dong, and L.V. Yang. Acidic tumor microenvironment and pH-sensing g protein-coupledreceptors.
Frontiers in Physiology , 4, 2013.[48] L.J. Kaufman, C.P. Brangwynne, K.E. Kasza, E. Filippidi, V.D. Gordon, T.S. Deisboeck, and D.A Weitz.Glioma expansion in collagen i matrices: analyzing collagen concentration-dependent growth and motilitypatterns.
Biophysical journal , 89(1):635–650, 2005.[49] L.D. Ke, Y.-X. Shi, S.-A. Im, X. Chen, and W.K.A. Yung. The relevance of cell proliferation, vascularendothelial growth factor, and basic fibroblast growth factor production to angiogenesis and tumorigenicityin human glioma cell lines.
Clinical Cancer Research , 6(6):2562–2572, 2000.[50] J. Kelkel and C. Surulescu. On some models for cancer cell migration through tissue networks.
MathematicalBiosciences and Engineering , 8(2):575–589, 2011.2751] J. Kelkel and C. Surulescu. A multiscale approach to cell migration in tissue networks.
MathematicalModels and Methods in Applied Sciences , 22(03):1150017, March 2012.[52] Y. Kim, S. Lawler, M.O. Nowicki, E.A. Chiocca, and A. Friedman. A mathematical model for patternformation of glioma cells outside the tumor spheroid core.
Journal of Theoretical Biology , 260(3):359–371,2009.[53] Y. Kim and S. Roh. A hybrid model for cell proliferation and migration in glioblastoma.
Discrete &Continuous Dynamical Systems - Series B , 18(4):969–1015, 2013.[54] L Ko, A Koestner, and W Wechsler. Characterization of cell cycle and biological parameters of trans-plantable glioma cell lines and clones.
Acta neuropathologica , 51(2):107–111, 1980.[55] N. Kolbe, N. Sfakianakis, C. Stinner, C. Surulescu, and J. Lenz. Modeling multiple taxis: tumor invasionwith phenotypic heterogeneity, haptotaxis, and unilateral interspecies repellence. arXiv:2005.01444v1.[56] E. Konukoglu, O. Clatz, P.-Y. Bondiau, H. Delingette, and N. Ayache. Extrapolating glioma invasionmargin in brain magnetic resonance images: Suggesting new irradiation margins.
Medical Image Analysis ,14(2):111–125, 2010.[57] P. Kumar, J. Li, and C. Surulescu. Multiscale modeling of glioma pseudopalisades: contributions from thetumor microenvironment. arXiv:2007.05297.[58] D.A. Lauffenburger and J.L. Lindermann.
Receptors. Models for binding, trafficing and signaling.
OxfordUniversity Press, 1993.[59] D.R. Lide (ed).
CRC handbook of chemistry and physics, vol. 85 . CRC Press, 2004.[60] T. Lorenz and C. Surulescu. On a class of multiscale cancer cell migration models: Well-posedness in lessregular function spaces.
Mathematical Models and Methods in Applied Sciences , 24(12):2383–2436, 2014.[61] S.Y. Lunt and M.G. Vander Heiden. Aerobic glycolysis: Meeting the metabolic requirements of cell prolif-eration.
Annual Review of Cell and Developmental Biology , 27(1):441–464, 2011. PMID: 21985671.[62] G.R. Martin and R.K. Jain. Noninvasive measurement of interstitial ph profiles in normal and neoplastictissue using fluorescence ratio imaging microscopy.
Cancer Research , 54(21):5670–5674, 1994.[63] N.K. Martin, E.A. Gaffney, R.A. Gatenby, R.J. Gillies, I.F. Robey, and P.K. Maini. A mathematical modelof tumour and blood pHe regulation: The buffering system.
Mathematical Biosciences , 230(1):1–11, 2011.[64] A. Mart´ınez-Gonz´alez, G.F. Calvo, L.A. P´erez Romasanta, and V.M. P´erez-Garc´ıa. Hypoxic cell wavesaround necrotic cores in glioblastoma: A biomathematical model and its therapeutic implications.
Bulletinof Mathematical Biology , 74(12):2875–2896, 2012.[65] R. Mart´ınez-Zaguil´an, E.A. Seftor, R.E.B. Seftor, Y.-W. Chu, R.J. Gillies, and M.J.C. Hendrix. Acidic pHenhances the invasive behavior of human melanoma cells.
Clinical & Experimental Metastasis , 14(2):176–186, 1996.[66] J.B. McGillen, E.A. Gaffney, N.K. Martin, and P.K. Maini. A general reaction–diffusion model of acidityin cancer invasion.
Journal of Mathematical Biology , 68(5):1199–1224, 2013.[67] G. Meral, C. Stinner, and C. Surulescu. A multiscale model for acid-mediated tumor invasion: Therapyapproaches.
Journal of Coupled Systems and Multiscale Dynamics , 3(2):135–142, 2015.[68] R.E. Moellering, K.C. Black, C. Krishnamurty, B.K. Baggett, P. Stafford, M. Rain, R.A. Gatenby, andR.J. Gillies. Acid treatment of melanoma cells selects for invasive phenotypes.
Clinical & ExperimentalMetastasis , 25(4):411–425, 2008.[69] J.D. Murray.
Mathematical Biology . Springer Berlin Heidelberg, 1989.[70] H.G. Othmer and T. Hillen. The diffusion limit of transport equations II: Chemotaxis equations.
SIAMJournal on Applied Mathematics , 62(4):1222–1250, 2002.[71] K.J. Painter and T. Hillen. Mathematical modelling of glioma growth: the use of diffusion tensor imaging(dti) data to predict the anisotropic pathways of cancer invasion.
Journal of Theoretical Biology , 323:25–39,2013. 2872] G. Pennarun, C. Granotier, L.R. Gauthier, D. Gomez, F. Hoffschir, E. Mandine, J.-F. Riou, J.-L. Mergny,P. Mailliet, and F.D. Boussin. Apoptosis related to telomere instability and cell cycle alterations in humanglioma cells treated by new highly selective g-quadruplex ligands.
Oncogene , 24(18):2917–2928, 2005.[73] B. Perthame, M. Tang, and N. Vauchelet. Derivation of the bacterial run-and-tumble kinetic equation froma model with biochemical pathway.
Journal of Mathematical Biology , 73(5):1161–1178, 2016.[74] K. Pham, A. Chauviere, H. Hatzikirou, X. Li, H.M. Byrne, V. Cristini, and J. Lowengrub. Density-dependent quiescence in glioma invasion: instability in a simple reaction–diffusion model for the migra-tion/proliferation dichotomy.
Journal of Biological Dynamics , 6(sup1):54–71, 2012.[75] F. Pignatti, M. van den Bent, D. Curran, C. Debruyne, R. Sylvester, P. Therasse, D. ´Afra, P. Cornu,M. Bolla, C. Vecht, and A.B.M.F. Karim. Prognostic factors for survival in adult patients with cerebrallow-grade glioma.
Journal of Clinical Oncology , 20(8):2076–2084, 2002.[76] S. Prag, E.A. Lepekhin, K. Kolkova, R. Hartmann-Petersen, A. Kawa, P.S. Walmod, V. Belman, H.C.Gallagher, V. Berezin, E. Bock, and N. Pedersen. Ncam regulates cell motility.
Journal of Cell Science ,115(2):283–292, 2002.[77] A.M. Rojiani and K. Dorovini-Zis. Glomeruloid vascular structures in glioblastoma multiforme: an im-munohistochemical and ultrastructural study.
Journal of neurosurgery , 85(6):1078–1084, 1996.[78] G.L. Semenza. Defining the role of hypoxia-inducible factor 1 in cancer biology and therapeutics.
Oncogene ,29(5):625–634, 2009.[79] M. Sidani, D. Wessels, G. Mouneimne, M. Ghosh, S. Goswami, C. Sarmiento, W. Wang, S. Kuhl, M. El-Sibai, J-M. Backer, R. Eddy, D. Soll, and J. Condeelis. Cofilin determines the migration behavior andturning frequency of metastatic cancer cells.
Journal of Cell Biology , 179(4):777–791, 2007.[80] K. Smallbone, R.A. Gatenby, and P.K. Maini. Mathematical modelling of tumour acidity.
Journal ofTheoretical Biology , 255(1):106–112, 2008.[81] C. Stinner, C. Surulescu, and G. Meral. A multiscale model for pH-tactic invasion with time-varyingcarrying capacities.
IMA Journal of Applied Mathematics , 80:1300–1321, 2015.[82] C. Stinner, C. Surulescu, and A. Uatay. Global existence for a go-or-grow multiscale model for tumorinvasion with therapy.
Mathematical Models and Methods in Applied Sciences , 26:2163–2201, 2016.[83] C.L. Stokes and D.A. Lauffenburger. Analysis of the roles of microvessel endothelial cell random motilityand chemotaxis in angiogenesis.
Journal of Theoretical Biology , 152(3):377–403, 1991.[84] C.L. Stokes, M.A. Rupnick, S.K. Williams, and D.A. Lauffenburger. Chemotaxis of human microvesselendothelial cells in response to acidic fibroblast growth factor.
Laboratory investigation; a journal oftechnical methods and pathology , 63(5):657668, 1990.[85] A. Swan, T. Hillen, J.C. Bowman, and A.D. Murtha. A patient-specific anisotropic diffusion model forbrain tumour spread.
Bulletin of Mathematical Biology , 80(5):1259–1291, 2017.[86] K.R. Swanson, R.C. Rockne, J. Claridge, M.A. Chaplain, E.C. Alvord, and A.R.A. Anderson. Quantifyingthe role of angiogenesis in malignant progression of gliomas: In silico modeling integrates imaging andhistology.
Cancer Research , 71(24):7366–7375, 2011.[87] A. Szab, R. nnep, E. Mhes, W. O. Twal, W. S. Argraves, Y. Cao, and A Czirk. Collective cell motion inendothelial monolayers.
Physical Biology , 7(4):046007, 2010.[88] P. Vaupel, F. Kallinowski, and P. Okunieff. Blood flow, oxygen and nutrient supply, and metabolic mi-croenvironment of human tumors: A review.
Cancer Research , 49(23):6449–6465, 1989.[89] Y.-L. Wang, J. Yao, A. Chakhoyan, C. Raymond, N. Salamon, L.M. Liau, P.L. Nghiemphu, A. Lai, W.B.Pope, N. Nguyen, M. Ji, T.F. Cloughesy, and B.M. Ellingson. Association between tumor acidity and hyper-vascularity in human gliomas using pH-weighted amine chemical exchange saturation transfer echo-planarimaging and dynamic susceptibility contrast perfusion MRI at 3t.
American Journal of Neuroradiology ,40(6):979–986, 2019.[90] B.A. Webb, M. Chimenti, M.P. Jacobson, and D.L. Barber. Dysregulated pH: a perfect storm for cancerprogression.
Nature Reviews Cancer , 11(9):671–677, 2011.2991] K.T. Weiß, M. Fante, G. K¨ohl, J. Schreml, F. Haubner, M. Kreutz, S. Haverkampf, M. Berneburg, andS. Schreml. Proton-sensing g protein-coupled receptors as regulators of cell proliferation and migrationduring tumor growth and wound healing.
Experimental Dermatology , 26(2):127–132, 2017.[92] F.J. Wippold, M. L¨ammle, F. Anatelli, J. Lennerz, and A. Perry. Neuropathology for the neuroradiologist:Palisades and pseudopalisades.
American Journal of Neuroradiology , 27(10):2037–2041, 2006.[93] M. Wrensch, Y. Minn, T. Chew, M. Bondy, and M.S. Berger. Epidemiology of primary brain tumors:Current concepts and review of the literature.
Neuro-Oncology , 4(4):278–299, 2002.[94] Q. Xie, S. Mittal, and M. E. Berens. Targeting adaptive glioblastoma: an overview of proliferation andinvasion.
Neuro-Oncology , 16(12):1575–1584, 2014.[95] L. Xu, D. Fukumura, and R.K. Jain. Acidic extracellular pH induces vascular endothelial growth fac-tor (VEGF) in human glioblastoma cells via ERK1/2 MAPK signaling pathway.
Journal of BiologicalChemistry , 277(13):11368–11374, 2001.[96] A. Zhigun, C. Surulescu, and A. Hunt. A strongly degenerate diffusion-haptotaxis model of tumour invasionunder the go-or-grow dichotomy hypothesis.