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 Louis et al. (2007), 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 Rong et al. (2006),Wrensch et al. (2002). It is characterized by fast, infiltrative spread and unchecked cell proliferation which triggers hypoxia and upregulation of glycolysis, usually accompanied locally by exuberant angiogenesis Brat et al. (2002),Fischer et al. (2006); one of the typical features of GBM is the development of a necrotic core Louis et al. (2007). Increased extracellular pressure from edema and expression of procoagulant factors putatively lead to vasoocclusion and thrombosis Brat and Van Meir (2004), 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 Brat et al. (2004), forming palisade-like structures exhibiting arrangements of elongated nuclei stacked in rows at the periphery of the hypocellular region around the occlusion site Wippold et al. (2006). Such histopathological patterns are typically observed in GBM and are used as an indicator of tumor aggressiveness Brat and Mapstone (2003),Kleihues et al. (1995). Pseudopalisades can be narrow, with a width less than and a fibrillar interior structure, medium-sized (200–400 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 in width and are surrounding extensive necrotic areas, most often containing central vessels Brat et al. (2004).

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. Böttger et al. (2012),Khain et al. (2011),Sander and Deisboeck (2002)) 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 Corbin et al. (2018),Engwer et al. (2015a),Engwer et al. (2015b),Engwer et al. (2016),Hunt and Surulescu (2016),Painter and Hillen (2013) upon starting on the mesoscale from the kinetic theory of active particles (KTAP) framework Bellomo (2008) 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 (Cai et al. 2016; Caiazzo and Ramis-Conde 2015) for agent-based approaches and to (Alfonso et al. 2016; Martínez-González et al. 2012) for continuous settings. Of the latter, Alfonso et al. (2016) 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 nonlinear 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 Martínez-González et al. (2012) 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, and 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 Sect. 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, correspondingly leading to diffusion and drift-dominated evolution, respectively, and depending on tissue properties (directed/undirected). In Sect. 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 a version of the macroscopic system obtained by parabolic scaling. A result concerning the asymptotic behavior of such solution is proved as well. Finally, Sect. 6 contains a discussion of the obtained results and an outlook on further problems of interest related to GBM pseudopalisades. In the Appendix we address an alternative modeling approach and its parabolic limit and provide a linear stability analysis with comments on pattern formation for a version of the macroscopic limit obtained via parabolic scaling.

Model set up on subcellular and mesoscopic scales

The approach in (Engwer et al. 2015a, b) 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 (Kelkel and Surulescu 2011, 2012; Lorenz and Surulescu 2014), 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 Sect. 1) the size of such samples is too small to allow a reliable assessment of underlying tissue distribution via DTI1, 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 Holzer (2009).

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 Bellomo (2008)) and by 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 be a threshold value, which, when exceeded, leads to cancer cell death. The corresponding binding/occupying kinetics are writtenso that we can write for the corresponding subcellular dynamics (upon rescaling )where and represent the reaction rates. We denote by the steady-state of the above ODE, thus we haveAs in (Engwer et al. 2015a, b) 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 (Engwer et al. 2015a, b) 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.2 Next, we consider the path of a single cell starting at position and moving with velocity in the acidic environment. Since the glioma cells are supposed to move away from the highly hypoxic site, we take:which leads toWe denote by the density function of glioma cells at time t , position , velocity , and with activity variable . We assume as in (Corbin et al. 2018; Engwer et al. 2015a, b, 2016; Hunt and Surulescu 2016; Painter and Hillen 2013) that the cells have a constant (average) speed , so that , i.e. only the cell orientation is varying on the unit sphere. In terms of deviations from the steady-state (we also call z activity variable) we consider for the evolution of p the kinetic transport equation (KTE)where denotes the turning operator modeling cell velocity adaptations due to tissue contact guidance and acidity sensing, with denoting the turning rate of cells. Thereby, is a turning kernel giving the likelihood of a cell with velocity to change its velocity regime into . We adopt the choice proposed in Hillen (2006), i.e. where is the (stationary) orientation distribution of tissue fibers with and . We take the turning rate as in Engwer et al. (2015a)where and 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 becomeswithWe also employ the notation 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 Engwer et al. (2015b), but now accounting for the effect of acidity, we consider a source term of the formwhere represents the likelihood of cells having activity state to go into activity state z under the influence of acidity : higher acid concentrations hinder proliferation and even lead to apoptosis. In particular, is a kernel with respect to z, i.e. . The acidity is reported again to the threshold value . The growth rate depends on the total amount 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 Sect. 4.1. Hence, the presence of tissue is supporting proliferation, which is maintained until the environment becomes too acidic even for tumor cells.

The micro-meso formulation for glioma dynamics including the KTE (2.4) with the turning and proliferation operators introduced in (2.6) and (2.8), respectively, is supplemented with the evolution of acidity described by the macroscopic PDEwhere 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 not3, 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 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 .

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 (Hillen 2006; Engwer et al. 2015a) bythe 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. . We rescale the time and space variables by , . Since proliferation is much slower than migration, we also rescale with the corresponding term, as in Engwer et al. (2015b). For notation simplification we will drop in the following the symbol from the scaled variables t and .

Thus, from (3.1) and (3.2) we get:Now, using Hilbert expansions for the moments:and identifying the equal powers of , we get

:::If we also expand around , (3.9) leads toIntegrating (3.6) w.r.t. we getThen from (3.5) we obtain . Integrating (3.8) w.r.t. givesThe assumption of undirected tissue gives , thus from the above equation we obtain , which in virtue of (3.8) implies The compact Hilbert-Schmidt operator is considered as in Hillen (2006) on the weighted space . It has kernel4, thus its pseudo-inverse can be taken on the orthogonal complement , to deduce from (3.7)We summarize our hitherto information about the moments:Now integrating (3.10) w.r.t. we obtainUsing (3.11)-(3.14), the previous equation becomeswhere: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 encoding information about the anisotropy of underlying tissue and the function g(S) which relates to the subcellular dynamics of proton sensing and transfer across cell membranes. The myopic diffusionis common to this and previous models (Engwer et al. 2015a, b; Hunt and Surulescu 2016; Painter and Hillen 2013) obtained by parabolic scaling from the KTAP framework.

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 is nonzero, as the orientation distribution q is unsymmetric.

Consider on the Chapman-Enskog expansion of the cell distribution function in the formwhere and . Then for the moments introduced at the beginning of this Sect. 3 and with the notations , , we haveNow we rescale the time and space variables by and drop again the symbol to simplify the notation. As before, the proliferation term is scaled by . With these, the Eqs (3.1) and (3.2) become, respectively:andUsing (3.17) we writeSince q is independent of time, these equations implyIntegrating (3.22) w.r.t. giveswhere we used the notation .

From (3.23) we get at leading orderPlugging this in (3.23), we obtain (again at leading order)whenceFrom (3.24):Plugging this into (3.22) we get (at leading order)Since the right hand side vanishes when integrated w.r.t. , we can pseudo-invert and use (3.25) to gethenceso that (3.24) becomesComparing 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 in the direction of the dominating drift.

Numerical simulations

Parameters and coefficient functions

We assume a logistic type growth of the tumor cells and choosewhere is the growth rate and 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:

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 . 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). For more details we refer e.g. to (Engwer et al. 2015a; Painter and Hillen 2013) and references therein. As mentioned in Sect. 1, pseudopalisades are comparatively small structures with a medium width of . 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 in the space points of such a narrow region. To this aim we proceed as in Painter and Hillen (2013) and consider the water diffusion tensorwhere . For the fiber distribution function, we consider a mixture between uniform and von Mises-Fisher distributions, as follows:Here, is a weighting coefficient, is the eigenvector corresponding to the leading eigenvalue of and is the modified Bessel function of first kind of order 0. Also, for , and ) = FA, where FA denotes the fractional anisotropy: in 2D it has the form Painter and Hillen (2013)with () denoting the eigenvalues of . The parameter characterizes the sensitivity of cells towards orientation of tissue fibers. For perfectly aligned tissue (i.e., maximum anisotropy), and . Taking means, however, that the cells are insensitive to even such alignment and the distribution in 4.5 becomes a uniform one. Taking 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 , and the rest of parameters are the same as in (4.5). The first summand, similar to the choice in Hillen and Painter (2013), generates an orientation along the diagonal , while the second leads to alignment along the positive x and y directions. Due to , 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 Engwer et al. (2015b) 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 , where is the characteristic (diffusion) time. The latter is determined by assuming the underlying stochastic process behind water diffusion tensor measurements to be a 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:wherewith being the largest eigenvalue of .

Numerical experiments

The system (4.2), (4.3) is solved numerically on a square domain (in ) 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 Mosayebi et al. (2010) 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 Weickert (1998) instead of the standard one. Thereby, the derivatives are calculated in newly chosen directions (diagonal directions of the -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 Fig. 1a, b: and 1c, d, respectively:

Fig. 1

Initial conditions. Upper row: set (4.8) for tumor cell density a and acidity distribution b, lower row: set (4.9) for tumor cell density c and acidity distribution d

Experiment 1

Fully isotropic tissue

We begin by considering a fully isotropic tissue, i.e. taking in (4.5). The corresponding fractional anisotropy is everywhere , and the macroscopic tissue density Q is shown in Fig. 2a.

Fig. 2

a: Macroscopic tissue density (Experiments 1, 2) and b: mesoscopic tissue distribution for Experiment 2, for a given fiber direction

The simulations show (see Figs. 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 Figs. 3 and 4.

Fig. 3

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

Fig. 4

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

Experiment 2

Anisotropic tissue

With the choice , we describe an underlying tissue with pronounced anisotropy (two crossing fibre bundles). The corresponding mesoscopic fiber distribution q is shown in Fig. 2b for a fixed fiber direction, while the macroscopic tissue density Q remains unchanged.

The results of this experiment are shown in Figs. 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 Figs. 5b and 6b). The patterns at later times (see Figs. 5c and 6c) 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 Figs. 5 and 6). The differences between the acidity distributions in Fig. 5(d–f) and those in Fig. 3(d–f) (and correspondingly Fig. 6(d–f) and, respectively, those in Fig. 4(d–f) for the set of initial conditions (4.9)) are less prominent, since the acidity concentration S obeys in both cases a PDE with linear diffusion, where the tissue anisotropy has minor influence.

Fig. 5

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

Fig. 6

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

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 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 Ramnani (2020) 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 glycolysis) 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 migration and growth. Notice the more extensive tumor spread in comparison with Fig. 6. In Appendix B 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.

Fig. 7

Tumor (upper row) and acidity (lower row) at several times for Experiment 2, initial conditions (4.9), and stronger proton buffering

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 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.

Fig. 8

Tumor (upper row) and acidity (lower row) at several times for Experiment 2, initial conditions (4.9), and source term instead of that in (4.2). All parameters as in Table 1

Table 1

Parameters (dimensional quantities)

Parameter	Meaning	Value	Reference
M_max	Glioma carrying capacity	10^5-10^8cells/mm^3	Banerjee et al. 2015; Hathout et al. 2014; Rockne et al. 2010)
S_max	Acidity threshold for cancer cell death	10^-6.4mol/l	Webb et al. (2011)
s	Speed of glioma cells	2.8·10^-6mm/s	estimated, Prag et al. (2002)
λ_0	Turning frequency coefficient	0.1s^-1	(Engwer et al. 2015a; Sidani et al. 2007)
λ_1	Turning frequency coefficient	0.01-0.2s^-1	(Engwer et al. 2015a; Sidani et al. 2007)
k^+	Interaction rate tumor cells-protons	0.004s^-1	Lauffenburger and Lindermann (1993)
k^-	Detachment rate	0.01s^-1	Lauffenburger and Lindermann (1993)
β	Proton production rate	10^-9 mol /(mm^3s)	Estimated, Martin and Jain (1994)
α	Proton removal rate	10^-11 /s	Estimated
D_s	Acidity diffusion coefficient	5·10^-8-10^-5 mm^2/s	Estimated
μ_0	Glioma growth rate	0.2/day	Stein et al. (2007) Eikenberry et al. (2009)

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 . 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).

Fig. 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)

To see the effect of drift dominance we also solve the macroscopic system (3.30), (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 . The initial conditions are those in set (4.9), as visualized in Figs.1(c, d). Here we consider an unsymmetric tissue with mesoscopic orientational distribution as in (4.6). Figure 10 shows the mean fiber orientation corresponding to along with a magnification to observe the directionality in the neighborhood of the crossing fiber strands, and with plotted for and a specific direction .

Fig. 10

Mean fiber orientation a and zoom near crossing of fiber strands b for as in (4.6) with . c: mesoscopic tissue density for direction

The results obtained by solving the system for the evolution of tumor cells and acidity are shown in Fig. 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 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 Fig. 6 obtained for the parabolic limit and q as in (4.5) shows the radically different behavior w.r.t. the two approaches.

Fig. 11

Tumor (upper row) and acidity (lower row) at several times for as in (4.6) with and initial conditions (4.9). Solutions of system (3.30), (2.9) obtained by hyperbolic scaling

Similar observations apply when the first von Mises distribution in (4.6) exerts full influence (for ). Figure 12 shows tissue characteristics for this case: fractional anisotropy FA, zoomed , and for . 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 Fig. 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.

Fig. 12

Zoomed mean fiber orientation (a), fractional anisotropy FA b for as in (4.6) with . Subfigure c: mesoscopic tissue density for direction

Fig. 13

Tumor (upper row) and acidity (lower row) at several times for as in (4.6) with and initial conditions (4.9). Solutions of system (3.30), (2.9) obtained by hyperbolic scaling

We also solved the system (3.30), (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 . 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.

Qualitative analysis of the macroscopic reaction-diffusion-taxis system

Main results

We consider system (4.2), (4.3) with a slight modification of the source term in (4.3):with and , where we use , , and to denote the corresponding constants occurring in the expression of g(S) as given in Sect. 4.1. is considered to be a bounded domain with sufficiently smooth boundary , all involved constants are positive, , , , and . The no-flux boundary conditions are obtained through the upscaling procedure (as done e.g., in Corbin et al. (2021) for a related problem).

For the tumor diffusion tensor we require

Assumption 5.1

, , is uniformly bounded in , and for ;

there exists such that for any and ,

Theorem 5.1

Let . Suppose that Assumption 5.1 holds. Then if and , system (5.1) admits a unique global bounded classical solution.

Theorem 5.2

Under the assumptions of Theorem 5.1, suppose moreover that for all . Then there exists , such that if , for any we haveMoreover, there exists and such that for all ,

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 Winkler (2010) or Cao (2014) for example).

Lemma 5.1

Let () be a bounded domain with smooth boundary. Suppose that the nonnegative functions are in . Then there exist and a unique pair of non-negative functions (M, S) satisfyingand solving (5.1) classically in . Moreover, if , then

Next we prove results relating to the global boundedness of solutions to (5.1).

Lemma 5.2

There exists such that

Proof

Taking () as a test function for the S-equation of (5.1), for any , we obtainfrom which we obtainand then by Gronwall’s inequalityfrom which we obtain that for any ,From the arbitrariness of we therefore obtainOn the other hand, from the - estimates for the Neumann heat semigroup on a bounded domain and the fact thatwe obtain for all ,where denotes the first nonzero eigenvalue of in under the Neumann boundary condition.

Lemma 5.3

Under the assumptions of Theorem 5.1, for any , there exists such that for , we have

Taking as a test function for the M-equation of (5.1) and denoting , , then from the no-flux boundary conditions we obtainwhereThus we obtain that for any ,

Proof of Theorem 5.1

From Lemma 5.3 and the standard Moser iteration process, there exists such that for all . 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 defined in (5.7) such that for , for all , the functionsatisfieswherewith D a constant defined in (5.8).

According to the strong maximum principle and the assumption , M is positive in . Testing the M-equation of (5.1) by , by Young’s inequality and the fact of for , for , using the no-flux boundary condition, we obtain that there exists such thatwith Testing the S-equation of (5.1) by , we obtainCombining (5.5) and (5.6), we obtainBy choosingwe obtain that leads to .

Proof of Theorem 5.2

The proof of Theorem 5.2 is very standard. We include the proof here for completeness. Denote . Noticing that and for all , we obtain that and F(t) is nonnegative. From Lemma 5.4, we have and thenfor all , from which we haveUsing a similar argument as in Lemma 3.10 of Tao and Winkler (2015), we can obtain the uniform convergence of solutions, namelyas . Then there exists such that for all , , which together with the fact thatimplies thatfor all . Hencefrom which we obtainSubstituting (5.10) into (5.9), we obtainwhich implies that there exists such that for all ,Furthermore, notice that there exists a constant such thatThus the Gagliardo–Nirenberg inequality yieldsfor all . Similarly, we can obtainThis 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 Sect. 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 with . 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 Fig. 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 (Alfonso et al. 2016; Martínez-González et al. 2012) of pseudopalisade formation and spread, which are accounting for vascularization and necrosis rather than for direct effects of acidity. Nevertheless, our 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. (Boyd et al. December 2017; Koltai et al. 2020) 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 McIntyre et al. (2012).

In our deduction of the macroscopic system from the KTAP framework we used for the turning rate . This could be made more general, e.g. upon considering any regular enough function , expanding it around the steady-state , and keeping the first two terms of the expansion: . 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 are no longer constants, but depend on the macroscopic variable S by way of .5 Consequently, 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 from (A.3).

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 (Hiremath and Surulescu 2015, 2016; Hiremath et al. 2018; Kloeden et al. 2016), 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 Hillen (2006) and Corbin et al. (2021) 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 B 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 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 Sect. 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 (Brat et al. 2002; Brat and Van Meir 2004; Wippold et al. 2006) 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.