A flux-limited model for glioma patterning with hypoxia-induced angiogenesis

We propose a model for glioma patterns in a microlocal tumor environment under the influence of acidity, angiogenesis, and tissue anisotropy. The bottom-up model deduction eventually leads to a system of reaction-diffusion-taxis equations for glioma and endothelial cell population densities, of which the former infers flux limitation both in the self-diffusion and taxis terms. The model extends a recently introduced [34] description of glioma pseudopalisade formation, with the aim of studying the effect of hypoxia-induced tumor vascularization on the establishment and maintenance of these histological patterns which are typical for high grade brain cancer. Numerical simulations of the population level dynamics are performed to investigate several model scenarios containing this and further effects.


Introduction
Classification of glioma, the most common type of primary brain tumors, typically comprises four grades, according to the degree of malignancy [32,39]. The highest grade IV (including the most aggressive type glioblastoma multiforme, shortly GBM) corresponds to characteristic histological patterns called pseudopalisades; they exhibit garland-like hypercellular structures surrounding necrotic regions usually centered around one or several sites of vasoocclusion [10,11,43,52], and are associated with poor patient survival prognosis [32].
More specifically, vascular occlusion and thrombosis [10,43] result into hypoxia, which promotes migration of glioma cells away from the highly acidic area, as well as tissue necrotization therein. Hypoxia around the pseudopalisading cells triggers the process of capillary formation, which involves a cascade of coordinated steps [37,47].
It typically starts with tumor cells overexpressing hypoxia-inducible regulators of angiogenesis, including vascular endothelial growth factor (VEGF) [11]. The latter acts as a chemoattractant and proliferation promotor for endothelial cells (ECs) lining the surrounding blood vessels. The chemotactic migration of ECs towards VEGF gradients is accompanied by formation of new capillaries off existing ones, thus leading to microvascular hyperplasia which is common in GBM [9,11] and which facilitates the growth of the neoplasm [11]. Understanding the mechanisms of tumor malignancy can contribute towards developing and tuning therapeutic strategies, e.g. by targeting angiogenesis [4] or by tumor alkalization [53], used in a (neo)adjuvant way to other approaches (such as surgery and radiotherapy).
Among the few mathematical models explicitly addressing glioma pseudopalisade formation (we refer to [12,13] for agent-based approaches and to [3,34,36] for continuous descriptions), the latter proposed a multiscale approach to deducing a reaction-(myopic)diffusion equation with repellent pH-taxis for glioma cell density coupled with a reaction-diffusion PDE for the acidity (expressed as proton concentration) produced by tumor cells. The derivation of this system started with describing the evolution of the glioma distribution function on the mesoscopic scale, in the KTAP (kinetic theory of active particles) framework [5]; thereby, the activity variable represented the amount of glioma transmembrane units occupied with protons.
A parabolic scaling led to the announced reaction-diffusion-taxis (RDT) equation on the so-called macroscopic level of the glioma population 1 , including in its terms influences of acidity and tissue coming from the lower modeling scales. The two equations for glioma and acidity were capable to qualitatively reproduce pseudopalisade-like patterns. In this work we aim to use a related approach and to extend the model in [34] such as to account for hypoxia-triggered angiogenesis and its effects on the pseudopalisade patterns. We thereby continue in the line of previous works [14, 16-18, 21, 23, 24, 29] dealing with cancer cell migration in heterogeneous tissues.
The rest of the paper is structured as follows: Section 2 is dedicated to setting up the micro-meso formulation in the KTAP framework and to deduce a macroscopic system of differential equations characterizing the dynamics of glioma density (RDT) with flux-limited diffusion and repellent pH-taxis, EC density (RDT) with linear diffusion and chemotaxis towards VEGF gradients, concentrations of acidity and VEGF (both with linear diffusion and tumor-and EC-dependent source terms). In Section 3 we perform numerical simulations for various scenarios of the macroscopic model obtained in Section 2 and interpret the results. Finally, Section 4 provides a discussion and an outlook on further problems of interest in this context.

Modeling
We begin by introducing some notations for the various variables and functions involved in our modeling process: we assume the cell speeds s, σ > 0 to be constant, S N−1 denotes the unit sphere in R N ; •v = v |v| ,θ = ϑ |ϑ | are unit vectors denoting the directions of vectors v ∈ V and ϑ ∈ Θ, respectively; • p(t, x, v): (mesoscopic) density function of glioma cells and M(t, x) = V p(t, x, v)dv: macroscopic density of glioma; • w(t, x, ϑ ): (mesoscopic) density function of ECs and W (t, x) = Θ w(t, x, ϑ )dϑ : macroscopic density of ECs.

Glioma cells
In the spirit of [17] we consider single cell (microscale) dynamics in the form where γ 1 , γ 2 ≥ 0 are some dimensionless weighting constants, K h , K M , X > 0 are further, dimensional constants (to be explained below 2 ), and the tensor B(v) := |v| 2 I N − v ⊗ v models biomechanical cell stress; thereby v ⊗ v represents active cell stress, while −|v| 2 I N is the isotropic part. Hence, similarly to [14,17,21], the velocity (and in particular the direction) of glioma cells is also influenced by a weighted sum of gradients (each inferring some limitation): increasing gradients of proton concentration have a repelling effect, and the cells also try to avoid too crowded regions. The constant α > 0 in (2.1) is a scaling parameter influenced by the gradients of h and M; for further details we refer to [17]. Notice that (2.1) actually implies d|v| 2 dt = 0, which is in line with our assumption of s = |v| being constant.
Correspondingly, on the mesoscale we have the following kinetic transport equation (KTE) for the evolution of glioma cell distribution function p(t, x, v): dv is a turning operator with constant turning rate λ > 0 and turning kernel K. As in [28] and many subsequent works [16-18, 23, 24, 34] we choose K(v, v ) := q(x,θ ) ω , meaning that the reorientation of cells is due to contact guidance by tissue, 2 with X to be chosen in Paragraph 3.1 with its orientation distribution q. Hence, x, v)). Moreover, as e.g., in [17] we assume p to be compactly supported in the velocity space.
The term P(M, h,W ) models glioma cell proliferation, which is supposed to depend on the amount of cells already present (irrespective of their orientation) and on favorable (nutrients provided by ECs of density W ) and unfavorable (acidity h) factors in the environment. More detailed descriptions account for proliferative interactions mediated by kernels depending on kinetic or activity variables, as in [17,24] or even differentiate between moving and proliferating cells [25,30]. Here we would rather focus on the macrolevel influences and choose a uniform velocity kernel to characterize such interactions, i.e.
where µ 1 , µ 2 > 0 are rates related to the cell growth or decay due to acidity and ECs. The positive constants K M and K W represent carrying capacities of glioma and endothelial cells, respectively, while K h > 0 denotes a threshold of proton concentration, which, when exceeded, leads to glioma cell death.
Similarly to [34] we consider a logistic type growth (as long as the environment is not too acidic), supplemented by a limited growth due to the interaction with ECs here representing the nutrient-supplying vasculature.
To deduce a macroscopic equation for the dynamics of macroscopic glioma density M from the KTE formulation, in [34] was used a parabolic upscaling 3 ; the microscale dynamics in that setting focused on the evolution of an activity variable quantifying the amount of transmembrane units occupied by protons. In [17], where velocity dynamics akin to (2.1) was considered (however without limitations as in the denominators of S), it was yet a parabolic scaling leading to the reaction-diffusion-taxis PDE for glioma evolution of the macroscale. Instead, we will consider here a moment closure approach similar to that in [14,29], ensuring the closure by way of an equilibrium distribution.
For that purpose we introduce the following quantities: commonly designated as momentum and pressure tensor, respectively. Thereby, U represents the so-called ensemble (macroscopic) cell velocity.
Multiplying (2.3) by v and integrating again w.r.t. v we get 3 with the observation that a hyperbolic one leads to a PDE which is too much drift-dominated and not able to reproduce pseudopalisade patterns hence due to (2.2) we can rewrite (2.8) as (2.9) System (2.8), (2.9) is not closed, since the pressure tensor P involves a second order moment of p and is therefore unknown. In order to close the system we proceed as e.g., in [29] and assume that it is close to equilibrium, which also dominates the second order moments. If we also assume that at equilibrium we can neglect cell proliferation (as it indeed happens much slower than the motility-related processes), then in virtue of (2.3) we can write which also leads to and hence at equilibrium the momentum of glioma population is driven by the cell ensemble aligning their motion to the average fiber orientationẼ q , and the 'population pressure' is generated by the variance-covariance matrix of the orientation distribution q/ω of tissue surrounding the cells.
Thus, we get from (2.9) The middle term on the right hand side of (2.13) stems from the description of forces acting on the cells; we suppose that the tensor MU ⊗ U therein is itself relating to the equilibrium situation, thus replace it by MẼ q ⊗Ẽ q . On the other hand, the left side in (2.13) can be interpreted as directional derivative of the momentum MU in the direction U of the 'cell flow'. As in [29] we consider that the cell flux relaxes quickly to its equilibrium, so that (2.14) If, furthermore, we assume that the tissue is undirected 4 , this translates intoẼ q = 0, which simplifies the above expression of the momentum: ω dv denoting as in previous works (see e.g. [23,34]) the so-called tumor diffusion tensor.
Plugging (2.15) into (2.7) then leads to the macroscopic PDE for glioma evolution (2.16) 4 as suggested by the simulation results in [34] which is of course merely an approximation, in view of the many assumptions made above.
The myopic diffusion term in (2.16) is made of a drift and a 'classical' diffusion term: hence using the expression of φ (h, M) we can rewrite (2. 16) in the form of a PDE in which both selfdiffusion (first term in (2.17)) and repellent pH-taxis (third term in (2.17)) are (at least partially) fluxlimited: (2.17) Notice that in order to have a genuine repellent pH-taxis the tensor s 2 λ I N − D T should be positive definite, which also ensures a nondegenerating, proper diffusion 5 . Under the hypothesis of undirected tissue, i.e. for E q = 0 this amounts to I N − V q = I N − S N−1 θ ⊗ θ q(θ )dθ being positive definite. This is, for instance, the case if q is the so-called peanut distribution [38]: where D water (x) denotes the (positive definite and usually diagonalized) water diffusion tensor. Same applies to our choice of q in Subsection 3.2. A much simplified situation is encountered when the stress tensor B(v) has only an isotropic part, while the active cell stress is neglected. The above deduction then directly yields instead of (2.17): For yet another choice of the right hand side in (2.1) leading in a different, but related way to a fluxlimited reaction-diffusion-taxis PDE for glioma evolving in a tissue (also under the influence of acidity and vascularization) we refer to [21].

Endothelial cells
We turn now to obtaining a macroscopic description for the evolution of ECs. For this purpose we can reproduce the derivation made in [16], the only difference being that ECs are supposed to respond here to (gradients of) VEGF concentration rather than proliferating glioma cells. 6 Therefore, we will only sketch the main steps and refer for details to [16].
We start with the KTE formulation for mesoscopic dynamics of EC distribution function w(t, x, ϑ ): with turning operator involving a uniform turning kernel and the turning rate depending upon the pathwise gradient D t g = g t + ϑ · ∇g of VEGF concentration. The coefficient function a(g) takes into account the way in which ECs perceive VEGF and respond to it, correspondingly adapting their turning. Following [16] we choose a(g) = χ a K R (K R +g/K g ) 2 , with χ a , K R > 0 constants. This corresponds to the rate of change of the expression y * (g) = g/K g g/K g +K R representing the equilibrium of EC receptor binding dynamics to g; the constant χ a scales changes in the turning rate per unit of change in dy * /dt.
The coefficient µ w (W, g) in the proliferation term is chosen as with µ W , K W , K g > 0 constants representing EC growth rate, EC carrying capacity, and maximum amount of VEGF, respectively. A parabolic scaling t ε 2 t, x εx performed in the usual way in combination with a Hilbert expansion w = w 0 + εw 1 + ε 2 w 2 + ..., a linearization of η, and the assumption that EC proliferation is much slower than migration then leads to the reaction-diffusion-chemotaxis PDE for the leading term in the (macroscopic) Hilbert expansion W = W 0 + εW 1 + ε 2 W 2 + ... (for details compare [16]): where D EC (x) := σ 2 Nη 0 I N and χ a (g) := σ 2 a(g) N I N = η 0 a(g)D EC (x). In the sequel we will (formally) approximate the evolution of W by that of its leading term W 0 and will use the notation W instead of W 0 .

Full macroscopic system
Equations (2.17) and (2.22) are coupled with the following equations set directly on the macroscopic scale and describing the dynamics of VEGF and acidity: where, D i , γ i and ζ i , for i = h, g, are positive constants representing diffusion coefficients, production, and uptake rates for the involved quantities.
The above PDE system needs to be supplied with initial conditions; they will be addressed in Section 3. Furthermore, since we are interested in the solution behavior inside a bounded region (representing e.g., a slice of a histologic sample), we restrict our equations (hitherto written for x ∈ R N ) to a bounded domain Ω ⊂ R N and consider no-flux boundary conditions for M,W, g, h. They are naturally assumed for g, h and can be obtained (at least for W ) during the upscaling process, as detailed in [17].

Numerical simulations
This section is dedicated to numerical simulations of the macroscopic system obtained in Section 2 for the quantities M,W, g, h.

Nondimensionalization and choice of parameters
We consider the following rescalings of the hitherto dimensional quantities: After dropping the tildes for simplicity of notation we therewith get the following nondimensionalized system: 1d) whose parameters are given by the above scalings and the values provided in Table 1.

Description of tissue
Typically the brain structure of each patient is assessed via diffusion tensor imaging (DTI) and the therewith provided water diffusion tensor D water already occurring in Subsection 2.1 above in the peanut distribution. The space scale on which glioma patterns are observed is, however, smaller than the best resolution of DTI, which does not go below the size of a voxel (ca. 1 mm 3 ): pseudopalisades are structures with a medium width of 200 − 400µm; for more details see the introduction of [34] and references therein. Therefore, as in that reference, we use an artificially created tissue with a corresponding water diffusion tensor as introduced in [38], in order to analyze the effect of tissue anisotropy. We consider the water diffusion tensor as where d(x, y) = 0.25e −0.005(x−450) 2 − 0.25e −0.005(y−450) 2 . A combination of uniform and von Mises-Fisher distributions is considered for the orientation distribution of tissue fibers: with δ ∈ [0, 1] being a weighting coefficient, ϕ 1 denotes the eigenvector corresponding to the leading eigenvalue of D water (x) and I 0 is the modified Bessel function of first kind of order 0. Furthermore, θ = (cos ξ , sin ξ ) for ξ ∈ [0, 2π], and k(x) = κFA(x), where κ ≥ 0 describes the sensitivity of cells towards orientation of tissue fibers and denotes the fractional anisotropy in 2D [38] with λ i (i = 1, 2) denoting the eigenvalues of D water (x).

Initial conditions
We consider the following initial densities and concentrations for the quantities involved:

Numerical method
We solve system (3.1) with no-flux boundary conditions and with the initial data given above on a square domain [0, 1000] × [0, 1000] (in µm) by using a finite difference method. The diffusion parts of acidity 7 , ECs and VEGF are discretized using a standard central difference scheme, while a nonnegativity-preserving discretization scheme, as proposed in [49], is used for the anisotropic diffusion of glioma cells. We apply a first order upwind scheme for the taxis terms in glioma and EC equations. For time discretization, an IMEX scheme is used for the acidity, EC, and VEFG equations, thereby treating the diffusion parts implicitly and discretizing the taxis and source terms with an explicit Euler method. We use an explicit Euler scheme for each term in the glioma equation.

Numerical experiments
We addressed in [34] the case with isotropic tissue (δ = 1) as well as that with anisotropy (δ = 0.2, k = 3); here similar observations can be made in this respect, thus we only consider the latter case. The parameters we use are given in Table 1, with changes correspondingly specified for each scenario.
Experiment 1 To start with, we solve the full system (3.1) with no-flux boundary conditions and with initial data as in Subsection 3.3. The results are shown in Figure 2. The garland-shaped glioma pattern forms (and is well visible after 60 simulated days) due to the high proton concentration in its central region, in spite of ECs being attracted by the glioma producing VEGF. In fact, ECs spread over almost the entire domain of simulation, coming from the left and right margins; the EC plot at t = 60 days exhibits large densities (except in the narrow central region, where the cells did not yet arrived), a phenomenon known as hyperplasia. Subsequently, ECs have accumulated enough near the tumor to uptake acidity and release nutrients, which leads among others to the tumor cells being able to move away from the acidic area and toward more favorable regions (due to γ 1 γ 2 the influence of acidity dominates the correction of self-diffusion), where they start producing acidity and VEGF. The pseudopalisades develop larger diameters and involve glioma garlands with lower densities. The cancer cells are followed by ECs, which again reduce the proton concentration, thus changing the glioma behavior. The pseudopalisades are disrupted and at later times both ECs and glioma cells can repopulate the previously too acidic regions.
Experiment 2 Same as Experiment 1, but with γ 1 = γ 2 , i.e. the flux-limited self-diffusion and pH-taxis being equally weighted. The results are illustrated in Figure 3. The behavior of solution components is qualitatively similar to the previous scenario, except the pseudopalisade-like patterns persisting for a longer time, being succeeded by a more uniform tumor-and EC spread occupying the whole domain (also migrating into the formerly most hypoxic areas) and exhibiting higher cell densities.

Experiment 3
Here we want to assess the effects of extending the model in [34]: , λ 0 , λ 1 , k D > 0 constants, by dynamics of quantities involved in angiogenesis, i.e. ECs and VEGF. Figure 4 shows differences between glioma density and proton concentration computed with model (3.1) and with model versions without vascularization, namely (3.5) (for the left part of Figure 4) and respectively for the right part of Figure 4.
The plots on the left column exhibit enhanced tumor growth and spread in the presence of ECs and VEGF. Earlier stages feature larger glioma densities (and therewith associated proton concentrations) distributed in a ring-shaped manner, which suggests more pronounced pseudopalisades. At later times, in the presence of vascularization, the tumor cells are not only able to proliferate, but also to reoccupy the formerly too acidic areas near the center of the simulation domain, while still responding to the spatial distribution of protons. This leads (as observed also in the previous numerical experiments) to patterns no longer having the pseudopalisade-like appearance. Similar effects can be observed in the plots on the right column, although less pronounced, due to the smaller structural difference between the two models compared therein. In fact, at later times there are only minimal differences between the plots on the left and on the right columns, which suggests that on the long run the effects of flux limitation in pH-taxis and additional self-diffusion from (3.6a) are increasingly diminished in comparison to the simple pH-taxis from (3.5), most probably because of the diffusion dominance in all these models. that the availability of vasculature permits an enhanced tumor spread in regions with higher diffusivity, but also enables invasion in areas where the cells modeled by (3.5) infer strongly degenerate diffusion and therewith associated annihilation of pH-taxis. Moreover, the model version (3.5) predicts high cell accumulation near the interface with highest diffusivity difference, whereas model (3.1) eludes such tendency potentially leading to singularities. This is presumably due to the flux-limitations in the taxis and additional self-diffusion terms (for more comments on this issue we refer to Section 4).

Experiment 4 Effects of flux limitations.
We compare the full model (3.1) with its counterpart involving in the motility terms only myopic diffusion and pH-taxis without flux saturation, i.e. characterizing glioma density dynamics by the rest of the equations in (3.1) remaining unchanged. The plots in Figure 5 illustrate the difference between solution components obtained with these two settings, in the above order. At a relatively early stage (30 days), with (3.1a) there are more glioma cells at the tumor edge, with an obvious preference for moving away from the more acidic, central area. The proton concentration is higher towards the domain margins, and with time passing (plots at 60 days) it becomes consistently larger all over the domain in the model version using (3.7), in which diffusion is dominating over taxis. In that setting the tumor cell density gets larger, too, which also applies to EC density; although (3.1) features (slightly) more glioma diffusion, however the flux-limited pH-taxis prevents a faster moving away from acidic regions. Instead, (3.7) ensures a wider spread of glioma, favored by the larger magnitude of the pH-tactic cell fluxes, which lets the tumor cells move faster and therewith produce the attracting VEGF at more locations. After a while (about further 30 days) the previous trend begins to get reversed, as far as glioma density is concerned: the flux limitation keeps the tumor more confined, but the acidity is still diffusing fast; the expression of VEGF is enhanced, which subsequently leads to hyperplasia, tumor enhancement, and revival of tiny glioma populations at more distant sites (in right half of the simulation domain). At the respective 'cancer spots' the cells stay more confined in the flux-limited case, while the model variant with (3.7) allows the cells to migrate faster (green-yellow margins in tumor difference plot at 90 days). The tumor cells accordingly produce acid and VEGF, attracting even more ECs at their main sites. Much later (300 days) (3.1) predicts larger densities of glioma cells in the center of the domain and near the boundary (the former because of the protons being increasingly removed by arriving ECs, and the latter due to the mutual proliferation support of the two cell types). The setting with (3.7) leads to pseudopalisades with a larger diameter, being preserved for a longer while. Still later (390 days and more) the differences become smaller and smaller in magnitude, while the areas with more or less tumor cells become mixed up, with the pseudopalisades disappearing and giving way to more or less heterogeneous patterns where the glioma cells and ECs are spread over the whole domain. The simulations (observe in particular the plots at 390 and 420 days) also suggest that the motion of ECs is dominated by diffusion (and its limitation near tumor cells) rather than chemotaxis towards VEGF. not; thus this model version should be less prone to singular structure formation. It also allows for sharper fronts of the cell mass, as can be seen in the right plot of the second row of Figure 6.

Discussion
The vast majority of continuous models for cell migration involve reaction-diffusion(-taxis) PDEs on the macroscopic scale, irrespective to whether those equations were obtained by a direct formulation or a more or less formal deduction from lower scales. Thereby, the diffusion is in some (most often constant) proportion to the density gradient, which is associated to infinite speed of propagation -an unrealistic feature, accompanied by smearing the sharp profile of the cell mass advancing into its surroundings; possible overcrowding effects cannot be properly captured either. These and other issues motivated the introduction of models with flux saturation also in this context. We refer e.g., to [6,40] for formal and respectively rigorous derivations of chemotaxis models from KTEs, and to [15,21,31] for settings specifically relating to glioma invasion. Among the latter, [21] contains a deduction of macroscopic equations with flux-saturated diffusion, haptotaxis, and chemotaxis; the passage from the mesoscopic KTAP framework to the RDT system is different from the one performed here, yet still formal. Very recently, [54] started from a KTE formulation similar to that in [21], hence also related to the one in the present work, and obtained by yet another upscaling method (with rigorous convergences) an RDT with myopic diffusion for the zeroth order approximation of the solution. That PDE did not involve flux-limitation, but the equation deduced for the first order correction did.
The model proposed here is still a bottom-up approach connecting subcellular and mesoscopic dynamics with the macroscopic evolution of cell populations in response to cues from the tumor microenvironment, leading to patterns which are qualitatively similar to those observed in histological imaging of glioblastoma. The setting is an extension of the previous one introduced in [34], since it accommodates the description of tumor vascularization in response to hypoxia. The main novelty here is the way of obtaining flux-limited pH-taxis and (supplementary) self-diffusion, namely upon considering on the single cell scale the joint effect of cell stress, chemical gradients, and population pressure to describe forces acting on the cells, cf. (2.1), (2.2). The approach is akin to that employed in [14,17,21], however it involves several modifications, as [14,17] did not contain any flux limitations, while [21] does, but models differently the single cell dynamics (allowing, among others, for non-constant speeds) and does not take into account cell reorientations (by way of a turning operator) in response to local tissue anisotropy and nonlocal orientation distribution of fibers and cells.
Our method also differs from [6,40], where the flux limitation stems from the choice of signal response function involved in the turning operator and depending on the directional derivative of the tactic signals (chemoattractants). We included in (2.2) only effects of pH-taxis and population pressure, but in modeling glioma migration through a highly anisotropic tissue one could also involve the influence of its gradients (via macroscopic tissue density, as in [21] or mesoscopic spatial and orientational distribution of tissue fibers, as in [17]). There are two 'sources' of diffusion in our model: the myopic diffusion comes from the description of cell reorientations via the turning operator, while the self-diffusion term with flux limitation originates from Newton's second law in (2.1), (2.2). The latter diffusion part has the effect of eluding diffusion degeneration when the tensor D T nullifies as a consequence of the mesoscopic tissue distribution q(x,v) vanishing locally. (Strongly) degenerating myopic diffusion, alone or in combination with taxis have been related to high cell aggregates and possible singularity formation even in lower dimensions [27,50,51], while flux limitation could alleviate such behavior, at least under certain conditions, see e.g., [7,8]. This tendency was also observed in Figure 7, in the case with strongly degenerating D T . Figure 6 showed a sharper tumor cell profile at the interface with large diffusivity drop, while the model without flux limitation generated a more uniform pattern with a tendency of smearing the interface and faster filling the areas of lower cell density. We have to stress, nevertheless, that our deduction of macroscopic PDEs is merely formal; as mentioned above, a rigorous one was done in [54], but for a simplified model and involving only limitation(s) of signal gradient(s) of the form ∇S 1+|∇S| instead of our choice encoded in φ (h, M). The obtained macroscopic model is able to reproduce, at least qualitatively, the histological patterns observed in patient biopsies. Indeed, the typical pseudopalisades due to severe hypoxia are formed even when (incipient) angiogenesis is present, followed by a disruption of the 'garlands' as consequence of acidity removal by vasculature (motile ECs) within the pseudopalisades and therewith associated repopulation with tumor cells of the formally acidic sites, which was not or very weakly noticed in simulations of the previous model [34]. The flux limitation does not seem to be necessary for the latter to happen, however it does influence the duration of pattern formation and maintenance, and (to a lesser extent) their appearance. We also noticed (cf. Figures 2 and 3) that the weights assigned to (flux-limited) pH-taxis and self-diffusion affect the patterns in terms of shape, duration, and persistence.
The macroscopic system (3.1) features two types of taxis: a repellent, flux-limited pH-taxis of glioma cells moving in the opposite direction of the proton concentration gradient ∇h, and chemotaxis of endothelial cells following the gradient of VEGF. As such, the setting does not fall into any class of models with multiple taxis as reviewed in [33], since each cell population is performing only one type of taxis, but it raises mathematical questions which are at least partially related to those models. The proton production by tumor cells gives the pH-taxis a direct character, while VEGF, the chemotactic cue of ECs, is produced by glioma cells in presence of acidity, while glioma proliferation is favorably influenced by ECs. The latter ensure uptake of protons as well as VEGF, hence exert an indirect and also a direct influence on their tactic signal. The mathematical analysis of system (3.1) with appropriate initial and boundary conditions is far from standard; this is not only due to the myopic diffusion and the nonlinearities arising from the source terms and from couplings, but especially to the flux-limited motility terms in the glioma reaction-diffusion-taxis PDE. Results about well-posedness and long time behavior of solutions to such systems are unknown. We merely provided a couple of examples of simulation-based, 1D pattern assessments in Figures 7 and 6 -with no pretension of completeness or rigor, but only for the purpose of giving a flavor of the patterns arising from such model. Beyond that, a traveling wave analysis, particularly of glioma and EC dynamics, would be of interest, both theoretically and from the application viewpoint.