Abstract
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 .
Similar content being viewed by others
1 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 \(100\ \mu m\) and a fibrillar interior structure, medium-sized (200–400 \(\mu 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\ \mu m\) 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.
2 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 DTIFootnote 1, 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 \(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
so that we can write for the corresponding subcellular dynamics (upon rescaling \(y\leadsto 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 (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.Footnote 2 Next, we consider the path of a single cell starting at position \({\mathbf {x}}_0\) and moving with velocity \({\mathbf {v}}\) in the acidic environment. Since the glioma cells are supposed to move away from the highly hypoxic site, we take:
which leads to
We denote by \(p(t, \mathbf { x, v},y)\) the density function of glioma cells at time t , position \({\mathbf {x}} \in {\mathbb {R}}^n\), velocity \({\mathbf {v}} \in V \subset {\mathbb {R}}^n \), and with activity variable \(y \in Y=(0,1)\). 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 \(s>0\), so that \(V=s{\mathbb {S}}^{n-1}\), i.e. only the cell orientation is varying on the unit sphere. In terms of deviations \(z\in Z\subset [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 \({\mathscr {L}}[\lambda (z)] p := -\lambda (z)p + \lambda (z) \int _{V} K(\mathbf {x,v,v^\prime })p({\mathbf {v}}^\prime )d{\mathbf {v}}^\prime \) denotes the turning operator modeling cell velocity adaptations due to tissue contact guidance and acidity sensing, with \(\lambda (z)\) denoting the turning rate of cells. Thereby, \(K(\mathbf { x, v, v^\prime })\) is a turning kernel giving the likelihood of a cell with velocity \({\mathbf {v}}^\prime \) to change its velocity regime into \({\mathbf {v}}\). We adopt the choice proposed in Hillen (2006), i.e. \(K(\mathbf {x,v,v^{\prime }})=\frac{q({\mathbf {x}},\hat{{\mathbf {v}}})}{\omega }\) where \(q({\mathbf {x}},\hat{{\mathbf {v}}})\) is the (stationary) orientation distribution of tissue fibers with \(\omega = \int _{V} q(\hat{{\mathbf {v}}})d{\mathbf {v}}= s^{n-1}\) and \( \hat{{\mathbf {v}}} = \frac{{\mathbf {v}}}{|{\mathbf {v}}|} \in {\mathbb {S}}^{n-1}\). We take the turning rate as in Engwer et al. (2015a)
where \(\lambda _0\) and \(\lambda _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 Engwer et al. (2015b), but now accounting for the effect of acidity, we consider a source term of the form
where \(\chi (z,z')\) represents the likelihood of cells having activity state \(z'\) to go into activity state z under the influence of acidity \(S(t,\mathbf{x})\): higher acid concentrations hinder proliferation and even lead to apoptosis. In particular, \(\chi \) is a kernel with respect to z, i.e. \(\int _Z\chi (z,z')dz=1\). The acidity is reported again to the threshold value \(S_{max}\). The growth rate \(\mu (M)\) depends on the total amount \(M(t,\mathbf{x})=\int _V\int _Z p(t,\mathbf{x},\mathbf{v},z)dzd\mathbf{v}\) 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 PDE
where \(D_s\) is the diffusion coefficient of protons, \(\beta \) is the proton production rate by tumor cells, and \(\alpha \) 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 notFootnote 3, 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.
3 Macroscopic limits
We consider the following moments with respect to \({\mathbf {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 \({\mathbb {R}}^n\times V\times 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 (Hillen 2006; Engwer et al. 2015a) by
the mean fiber orientation and the variance-covariance matrix for the orientation distribution of tissue fibers, respectively.
3.1 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. \(\int _{{\mathbb {S}}^{n-1}}q(\mathbf{x}, \varvec{\theta })d\varvec{\theta }=\int _{{\mathbb {S}}^{n-1}}q(\mathbf{x}, -\varvec{\theta })d\varvec{\theta }\). We rescale the time and space variables by \({\tilde{t}}:= \epsilon ^2 t\), \({\tilde{\mathbf{x}}}:= \epsilon \mathbf{x}\). Since proliferation is much slower than migration, we also rescale with \(\epsilon ^2\) the corresponding term, as in Engwer et al. (2015b). For notation simplification we will drop in the following the \(\tilde{}\) symbol from the scaled variables t and \(\mathbf{x}\).
Thus, from (3.1) and (3.2) we get:
Now, using Hilbert expansions for the moments:
and identifying the equal powers of \(\epsilon \), we get
\({\epsilon ^0}\):
\({\epsilon ^1}\):
\({\epsilon ^2}\):
If we also expand \(\mu \) around \(M_0\), (3.9) leads to
Integrating (3.6) w.r.t. \({\mathbf {v}}\) we get
Then from (3.5) we obtain \(m_0 = \frac{q}{\omega }M_0\). Integrating (3.8) w.r.t. \({\mathbf {v}}\) gives
The assumption of undirected tissue gives \({\mathbb {E}}_q=\mathbf{0}\), thus from the above equation we obtain \(M_1^z = 0\), which in virtue of (3.8) implies
The compact Hilbert-Schmidt operator \({\mathscr {L}}[\lambda _0]m_1 = -\lambda _0 m_1 + \lambda _0 \frac{q}{\omega } M_1\) is considered as in Hillen (2006) on the weighted space \(L^2_{\frac{q}{\omega }}(V):=\{\zeta \ :\ \int _V|\zeta (\mathbf{v})|^2\ \frac{d\mathbf{v}}{\frac{q(\hat{{\mathbf {v}}})}{\omega }}<\infty \}\). It has kernelFootnote 4\(\langle \frac{q}{\omega } \rangle \), thus its pseudo-inverse can be taken on the orthogonal complement \(\langle \frac{q}{\omega } \rangle ^{\perp } \), to deduce from (3.7)
We summarize our hitherto information about the moments:
Now integrating (3.10) w.r.t. \({\mathbf {v}}\) we obtain
Using (3.11)-(3.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 \({\mathbb {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 across cell membranes. The myopic diffusion
is 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.
3.2 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 \({\mathbb {E}}_q\) is nonzero, as the orientation distribution q is unsymmetric.
Consider on \(L^2_{\frac{q}{\omega }}(V)=<q/\omega>\oplus <q/\omega >^\perp \) the Chapman-Enskog expansion of the cell distribution function \(p(t,\mathbf{x},\mathbf{v},z)\) in the form
where \(\int _V p^{\perp }(t,\mathbf{x},\mathbf{v},z)d\mathbf{v}=0\) and \({\bar{p}}(t,\mathbf{x},z):= \int _V p(t,\mathbf{x},{\mathbf {v}},z)d{\mathbf {v}}\). Then for the moments introduced at the beginning of this Sect. 3 and with the notations \(m^{\perp }(t,\mathbf{x},{\mathbf {v}}):=\int _Zp^{\perp } (t,\mathbf{x},{\mathbf {v}},z)dz\), \(m_\perp ^z(t,\mathbf{x},\mathbf{v}):=\int _Z z p^{\perp }(t,\mathbf{x},{\mathbf {v}},z)dz\), we have
Now we rescale the time and space variables by \({{\tilde{t}}}:=\varepsilon t,\ {{\tilde{\mathbf{x}}}} :=\varepsilon \mathbf{x}\) and drop again the \(\tilde{}\) symbol to simplify the notation. As before, the proliferation term is scaled by \(\varepsilon ^2\). With these, the Eqs (3.1) and (3.2) become, respectively:
and
Using (3.17) we write
Since q is independent of time, these equations imply
Integrating (3.22) w.r.t. \({\mathbf {v}}\) gives
where we used the notation \(\tilde{{\mathbb {E}}}_q(\mathbf{x}):=\int _V\mathbf{v}\frac{q}{\omega }(\mathbf{x},{{\hat{\mathbf{v}}}} )d\mathbf{v}=s{\mathbb {E}}_q\).
From (3.23) we get at leading order
Plugging this in (3.23), we obtain (again at leading order)
whence
From (3.24):
Plugging this into (3.22) we get (at leading order)
Since the right hand side vanishes when integrated w.r.t. \(\mathbf{v}\), we can pseudo-invert \({\mathscr {L}}[\lambda _0]\) and use (3.25) to get
hence
so that (3.24) 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 \(\varepsilon \)-corrections of the leading transport terms - together with the new advection which drives cells with velocity \(\frac{\varepsilon }{\lambda _0}{\mathbb {E}}_q\nabla \cdot {\mathbb {E}}_q\) in the direction of the dominating drift.
4 Numerical simulations
4.1 Parameters and coefficient functions
We assume a logistic type growth of the tumor cells and choose
where \(\mu _0\) is the growth rate and \(M_{max}\) is the carrying capacity of tumor cells.
4.2 Nondimensionalization
Considering the following nondimensional quantities:
and dropping the tildes for simplicity of notation, we get the following nondimensionalized system:
4.3 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 \({\mathbb {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 (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 \(200-400\ \mu 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 \({\mathbb {D}}_T\) 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 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, \(\delta \in \left[ 0,1\right] \) is a weighting coefficient, \(\varphi _1\) is the eigenvector corresponding to the leading eigenvalue of \({\mathbb {D}}_w({\mathbf {x}})\) and \(I_0\) is the modified Bessel function of first kind of order 0. Also, \(\varvec{\theta }= \left( \cos \xi , \sin \xi \right) \) for \(\xi \in \left[ 0,2\pi \right] \), and \(k({\mathbf {x}}\)) = \(\kappa \)FA\(({\mathbf {x}})\), where FA\(({\mathbf {x}})\) denotes the fractional anisotropy: in 2D it has the form Painter and Hillen (2013)
with \(\lambda _i\) (\(i=1,2\)) denoting the eigenvalues of \({\mathbb {D}}_w(\mathbf{x})\). The parameter \(\kappa \ge 0\) characterizes the sensitivity of cells towards orientation of tissue fibers. For perfectly aligned tissue (i.e., maximum anisotropy), \(FA(\mathbf{x})=1\) and \(k({\mathbf {x}}) = \kappa \). Taking \(\kappa =0\) means, however, that the cells are insensitive to even such alignment and the distribution in 4.5 becomes a uniform one. Taking \(\delta =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 \(k_h({\mathbf {x}}) = 0.05e^{-10^{-6} \left( (x-450)^2 + (y-450)^2 \right) }\), \(\gamma = \left( 1/\sqrt{2}, 1/\sqrt{2} \right) ^T\) 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 \(\gamma \), while the second leads to alignment along the positive x and y directions. Due to \(k_h({\mathbf {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 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 \(l_c=\sqrt{tr ({\mathbb {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 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:
where
with \(l_1\) being the largest eigenvalue of \({\mathbb {D}}_w\).
4.4 Numerical experiments
The system (4.2), (4.3) is solved numerically on a square domain \([0,1000]\times [0,1000]\) (in \(\mu 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 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 \(3 \times 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 Fig. 1a, b:
and 1c, d, respectively:
Experiment 1
Fully isotropic tissue
We begin by considering a fully isotropic tissue, i.e. taking \(\delta =1\) in (4.5). The corresponding fractional anisotropy is everywhere \(FA=0\), and the macroscopic tissue density Q is shown in Fig. 2a.
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.
Experiment 2
Anisotropic tissue
With the choice \(\delta =0.2\), \(\kappa =3\) 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.
Running numerical simulations for several different parameter sets led to the following observations:
-
The decisive parameter in this system seems to be \(\alpha \), 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 \(\mu _0\) plays a role, too, but a less prominent one. Concretely, pseudopalisade patterns form for very low values of \(\alpha \) (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 \(\alpha \), 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.
-
If \(\alpha \) 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 \(\mu _0 (1-M)\frac{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.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 \(\varepsilon =10^{-5}\). 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 \(q_h\) as in (4.6). Figure 10 shows the mean fiber orientation \({\mathbb {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 \(\delta =0.2\) and a specific direction \(\xi =\pi /2\).
Mean fiber orientation \({\mathbb {E}}_q\) a and zoom near crossing of fiber strands b for \(q_h\) as in (4.6) with \(\delta =0.2\). c: mesoscopic tissue density \(q_h\) for direction \(\xi =\pi /2\)
Zoomed mean fiber orientation \({\mathbb {E}}_q\) (a), fractional anisotropy FA b for \(q_h\) as in (4.6) with \(\delta =1\). Subfigure c: mesoscopic tissue density \(q_h\) for direction \(\xi =\pi /2\)
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 \(\delta =0.2\) makes the second term in (4.6) dominant). The cells ’escaping’ that influence move fast along the diagonal \(\gamma \) 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.
Similar observations apply when the first von Mises distribution in (4.6) exerts full influence (for \(\delta =1\)). Figure 12 shows tissue characteristics for this case: fractional anisotropy FA, zoomed \({\mathbb {E}}_q\), and \(q_h\) for \(\xi =\pi /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 \(\gamma \), 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.
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 \(\varepsilon \in [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.
5 Qualitative analysis of the macroscopic reaction-diffusion-taxis system
5.1 Main results
We consider system (4.2), (4.3) with a slight modification of the source term in (4.3):
with \(g(S):=\frac{\Lambda }{(S+K)^2(S+K+B)}\) and \(f(M,S)=\mu _0M(1-M)(1-S)\), where we use \(\Lambda =\lambda _1k_D\), \(K=k_D\), and \(B=\lambda _0\) to denote the corresponding constants occurring in the expression of g(S) as given in Sect. 4.1. \(\Omega \subset {\mathbb {R}}^N\) is considered to be a bounded domain with sufficiently smooth boundary \(\partial \Omega \), all involved constants are positive, \(M_0\in L^{\infty }(\Omega )\), \(M_0\not \equiv 0\), \(M_0, S_0\ge 0\), and \(S_0\in W^{1,\infty }(\Omega )\). 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 \({\mathbb {D}}_T\) we require
Assumption 5.1
-
(A)
\({\mathbb {D}}_T(\mathbf{x})\in \Big (C^{2,\gamma }(\Omega )\cap C({{\overline{\Omega }}})\Big )^{N\times N}\), \(\gamma \in (0,1)\), \(\mathbf{u}=\nabla \cdot {\mathbb {D}}_T\) is uniformly bounded in \(\Omega \), and \(\mathbf{u}(\mathbf{x})=0\) for \(\mathbf{x}\in \partial \Omega \);
-
(B)
there exists \(\vartheta >0\) such that for any \(\varvec{\xi }\in {\mathbb {R}}^N\) and \(\mathbf{x}\in \Omega \),
$$\begin{aligned} \varvec{\xi }^{\top }\cdot {\mathbb {D}}_T(\mathbf{x})\cdot \varvec{\xi }\ge \vartheta |\varvec{\xi }|^2. \end{aligned}$$
Theorem 5.1
Let \(N\ge 1\). Suppose that Assumption 5.1 holds. Then if \(\zeta <\alpha \) and \(\Vert S_0\Vert _{L^\infty (\Omega )}<1\), system (5.1) admits a unique global bounded classical solution.
Theorem 5.2
Under the assumptions of Theorem 5.1, suppose moreover that \(\nabla \cdot \mathbf{u}(\mathbf{x})=0\) for all \(\mathbf{x}\in \Omega \). Then there exists \(\mu ^*>0\), such that if \(\mu _0>\mu ^*\), for any \(\mathbf{x}\in \Omega \) we have
Moreover, there exists \(C>0\) and \(D>0\) such that for all \(t\in [0,+\infty )\),
5.2 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 \(\Omega \subset {\mathbb {R}}^N\) (\(N\ge 1\)) be a bounded domain with smooth boundary. Suppose that the nonnegative functions \(M_0, S_0\) are in \(W^{1,\infty }(\Omega )\). Then there exist \(T_{max}\in (0,\infty ]\) and a unique pair of non-negative functions (M, S) satisfying
and solving (5.1) classically in \(\Omega \times (0,T_{max})\). Moreover, if \(T_{max}<\infty \), then
Next we prove results relating to the global boundedness of solutions to (5.1).
Lemma 5.2
There exists \(C_S>0\) such that
Proof
Taking \(pS^{p-1}\) (\(p>1\)) as a test function for the S-equation of (5.1), for any \(\varepsilon \in (0,1)\), we obtain
from which we obtain
and then by Gronwall’s inequality
from which we obtain that for any \(t\in [0,T_{max})\),
From the arbitrariness of \(\varepsilon \in (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\in (0,T_{max})\),
where \(\lambda _1>0\) denotes the first nonzero eigenvalue of \(-\Delta \) in \(\Omega \subset {\mathbb {R}}^N\) under the Neumann boundary condition. \(\square \)
Lemma 5.3
Under the assumptions of Theorem 5.1, for any \(p>1\), there exists \(C(p)>0\) such that for \(t\in (0, T_{max})\), we have
Proof
Taking \(pM^{p-1}\) as a test function for the M-equation of (5.1) and denoting \(D_0:=\Vert {\mathbb {D}}_T(\cdot )\Vert _{L^\infty (\Omega )}\), \(D_1:=\Vert \mathbf{u}(\cdot )\Vert _{L^{\infty }(\Omega )}\), then from the no-flux boundary conditions we obtain
where
Thus we obtain that for any \(t\in (0,T_{max})\),
\(\square \)
Proof of Theorem 5.1
From Lemma 5.3 and the standard Moser iteration process, there exists \(C>0\) such that \(\Vert M(t,\cdot )\Vert _{L^\infty (\Omega )}\le C\) for all \(t\in (0,T_{max})\). Then in view of Lemma 5.1, Theorem 5.1 is a direct consequence of Lemma 5.2. \(\square \)
5.3 Long time behavior
Lemma 5.4
Under the assumptions of Theorem 5.2, there exists \(\mu ^*>0\) defined in (5.7) such that for \(\mu _0>\mu ^*\), for all \(t>0\), the function
satisfies
where
with D a constant defined in (5.8).
Proof
According to the strong maximum principle and the assumption \(M_0\not \equiv 0\), M is positive in \((0,\infty )\times \Omega \). Testing the M-equation of (5.1) by \(1-\frac{1}{M}\), by Young’s inequality and the fact of \(\nabla \cdot \mathbf{u}(\mathbf{x})=0\) for \(\mathbf{x}\in \Omega \), \(\mathbf{u}(\mathbf{x})=0\) for \(\mathbf{x}\in \partial \Omega \), using the no-flux boundary condition, we obtain that there exists \(C_M>0\) such that
with \(C_M:=\frac{D_0^2\Lambda ^2}{4\vartheta K^4(K+B)^2}.\) Testing the S-equation of (5.1) by \(S-\frac{\zeta }{2\alpha }\), we obtain
Combining (5.5) and (5.6), we obtain
By choosing
we obtain that \(\mu _0>\mu ^*\) leads to \(F'(t)\le -D(t)\). \(\square \)
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-\frac{1}{s}\) and \(h''(s)=\frac{1}{s^2}>0\) for all \(s>0\), we obtain that \(h(s)\ge h(1)=0\) and F(t) is nonnegative. From Lemma 5.4, we have \(F'(t)\le -D(t)\) and then
for all \(t>0\), from which we have
Using a similar argument as in Lemma 3.10 of Tao and Winkler (2015), we can obtain the uniform convergence of solutions, namely
as \(t\rightarrow \infty \). Then there exists \(t_0>0\) such that for all \(t>t_0\), \(\Vert M-1\Vert _{L^\infty (\Omega )}\le \frac{1}{2}\), which together with the fact that
implies that
for all \(t>t_0\). Hence
from which we obtain
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. \(\square \)
Remark 5.1
For the above rigorous results to hold we required among others that \(\zeta <\alpha \), 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 \(\zeta \) substantially exceeds \(\alpha \).
In fact, in Sect. 4.3 we already observed that \(\alpha \) (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 \(\zeta \); 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).
6 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 \(\lambda (z)=\lambda _0-\lambda _1z>0\). This could be made more general, e.g. upon considering any regular enough function \(\lambda \), expanding it around the steady-state \(y^*\), and keeping the first two terms of the expansion: \(\lambda (z)\simeq \lambda (y^*)-\lambda '(y^*)z:=\lambda _0(S)-\lambda _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 \(\lambda _0, \lambda _1\) are no longer constants, but depend on the macroscopic variable S by way of \(y^*\).Footnote 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 \(\lambda (\mathbf{v},S)\) 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 \(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 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.
Notes
The typical size of a voxel is ca. \(1\ mm^3\)
In fact, keeping this dependency and accounting for the correct scales during the macroscopic limit leads to omission of the arising supplementary term, so that the outcome is the same whether we do such assumption at this step or later on.
By ’undirected’ we mean as in Hillen (2006) that the fibers are symmetric along their axis and both fiber directions are identical, while ’directed’ means unsymmetric fibers; since the common medical imaging techniques are not providing the necessary resolution, it is actually not known whether brain tissue is directed, but such feature might play a role in the formation of glioma patterns
The brackets \(\langle \cdot \rangle \) denote the span set
An essential requirement on \(\lambda \) is thereby to ensure that \(\lambda (y^*)>0\) and \(\lambda '(y^*)>0\).
Notice the difference of sign in the exponent: the cells are supposed to follow the decreasing gradient of signal S
The latter could serve, however, for a model extension including angiogenesis.
References
Alfonso JCL, Köhn-Luque A, Stylianopoulos T, Feuerhake F, Deutsch A, Hatzikirou H (2016) Why one-size-fits-all vaso-modulatory interventions fail to control glioma invasion: in silico insights. Sci Rep 6:37283
Banerjee S, Khajanchi S, Chaudhuri S (2015) A mathematical model to elucidate brain tumor abrogation by immunotherapy with t11 target structure. PLOS ONE 10(5):e0123611
Bartel P, Ludwig FT, Schwab A, Stock C (2012) pH-taxis: directional tumor cell migration along pH-gradients. Acta Physiologica 204(Suppl. 689):113
Bellomo N (2008) Modeling complex living systems. Birkhäuser, Boston
Böttger K, Hatzikirou H, Chauviere A, Deutsch A (2012) Investigation of the migration/proliferation dichotomy and its impact on avascular glioma invasion. Math Model Nat Phenom 7:105–135
Boyd NH, Walker K, Fried J, Hackney JR, McDonald PC, Benavides GA, Spina R, Audia A, Scott SE, Libby CJ, Tran AN, Bevensee MO, Griguer C, Nozell S, Gillespie GY, Nabors B, Bhat KP, Bar EE, Darley-Usmar V, Xu B, Gordon E, Cooper SJ, Dedhar S, Hjelmeland AB (2017) Addition of carbonic anhydrase 9 inhibitor SLC-0111 to temozolomide treatment delays glioblastoma growth in vivo. JCI Insight 2(24):256
Brat DJ, Mapstone TB (2003) Malignant glioma physiology: cellular response to hypoxia and its role in tumor progression. Ann Intern Med 138(8):659–668
Brat DJ, Van Meir EG (2004) Vaso-occlusive and prothrombotic mechanisms associated with tumor hypoxia, necrosis, and accelerated growth in glioblastoma. Lab Investig 84(4):397–405
Brat DJ, Castellano-Sanchez A, Kaur B, Van Meir EG (2002) Genetic and biologic progression in astrocytomas and their relation to angiogenic dysregulation. Adv Anat Pathol 9(1):24–36
Brat DJ, Castellano-Sanchez AA, Hunter SB, Pecot M, Cohen C, Hammond EH, Devi SN, Kaur B, Van Meir EG (2004) Pseudopalisades in glioblastoma are hypoxic, express extracellular matrix proteases, and are formed by an actively migrating cell population. Cancer Res 64(3):920–927
Cai Y, Wu J, Li Z, Long Q (2016) Mathematical modelling of a brain tumour initiation and early development: a coupled model of glioblastoma growth, pre-existing vessel co-option, angiogenesis and blood perfusion. PloS one 11(3):e0150296
Caiazzo A, Ramis-Conde I (2015) Multiscale modelling of palisade formation in gliobastoma multiforme. J Theor Biol 383:145–156
Cao X (2014) Boundedness in a quasilinear parabolic-parabolic keller-segel system with logistic source. J Math Anal Appl 412(1):181–188
Corbin G, Hunt A, Klar A, Schneider F, Surulescu C (2018) Higher-order models for glioma invasion: from a two-scale description to effective equations for mass density and momentum. Math Models Methods Appl Sci 28(09):1771–1800
Corbin G, Engwer C, Klar A, Nieto J, Soler J, Surulescu C, Wenske M (2021) Modeling glioma invasion with anisotropy- and hypoxia-triggered motility enhancement: from subcellular dynamics to macroscopic pdes with multiple taxis. Math Models Methods Appl Sci 31(01):177–222
Eikenberry SE, Sankar T, Preul MC, Kostelich EJ, Thalhauser CJ, Kuang Y (2009) Virtual glioblastoma: growth, migration and treatment in a three-dimensional mathematical model. Cell Prolif 42(4):511–528
Engwer C, Hillen T, Knappitsch M, Surulescu C (2015a) Glioma follow white matter tracts: a multiscale DTI-based model. J Math Biol 71(3):551–582
Engwer C, Hunt A, Surulescu C (2015b) Effective equations for anisotropic glioma spread with proliferation: a multiscale approach and comparisons with previous settings. Math Med Biol: an IMA J 33(4):435–459
Engwer C, Knappitsch M, Surulescu C (2016) A multiscale model for glioma spread including cell-tissue interactions and proliferation. J Math Biosci Eng 13:443–460
Fischer I, Gagner J-P, Law M, Newcomb EW, Zagzag D (2006) Angiogenesis in gliomas: biology and molecular pathophysiology. Brain Pathol 15(4):297–310
Hathout L, Ellingson BM, Cloughesy T, Pope WB (2014) A novel bicompartmental mathematical model of glioblastoma multiforme. Int J Oncol 46(2):825–832
Hillen T (2006) \(M^5\) mesoscopic and macroscopic models for mesenchymal motion. J Math Biol 53:585–616
Hillen T, Painter KJ (2013) Transport and anisotropic diffusion models for movement in oriented habitats In Dispersal, individual movement and spatial ecology. Springer, Heidelberg
Hiremath S, Surulescu C (2015) A stochastic multiscale model for acid mediated cancer invasion. Nonlinear Anal: Real World Appl 22:176–205
Hiremath S, Surulescu C (2016) A stochastic model featuring acid-induced gaps during tumor progression. Nonlinearity 29(3):851–914
Hiremath S, Surulescu C, Zhigun A, Sonner S (2018) On a coupled SDE-PDE system modeling acid-mediated tumor invasion. Discret Contin Dyn Syst - B 23(6):2339–2369
Holzer P (2009) Acid-sensitive ion channels and receptors. Sensory nerves. Springer, Berlin Heidelberg, pp 283–332
Hunt A, Surulescu C (2016) A multiscale modeling approach to glioma invasion with therapy. Vietnam J Math 45(1–2):221–240
Kelkel J, Surulescu C (2011) On some models for cancer cell migration through tissue networks. Math Biosci Eng 8(2):575–589
Kelkel J, Surulescu C (2012) A multiscale approach to cell migration in tissue networks. Math Models Methods Appl Sci 22(03):1150017
Khain E, Katakowski M, Hopkins S, Szalad A, Zheng X, Jiang F, Chopp M (2011) Collective behavior of brain tumor cells: the role of hypoxia. Phys Rev E 83:031920
Kleihues P, Soylemezoglu F, Schäuble B, Scheithauer BW, Burger PC (1995) Histopathology, classification and grading of gliomas. Glia 5:211–221
Kloeden PE, Sonner S, Surulescu C (2016) A nonlocal sample dependence SDE-PDE system modeling proton dynamics in a tumor. Dis Contin Dyn Syst - Series B 21(7):2233–2254
Koltai T, Reshkin SJ, Harguindey S (2020) The pH-centered paradigm in cancer. in an innovative approach to understanding and treating cancer targeting pH. Elseiver, Amsterdam, pp 53–97
Lauffenburger DA, Lindermann JL (1993) Receptors. models for binding, trafficing and signaling. Oxford University Press, Oxford
Lorenz T, Surulescu C (2014) On a class of multiscale cancer cell migration models: Well-posedness in less regular function spaces. Math Models Methods Appl Sci 24(12):2383–2436
Louis DN, Ohgaki H, Wiestler OD, Cavenee WK, Burger PC, Jouvet A, Scheithauer BW, Kleihues P (2007) The 2007 WHO classification of tumours of the central nervous system. Acta Neuropathologica 114(2):97–109
Loy N, Preziosi L (2020) Kinetic models with non-local sensing determining cell polarization and speed according to independent cues. J Math Biol 80:374–421
Martin GR, Jain RK (1994) Noninvasive measurement of interstitial pH profiles in normal and neoplastic tissue using fluorescence ratio imaging microscopy. Cancer Res 54(21):5670–5674
Martínez-González A, Calvo GF, Pérez Romasanta LA, Pérez-García VM (2012) Hypoxic cell waves around necrotic cores in glioblastoma: a biomathematical model and its therapeutic implications. Bull Math Biol 74(12):2875–2896
McIntyre A, Patiar S, Wigfield S, Li JI, Ledaki I, Turley H, Leek R (2012) Carbonic anhydrase IX promotes tumor growth and necrosis in vivo and inhibition enhances anti-VEGF therapy. Clin Cancer Res 18(11):3100–3111
Mosayebi P, Cobzas D, Jagersand M, Murtha A(2010) Stability effects of finite difference methods on a mathematical tumor growth model. In 2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition-Workshops, pages 125–132. IEEE,
Othmer HG, Hillen T (2002) The diffusion limit of transport equations II: chemotaxis equations. SIAM J Appl Math 62(4):1222–1250
Painter K, Hillen T (2013) Mathematical modelling of glioma growth: the use of diffusion tensor imaging (DTI) data to predict the anisotropic pathways of cancer invasion. J Theor Biol 323:25–39
Paradise RK, Whitfield MJ, Lauffenburger DA, Van Vliet KJ (2013) Directional cell migration in an extracellular pH gradient: a model study with an engineered cell line and primary microvascular endothelial cells. Exper Cell Res 319(4):487–497
Perthame B, Tang M, Vauchelet N (2016) Derivation of the bacterial run-and-tumble kinetic equation from a model with biochemical pathway. J Math Biol 73(5):1161–1178
Prag S, Lepekhin EA, Kolkova K, Hartmann-Petersen R, Kawa A, Walmod PS, Belman V, Gallagher HC, Berezin V, Bock E, Pedersen N (2002) Ncam regulates cell motility. J Cell Sci 115(2):283–292
Ramnani D WebPathology - visual survey of surgical pathology. https://www.webpathology.com
Rockne R, Rockhill JK, Mrugala M, Spence AM, Kalet I, Hendrickson K, Lai A, Cloughesy T, Alvord EC, Swanson KR (2010) Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach. Phys Med Biol 55(12):3271–3285
Rong Y, Durden DL, Van Meir EG, Brat DJ (2006) ‘Pseudopalisading’ necrosis in glioblastoma: a familiar morphologic feature that links vascular pathology, hypoxia, and angiogenesis. J Neuropathol Exper Neurol 65(6):529–539
Sander LM, Deisboeck TS (2002) Growth patterns of microscopic brain tumors. Phys Rev E 66(5):051901
Sidani M, Wessels D, Mouneimne G, Ghosh M, Goswami S, Sarmiento C, Wang W, Kuhl S, El-Sibai M, Backer JM et al (2007) Cofilin determines the migration behavior and turning frequency of metastatic cancer cells. J Cell Biol 179(4):777–791
Stein AM, Demuth T, Mobley D, Berens M, Sander LM (2007) A mathematical model of glioblastoma tumor spheroid invasion in a three-dimensional in vitro experiment. Biophys J 92(1):356–365
Tao Y, Winkler M (2015) Large time behavior in a multidimensional chemotaxis-haptotaxis model with slow signal diffusion. SIAM J Math Anal 47(6):4229–4250
Webb BA, Chimenti M, Jacobson MP, Barber DL (2011) Dysregulated pH: a perfect storm for cancer progression. Nat Rev Cancer 11(9):671–677
Weickert J (1998) Anisotropic diffusion in image processing. Teubner Stuttgart, Germany
Winkler M (2010) Global solutions in a fully parabolic chemotaxis system with singular sensitivity. Math Methods Appl Sci 34(2):176–190
Winkler M, Surulescu C (2017) Global weak solutions to a strongly degenerate haptotaxis model. Communi Math Sci 15(6):1581–1616
Wippold FJ, Lämmle M, Anatelli F, Lennerz J, Perry A (2006) Neuropathology for the neuroradiologist: palisades and pseudopalisades. Am J Neuroradiol 27(10):2037–2041
Wrensch M, Minn Y, Chew T, Bondy M, Berger MS (2002) Epidemiology of primary brain tumors: current concepts and review of the literature. Neuro-Oncol 4:278–299
Acknowledgements
PK acknowledges funding by DAAD in form of a PhD scholarship. CS was funded by BMBF in the project GlioMaTh 05M2016.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix
A An alternative approach to including acidity effects
In the derivation performed via parabolic scaling in Sect. 3.1 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 Othmer and Hillen (2002) and re-employed in Loy and Preziosi (2020) 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 Perthame et al. (2016), 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 \(\rho (t,\mathbf{x}, \mathbf{v})\):
and its parabolic limit. Here the proliferation operator is defined with the same \(\mu (M,S)\) as previously used in Sect. 3 and takes the form
while for the turning rate we setFootnote 6
where \(D_tS:=S_t+\mathbf{v}\cdot \nabla _\mathbf{x}S\) is the pathwise gradient of S. The coefficient function h(S) is to be chosen later.
We use again a parabolic scaling \({\tilde{t}}:= \epsilon ^2 t\), \({\tilde{x}}:= \epsilon x\) and rescale as before the proliferation term by \(\epsilon ^2\), thus
Performing a Hilbert expansion \(\rho =\rho _0+\varepsilon \rho _1+\varepsilon ^2 \rho _2+\dots \) and equating the powers of \(\epsilon \) yields
\(\varvec{\varepsilon ^0}\):
thus
\(\varvec{\varepsilon ^1}\):
thus by (A.6) and the assumption of undirected tissue
This can be rewritten as
Since the integral of the right hand side w.r.t. \(\mathbf{v}\) vanishes we can pseudo-invert \({\mathscr {L}}[\lambda _0]\) as before, to get
\(\varvec{\varepsilon ^2}\):
Integrating (A.10) w.r.t \({\mathbf {v}}\) gives
Using (A.9) leads to
which differs from (3.15) only by the tactic sensitivity coefficient. This suggests to choose
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 \(\lambda _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^-+\lambda _0)\). While a factor \(\text {const}\cdot f'(S)\) is also encountered in Othmer and Hillen (2002), the denominator obtained here for the tactic sensitivity appears due to the specific choice of the turning rate \(\lambda (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).
B Linear stability analysis for a version of (4.2), (4.3) with constant tumor diffusion coefficient
For simplicity we perform a 1D analysis; the extension to a higher dimensional case involving a constant tumor diffusion tensor in diagonal form is straightforward.
The uniform steady-states are \(P_1=(0,0)\), \(P_2=(1,\frac{1}{\alpha })\), and \(P_3=(\alpha ,1)\). In the absence of diffusion and taxis, \(P_1\) and \(P_2\) are saddles, while \(P_3\) is stable for \(0<\alpha <1\). Thus, we investigate the possibility of Turing-like patterns only around \(P_3\) for such biologically relevant \(\alpha \).
Let \(P_*:=(M_*,S_*)\) be a steady-state and consider the perturbations \(u:=M-M_*\), \(\sigma :=S-S_*\). Linearizing (4.2), (4.3) (with \(D_T=D\) constant) around \(P_*\) leads to
where \({\mathbb {A}}=\begin{pmatrix} D&{}D\alpha g(S_*)S_*\\ 0&{}1\end{pmatrix}\) and \({\mathbb {B}}=\begin{pmatrix}\mu _0(1-2\alpha S_*)(1-S_*)&{}-\alpha \mu _0S_*(1-\alpha S_*)\\ 1&{}-\alpha \end{pmatrix}\).
We look for solutions of the form \(\sum _{k\ne 0}T_k(t)\mathbf{X}_k(x)\), with \(\Delta \mathbf{X}_k(x)+k^2\mathbf{X}_k(x)=\mathbf{0}\) in \(\Omega \) and \(\nabla \mathbf{X}_k\cdot \varvec{\nu }=\mathbf{0}\) on \(\partial \Omega \), thus \(\mathbf{X}_k\) are eigenfunctions of \(-\Delta \), each corresponding to the wavenumber k.
The terms making up the solution \((u,\sigma )\) will involve exponents of the eigenvalues \(\lambda _k\) of the matrix \({\mathbb {B}}-k^2{\mathbb {A}}\). It holds that
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
For Turing-like patterns around \(P_*\) we need \(\text{ det }({\mathbb {B}}-k^2{\mathbb {A}})<0\), which means that there is a positive eigenvalue. This condition takes for \(P_*=P_3\) the concrete form
This cannot be satisfied for \(\alpha < 1\), which is typical for higher grade tumors (especially for GBM). In fact, in view of the nondimensionalization done in Sect. 4.1, it is very improbable to have \(\alpha > 1\) for this problem.
To see the effect of pH-taxis (repellence by acidity) we set \(g(S)=0\), which only enhances the chances of (B.4) to hold (still for \(\alpha <1\)), thus the presence of pH-taxis does not dramatically change the patterning behavior; this is due to the tactic bias being repellent; an attractive pH-taxis (as proposed in Bartel et al. (2012) for melanoma cells or in Paradise et al. (2013) for endothelial cellsFootnote 7) could render (B.4) valid even for \(\alpha <1\), but 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 Fig. 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 Fig. 14 shows the case where \(D_T\) degenerates on a 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 Winkler and Surulescu (2017) for a closely related problem, however with haptotaxis).
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Kumar, P., Li, J. & Surulescu, C. Multiscale modeling of glioma pseudopalisades: contributions from the tumor microenvironment. J. Math. Biol. 82, 49 (2021). https://doi.org/10.1007/s00285-021-01599-x
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s00285-021-01599-x
Keywords
- Glioblastoma
- Pseudopalisade patterns
- Hypoxia-induced tumor behavior
- Kinetic transport equations
- Upscaling
- Reaction-diffusion-taxis equations
- Global existence
- Uniqueness
- Long time behavior
- Multiscale modeling
- Directed/undirected tissue