Multiscale modeling of glioma pseudopalisades: contributions from the tumor microenvironment

Gliomas are primary brain tumors with a high invasive potential and infiltrative spread. Among them, glioblastoma multiforme (GBM) exhibits microvascular hyperplasia and pronounced necrosis triggered by hypoxia. Histological samples showing garland-like hypercellular structures (so-called pseudopalisades) centered around the occlusion site of a capillary are typical for GBM and hint on poor prognosis of patient survival. We propose a multiscale modeling approach in the kinetic theory of active particles framework and deduce by an upscaling process a reaction-diffusion model with repellent pH-taxis. We prove existence of a unique global bounded classical solution for a version of the obtained macroscopic system and investigate the asymptotic behavior of the solution. Moreover, we study two different types of scaling and compare the behavior of the obtained macroscopic PDEs by way of simulations. These show that patterns (not necessarily of Turing type), including pseudopalisades, can be formed for some parameter ranges, in accordance with the tumor grade. This is true when the PDEs are obtained via parabolic scaling (undirected tissue), while no such patterns are observed for the PDEs arising by a hyperbolic limit (directed tissue). This suggests that brain tissue might be undirected - at least as far as glioma migration is concerned. We also investigate two different ways of including cell level descriptions of response to hypoxia and the way they are related .


Introduction
Classified as grade IV astrocytoma by WHO [38], glioblastoma multiforme (GBM) is considered to be the most aggressive type of glioma, with a median overall survival time of 60 weeks, in spite of state-of-the-art treatment [51,62].It is characterized by fast, infiltrative spread and unchecked cell proliferation which triggers hypoxia and upregulation of glycolysis, usually accompanied locally by exuberant angiogenesis [7,20]; one of the typical features of GBM is the development of a necrotic core [38].Increased extracellular pressure from edema and expression of procoagulant factors putatively lead to vasoocclusion and thrombosis [10], hence impairing oxygen supply at the affected site, which becomes hypoxic and induces tissue necrotization.As a consequence, glioma cells actively and radially migrate away from the acidic area [8], forming palisade-like structures exhibiting arrangements of elongated nuclei stacked in rows at the periphery of the hypocellular region around the occlusion site [61].Such histopathological patterns are typically observed in GBM and are used as an indicator of tumor aggressiveness [9,33].Pseudopalisades can be narrow, with a width less than 100 µm and a fibrillar interior structure, medium-sized (200 − 400 µm wide) with central necrosis and vacuolization, but with a fibrillar zone in the immediate interior proximity of the hypercellular garland-like formation.Finally, the largest ones exceed 500 µm in width and are surrounding extensive necrotic areas, most often containing central vessels [8].
Mathematical modeling has become a useful means for supporting the investigation of glioma dynamics in interaction with the tumor microenvironment.Over the years several modeling approaches have been proposed.While discrete and hybrid models (see e.g.[5,32,52]) use computing power to assess the rather detailed interplay between glioma cells and their surroundings, the continuum settings enable less expensive simulations and mathematical analysis of the resulting systems of differential equations.Since the structure of brain tissue with its patient-specific anisotropy is (among other factors) essential for the irregular spread of glioma, the mathematical models should be able to include in an appropriate way such information, which is available from diffusion tensor imaging (DTI) data.The macroscopic evolution of a tumor is actually determined by processes taking places on lower scales, thus it is important to deduce the corresponding population dynamics from descriptions of cell behavior on mesoscopic or even subcellular levels, thereby taking into account the interactions with the underlying anisotropic tissue and possibly further biochemical and/or biophysical traits of the extracellular space.This has been done e.g. in [15, 17-19, 29, 45] upon starting on the mesoscale from the kinetic theory of active particles (KTAP) framework [4] and obtaining with an adequate upscaling the macroscopic PDEs of reaction-(myopic) diffusion(-taxis) type for the tumor dynamics.Depending on whether the modeling process included subcellular events, these PDEs contain in their coefficients information from that modeling level, thus receiving a multiscale character.It is this approach that we plan to follow here, however with the aim of obtaining cell population descriptions for the pseudopalisade formation rather than for the behavior of the whole tumor.
Mathematical models addressing glioma pseudopalisade formation are scarce; we refer to [11,12] for agent-based approaches and to [1,41] for continuous settings.Of the latter, [1] investigated the impact of blood vessel collapse on glioma invasion and the phenotypic switch in the migration/proliferation dichotomy.It involves a system of PDEs coupling the nonnlinear dynamics of glioma population with that of nutrient concentration and vasculature, thus not explicitly including acidity.The PDE for the evolution of tumor cell density was obtained upon starting from a PDE/ODE system for migrating/proliferating glioma densities and performing transformations relying on several assumptions.The work in [41] describes interactions between normoxic/hypoxic glioma, necrotic tissue, and oxygen concentration.The model confirms the histological pattern behavior and shows by simulations a traveling wave concentrically moving away from the highly hypoxic site toward less acidic areas.The PDE system therein was set up in a heuristic manner directly on the macroscale; it features reaction-diffusion equations without taxis or other drift.
In the present work we are interested in deducing effective equations for the space-time evolution of glioma cell density in interaction with extracellular acidity (concentration of protons), thereby accounting for the multiscality of the involved processes and for the anisotropy of brain tissue.The deduced model should be able to reproduce pseudopalisade-like patterns and to investigate the influence of acidity and tissue on their behavior.The rest of this paper is organized as follows: in Section 2 we formulate our model upon starting from descriptions of cell dynamics on the microscopic and mesoscopic scales.Section 3 is concerned with obtaining the macroscopic limits of that setting; we will investigate parabolic as well as hyperbolic upscalings, corresondingly leading to diffusion and drift-dominated evolution, respectively, and depending on tissue properties (directed/undirected).An alternative modeling approach and its parabolic limit will be addressed as well.In Section 4 we provide an assessment of parameters and functions involved in the deduced macroscopic PDE systems and perform numerical simulations, also providing a comparison between the studied modeling approaches.Section 5 is dedicated to establishing the existence and uniqueness of a global bounded classical solution to the macroscopic system obtained by parabolic scaling.A result concerning the asymptotic behavior of such solution is proved as well.Finally, Section 6 contains a discussion of the obtained results and an outlook on further problems of interest related to GBM pseudopalisades.

Model set up on subcellular and mesoscopic scales
The approach in [17,18] led to (hapto)taxis of glioma cell population on the macroscale upon taking into account receptor binding of cells to the surrounding tissue.As such, it was a simplification of subcellular dynamics as considered in [30,31,37], where the cancer cells were supposed to interact with the tissue and with a soluble ligand acting as a chemoattractant.The latter works, however, were concerned with the micro-meso-macro formulation of cancer cell evolution in dynamic interaction with the tissue (as a mesoscopic quantity) and the ligand (obeying a nonlocal macroscopic PDE), along with the analysis therewith, whereas here we intend to obtain a system of effective macroscopic PDEs for glioma population density in interaction with space-time dependent acidity.Here the macroscopic scale is smaller than in the mentioned previous works: it is not the scale of the whole tumor, but that of a subpopulation, localized around one or several vasoocclusion sites in a comparatively small area of the tumor -corresponding to a histological sample.Since (see end of first paragraph in Section 1) the size of such samples is too small to allow a reliable assessment of underlying tissue distribution via DTI 2 , we do not describe a detailed cell-tissue interaction via cell activity variables as in the mentioned works.However, tissue anisotropy might be relevant even on such lower scale, therefore we consider, instead, an artificial structure by way of some given water diffusion tensor, in order to be able to test such influences on the glioma pattern formation.Therefore, on the subcellular level we only account explicitly for interactions between extracellular acidity and glioma transmembrane units mediating them.The latter can be ion channels and membrane transporters ensuring proton exchange, or even proton-sensing receptors [28].
We denote by y(t) the amount of transmembrane units occupied with protons (in the following we will call this the activity variable, in line with the KTAP framework in [4]) and by R 0 the total amount of such units (ion channels, receptors, etc), which for simplicity we assume to be constant.Let S denote the concentration of (extracellular) protons and S max be a threshold value, which, when exceeded, leads to cancer cell death.The corresponding binding/occupying kinetics are written 2 the typical size of a voxel is ca. 1 mm3 so that we can write for the corresponding subcellular dynamics (upon rescaling y y/R 0 ) where k + and k − represent the reaction rates.We denote by y * the steady-state of the above ODE, thus we have As in [17,18] we will consider deviations from the equilibrium of subcellular dynamics: Since the events on this scale are much faster than those on the mesoscopic and especially macroscopic levels, the equilibrium is supposed to be quickly attained, so z is very small.We will use this assumption in the subsequent calculations; as in [17,18] it will allow us to get rid of higher order moments during the upscaling process, thus to close the system of moments leading to the macroscopic formulation.
This assumption also allows us to ignore on this microscopic scale the time dependency of S. 3 Next, we consider the path of a single cell starting at position x 0 and moving with velocity v in the acidic environment.Since the glioma cells are supposed to move away from the highly hypoxic site, we take: We denote by p(t, x, v, y) the density function of glioma cells at time t , position x ∈ R n , velocity v ∈ V ⊂ R n , and with activity variable y ∈ Y = (0, 1).We assume as in [15, 17-19, 29, 45] that the cells have a constant (average) speed s > 0, so that V = sS n−1 , i.e.only the cell orientation is varying on the unit sphere.In terms of deviations z ∈ Z ⊂ [y * − 1, y * ] from the steady-state (we also call z activity variable) we consider for the evolution of p the kinetic transport equation (KTE) where )dv denotes the turning operator modeling cell velocity innovations due to tissue contact guidance and acidity sensing, with λ (z) denoting the turning rate of cells.Thereby, K(x, v, v ) is a turning kernel giving the likelihood of a cell with velocity v to change its velocity regime into v.We adopt the choice proposed in [23], i.e.K(x, v, v ) = q(x,v) ω where q(x, v) is the (stationary) orientation distribution of tissue fibers with ω = V q(v)dv = s n−1 and v = v |v| ∈ S n−1 .We take the turning rate as in [17] where λ 0 and λ 1 are positive constants.The choice means that the turning rate is increasing with the amount of proton-occupied transmembrane units.The turning operator in (2.4) thus becomes with We also employ the notation f (S) = y * to emphasize that the steady-state of subcellular dynamics depends on the proton concentration S. The last term in (2.4) represents growth or depletion, according to the acidity level in the tumor microenvironment.Similarly to [18], but now accounting for the effect of acidity, we consider a source term of the form where χ(z, z ) represents the likelihood of cells having activity state z to go into activity state z upon interacting with acidity S(t, x).In particular, χ is a kernel with respect to z, i.e.Z χ(z, z )dz = 1.The acidity is reported again to the threshold value S max .The growth rate µ(M) depends on the total amount M(t, x) = V Z p(t, x, v, z)dzdv of glioma cells, irrespective of their orientation or activity state and takes into account limitations by overcrowding.We will provide a concrete choice later in Subsection 4.1.Hence, the presence of tissue is supporting proliferation, which is maintained until the environment becomes too acidic even for tumor cells.
The above micro-meso formulation for glioma dynamics is supplemented with the evolution of acidity, described by the macroscopic PDE where D s is the diffusion coefficient of protons, β is the proton production rate by tumor cells, and α denotes the rate of acidity decay.The high dimensionality of the above setting makes the numerics too expensive, thus we aim to deduce macroscopic equations which can be solved more efficiently and, moreover, facilitate the observation of the glioma cell population and its patterning behavior.In order to investigate the possible effects of the tissue being directed or not 4 , we will perform two kinds of macroscopic limit: the parabolic one, for the diffusion-dominated case of undirected tissue, and the hyperbolic limit for directed tissue, which should be drift-dominated.Both types of limits are performed in a formal way, as the rigorous processes would require analytical challenges which go beyond the aims of this note.

Macroscopic limits
We consider the following moments with respect to v and z: and neglect higher order moments w.r.t.z due to the assumption of the steady-state of subcellular dynamics being rapidly reached.Moreover, we assume p to be compactly supported in the phase space R n ×V × Z.
Integrating (2.4) w.r.t z, we get: Multiplying (2.4) by z and integrating w.r.t.z we get: In the following we denote as e.g., in [17,23] by the mean fiber orientation and the variance-covariance matrix for the orientation distribution of tissue fibers, respectively.

Parabolic limit
In this subsection we consider the tissue to be undirected, which translates into the directional distribution function for tissue fibers being symmetric, i.e. S n−1 q(x, θ )dθ = S n−1 q(x, −θ )dθ .We rescale the time and space variables by t := ε 2 t, x := εx.Since proliferation is much slower than migration, we also rescale with ε 2 the corresponding term, as in [18].For notation simplification we will drop in the following the˜symbol from the scaled variables t and x.
Thus, from (3.1) and (3.2) we get: Now, using Hilbert expansions for the moments: and identifying the equal powers of ε, we get ε 0 : If we also expand µ around M 0 , (3.9) leads to Then from (3.5) we obtain The assumption of undirected tissue gives E q = 0, thus from the above equation we obtain M z 1 = 0, which in virtue of (3.8) implies .
We summarize our hitherto information about the moments: 14), the previous equation becomes where: This macroscopic PDE forms together with (2.9) the system characterizing glioma evolution under the influence of acidity.It involves a term describing repellent pH-taxis (the glioma cells move away from large acidity gradients), in which the tactic sensitivity function contains the tumor diffusion tensor D T encoding information about the anisotropy of underlying tissue and the function g(S) which relates to the subcellular dynamics of proton sensing and transfer accros cell membranes.The myopic diffusion is common to this and previous models [17,18,29,45] obtained by parabolic scaling from the KTAP framework.

An alternative approach to including acidity effects
In the previous derivation the term with repellent pH-taxis was obtained as a consequence of including subcellular level dynamics in the mesoscopic KTE (2.4) by way of the transport term w.r.t. the activity variable z and by letting the turning rate depend on it.In the context of bacteria motion an alternative approach was proposed in [44] and re-employed in [39] also for eukaryotes having a more complex motility behavior.It does not explicitly include subcellular dynamics (thus no activity variables and corresponding transport terms are considered), but lets instead the cell turning rate depend on the pathwise gradient of some chemoattractant concentration which is supposed to bias the cell motion.
The relationship between the two mesoscopic modeling approaches was studied for bacteria dispersal in [47], where it was rigorously shown that the alternative approach follows from the former one, under certain assumptions made on the receptor binding dynamics on the subcellular level (along with fast relaxation towards equilibrium of external signal transduction and stiff response of the activity variables), on the turning rate, and on the initial data.
Here we intend to investigate two such approaches for the problem at hand (glioma cells moving away from acidity) from a less rigorous perspective, namely looking into their (formal) macroscopic limits and comparing the numerical results obtained therewith.Concretely, we want to compare our KTE (2.4) and its parabolic limit with the following simpler KTE for the cell density function ρ(t, x, v): and its parabolic limit.Here the proliferation operator is defined with the same µ(M, S) as previously used in this section and takes the form while for the turning rate we set 5 λ (v, S) := λ 0 exp (h(S)D t S) (3.18) where D t S := S t + v • ∇ x S is the pathwise gradient of S. The coefficient function h(S) is to be chosen later.
We use again a parabolic scaling t := ε 2 t, x := εx and rescale as before the proliferation term by ε 2 , thus Performing a Hilbert expansion ρ = ρ 0 + ερ 1 + ε 2 ρ 2 + . . .and equating the powers of ε yields thus by (3.21) and the assumption of undirected tissue This can be rewritten as Since the integral of the right hand side w.r.t.v vanishes we can pseudo-invert L[λ 0 ] as before, to get ε 2 : 5 notice the difference of sign in the exponent: the cells are supposed to follow the decreasing gradient of signal S Using (3.24) leads to which differs from (3.15) only by the tactic sensitivity coefficient.This suggests to choose h(S) = g(S) Notice that this function corresponds to the rate of change of f (S) (equilibrium of transmembrane interactions between glioma cells and acidity), scaled by a constant λ 1 to account for changes in the turning rate per unit of change in dy * /dt, and multiplied by 1/(k + S/S max + k − + λ 0 ).While a factor const • f (S) is also encountered in [44], the denominator obtained here for the tactic sensitivity appears due to the specific choice of the turning rate λ (z) in (2.5), which also facilitated the upscaling.The whole sensitivity function is tightly related to the first order moment w.r.t. the receptor binding ('activity') variable z, recall (3.12) and (3.14).

Hyperbolic scaling
In this subsection we investigate the macroscopic limit of (2.4) in the case where the tissue is directed.
In particular, this means that the mean fiber orientation E q is nonzero, as the orientation distribution q is unsymmetric.We will only treat the case explicitly accounting for the dynamics of activity variables.
where V p ⊥ (t, x, v, z)dv = 0 and p(t, x, z) := V p(t, x, v, z)dv.Then for the moments introduced at the beginning of this Section 3 we have (3.28) Now we rescale the time and space variables by t := εt, x := εx and drop again the˜symbol to simplify the notation.As before, the proliferation term is scaled by ε 2 .With these, the equations (3.1) and (3.2) become, respectively: and Using (3.28) we write Since q is independent of time, these equations imply where we used the notation Ẽq (x) := V v q ω (x, v)dv = sE q .From (3.34) we get at leading order Plugging this in (3.34), we obtain (again at leading order) From (3.35): Plugging this into (3.33)we get (at leading order) Since the right hand side vanishes when integrated w.r.t.v, we can pseudo-invert L[λ 0 ] and use (3.36) to get so that (3.35) becomes Comparing this with the parabolic limit obtained in (3.15) we observe that we obtain the same form for the (myopic) diffusion, repellent pH-taxis, and proliferation terms, but here they are ε-corrections of the leading transport terms -together with the new advection which drives cells with velocity ε λ 0 E q ∇ • E q in the direction of the dominating drift.

Parameters and coefficient functions
We assume a logistic type growth of the tumor cells and choose where µ 0 is the growth rate and M max is the carrying capacity of tumor cells.

Nondimensionalization
Considering the following nondimensional quantities: and dropping the tildes for simplicity of notation, we get the following nondimensionalized system: 2)

Description of tissue
The structure of brain tissue can be assessed by way of biomedical imaging, e.g.diffusion tensor imaging (DTI) which provides for each voxel the water diffusion tensor D w .The corresponding resolution is, however, too low and does not deliver information about the (orientation) distribution of tissue fibers below the size of a voxel (ca. 1 mm 3 ).For more details we refer e.g. to [17,45] and references therein.
As mentioned in Section 1, pseudopalisades are comparatively small structures with a medium width of 200 − 400 µm.Thus, in order to investigate the possible effect of (local) tissue anisotropy on these patterns we will create a synthetic DTI data set which will allow to compute the tumor diffusion tensor D T in the space points of such a narrow region.To this aim we proceed as in [45] and consider the water diffusion tensor where d(x, y) = 0.25e −0.005(x−450) 2 − 0.25e −0.005(y−450) 2 .For the fiber distribution function, we consider a mixture between uniform and von Mises-Fisher distributions, as follows: Here, δ ∈ [0, 1] is a weighting coefficient, ϕ 1 is the eigenvector corresponding to the leading eigenvalue of D w (x) and I 0 is the modified Bessel function of first kind of order 0. Also, θ = (cos ξ , sin ξ ) for ξ ∈ [0, 2π], and k(x) = κFA(x), where FA(x) denotes the fractional anisotropy: in 2D it has the form [45] FA with λ i (i = 1, 2) denoting the eigenvalues of D w (x).The parameter κ ≥ 0 characterizes the sensitivity of cells towards orientation of tissue fibers.For perfectly aligned tissue (i.e., maximum anisotropy), FA(x) = 1 and k(x) = κ.Taking κ = 0 means, however, that the cells are insensitive to even such alignment and the distribution in 4.5 becomes a uniform one.Taking δ = 1 has the same effect.
For the model deduced by hyperbolic scaling in Subsection 3.2, we consider for the orientation distribution of tissue fibers the following combination of two unsymmetric unimodal von Mises distributions: where T and the rest of parameters are the same as in (4.5).The first summand, similar to the choice in [24], generates an orientation along the diagonal γ, while the second leads to alignment along the positive x and y directions.Due to k h (x), the strength of diagonal orientation of tissues decreases from the chosen center (450, 450).
The macroscopic tissue density Q is obtained in the same way as in [18] by using the free path length from the diffusivity obtained from the data, more precisely from the water diffusion tensor.In that approach the occupied volume is obtained upon computing a characteristic (diffusion) length l c = tr(D w )t c , where t c is the characteristic (diffusion) time.The latter is determined by assuming the underlying stochastic process behind water diffusion tensor measurements to be the Brownian motion and considering the expected exit time from the minimal ball with radius r containing a square with side length h as smallest unit in our grid.Therefore, the tissue density Q (area fraction occupied by tissue) is: where with l 1 being the largest eigenvalue of D w .

Numerical experiments
The system (4.2),(4.3) is solved numerically on a square domain [0, 1000] × [0, 1000] (in µm) using appropriate finite difference methods for spatial discretization and an IMEX method for time discretization, where the diffusion part is handled implicitly, while the advection and source terms are treated explicitly.We use a standard central difference scheme (5-point stencil) for the acidity diffusion.
To avoid numerical instability [43] due to negative values in the stencil obtained from discretization of mixed derivative terms in the myopic tumor diffusion, we use the non-negative discretization scheme proposed in [58] instead of the standard one.Thereby, the derivatives are calculated in newly chosen directions (diagonal directions of the 3 × 3-stencil in 2D) in addition to the standard x,y-directions and mixed term derivatives are replaced by directional derivatives.To discretise the advection terms, we use a first order upwind scheme for the parabolic scaling model, while for the system obtained via hyperbolic scaling we employ a second order upwind scheme with Van Leer flux limiter.Implicit and explicit Euler methods are used for IMEX time discretization.The systems are solved with no-flux boundary conditions and the following sets of initial conditions as illustrated in Figures 1(a Experiment 1 -Fully isotropic tissue We begin by considering a fully isotropic tissue, i.e. taking δ = 1 in (4.5).The corresponding fractional anisotropy is everywhere FA = 0, and the macroscopic tissue density Q is shown in Figure 2(a).The simulations show (see Figures 3 and 4) the formation of a pseudopalisade-like pattern, with a very acidic, hypocellular center region surrounded by relatively high glioma cell densities.Thereby, the initial distribution of the tumor cell aggregates and their corresponding pH distribution decisively influence the shape and size of the pattern and the space-time acidity distribution; compare Figures 3 and 4.

Experiment 2 -Anisotropic tissue
With the choice δ = 0.2, κ = 3 we describe an underlying tissue with pronounced anisotropy (two crossing fibre bundles).The corresponding mesoscopic fiber distribution q is shown in Figure 2(b) for a fixed fiber direction, while the macroscopic tissue density Q remains unchanged.The results of this experiment are shown in Figures 5 and 6.The simulated patterns have similar shapes with those in Experiment 1, but here the tissue anisotropy determines the cells to follow the main orientation of the fiber bundles, which leads to a longer persistence of (small amounts of) cells in the central region with more localized cell aggregates exhibiting higher maxima (see Subfigures 5(b) and 6(b)).The patterns at later times (see Subfigures 5(c) and 6(c)) still bear traits of the degraded tissue; the cells are still forming garland-like structures around the hypoxic centers, with the highest cell density located at one or several peripheral sites with highly aligned tissue, farthest away from the main     sources of (initial) acidity.As before, the initial distributions of tumor and acidity influence the shape of the patterns (compare Figures 5 and 6).The differences between the acidity distributions in Subfigures 5(d)-5(f) and those in Subfigures 3(d)-3(f) (and correspondingly Subfigures 6(d)-6(f) and, respectively, those in Subfigures 4(d)-4(f) for the set of initial conditions (4.9)) are less prominent, since the acidity concentration h obeys in both cases a PDE with linear diffusion, where the tissue anisotropy has minor influence.Running numerical simulations for several different parameter sets led to the following observations: • The decisive parameter in this system seems to be α, which relates to the proton buffering efficacy (in the nondimensionalized form (4.2) it is basically the ratio between the acidity removal and proton production rates).The tumor growth rate µ 0 plays a role, too, but a less prominent one.Concretely, pseudopalisade patterns form for very low values of α (weak buffering).If the system is able to remove protons more efficiently (e.g., because there is a functioning capillary network), then these garland-like patterns typical for GBM do not form in a time span which is relevant for this cancer (less than a year); instead, there are rather homogeneous structures with dense cellular areas and no necrosis -which corresponds to a lower tumor grade, without (local) occlusions of capillaries and corresponding necrotization (anaplastic astrocytoma), some with partially preserving the underlying tissue structure (fibrillar astrocytoma), see [49] for WHO-grading on the basis of histopathological samples.Figure 7 shows the evolution of glioma and acidity at several times for the system with initial conditions (4.9) and the same parameters as in the simulations of Experiment 2, with the exception of α, which is now still very small, but four orders of magnitude higher.The tumor cells are producing acidity (by glycolyis) and the inner region begins to degrade, as in the previous simulations.However, due to the stronger acidity removal ratio, it does not become severely hypoxic, which allows the tumor cells to repopulate it, while the rest of the neoplasm is expanding outwards.The underlying tissue structure is thereby supporting both migraton and growth.Notice the more extensive tumor spread in comparison with Figure 6.In 6 we do a short linear stability analysis   of system (4.2),(4.3) with a constant tumor diffusion coefficient; it turns out that no Turing patterns are formed -which does not mean, however, that other types of patterns are not possible.
• If α exceeds a certain threshold value (in our simulations it was one order of magnitude higher than in the computations for Figure 7) then the solution blows up already in 1D.
• The shape of the source term in the equation for tumor cell density has itself a substantial influence on the pattern.It should be chosen in such a way that proliferation is reduced for higher acidity levels.This is, however, not enough for pseudopalisade formation: for instance, a source term of the form µ 0 (1 − M) M 1+S instead of that in (4.2) does not lead to such patterns, as there is no decay of glioma cells due to hypoxia. Figure 8 shows the behavior of tumor and acidity for this alternative choice of the source term, in the framework of Experiment 2.  Including repellent pH-taxis leads to the formation of wider pseudopalisades, with thinner cell 'garlands' than in the case where the glioma cells are only performing myopic diffusion and being degraded by excessive acidity.Figure 9 illustrates the differences between tumor density and acidity in the two cases, i.e. between solutions of (4.2), (4.3) and those of the same system with g(S) = 0.The differences are more pronounced in earlier stages of pattern formation and become smaller with advancing time.The plots also show that pseudopalisades are formed even if there is no pH-taxis, suggesting that the latter merely enhances the effect of the source/decay term in (4.2) who is actually driving the pattern -together with an opportune parameter combination (in particular, adequate proton buffering).
To see the effect of drift dominance we also solve the macroscopic system (3.41),(2.9) obtained by hyperbolic scaling.Thereby we use (where applicable) the same set of parameters and boundary conditions as before for the parabolic scaling (Table 1).For the scaling parameter we take ε = 10 −5 .The initial conditions are those in set (4.9), as visualized in Subfigures 1(c), 1(d).Here we consider an unsymmetric tissue with mesoscopic orientational distribution q h as in (4.6). Figure 10 shows the mean fiber orientation E q corresponding to q h along with a magnification to observe the directionality in the neighborhood of the crossing fiber strands, and with q h plotted for δ = 0.2 and a specific direction instead of that in (4.2).All parameters as in Table 1.The results obtained by solving the system for the evolution of tumor cells and acidity are shown in Figure 11.Although we ran the simulations for a longer time than we did for the system obtained via parabolic scaling no pseudopalisade patterns are formed.Rather, the drift-dominated PDE for glioma cell density drives the cells along the positive x and y directions (as δ = 0.2 makes the second term in (4.6) dominant).The cells 'escaping' that influence move fast along the diagonal γ towards the right upper corner and cannot form the pattern in due time.A quick comparison with Figure 6 obtained for the parabolic limit and q as in (4.5) shows the radically different behavior w.r.t. the two approaches.
Similar observations apply when the first von Mises distribution in (4.6) exerts full influence (for δ = 1).Figure 12 shows tissue characteristics for this case: fractional anisotropy FA, zoomed E q , and q h for ξ = π/2.The very low FA values indicate a highly isotropic tissue.Figure 13 illustrates the behavior of tumor cell density and acidity for this case, in which the glioma cells are migrating very fast along the diagonal γ, accompanied by acidity they produce.When reaching the right uppermost corner of the domain they remain there (due to the no-flux boundary conditions) and further express acidity, eventually both solution components getting depleted.This is again in striking contrast to the solution behavior obtained by parabolic scaling for isotropic and undirected tissue (compare with Figure 4) and, since such evolution is not seen in histologic patterns, it endorses the suspicion of the underlying tissue being undirected, at the same time speaking against hyperbolic scaling.As a casual observation there can be blow-up also in this case, however for a much (three orders of magnitude) stronger proton buffering than in the parabolic case.
We also solved the system (3.41),(2.9) upon using several other initial conditions and parameter sets, none of which led to the formation of pseudopalisades.The observed behavior does not significantly change for any choice of the scaling parameter ε ∈ [10 −6 , 10 −2 ].Thus, since such patterns are actually observed in histologic samples of glioblastoma, the simulations suggest that the fibers of brain tissue do not seem to be directed.This endorses the parabolic upscaling approach and goes along with a diffusion dominated motion, correspondingly biased by acidity gradients.With these, the cells are primarily driven by acidity, but also influenced by the underlying, undirected tissue.The interplay between these actors leads to various types of patterns, depending on the parameter range and the relationship between the parameters.

Main results
We consider system (4.2),(4.3) with a slight modification of the source term in (4.3): with g(S) := , where we use Λ = λ 1 k D , K = k D , and B = λ 0 to denote the corresponding constants occurring in the expression of g(S) as given in Subsection 4.1.Ω ⊂ R N is considered to be a bounded domain with sufficiently smooth boundary ∂ Ω, all involved constants are positive, M 0 ∈ L ∞ (Ω), M 0 ≡ 0, M 0 , S 0 ≥ 0, and S 0 ∈ W 1,∞ (Ω).The no-flux boundary conditions are obtained through the upscaling procedure (as done e.g., in [14] for a related problem).For the tumor diffusion tensor D T we require Assumption 5.1 Theorem 5.1 Let N ≥ 1. Suppose that Assumption 5.1 holds.Then if ζ < α and S 0 L ∞ (Ω) < 1, system (5.1)admits a unique global bounded classical solution.

Global existence, uniqueness, and boundedness of solutions
Firstly, we state a result concerning local existence of classical solutions, which can be proved by well-established methods involving standard parabolic regularity theory and an appropriate fixed point framework.Moreover, one can thereby derive a sufficient condition for extensibility of a given local-in-time solution(see [59] or [13] for example).

PROOF.
Taking pS p−1 (p > 1) as a test function for the S-equation of (5.1), for any ε ∈ (0, 1), we obtain from which we obtain and then by Gronwall's inequality from which we obtain that for any t ∈ [0, T max ), From the arbitrariness of ε ∈ (0, 1) we therefore obtain On the other hand, from the L p -L q estimates for the Neumann heat semigroup on a bounded domain and the fact that we obtain for all t ∈ (0, T max ), where λ 1 > 0 denotes the first nonzero eigenvalue of −∆ in Ω ⊂ R N under the Neumann boundary condition.
PROOF.Taking pM p−1 as a test function for the M-equation of (5.1) and denoting D 0 := D T (•) L ∞ (Ω) , from the no-flux boundary, we obtain where Thus we obtain that for any t ∈ (0, T max ), Proof of Theorem 5.1.From Lemma 5.3 and the standard Moser iteration process, there exists C > 0 such that M(t, •) L ∞ (Ω) ≤ C for all t ∈ (0, T max ).Then in view of Lemma 5.1, Theorem 5.1 is a direct consequence of Lemma 5.2.

Long time behavior
Lemma 5.4 Under the assumptions of Theorem 5.2, there exists µ * > 0 defined in (5.7) such that for µ 0 > µ * , for all t > 0, the function where with D a constant defined in (5.8).

PROOF.
According to the strong maximum principle and the assumption M 0 ≡ 0, M is positive in (0, ∞) × Ω. Testing the M-equation of (5.1) by 1 − 1 M , by Young's inequality and the fact of ∇ • u(x) = 0 for x ∈ Ω, u(x) = 0 for x ∈ ∂ Ω, using the no-flux boundary condition, we obtain that there exists Testing the S-equation of (5.1) by S − ζ 2α , we obtain 1 2 Combining (5.5) and (5.6), we obtain By choosing we obtain that µ 0 > µ * leads to F (t) ≤ −D(t).
Proof of Theorem 5.2.The proof of Theorem 5.2 is very standard.We include the proof here for completeness.Denote h(s) := s − 1 − ln s.Noticing that h (s) = 1 − 1 s and h (s) = 1 s 2 > 0 for all s > 0, we obtain that h(s) ≥ h(1) = 0 and F(t) is nonnegative.From Lemma 5.4, we have F (t) ≤ −D(t) and then t 0 D(τ)dτ ≤ F(0) for all t > 0, from which we have Using a similar argument as in Lemma 3.10 of [56], we can obtain the uniform convergence of solutions, namely Then there exists t 0 > 0 such that for all t > t 0 , M − 1 L ∞ (Ω) ≤ 1 2 , which together with the fact that 1 3 for all t > t 0 .Hence F (t) ≤ −D(t) ≤ −DF(t), from which we obtain F(t) ≤ F(t 0 )e −D(t−t 0 ) .
(5.10) Substituting (5.10) into (5.9),we obtain which implies that there exists C > 0 such that for all t > t 0 , Furthermore, notice that there exists a constant C 1 > 0 such that Thus the Gagliardo-Nirenberg inequality yields for all t > 0. Similarly, we can obtain This concludes the proof of Theorem 5.2.
Remark 5.1 For the above rigorous results to hold we required among others that ζ < α, which means that the acidity buffering by the tumor environment is stronger than the production of protons by the cancer cells.While this is true for lower grade tumors, it no longer holds for more advanced neoplasms like GBM.Numerical simulations show that no pseudopalisades are forming, unless ζ substantially exceeds α.
In fact, in Subsection 4.3 we already observed that α (which controls proton buffering) was the decisive parameter for the fate of the patterns and even for singularity formation.The weakening of proton production considered in this Section enhances the influence of acidity depletion, which due to the form of f (M, S) contributes to keeping the glioma density bounded by its carrying capacity.A similar result can be obtained by replacing in f (M, S) the factor 1 − S with 1/(1 + S).In that case there is no smallness requirement for ζ ; moreover, the results hold even if (4.3) is considered instead of the S-equation in (5.1).
No pseudopalisades are forming in this case either (recall Figure 8).

Discussion
The multiscale approach employed in this work allows to obtain a macroscopic description for the evolution of glioma cell density featuring repellent pH-taxis and providing the concrete forms of involved diffusion, transport, and taxis coefficients, upon starting from modeling on the microscopic level of cell-acidity interactions.This fully continuous setting is quite different from previous models [1,41] of pseudopalisade formation and spread, which are rather accounting for vascularization and necrosis than for direct effects of acidity.Nevertheless, the system of two PDEs of reaction-(myopic) diffusion-advection type obtained by parabolic upscaling from lower levels of description is able to reproduce biologically observed patterns, whereby repellent pH-taxis does not seem to effectively trigger, but merely to enlarge such structures; depending on the acidity buffering potential of the tumor cells and their environment in relationship to their ability to proliferate, the resulting patterns can be assigned to lower or higher tumor grades, with pseudopalisades corresponding to the latter.This endorses the idea that proton buffering might be beneficial for decelerating progression towards GBM, see e.g.[6,35] and references therein.For instance, genetic targeting of carbonic anhydrase 9 (a common hypoxia marker catalyzing the conversion of carbon dioxide to bicarbonate and protons) provided evidence of delayed tumor growth in the GBM cell line U87MG [42].
In our deduction of the macroscopic system from the KTAP framework we used for the turning rate λ (z) = λ 0 − λ 1 z > 0. This could be made more general, e.g.upon considering any regular enough function λ , expanding it around the steady-state y * , and keeping the first two terms of the expansion: λ (z) λ (y * ) − λ (y * )z := λ 0 (S) − λ 1 (S)z.The higher order terms will get anyway lost during the scaling process, due to ignoring the higher order moments w.r.t.z.The new coefficients λ 0 , λ 1 are no longer constants, but depend on the macroscopic variable S by way of y * . 6Consequently, the obtained macroscopic PDE for the glioma population density will have diffusion and taxis coefficients depending on S, thus leading to a more intricate coupling of the PDE system for M and S.
Beside including subcellular level information via a transport term w.r.t. the activity variable(s) and a turning rate depending therewith, we also considered an alternative way to account for cell reorientations in response to acidity levels.Trying to recover the same macroscopic limit led to a well-determined choice of the acidity-dependent function h involved in the turning rate λ (v, S) from (3.18).
For the sake of simplicity we considered in (2.9) a genuinely macroscopic PDE of reaction-diffusion type for the evolution of acidity.More detailed models involving intra-and extracellular proton dynamics with randomness have been introduced in [25][26][27]34], some of them connecting it to the dynamics of tumor cells.The latter inferred, however, a rather heuristic, mainly macroscopic description, with coefficients possibly depending on such microscopic quantities like concentration of intracellular protons.Connecting multiscale formulations of proton and cell dynamics and identifying an appropriate way of upscaling to deduce the corresponding macroscopic equations would be a first step towards accounting for subcellular processes in a manner which is detailed enough to capture such low-scale events, but also eventually simplified enough to still enable efficient computations.
The observation that no pseudopalisades seem to emerge for a transport-dominated system as obtained by hyperbolic scaling of the micro-meso setting suggests that the microscopic brain tissue is undirected, at least w.r.t.glioma migration along its constituent fibers.This is a relevant information for the existing models of glioma invasion built upon ideas commonly employed within the KTAP framework and which take into account the underlying brain structure and its properties in trying to predict the tumor extent and its aggressiveness; we refer to [23] and [14] for two works where such issue is explicitly addressed.On the other hand, this could also be relevant from a biological viewpoint; indeed, to our knowledge such information is not available in the biological literature.We are far from claiming to have a watertight evidence; it is rather a cue to motivate such speculation which should of course be properly verified by appropriately designed biological experiments.
The linear stability analysis performed in Appendix A for constant diffusion coefficients suggests that pseudopalisades are a rather nonstandard type of patterns -at least as far as this model is concerned.The pH-chemorepellence is enhancing the diffusive effect, driving the tumor cells away from the strongly hypoxic site(s).Thereby, the form of the space-dependent tumor diffusion coefficient seems to play a decisive role for the shape of the tumor pattern, as simulations show.The formation of garland-like structures can be observed during the first half of the simulation time, after which there is no 'ring-like closure' of the cell aggregates, although these seem to develop on each side of a hypocellular, acidic region.A rigorous analysis has still to be done, even in the case D T (x) > 0 for all x.
To acquire more qualitative information about the solutions of the macroscopic system deduced via parabolic scaling, we also performed a well-posedness analysis.As the global behavior of solutions to (4.2), (4.3) seems out of reach, we assumed the production of protons by tumor cells to infer saturation and proved that the corresponding system has a unique global bounded nonnegative solution in the classical sense -for which certain assumptions on the tumor diffusion tensor were needed.In the case of solenoidal drift velocity and sufficiently large tumor growth, we proved that the solution approaches asymptotically the steady-state in which the tumor is at its carrying capacity, with a corresponding acidity concentration.The patterning behavior for the system with saturated, but sufficiently high net proton production is the same as for system (4.2),(4.3) and numerical simulations show, too, the same qualitative behavior of solutions.The rigorous qualitative study of system (4.2),(4.3) (without the modifications and assumptions made in Section 5) in terms of global well-posedness and singularity formation remains open.
The model could be extended to include effects of vascularization and necrosis.Indeed, it is largely accepted [7,10,61] that the hypoxic glioma cells induced to migrate away from sites with very low pH express, among others, proteases and vascular endothelial growth factors (VEGF) initiating and sustaining angiogenesis.Endothelial cells (ECs) are attracted chemotactically towards the garland-like structures of high glioma density surrounding the hypoxic area, which leads to further invasion and overall tumor expansion.A corresponding model should contain an adequate description of macroscopic EC dynamics, which could be obtained as well from an originally multiscale setting, similarly to that for glioma cells but taking into account the features specific to EC migration.does not seem to be appropriate for describing pseudopalisade formation.
The above suggests that pseudopalisades are not Turing-like patterns -at least as far as our model with constant diffusion coefficient is used for their description.To see the effect of the diffusion coefficient D T (x) we plot in Figure 14 the patterns obtained in 1D for the same parameter combination and several choices of D T (x), including various kinds of degeneracy.Thus, the second row in Figure 14 shows the case where D T degenerates on countable set, while the last row illustrates the situation with a strong degeneracy, i.e. on whole subintervals of the space domain (D T is of the type considered in [60] for a closely related problem, however with haptotaxis).

Figure 3 :
Figure 3: Tumor (upper row) and acidity (lower row) at several times for Experiment 1 and initial conditions (4.8).

Figure 4 :
Figure 4: Tumor (upper row) and acidity (lower row) at several times for Experiment 1 and initial conditions (4.9).

Figure 5 :
Figure 5: Tumor (upper row) and acidity (lower row) at several times for Experiment 2 and initial conditions (4.8)

Figure 6 :
Figure 6: Tumor (upper row) and acidity (lower row) at several times for Experiment 2 and initial conditions (4.9).

Figure 9 :
Figure 9: Difference between tumor (upper row) and acidity (lower row) at several times computed for System (4.2), (4.3) with and without pH-taxis in the framework of Experiment 2, initial conditions (4.9).

Figure 14 :
Figure 14: Patterns (in 1D) for tumor (left column) and acidity (middle column) for several choices of the tumor diffusion coefficient D T (x) (plots of the latter shown in the right column).Simulations are done for the dimension-free system (4.2),(4.3), the maximum simulation time corresponds to 10 months in dimensional framework.