Multiphase modelling of glioma pseudopalisading under acidosis

1.   Introduction

Glioblastoma is the most common type of primary brain tumors in adults, with a dismal prognosis. The histological features include characteristic patterns called pseudopalisades, which exhibit garland-like structures made of aggregates of glioma cells stacked in rows at the periphery of regions with low pH and high necrosis surrounding the occlusion site(s) of one or several capillaries ^[39]. Such patterns are used to grade the tumor and are essential for diagnosis ^[6,23].

The few continuous mathematical models proposed for the description of pseudopalisade patterns ^[1,26,27,30], involve systems of ODEs and PDEs which were set up in a heuristic manner or obtained from lower scale dynamics. They account for various aspects like phenotypic switch between proliferating and migrating glioma as a consequence of vasoocclusion and nutrient suppression associated therewith ^[1], interplay between normoxic/hypoxic glioma, necrotic matter, and oxygen supply ^[30], or tissue anisotropy and repellent pH-taxis, with ^[27] or without ^[26] vascularisation. The latter two works deduced effective equations for glioma dynamics from models set on the lower, microscopic and mesoscopic scales in the kinetic theory of active particles (KTAP) framework explained e.g., in ^[4]. Those deductions are in line with previous works concerning glioma invasion in anisotropic tissue ^{[9,10,11,14,15,16,18,31]}, which use parabolic scaling to obtain the equations for glioma evolution on the macroscale from descriptions of subcellular and mesoscopic dynamics. Here we propose yet another approach, relying on the interpretation of the relevant components (glioma, normal tissue, necrotic matter, and acidity) as phases in a mixture - with the exception of acidity, which is characterised by proton concentration in the volume occupied by the phases.

Multiphase models in the framework of mixture theory ^[2,12] have been considered in connection with cancer growth and invasion, (arguably) starting with ^[7,8] and followed by many others, involving different mathematical and biological aspects; see e.g., ^{[17,19,32,33,34,36,37,40]} and references therein. Particularly ^[34] provides a comprehensive discussion of the classes of multiphase models, along with their strengths and drawbacks. Such models employ a population level description with conservation laws analogous to balance equations for single cells, with supplementary terms accounting for interphase effects. Thereby, the living cells and tissues are most often seen as viscous fluids, while the interstitial fluid is inviscid. Few of the existing models (e.g., ^[17,34,40]) are explicitly handling necrosis, which is, however, an important component of advanced tumors. Likewise, there are relatively few models accounting for chemoattracting or chemorepellent agents (e.g., ^[17,37]), which are known to bias in a decisive way the expansion of the neoplasm and therewith associated dynamics of tumor cells and surrounding tissue. To our knowledge there are no multiphase models for glioma pseudopalisade development. In this context we propose such a model comprising necrotic matter as one of the phases in the mixture (although there are a few issues related to this approach, as mentioned, e.g., in ^[34]). As our main aim is to describe glioma behavior in an acidic environment fastly leading to extensive necrosis, we explicitly include these influences in our model. From the mass and momentum balance equations written for the phases composing the neoplasm and supplemented with appropriate constitutive relations, we then deduce, under some simplifying assumptions, a system of reaction-cross diffusion equations with repellent pH-taxis for the glioma cells and normal and necrotic tissues, also including a solenoidality constraint on the total flux. Our deduction has some similarities with that in ^[19], however is done in N dimensions instead of 1D, it accounts explicitly for the evolution of necrotic matter and the effects of acidity, and does not require the drag coefficients between the phases to be equal. The obtained equations are able to reproduce qualitatively the typical pseudopalisade patterns with all their aspects related to the glioma aggregates, acidity, necrotic inner region, normal tissue dynamics.

The rest of this paper is organized as follows: Section 2 contains some notations and conventions. Section 3.1 contains the model setup with the considered phases (glioma cells, normal tissue, necrotic matter, and interstitial fluid) and the corresponding mass and momentum balance equations, along with their associated constitutive relations. That description is used to deduce the announced reaction-cross diffusion equations. Thereby, we investigate a model with an immovable component and observe that enforcing one of the phases to be fixed prevents the simultaneous validity of all basic conservation laws needed in the setting. Then we neglect the interstitial fluid, thus reducing the number of phases, and obtain the PDEs with reaction, nonlinear diffusion and taxis terms, while still ensuring all balance equations. The total flux needs to be given, satisfying a solenoidality constraint. Section 4 is dedicated to the existence of weak solutions to the obtained cross diffusion system coupled with the reaction-diffusion equation for acidity dynamics in the volume of interest. Finally, in Section 5 we provide numerical simulations for that system, studying several parameter scenarios in order to put in evidence the effect of acidity and of different drag coefficients on the obtained patterns exhibiting pseudopalisade formation.

2.   Preliminaries

We will use the following notation throughout this paper:

Notation 2.1. $1)$ By vector we always mean a column vector.

$2)$ We denote by $e$ the vector of length $4$ and all components equal to one. As usual, $I_l$ stands for the identity matrix of size $l$.

$3)$ For two vectors $w^{(1)}$ and $w^{(2)}$ of the same length we denote by $w^{(1)}.w^{(2)}$ the vector with elements $w^{(1)}_iw^{(2)}_i$. Similarly, for a vector $w$ with nonzero elements we write $w.^{-1}$ meaning the vector with coordinates $w_i^{-1}$.

$4)$ For a vector $w$ we denote by ${\operatorname{diag}}(w)$ the diagonal matrix with $({\operatorname{diag}}(w))_{ii} = w_i$.

$5)$ As usual, $\cdot$ denotes the scalar product.

$6)$ As usual, $\nabla$ refers to the gradient with respect to the spatial variable $x$.

Notation 2.2. Throughout the paper we often skip the arguments of coefficient functions.

3.   Modelling

3.1.   Model setup

Model variables Motivated by the multiphase approach developed in ^[19] we view the tumour and its environment as a saturated mixture of several components. We assume these components to be: tumour cells, normal tissue (mainly the extracellular matrix, but also normal cells), necrotic tissue (was not included in ^[19]), and interstitial fluid. The latter, in turn, has several constituents, among which are protons. Depending on their concentration, the environment can be more or less acidic. pH levels drop due to enhanced glycolytic activity of neoplastic cells. We assume that neither cells nor living tissue are produced due to an already very acidic, and hence very unfavourable, environment. Thus we assume that the total volume of the mixture is preserved, the phases only transferring from one into another. The main variables in our models, all depending on time $t\geq0$ and position in space $x\in\Omega\subset \mathbb{R}^N$, $\Omega$ a domain, are:

● vector $u:[0, \infty)\times\Omega\rightarrow [0, 1]^{ 4}$ of volume fractions of

– tumour cells, $u_c$,

– normal tissue, $u_m$,

– necrotic tissue, $u_n$,

– interstitial fluid, $u_w$.

● acidity (concentration of protons), $h:[0, \infty)\times\Omega\rightarrow [0, \infty)$.

Unlike ^[19] we do not require the space dimension, $N\in \mathbb{N}$, to be one. The main goal of this section is to derive equations for the volume fractions based primarily on the conservation laws for mass and momentum as well as additional assumptions on some of the phases. To write down the physical laws, we introduce additional variables that are subsequently eliminated. These are:

● fluxes of the components, $J_i:[0, \infty)\times\Omega\rightarrow \mathbb{R}^{N}$, $i\in \{c, m, n, w\}$,

● common pressure, $p:[0, \infty)\times\Omega\rightarrow \mathbb{R}$.

Model parameters The equations we develop below involve the following set of parameters:

● matrix of drag coefficients associated with each pair of components, $K\in \mathbb{R}^{ 4\times 4}$;

● additional, isotropic pressures by the components, $\tau_i:[0, 1]^{ 4}\times \mathbb{R}\rightarrow \mathbb{R}$, $i\in \{c, m, n, w\}$;

● reaction terms, $f_i: [0, 1]^{ 4}\times \mathbb{R}\rightarrow \mathbb{R}$, $i\in \{c, m, n, w, h\}$;

● diffusion coefficient of the protons, $D_h > 0$.

Assumptions on $K$. We assume throughout that

Assumptions on $f_i$'s. In order to ensure that the sum of all volume fractions in the mixture is always one, we require

Since $u_i$'s and $u_h$ should be nonnegative, we require further that

Possible choices for the reaction terms are:

Example 3.1. $f_c(u_c, u_n, h) = -c_1u_cu_n(h-h_{max})$, $f_m(u_m, h) = -\frac{c_2u_m(h-h_{max})_+}{1+u_m+h}-c_3u_m$, $f_n = -f_c-f_m$, $f_w = 0$. This accounts for the fact that the amount of both tissue and viable cancer cells (the latter being in interaction with the necrotic matter embedded in the acidic environment) decreases when the proton concentration $h$ exceeds a certain maximum threshold, leading to acidosis and hypoxia. The tissue infers degradation due to causes other than direct influence of acidity or tumor cells, and, on the other hand, it experiences a certain amount of self-regeneration, which is limited by the already available tissue and low pH. Thereby, $c_1, c_2, c_3 > 0$ are constant rates; they could, however, also include further dependencies. For the reaction term in the acidity equation we choose e.g., $f_h = -a_1u_wh/(1 + u_wh/h_{max}) + a_2u_c - a_3h$, with $a_i\ge 0$ ($i = 1, 2, 3$) constants. This choice ensures proton uptake by the interstitial fluid (with saturation), production by glioma cells, and decay with rate $a_3$.

Assumptions on $\tau_i$'s. We assume that the non-living components (necrotic matter and extracellular fluid) acquire a tendency to reach equilibrium, so that only living matter exerts additional pressure, thus

The mechanical properties of living matter (tumour cells and normal tissue) are different; among others, they are able to generate both intra- and interphase forces in response to changes in the local environment. The influences include variations in the volume fraction of the cell phase and the presence of chemical cues, of which we account for the local acidity (via proton concentration). The corresponding forces manifest themselves as an additional intraphase pressure. It is assumed that the pressure in the tumour and normal tissue phases increases with their respective densities; moreover, glioma cells exert supplementary (isotropic) pressure on the normal tissue. Taking these effects into account, we choose for the additional pressure terms

where $\alpha_c, \alpha_m, \theta, h_{max}, \chi > 0$ are constants. In (3.3) we combine the influence of local cell mass and acidity, both adding to cancer cell pressure which, in turn, enhances glioma motility. Cell stress increases with growing cell mass, pushing the cells away from overcrowded regions. Stress due to acidity saturates for large acidity levels. Indeed, in highly acidic regions cell ion channels and pH-sensing receptors are mostly occupied, thus making cells insensitive to the presence of protons in their environment. As in ^{[11,24,26,27]}, we call this a repellent pH-tactic behaviour. Similarly, the choice in (3.4) means that there is some intraspecific tissue stress (compression) further accentuated by the interaction between glioma cells and their fibrous environment. The forms proposed in (3.3), (3.4) are reminiscent of those chosen in ^[19].

Main equations To derive our models we will rely on a set of equations describing physical laws. They are as follows.

● Since the the components of $u$ are volume fractions, we have

This is the so-called 'no void' condition.

● Mass conservation for $i$th component is given by

It would be reasonable to assume that this equation should hold for all phases. However, this may come into conflict with additional assumptions, see subsequent Subsection 3.2.

● Conservation of momentum for $i$th component is given by

where

is the stress tensor involving the common pressure $p$ and the additional pressure $\tau_i$. It is assumed that the response of the tumour to stress is elastic and isotropic, which implies that the deformation induced by the applied force on each of the considered cellular and tissue components is limited. An unlimited deformation would correspond to a viscoelastic material ^[20], which is not considered here. The forms of the stress tensors depend on the material properties of each phase and their response to mechanical and chemical cues in the environment. As in ^[19,28] the interstitial fluid is considered to be inert and isotropic and viscous effects within each phase are neglected.

The right hand side in (3.7) represents the force acting on the $i$th phase and, neglecting any inertial effects and exterior body forces, it takes the following form:

The first term on the right hand side in (3.9) accounts as usual (see, e.g., ^[7,19,28]) for the pressure distribution at the interface between phases. The remaining term represents viscous drag between the phases, with drag coefficients $K_{ij}$.

Plugging (3.8) and (3.9) into (3.7) yields an equation for the variables we introduced above:

Equation (3.10) is assumed to hold for all four phases.

● The proton concentration satisfies the reaction-diffusion equation

In the remainder of this section we derive two models based on these laws. To begin with, we simplify Eq (3.10). Given that $p$ and $\tau_i$ are scalar functions, we have for all $i\in \{c, m, n, w\}$

Since $K$ is symmetric, adding together the left-hand sides of (3.12) for all $i\in \{c, m, n, w\}$, we find that

Combining (3.12) with (3.5) and (3.13), we obtain

so that

Plugging (3.14) into (3.12) we arrive at a system which no longer involves $p$:

and it holds that

Next, we introduce a matrix function

The symmetry of $K$ once again implies that

System (3.15) can now be written in the form

where we denote by $J^{(l)}$ the vector made up of the $l$th coordinates of each $J_i$. Due to (3.16) and (3.18) for each $l$ (at least) one equation is redundant. We exploit this in more detail in the next Subsections.

3.2.   A model with an immovable component

Recall that $u_m$ and $u_n$ correspond to normal and necrotic tissues, respectively. A standard modelling assumption is that any tissue is completely immovable. This means that its flux is a zero function, turning the mass conservation law (3.6) into an ODE. However, as we show in this Section, already presupposing, e.g.,

is problematic.

Singling out the $n$th phase, we use the following convenient notation.

Notation 3.2. For a vector $y$ with components corresponding to the four phases we denote by $\tilde y$ the vector which is obtained by removing the component which corresponds to $n$th phase. Similarly, for a matrix $A$ we denote by $\tilde A$ the submatrix which results from removing the row and column of $A$ corresponding to that phase.

Recall that fluxes $J_i$, $i\in \{c, m, n, w\}$, satisfy system (3.19) which is underdetermined. We notice that the submatrix $\tilde A$ is diagonal-dominant due to (3.18) and the nonnegativity of $u_i$'s and $K_{ij}$'s. Using the Gershgorin circle theorem, we infer that the real parts of the eigenvalues of $\tilde A$ do not exceed $-\rho_A(u_i)$, where

Therefore, for $u_n > 0$ the submatrix $\tilde A$ of $A$ is invertible. Consequently, system (3.19)–(3.20) is uniquely solvable.

To simplify the calculations, let us now assume that the following technical assumption holds:

In other words, the symmetric matrix $K$ has rank one. One readily verifies the following formulas:

Lemma 3.3. Let assumption (3.22) hold. Then:

and

so that

Using (3.25) we can resolve system (3.19)–(3.20) with respect to the components of the fluxes. For $l\in\{1, \dots, N\}$ and $i\in\{c, m, w\}$ we obtain

We used (3.5) and (3.2) in the third and fourth equalities, respectively. Plugging (3.26) into (3.6) and using (3.2), we obtain

For variable $u_n$ we have due to (3.6) and (3.20) that it satisfies an ODE:

However, the structure of the fluxes in (3.27) would ensure (3.5) only if

This condition fails to hold in general, unless for $\tau_c$ and $\tau_m$ such that

It is not fulfilled by the coefficients we have in mind, see Subsection 3.1.

Let us assume that (3.27) holds only for $i\in\{c, w\}$. Assume further that

i.e., that the liquid phase is negligible. In this case we obtain a haptotaxis model:

Unless $\tau_c$ is zero for $u_c+u_n = 1$, this model cannot ensure $u_m\geq0$.

Altogether we see that trying to explicitly enforce an immovable phase leads to a situation where not all basic conservation laws can hold at the same time.

3.3.   A model with total flux control

Let us now assume that the number of components in the mixture is three since

thus, unlike previous multiphase models in a related context (see e.g., ^[19,37]), but compare also ^[36], we neglect the interstitial fluid and only take into account the more 'solid' components, namely cancer cells, necrotic matter, and normal tissue, all of which are assumed to be more or less heterogeneously interspersed within the volume of interest. In this Subsection we assume all physical laws from Subsection 3.1 to hold in full, so that, in particular, the mass conservation law (3.6) holds for each $i\in\{c, m, n\}$. Due to (3.1) and (3.5) we can replace (3.6) for $i = n$ by

where

As we have seen previously, system (3.19) is underdetermined. However, if $J_{sum}$ is given, then the system is equivalent to the following matrix equation:

where

and

One can readily verify that

where

Using (3.5), we can rewrite (3.35) as

where

Resolving (3.33) with respect to $J_i$'s and using (3.37), we get

Combining (3.38)–(3.40), we obtain

and

Recalling that $\tau_{n} = 0$ due to (3.2), we conclude from (3.41)–(3.42) that

and, similarly,

Set

Then (3.43) and (3.44) take the form

Finally, combining (3.2)–(3.6), (3.32), (3.49), and (3.50) we arrive at the system

Our construction guarantees that $u_n$ satisfies an equation similar to (3.51a) and (3.51b):

where

and

This a priori ensures that $u_n\geq0$ is satisfied if the solution components are smooth.

Remark 3.4. Model (3.51) is a generalisation of the model derived in ^[19]. Unlike that work we require neither the constants $K_{ij}$ to be equal nor the space dimension $N$ to be one.

Remark 3.5. Our allowing the tissues to be displaced (thus giving up the hypothesis of their immovability) led to nonlinear diffusion and drift in (3.51d) (or, for that matter, in (3.52)). This might seem unusual, as most reaction-difusion-transport systems describing cell motility in tissues assume the latter fixed and use ODEs for their evolution. We emphasize that the obtained terms do not state a self-driven motion of these components, but rather the effect of population- and biochemical pressure exerted thereon.

Remark 3.6. To close model (3.51), one needs to choose some divergence-free $J_{sum}$. In dimension one and for no-flux boundary conditions $J_{sum}\equiv0$ is the only option. In higher dimensions $J_{sum}$ can be a curl field (in $3D$) or, more generally in $N$ dimensions (cf. ^[3]), an exterior product of gradients: $J_{sum} = \nabla \gamma_1\wedge\nabla \gamma _2\wedge \dots \wedge \nabla \gamma _{N-1}$. In the physically relevant $3D$ case this means that $J_{sum} = \nabla \gamma _1\times \nabla \gamma _2$, thus there exist some scalar quantities $\gamma_1, \gamma _2$ (e.g., densities/volume fractions/concentrations) such that the total flux is orthogonal to each of their gradients, or, put in another way, $J_{sum}$ is tangential to a curve which lies in the intersection of the surfaces described by $\gamma_1$ and $\gamma_2$. This would mean that the total flux is following a direction which is equally biased by those two species. The solenoidality of $J_{sum}$ is a reasonable assumption in view of our previous requirement that there were no sources or sinks of material, either.

Diffusion matrix Without loss of generality we assume that

One readily verifies that

where

Due to assumptions (3.53) we have for $0\leq u_c, u_m, u_c+u_m\leq1$ that

In particular, for small $\varepsilon_i$ matrix (3.54) can be regarded as a perturbation of matrix

corresponding to

i.e., to

This case was addressed in ^[22]. Further, we compute

Overall, the diffusion matrix of equations (3.51a) and (3.51b) takes the form

The complete diffusion matrix of (3.51a), (3.51b), and (3.11) includes the third line

4.   Existence of weak solutions to (3.51a)–(3.51c), (3.11)

In this section we use the method that was presented in ^[21] in order to establish an existence result for our cross diffusion system (3.51a)–(3.51c), (3.11). The key to applying the method is finding a suitable so-called entropy density. Motivated by the study in ^[22], where the case of equal $K_{ij}$'s and no acidity was treated, we consider the following entropy density:

Here $L$ is the well-known logarithmic entropy and $a > 0$ is a sufficiently large constant yet to be fixed. For the subsequent computations we need the matrix of second-order partial derivatives of $E$:

{In order to be able to apply} the method from ^[21]{, we need to ensure positive (semi-)definiteness of} matrix $(D^2E)M$ in ${\cal D}$ where $M$ is the diffusion matrix of (3.51a), (3.15b), and (3.11). In this Subsection we verify this property for the parameter values satisfying the following conditions:

and

where $\varepsilon_c$ and $\varepsilon_m$ are constants defined in (3.55).

Remark 4.1. Since both functions which need to be minimised in (4.3) and (4.4) are fourth degree polynomials in $y$, these conditions can be readily checked numerically for a set of given parameters.

Lemma 4.2 (Uniform ellipticity). Let (4.3) and (4.4) hold. Let $E$ be as defined in (4.1). Then there exist some constants $a > 0$ and $0 < {C_1}\leq{C_2} < \infty$, such that:

Proof. To begin with, we compute

where

so that

Combining (3.61), (3.62), (4.2), and (4.6), we obtain (arguments of functions are omitted)

We introduce the symmetric matrix

where

Next, we study $\det(\tilde P)$ in ${\cal D}_{cm}$. We compute

Observe that $d_P$ is quadratic with respect to $u_c$, and the coefficient of $u_c^2$ is negative. Consequently, $d_P$ cannot attain its minimum inside ${\cal D}_{cm}\cup((0, 1)\times\{0\})$. It remains to ensure that $d_P$ is positive on the sets $\{0\}\times[0, 1]$ and $\{(u_c, u_m)\in(0, 1]\times[0, 1):\ u_c+u_m = 1\}$. This the case if

and

By Sylvester's criterion, $P$ is positive definite if and only if $\tilde P$ is positive definite and $\det(P) > 0$. Due to (3.58)–(3.59), (4.7), and $g$ Lipschitz all functions involved in (4.9) are bounded and functions ${S_1}, {S_2}, \beta_c, \beta_m$ have positive lower bounds in ${\cal D}$. In particular, $\tilde P_{22}\geq 2\alpha_m > 0$, so that $\tilde P$ is positive definite if and only if $\det(\tilde P) > 0$. Further, we have

where $\varphi:\overline{{\cal D}}\rightarrow \mathbb{R}$ is a bounded function ({recall that $g$ is Lipschitz}). Thus, for sufficiently large $a$, $\det(P) > 0$ holds provided that $\det(\tilde P) > 0$. Altogether, we conclude that conditions (4.10) and (4.11) imply that matrix $P$ is positive for all triples {in} $\overline{{\cal D}}$. Assuming these conditions to be satisfied, let $\lambda(P)$ denote an eigenvalue of $P$. Then

In $\overline{{\cal D}}$, functions ${\operatorname{trace}}(P)$ and $\det(P)$ are bounded from above and from below, respectively, by some positive constants. Hence, we have positive lower and upper bounds for the eigenvalues of $P$.

Now we are ready to state our existence result.

Theorem 4.3. Let $0 < K_{cm}\leq K_{cn}\leq K_{mn}$ and $\alpha_c, \alpha_m, \theta, h_{max}, \chi > 0$ be some constants which satisfy (3.55), (4.3), and (4.4). Let coefficients $\tau_c$ and $\tau_m$ be as defined in (3.3) and (3.4) and let Lipschitz functions $f_c, f_m, f_h: {\overline{{\cal D}}}\rightarrow \mathbb{R}$, with domain ${\cal D}$ as in (4.1a), be such that

Then for every given $J_{sum}\in L^{2}_{loc}([0, \infty); (L^2(\Omega))^n)$ and $(u_{c0}, u_{m0}, h)\in (L^{\infty}(\Omega))^2\times L^2(\Omega)$ such that $(u_{c0}, u_{m0}, h)(x)\in\overline{{\cal D}}$ for a.a. $x\in\Omega$ there exists a weak solution $(u_c, u_m, h):[0, \infty)\times\Omega\rightarrow {\overline{{\cal D}}}$ to system (3.51a)–(3.51c), (3.11) under no-flux boundary conditions. This means that:

a.e. in $(0, \infty)$ for all $\varphi\in H^1(\Omega)^3$, and the initial conditions for each variable are satisfied in $L^2(\Omega)$-sense.

Remark 4.4. Using the assumptions on the coefficients of system (3.51a)–(3.51c), (3.11), one readily verifies that if $(u_c, u_m, h)$ satisfies (4.13) and $(u_c, u_m, h)\in{\cal D}$ a.e., then the gradients $\nabla(u_c\tau_c(u_c, u_m))$ and $\nabla(u_m\tau_m(u_c, u_m))$ as well as the fluxes in (4.14a) and (4.14b) belong to $L^2_{loc}([0, \infty); ({L^2(\Omega)})^N)$.

Proof of 4.3 (sketch). We rely on the theory of weak solvability for cross diffusion systems which was developed in ^[21].

We recall that the main tool of the method in ^[21] is a suitable entropy density. Here we use function $E:{\cal D}\rightarrow \mathbb{R}$ previously defined in (4.1). This function is the sum of the standard logarithmic entropy $L$ for variables $u_1$ and $u_2$ and the quadratic $ah^2$. Combining the corresponding properties of the logarithmic entropy and our previous calculations in this section, we obtain:

${\rm{H1}}:$ $E\in C^2({\cal D})$ and is convex and bounded below by

The derivative

is invertible, the inverse being

${\rm{H2':}}$ For sufficiently large $a$, left multiplying the diffusion matrix $M$ by $D^2E$ yields a uniformly positive definite matrix (see Lemma 4.2).

${\rm{H2'':}}$ The diffusion matrix $M$ is uniformly bounded in ${\cal D}$ (compare (3.61) and (3.62), recall that $g$ is Lipschitz).

${\rm{H3:}}$ $M$ and $f_i$, $i\in\{c, m, h\}$, are continuous mappings and there exists a constant ${C_3}$ such that

Estimate (4.16) is a standard consequence of the assumptions (4.12). Indeed, since $z\mapsto z\ln z$ is negative and bounded below in $(0, 1)$, the part originating from the logarithmic entropy satisfies

for some constant ${C_4} > 0$. Further, since $f_h$ is Lipschitz, we also have

for some constant ${C_5} > 0$. Adding (4.17) and (4.18) together, we obtain

since obviously $\min E = \min L$.

The main result of ^[21], Theorem 2 on existence of bounded weak solutions, cannot be directly applied in our case for two reasons: firstly, apart from diffusion and reaction, equations for $u_c$ and $u_m$ also involve transport in the direction of a given vector-valued function $J_{sum}$; secondly, the domain ${\cal D}$ is not bounded. Were it not for these differences, the above properties H1–H3 would correspond to the hypotheses H1–H3 in ^[21]. Still, the proof of existence of weak solutions can be carried out very similar to the proofs presented in ^[21]. The latter go through the following steps: approximation of the time derivative using the implicit Euler scheme, regularisation of the diffusion operator by adding a higher order differential operator such as $\epsilon(I+(-\Delta)^m)$ for $m > N/2$ and small $\epsilon$, solving the linearised approximation problem using the Lax-Milgram lemma, solving the nonlinear approximation problem using the Leray-Schauder theorem, establishing uniform estimates, and, finally, using the compactness method in order to pass to the limit and solve the original problem. Our case can be handled in the very same way. Indeed, on the one hand the linear transport along $J_{sum}$ is subordinate to diffusion, so that it, e.g., does not hinder the derivation of estimates. On the other hand even though ${\cal D}$ is not bounded, we actually know that for bounded $h_0$ the $h$-component of a solution is a priori bounded on all finite time cylinders. This is a consequence of the standard theory of semilinear parabolic PDEs. General $h_0\in L^2(\Omega)$ can be regularised and a corresponding solution obtained by means of yet one more limit procedure.

We omit further details and refer the interested reader to ^[21] where the complete proofs for very similar cases can be found.

5.   Numerical study

We perform numerical simulations of the system (3.51) supplemented with the PDE (3.11) for the evolution of proton concentration $h$, endowed with no-flux boundary conditions. For simplicity we choose $J_{sum} = 0$. We also perform a nondimensionalisation of the model and use it in our simulations (for its concrete form and for the employed parameters refer to the Appendix). The initial conditions are as follows:

These are illustrated in Figure 1. The problem is set in a square $[0, 1000]\times [0, 1000]$ (in $\mu m$), corresponding to the size of large pseudopalisades ^[5]. For the discretisation we use a method of lines approach. Thereby, the diffusion terms in all involved PDEs are computed by using a standard central difference scheme, while a first order upwind scheme is employed for the advection terms occurring in the equations for glioma and normal tissue. For the acidity equation we discretise time upon using an implicit-explicit (IMEX) method, with forward and backward Euler schemes for the diffusion and reaction terms, respectively. The glioma and normal tissue PDEs are discretised in time by an explicit Euler method.

Figure 2 shows the computed volume fractions of glioma (first column), normal tissue (3rd column), necrotic matter (last column), and acidity concentration (second column) at several times, in a total time span which is relevant for pseudopalisade formation. The typical garland-like structure of glioma pseudopalisades is clearly visible. The tumor cells encircle a highly acidic and necrotic region, while outwards, beyond the glioma ring, the normal tissue remains non-depleted and the acidity decreases.

To investigate the effect of acidity we compare the previous results with those obtained for $\chi = 0$, i.e., in the absence of additional pressure on glioma cells due to acidosis, recall (3.3). The corresponding plots are shown in Figure 3. The effect of repellent pH-taxis can be seen in the difference between the solution components in two situations (with $\chi > 0$, with pH-taxis and acid-influenced motility, respectively with $\chi = 0$). The pseudopalisades develop faster and get thicker and more extensive in the former case, accompanied by enhanced acidity in the proximity of higher cancer cell densities and enhanced normal tissue depletion inside the ring-shaped structures and around the glioma aggregates, which also triggers a higher volume fraction for the necrotic matter. These observations are in line with those obtained in ^[26] by another modeling approach: the repellent pH-taxis is not the driving factor of pseudopalisade formation, but it leads to wider such structures.

Finally, we also study the effect of drag coefficients $K_{cm}, K_{mn}, K_{cn}$ in (3.9). Recall that they represent the resistance inferred by the phases in contact when passing over each other. They are contained in the diffusion and drift coefficients of system (3.51). In the computations for Figure 2 we had considered them to be different, more precisely $K_{cm} < K_{mn} < K_{cn}$. This ensured a relatively easy movement of glioma through normal tissue when compared to the shift over necrotic matter of both normal tissue and cancer cells. Now let all these drag coefficients be equal (as was assumed in ^[19]) and take the previously lowermost value $K_{cm}$, thus reducing the drag between necrotic matter and the other two phases. The differences between the two cases are plotted in Figure 4 and show in the second case an enhanced outward migration of glioma cells, away from the highly acidic area at the core of the pseudopalisade and with a pronounced suppression of normal cells, and emergence of larger necrotic matter. An opposite effect is noticed when letting $K_{cm} = K_{mn} = K_{cn}$, all at the previously highest value $K_{cn}$; we do not show here those results (we refer for them to ^[25]), but rather illustrate in Figure 5 the situation when all drag coefficients equal the previously intermediate value $K_{mn}$. This means that the drag between cancer cells and normal tissue is increased, while they can easier shift over the necrotic matter. As a consequence, the extent of the pseudopalisades is reduced, but the glioma aggregates in the garland-like structures are larger. This is due not only to the higher drag between cancer cells and normal tissue, but also to the lower $K_{cn}$ value, which ensures that the glioma cells can leave faster than previously the acidic area inside the 'ring'. Around the site where the cancer cells were initially more concentrated the necrosis is reduced and the normal tissue is better preserved.

Acknowledgement

P. K. acknowledges funding by the German Academic Exchange Service (DAAD) in form of a PhD scholarship.

Conflict of interest

The authors declare no conflict of interest.

Appendix

Non-dimensionalisation

The following rescalings of the dimensional quantities involved in the obtain equations are considered:

Dropping the tildes for simplicity, system (3.51), (3.11) writes in the nondimensional form as

Parameters

The following table (see Table 1) summarises the parameters used in the simulations of Section 5.