On a mathematical model for cancer invasion with repellent pH-taxis and nonlocal intraspecific interaction

Starting from a mesoscopic description of cell migration and intraspecific interactions we obtain by upscaling an effective reaction-difusion-taxis equation for the cell population density involving spatial nonlocalities in the source term and biasing its motility and growth behavior according to environmental acidity. We prove global existence, uniqueness, and boundedness of a nonnegative solution to a simplified version of the coupled system describing cell and acidity dynamics. A 1D study of pattern formation is performed. Numerical simulations illustrate the qualitative behavior of solutions.


Introduction
Migration, proliferation, and differentiation of cells are influenced by biochemical and biophysical characteristics of their surroundings, which they perceive by way of transmembrane units like ion channels, receptors, etc. Increasing experimental evidence suggests that cells are able to sense such cues not only where they are, but also at larger distances, up to several cell diameters around their current position [24,30,47].This led to mathematical models accounting for various types of nonlocalities, most of them addressing cell-cell and/or cell-matrix adhesions; we refer to the review article [9] and references therein.The settings typically involve reaction, diffusion and drift terms, whereby the latter contain an integral operator to characterize the so-called adhesion velocity over the interaction range.In [18] was performed a rigorous passage from a cell-matrix adhesion model to a reaction-diffusion-haptotaxis equation when the sensing radius is becoming infinitesimally small, thus recovering the local PDE formulation from that featuring the mentioned nonlocality.The remote sensing of signals by cells affects, however, not only motility, but also proliferation, growth, and phenotypic switch, either directly -by occupancy of transmembrane units on cellular extensions like cytonemes and folopodia and subsequently initiated signaling pathways, or in an indirect manner -as effects of altered migratory and aggregation behavior.Models involving reaction-diffusion equations with nonlocal source terms have been proposed in various contexts, including biological and ecological ones, see e.g.[28,54] and references therein for rather generic settings, [4,5,46] for chemotaxis systems, and [41,49,50] for equations dedicated to tumor growth.We refer to [9,28,54] for some reviews of model classes addressing this type of nonlocality.
As far as growth and migration of cell populations are concerned, the reaction-diffusion models with nonlocal source terms u t " ∇ ¨pD∇uq `F puq (1.1) typically feature F puq " µJ ˚up1 ´uq to describe nonlocal stimulation of growth (see e.g.[50,54]), or F puq " µu α p1´J ˚uβ q, which characterizes competition between (bunches of) cells for available resources in their surroundings, attempting, e.g., to prevent overcrowding.In the context of (tumor) cell migration such models have been handled e.g. in [49], where intra-and interspecific nonlocal interactions led to an ODE-PDE system for the interplay between cancer cells peforming linear diffusion and haptotaxis with the extracellular matrix being (nonlocally) degraded by the cells and remodeled with the mentioned growth limitation.We also refer to [9,35] for short reviews of models with source terms of this type and therewith associated mathematical challenges.
In this note we propose and analyze a model for tumor cell migration involving myopic diffusion, repellent pH-taxis, and a nonlocal source term of the competition type mentioned above.The cross-diffusion system is obtained upon starting from the mesoscopic description of cell migration via a kinetic transport equation for the space-time distribution function of cells sharing some velocity regime.An appropriate upscaling relying on diffusion dominance then leads to the effective macroscopic equation for the cancer cell density, with precisely specified diffusion and drift coefficients.The remaining of this paper is structured as follows: Section 2 contains the model deduction with the mentioned upscaling.Section 3 is dedicated to the mathematical analysis of the obtained nonlocal macroscopic system, in terms of global existence, uniqueness, and boundedness of a solution to a simplified version of the problem.In section 4 we study the asymptotic behavior.Section 5 offers a 1D study of pattern formation for the equations handled in Section 3, but only involving constant motility coefficients.In Section 6 we provide numerical simulations to illustrate the qualitative behavior of solutions to the investigated nonlocal problem.Section 7 contains a discussion of the results.

Modeling
In this section we start from a mesoscopic description of cell migration and intrapopulation interactions and deduce (in a non-rigorous way) effective equations on the macroscopic scale of cell population dynamics.The deduction closely follows that in [35], however extends it, by accounting here for the repellent effects of acidity eventually leading on the population scale to chemorepellent pH-taxis.
Tumor migration and spread are typically assessed on the macroscopic scale of the cancer cell population via biomedical imaging.The involved processes are, however, highly complex and originate at the lower levels of cell aggregates sharing -beside time-space dynamics-one or several further traits (e.g., velocity, phenotypic state or other so called 'activity variables'), down to microscopic events on individual cells.This multiscale character of cell migration can be captured (at least partially) by models within the kinetic theory of active particles (KTAP) framework formulated by Bellomo et al. (see e.g., [1,3] and references therein).Starting from kinetic transport equations (KTEs), a large variety of (spatially) local and nonlocal models have been proposed and various kinds of upscaling and moment closure methods have been performed in order to deduce their macroscopic limits which enable a mathematically more efficient handling, see e.g.[7, 8, 10-14, 17, 20-22, 26, 31, 32, 37, 59].The obtained macroscopic equations carry in the coefficients of their motility and source terms some of the traits from the mesoscale on which KTEs were formulated.Those coefficients are no longer 'guessed' as in the case of stating reactiondiffusion-taxis directly on the population level and the diffusion is often of the 'myopic' type, involving a drift correction.We will perform here a diffusion-dominated upscaling of mesoscale dynamics.We will use the following notations: • p " ppt, x, vq: distribution function of cells having at time t and position x P R n the velocity v P V ; • V " rs 1 , s 2 sˆS n´1 : velocity space.Thereby, s 1 , s 2 denote the minimum, respectively the maximum speed of a cell, θ P S n´1 represents the cell direction; • upt, xq " ş V ppt, x, vq dv: macroscopic cell density; • hpt, xq: concentration of protons.This is a macroscopic quantity throughout this note.

The kinetic transport equation (KTE)
characterizes the mesoscopic dynamics of the considered cell population.This is the framework set in [43], which assumes that changes in p are due to velocity jumps accompanied by reorientations dictated by a turning kernel contained in the operator Lrhs.
The first term on the right hand side of (2.1) represents the so-called turning operator.The second term describes growth/decay of cells due to intraspecific proliferative/competitive interactions, while μ ą 0 is the constant interaction rate. 1 With a small constant ϵ ą 0 relating to the cell size and to the distance at which cells can sense signals in their proximity, we will assume that μ " ϵ 2 µ.This means that cells have a much higher preference to motility (in particular, to changing direction) than to interaction and crowding.
We assume that the turning operator is of the form with the turning rate T rhspv, v 1 q ě 0 chosen such that the reorientation is a Poisson process with rate λrhs " hence such that T rhs{λrhs is a kernel giving the probability density for a change of the velocity regime of a cell from v 1 to v.In particular, this means that Lrhs is preserving mass.The reorientation of cells depends on the acidity of their environment (expressed by the concentration h of protons).
In the following we assume that the turning rate has an asymptotic expansion of the form thus the turning operator admits itself an expansion where L 0 rhs and L 1 rhs are linear operators, For I we consider as in [35] the form where: α, β ą 0 are constants, Jpx, x 1 q is a function weighting the interactions between (bunches of) cells sharing the same velocity regime within a bounded domain Ω Ă R n .We assume that J depends on the distance between interacting (clusters of) cells and take Jpx, x 1 q " Jpx ´x1 q, also requiring J to satisfy J ě η for some η ą 0. (2.8) We also assume that there exists a bounded velocity distribution M px, vq ą 0 such that: ş V vM px, vq dv " 0, i.e. the flow produced by the equilibrium distribution M pvq vanishes.3. The rate T 0 rhspv, v 1 q satisfies the detailed balance equation 4. The turning rate T 0 rhspv, v 1 q is bounded and there exists σ ą 0 such that T 0 rhspv, v 1 q ě σM px, vq, for all pv, v 1 q P V ˆV, x P R n , t ą 0.
The following lemma summarizing the properties of the operator ´L0 can be easily verified (see e.g.[2,7]).
Lemma 2.1.Let L 0 rhs be the operator defined in (2.5).Then ´L0 rhs has the following properties: (i) ´L0 rhs is positive definite w.r.t. the scalar product and the associated norm in the weighted space L 2 pV, dv M px,vq q, and self-adjoint: for all p, ζ P L 2 pV, dv M px,vq q it holds that (ii) For ϕ P L2 pV, dv M px,vq q, the equation L 0 rhspζq " ϕ has a unique solution ζ P L 2 pV, dv M px,vq q satisfying 2 ζ " 0 iff φ " 0.
Example 2.2.Consider T 0 rhspv, v 1 q :" λ 0 rhsM pvq, with λrhs ě λ 0 rhs ą 0 for any h.This obviously satisfies the properties 3. and 4. in our above assumption.With this choice, L 0 rhsppq " λ 0 rhspM pvqu ´pq (2.9) and it is straightforward to see that this operator satisfies the properties in Lemma 2.1 and the function ψ in (iv) becomes ψpvq " ´vM pvq{λ 0 rhs if ψ P pspan pM pvqqq K .3 Equation (2.1) is supplemented with the macroscopic PDE for proton concentration: where D H ą 0 is the diffusion constant and gpu, hq represents production by tumor cells and uptake (e.g., by blood capillaries -not explicitly modeled in this note) or decay.We also consider initial conditions for p and h: pp0, x, vq " p 0 px, vq, hp0, xq " h 0 pxq, x P Ω Ď R n , v P V. (2.11) Together with these, equations (2.1),(2.10)form a meso-macro system describing the dynamics of the (mesoscopic) cell distribution in response to acidity in the extracellular space.
We perform a parabolic scaling to obtain the diffusion limit of the KTE (2.1).This means that we rescale the time and space variables as follows: t " ϵ 2 t, x " ϵx.
Subsequently we will drop the 'ˆ' symbol and the ϵ-dependency of the solution p ϵ to the resulting KTE, in oder not to complicate the writing.Then, (2.1) becomes ϵp t `v ¨∇x p " 1 ϵ Lrhsp `µϵIrp, ps. (2.12) Now consider the decomposition (Chapman-Enskog expansion) ppt, x, vq " F puqpt, x, vq `ϵp K pt, x, vq, (2.13) with ş V p K pt, x, vq dv " 0, thus p K P pspan pM pvqqq K , and F puq P span pM pvqq such that ş V F puq dv " u.A natural choice is F puqpt, x, vq :" M px, vqupt, xq, which we will subsequently adopt.

14)
Let P : L 2 pV, dv M pvq q Ñ Ker L 0 rhs be the projection operator.Then It is easy to verify that the following lemma holds (see, e.g., [2]).
Lemma 2.3.The projection operator P has the following properties: (i) pI ´P qpM pvquq " P pp K q " 0.
If we now apply I ´P to (2.14) we get For the first transport term on the left hand side we compute where we applied the observations made at the end of Example 2.2 and denoted by the diffusion tensor of tumor cells.
The first term on the right hand side of (2.21) represents (myopic) diffusion, the second one characterizes repellent chemotaxis, away from increasing gradients of proton concentration6 , while the last is a source term accounting for tumor cell growth enhanced or limited by intraspecific interactions.
The above deduction of a macroscopic reaction-diffusion-taxis is merely formal; the nonlinear source term prevents applying the proof of the rigorous derivation from [7].The following section will be dedicated to proving global existence and boundedness of nonnegative solutions to the coupled PDE system for u and h obtained on the macrolevel by considering the above much simplified forms of the coefficient functions λ 0 , a, b.The previous calculations were made for x P R n , however we can restrict to a bounded domain Ω Ă R n upon proceeding as in [14,17,45] and assuming no flux of cells or protons through the boundary.

Mathematical analysis
Let Ω Ă R n be a bounded domain with smooth enough boundary and outer unit normal ν.We consider the model where u denotes the cell density and h the acid concentration.Here, the convolution over Ω is as usually given by J ˚uβ px, tq " ż Ω Jpx ´yqu β py, tq dy.
For our diffusion tensor D " pd ij q i,j"1,...,n , we assume that d ij P C 1 p Ωq.Moreover, D satisfies the uniform parabolicity and boundedness condition, i.e. there are B 1 , B 2 ą 0 such that for all ξ P R n and x P Ω it holds that Additionally, we assume that for x P BΩ and ξ P R n with ξ ¨ν " 0 on BΩ and |ξ| ‰ 0 it holds that The exponents α, β ě 1 satisfy (as in [35]) On the remaining functions and parameters we make the subsequent assumptions: • u 0 P Cp Ωq and u 0 ě 0, • h 0 P W 1,8 pΩq and 0 ď h 0 ď H, h 0 ı H, where H is a positive constant, • µ is Lipschitz-continuous with constant L µ , satisfying 0 ď µ and µphq ě δ ą 0 for h ď H, • g P C 1 pR 0 ˆR0 q with ∇g P pL 8 pR 0 ˆR0 qq 2 , 0 ď gpu, 0q ď G and gpu, Hq ď 0 for u P R 0 , • J P L p pB diampΩq p0qq for some p P p1, 8q, 0 ă η ď J, • D H ą 0.
By convention the term K i ą 0 denotes a positive constant for all i P N (or, respectively, a positive function of its arguments).

Local existence in an approximate problem
The Stone-Weierstraß theorem implies that there is a sequence of diffusion tensors pD l q lPN with D l " pd lij q i,j"1,...,n s.t.d lij P C 2`ϑ p Ωq for ϑ P p0, 1q and D l Ñ D in C 1 p Ωq nˆn for l Ñ 8.Moreover, D l satisfies (3.3) and the uniform parabolicity condition for all l P N, i.e. there are 0 ă such that for all ξ P R n , x P Ω and l P N it holds that For l P N we consider the approximate problem x P Ω, t ą 0, pD l pxq∇u l `∇ ¨Dl pxqu l `Dl pxqu l ∇h l q ¨ν " ∇h l ¨ν " 0, x P BΩ, t ą 0, u l px, 0q " u 0 pxq, hpx, 0q " h 0 pxq, x P Ω.
(3.6) Lemma 3.1.For all l P N there are T max ą 0 and a weak solution pu l , h l q of (3.6) such that for all T P p0, T max q it holds that u l P Cp Ωˆr0, T sqqXL 2 p0, T ; H 1 pΩqq and 7 h l P Cp Ωˆr0, T sqXL 8 p0, T ; W 1,8 pΩqqX W 2,1 2 pΩ ˆp0, T qq and pu, hq satisfies for a.e.t P p0, T q and all η P W 1,1 2 pΩ ˆp0, T qq it holds that p∇ ¨Dl u l `Dl ∇u l `Dl u l ∇h l q ¨∇η dx ds 7 W k,l p pΩ ˆp0, T qq denotes the Sobolev space of functions having weak derivatives in L p pΩ ˆp0, T qq, namely up to order k w.r.t.space and up to order l w.r.t.time " and B t h l " D H ∆h l `gpu l , h l q a.e. in Ω ˆp0, T max q (3.9) ∇h l ¨ν " 0 a.e. in BΩ ˆp0, T max q, (3.10) It holds either T max " 8 or T max ă 8 and Proof.Fix l P N. Due to the Stone-Weierstrass theorem there is a sequence pu 0k q kPN Ă C 0,1 p Ωq, u 0k ě 0 with limit u 0 in Cp Ωq.We set M :" sup kPN }u 0k } L 8 pΩq ă 8.For h ă 0 and ū ě 0 extend the coefficients by gpū, hq :" 2gpū, 0q ´gpū, ´hq and µphq :" µp´hq.
We show the existence of a solution pu lk , h lk q of (3.6) with initial value u 0k instead of u 0 in the sense of (3.7) and (3.9) for k P N by showing the existence of a fixed point of the operator F introduced below similarly to [51].Namely, we define for some small enough T ą 0 the set S :" tū P L 8 pΩ ˆp0, T qq : 0 ď ū ď M `1 a.e. in Ω ˆp0, T qu.

$ ' & ' %
B t h lk " D H ∆h lk `gpū, h lk q, x P Ω, t P p0, T q, ∇h lk ¨ν " 0, x P BΩ, t P p0, T q, h lk px, 0q " h 0 pxq, x P Ω. (3.14) Here, T can be chosen independent of ū and k.Through a fixed point argument similar to [27], we conclude that there is a unique function h lk P L 8 p0, T ; W 1,8 pΩqq that satisfies Moreover, h lk is the unique weak solution of (3.14) in the sense that for a.e.t P p0, T q and all η P W 1,1 2 pΩ ˆp0, T qq it holds that As in Lemma 3.3 below and due to Theorem IV.9.1 (and the remark at the end of that section) in [33], it follows that and Hence, h lk solves (3.14) in the sense of (3.9).Moreover, the continuity of h lk follows from Theorem 4 in [16] due to W 2,1 2 pΩ ˆp0, T qq Ă Cp0, T ; L 2 pΩqq and the embedding of W 1,8 pΩq into some Hölder space on Ω.Now, Theorems III.5.1 and 7.1 in [33] (that also hold for our no-flux boundary condition), Theorem 4 in [16], and Gronwall's inequality imply that there is a unique u lk in the space S X C κ, κ 2 p Ω ˆr0, T sq, such that it solves (3.13) in the sense of (3.7) (with u 0k instead of u 0 ) and satisfies for some κ P p0, 1q.Note that C 1 , C 2 and C 3 are independent from ū and k.Hence, the operator where u lk solves (3.13) for ū in the sense of (3.7), is well-defined.Moreover, as C κ, κ 2 p Ω ˆr0, T sq ãÑãÑ Cp Ω ˆr0, T sq due to the Arzelà-Ascoli theorem, F maps bounded sets on precompact ones.To apply the Leray-Schauder theorem it remains to show that F is closed and, consequently, a compact operator.Let We want to show that F pūq " u lk .
Let h lkm be the solution of (3.14) that corresponds to ūm for m P N. Due to (3.17) we conclude from the Lions-Aubin lemma and the Banach-Alaoglu theorem that there are h lk P W 2,1 2 pΩ ˆp0, T qq and subsequences Therefore, due to (3.21) and the Lipschitz-continuity of g for a.e.t P p0, T q and all η P W 1,1 2 pΩ ˆp0, T qq it holds that h lk is a solution of (3.14) in the sense of (3.9).From (3.7) we conclude as in [33] that for a.e.t P p0, T q and all η P W Inserting this into (3.22) and using (3.16),Young's inequality, the continuity of h lk and the Lipschitzcontinuity of µ, we conclude that for a.e.t P p0, T q it holds 1 2 }u lkm p¨, tq} 2 From Gronwall's inequality we obtain a constant C 6 pl, k, T q ą 0 such that }∇u lkm } L 2 p0,T ;Ωq ď C 6 pl, k, T q for all m P N. 8 Hence, the Banach-Alaoglu theorem implies that (by switching to a subsequence, if necessary) Then, (3.19) -(3.21), (3.24) and the dominated convergence theorem imply that u lk is a solution of (3.13) in the sense of (3.7), and therefore F pūq " u lk and F is a compact operator.Consequently, by a 8 The majority of subsequent constants will depend on T , but we will omit it in the writing.
Leray-Schauder argument we obtain the existence of a fixed point u lk of F , that satisfies for a.e.t P p0, T q and all η P W 1,1 2 pΩq the weak formulation (3.7) for u 0 replaced by u 0k .Now, (3.18) and the compact embedding of C κ, κ 2 p Ω ˆr0, T sq in Cp Ω ˆr0, T sq imply that there is a convergent subsequence of pu lk q k such that Then, with the same arguments as before we obtain the desired weak solution pu l , h l q of (3.6).
Finally, for such pair property (3.12) follows from a standard extensibility argument.Theorem 3.2.There is T max P p0, 8s and a unique solution pu l , h l q of (3.6) with 0 ď u l and 0 ď h l ăH, u l , h l P Cp Ω ˆr0, T max qq X C 2,1 p Ω ˆp0, T max qq. (3.25) Proof. 1. Regularity: Let l P N, 0 ă T 1 ă T max and consider the weak solution pu l , h l q from Lemma 3.1.
The boundedness and nonnegativity of h l (0 ď h l ă H) follow from the comparison principle of the semilinear heat equation with Neumann boundary condition and our assumptions on g.
3.2 Global existence and boundedness of u in the approximate problem Lemma 3.3.It holds that Proof.With Lemma 1.3 from [55] and (3.15) we estimate for t P p0, T max q that ˙e´λ1D H pt´sq }gpu l , h l q} L 8 ˙ds holds for q P p1, 8q, where λ 1 is the first eigenvalue of ´∆ on Ω with Neumann boundary condition.
Using the properties of g and the boundedness of h, we obtain that for t P pT 0 , T max q it holds that }∇h l } L q pΩq ďC 13 e ´λ1DH t }∇h 0 } L q pΩq `C13 p}B h g} L 8 H `Gq ˙e´λ1D H pt´sq ds Consequently, We will show the global boundedness of u as in the proof of Theorem 1.1 in [35].
Lemma 3.4.It holds that u l P L 8 p0, T max ; L q pΩqq for q P r1, 8q.
Proof.Let q ě maxt1, β `α ´1u.Due to (3.25) the terms in the estimates below are well-defined for a.e.t P p0, T max q.Multiplying the first equation of (3.6) by qu q´1 l , integrating over Ω and using partial integration, we obtain Using the uniform parabolicity of D, we estimate Further, due to Young's inequality we obtain the estimate Inserting these estimates into (3.27) and using our assumptions on µ and J, the boundedness of h, and it follows that where Adding qC 14 pq, lq}u l } q L q on both sides of (3.29) and using Young's inequality one more time, we obtain d dt Similarly to Step 1 in the proof of Theorem 1.1 in [35] it follows that 2qC 14 pq, lq Here, C S ą 0 denotes the Sobolev embedding constant from page 8 in [35], C P ą 0 the constant from the Poincaré inequality, and (3.31) Hence, for t P p0, T max q we conclude that d dt }u l } q L q pΩq `qC 14 pq, lq}u} q L q pΩq ď 2qC 14 pq, lqpC 15 pq, lq `|Ω|q. (3.32) Hence, for t P p0, T max q and q ě maxtβ `α ´1, 1u we obtain from (3.32) the upper bound }u l p¨, tq} L q pΩq ď q b 2C 15 pq, lq `|Ω|p2 `}u 0 } q L 8 pΩq q. (3.33) Remark 3.5.As in [35] we cannot directly conclude from Lemma 3.4 that u l is bounded on Ωˆp0, T max q as lim qÑ8 q b 2C 15 pq, lq `|Ω|p2 `}u 0 } q L 8 pΩq q " 8.
Theorem 3.6.For all l P N there is a unique bounded and nonnegative solution pu l , h l q of (3.6) consisting of nonnegative functions u l , h l P Cp Ω ˆr0, 8qq X C 2,1 p Ω ˆp0, 8qq and h l ă H. Thereby, there is some C 18 ą 0 that does not depend on l s.t.u l ď C 18 .Moreover, if Ω is convex for some 'small' enough choice of parameters specified in (3.38) and (3.39) below, then for any K ą 1 where for some function G to be adequately chosen and we set s from (3.31) equal to 8 for n " 2.
Proof.Let l P N. We proceed with a Moser iteration as in Step 2 in Theorem 1.1.in [35].Set q k :" 2 k `a with a :" 2ps´1qpα´1q s´2 for k P N and s as in (3.31).Analogously to [35] we obtain for t P p0, T max q the estimate where As D l tends to D in C 1 p Ωq nˆn , there is C 21 ą 0 s.t.}D l } C 1 p Ωq nˆn ď C 21 for all l P Nq.We can further estimate that For k P N and t P p0, T max q we set y k ptq :" }u l p¨, tq} q k L q k .
Inserting this into (3.35)we obtain Moreover, we estimate that Hence, from Lemma 2.1 in [35] it follows that for m ě 1 and t P p0, T max q it holds that ż Ω u q k l dx ď p4C 22 q 2 k´m`1 2 2s s´2 p2p2 k´m ´1q`m2 k´m`1 ´kq max Consequently, for t P p0, T max q and m ě 1 it holds that }u l } L 8 pΩq " lim kÑ8 }u l } L q k pΩq ď p4C Consequently, u l is bounded on Ω ˆr0, T max q.Combining this with the boundedness of h l , Lemma 3.3 and (3.12) in Theorem 3.2, T max " 8 follows.
If Ω is convex we proceed as in Step 3 of [35].First, we fix some m and choose our parameters sufficiently 'small' such that it holds that and This depends on our choice of D, µ, g, }∇h 0 } L 8 , and D H . Consequently, we conclude as in [35] that for any K ą 1 we find 'small' enough parameters (satisfying (3.38)
Proof.Let φ P H 1 pΩq.Obviously, for a.e.T ą 0 and each l P N the function u l satisfies ż Ω B t u l φ dx " ´żΩ pD l ∇u l `∇ ¨Dl u l `Dl u l ∇h l q ¨∇φ dx `żΩ µph l qu α l p1 ´J ˚uβ l q dx, (3.40) as it is a classical solution.Due to Theorem IV.9.1 (and the remark at the end of that chapter) in [33], h l satisfies where C 25 ą 0 is independent from l due to the properties of g and h l ăH for all l P N.
Setting φ " u l in (3.40) and using Hölder's and Young's inequalities, the facts that }D l } C 1 p Ωq nˆn ď C 21 , and the uniform boundedness of u l from Lemma 3.6, we can estimate that 1 2 Consequently, from Gronwall's inequality follows Similarly it follows that Putting this together with the uniform boundedness of u l and (3.41), the Lions-Aubin and Banach-Alaoglu theorems imply that there are u P Cp0, T ; L 2 pΩqq X L 2 p0, T ; H 1 pΩqq with B t u P L 2 p0, T ; pH 1 pΩqq ˚q and h P Cp0, T ; H 1 pΩqq X W 2,1 2 pΩ ˆp0, T qq s.t.(after switching to a subsequence if necessary) in L 2 pΩ ˆp0, T qq and pointwise a.e., (3.42) , T ; H 1 pΩqq and pointwise a.e., From this, the dominated convergence theorem, the uniform boundedness of u l , h l and ∇h l from Lemma 3.3 and the Lipshitz-continuity of µ and g, it follows that pu, hq solves (3.1) in the required sense.
The a.e.boundedness and nonnegativity of u and h follow from the pointwise convergence and the uniform boundedness and nonnegativity of u l and h l .Uniqueness follows similarly to Theorem 3.2.

Long time behavior
We consider the long time behavior of our solution under the additional assumptions that we make from now on: • the domain Ω is convex, • the parameters satisfy (3.38) and (3.39), • u 0 ı 0, • we extend J by 0 to R n zB diampΩq p0q and assume }J} L 1 pR n q " 1, • there are h ˚P r0, Hs and constants C H ą 0 and C U ě 0 s.t.
We assume that D is in the closure of M in the pC 1 p Ωqq nˆn -norm and that the sequence pD l q lPN from above is the sequence in M that approaches D.
We proceed by combining the methods from [35] and [31].Proof.Let l P N. We conclude from the strong maximum principle and the assumption u 0 ı 0 that u l ą 0 holds in Ω ˆp0, 8q.
As in [35] we define apsq :" s β ´1 β lnpsq ´1 β with apsq ě 0 for s P p0, 8q.By multiplying the equation for u l in (3.6) by u β´1 l ´u´1 l , integrating over Ω and using partial integration and our additional assumptions on D l , we obtain where using again partial integration Hence, we can estimate using the positivity of u l that Then, we extend u l by 0 to R n zΩ and proceed similarly to the proof of Proposition 3.1 in [35] to obtain using Hölder's inequality that where the interval on the right hand side is nonempty due to our assumptions on C A and C B .
Moreover, due to the convexity of Ω and using the uniform boundedness of pu l q l by U we can estimate Now, inserting this in (4.4), using our assumptions on J and µ, and the uniform boundedness of pu l q l , we conclude Inserting this into (4.3), it follows that Multiplying the equation for h l by h l ´h˚a nd using (4.1), we obtain 1 2 Further, we multiply this by C 30 :" and add it to (4.5) to obtain Due to our assumptions on C A , C B and our choice of ε it holds that p1 ´εqδη|Ω| ´CU C 30 ą 0. Now, Gronwall's inequality implies (4.2).
Combining this with the uniform continuity of u l for all l P N we get the convergence stated in (4.8).
Proof.For h and in the case α " 1 convergence follows directly from Theorem 4.3.We consider the case α ą 1.From (3.42) we know that (after switching to a subsequence if necessary) u l converges to u pointwise a.e. in Ω ˆp0, 8q.Hence, for all px, tq the sequence pu l px, tqq l is Cauchy and we can conclude that for large enough l the sequence u l ptq converges to the same c P t0, 1u.Hence, for a.e.px, tq and large enough l it holds that |upx, tq ´c| ď |upx, tq ´ul px, tq| `|u l px, tq ´c| Ñ 0 for t Ñ 8 and l Ñ 8 due to (3.42) and Theorem 4.3.

Pattern formation: a 1D study
We want to investigate pattern formation in our model (see [39]).For this aim we adapt some of the assumptions on our functions and parameters: • J P L 1 pRq, Jpxq " Jp´xq for x P R and ş R Jpxq dx " 1 and J f upxq :" ş R Jpx ´yqupyq dy, whereas we drop the condition that 0 ă η ă J; • there is exactly one h ˚ą 0 with gp1, h ˚q " 0, moreover, for this h ˚it holds that µph ˚q ą 0, B u gp1, h ˚q ě 0 and B h gp1, h ˚q ă 0. This means that when the cancer cells are at their carrying capacity (corresponding to an acidity level h ˚), the production of protons is increasing with the cell mass and decreasing with enhancing proton concentration.Indeed, crowded tumor cells are highly hypoxic, and a too acidic environment leads to quiescence or necrosis, thus reducing proton expression.Moreover, we assume that µ 1 ph ˚q ă 0, thus the growth rate is decreasing with the proton concentration in the neighborhood of the critical value h ˚.
• w.l.o.g.we consider Ω " r´a, as for a P R.

Stability in the local model without diffusion and taxis
We start by establishing the equilibria of the non-spatial local model that corresponds to (5.1), i.e. # B t u " µphqu α p1 ´uβ q, B t h " gpu, hq. (5. 2) The biologically more interesting one is given by pu ˚, h ˚q " p1, h ˚q, where h ˚is the unique solution of gp1, hq " 0. The corresponding characteristic equation of the Jacobian in p1, h ˚q is given by λ 2 `pβµph ˚q ´Bh gp1, h ˚qqλ ´βµph ˚qB h gp1, h ˚q " 0.
The corresponding eigenvalues are λ 1 " ´βµph ˚q and λ 2 " B h gp1, h ˚q and both have negative real parts due to the assumption B h gp1, h ˚q ă 0. Hence, the steady state p1, h ˚q is stable in this case.

Stability in the local model with diffusion and taxis
We continue by adding again the diffusion and taxis term to the local model (5.2).Adapting the ansatz from [44] we consider perturbations of p1, h ˚q of the form u " 1 `εūpkq and h " h ˚`ε hpkq, where ūpkq " ũe λpkqt cospkxq and hpkq " he λpkqt cospkxq for ũ, h P R and wavenumber k P N and |ε| ăă 1.
Here, λpkq denotes some eigenvalue of the corresponding characteristic equation.As in [42] we use the fact that e ikx `e´ikx 2 " cospkxq to ensure that our perturbations are real.Inserting these u and h into our model and linearizing about the steady state p1, h ˚q, we obtain # λpkqū " ´k2 dū ´dk 2h ´βµph ˚qū, λpkq h " ´DH k 2h `Bu gp1, h ˚qū `Bh gp1, h ˚qh . (5. 3) The corresponding eigenvalues are given by λ 1,2 pkq " trpJ u,h pkqq ˘atrpJ u,h pkqq 2 ´4 detpJ u,h pkqq 2 , where we denote by J u,h the Jacobian of the right hand side in system (5.3) at p1, h˚q and its determinant and trace are, respectively, given by trpJ u,h pkqq " ´pd `DH qk ´βµph ˚qB h gp1, h ˚q ą 0.
Hence, the equilibrium p1, h ˚q is stable.The local model does not lead to any Turing type patterns.

Stability in the nonlocal model
We consider u and h as in the previous section and linearize the convolution term about p1, h ˚q similarly to [44].Hence, inserting u in the convolution term and using the symmetry of J, we compute that Here, Ĵ denotes the Fourier transform of J. Hence, linearizing system (5.1),we obtain # λpkqū " ´dk 2 ū ´dk 2h ´βµph ˚qp2πq 1 2 Ĵpkqū, λpkq h " ´DH k 2h `Bu gp1, h ˚qū `Bh gp1, h ˚qh . (5.4) The corresponding eigenvalues are as above given by λ 1,2 pkq " trpJ u,h pkqq ˘atrpJ u,h pkqq 2 ´4 detpJ u,h pkqq 2 , where we denote by J u,h the Jacobian of the right hand side in (5.4) at p1, h ˚q and its trace and determinant are given by trpJ u,h qpkq " ´pd `DH qk 2 ´βµph ˚qp2πq The sign of the real part of the eigenvalues is ambiguous here and depends especially on the sign of Ĵpkq, which depends on k.As above, we have stability here if trJ u,h pkq ă 0 and det J u,h pkq ą 0 (5.5) for all k " π a z, where z P Z.We make this restriction due to our boundary condition u x " h x " 0. Now, we are looking for a critical k c (that is not necessarily of the form π a z) depending on our choice of parameters, where we distinguish as in [42] the occurrence of Turing instabilities in the case Impλpk c qq " 0 for some arbitrary critical k c , Hopf instabilities in the case Impλp0qq ‰ 0, and wave instabilities in the case Impλpk c qq ‰ 0 for some critical k c ‰ 0. If Ĵ is symmetric it suffices to consider only positive k c .
A Turing bifurcation can occur if we find k c such that detpJ u,h qpk c q " 0 and trpJ u,h qpk c q ă 0. Now, rewriting these conditions we conclude that the equality Ĵpk c q " ´d k 2 c βp2πq ă Ĵpk c q (5.7) have to hold for one or several critical k c in a set K c , whereas (5.5) holds for all k R K c that are of the form π n z.Such k c exist depending on the choice of parameters, on the functions µ and g, and especially on the sign of the Fourier transform of J.Moreover, due to our assumptions the terms on the right-hand side of (5.6) and on the left-hand side of (5.7) are negative and tend to ´8 for k Ñ ˘8.On the other hand, a Hopf or a wave instability can occur if we find k c such that trpJ u,h qpk c q " 0 and det J u,h pk c q ą 0, whereas (5.5) holds for all k that do not satisfy this and are of the form π n z.Hence, a Hopf instability occurs if " Ĵp0q and Ĵp0q ą 0, whereas (5.5) holds for all k ‰ 0. On the other hand, a wave instability occurs if holds for one or several k c ‰ 0, whereas (5.5) holds for all other k that are of the form π n z and do not satisfy the above equality and inequality.
From the above considerations we conclude that the occurrence of a Turing, Hopf or wave instability depends on the concrete choice of J, as we need to find suitable k of the form π a z.If the Fourier transform Ĵ is nonnegative, no Turing patterns occur.
Remark 5.1.If there is a steady state of the form p0, h ˚˚q for some h ˚˚ą 0 and B h gp0, h ˚˚q ď 0, then this equilibrium is stable in the case with diffusion, taxis and nonlocal term.If, on the other hand, B h gp0, h ˚˚q ą 0, this steady state is unstable already in the case without diffusion and taxis.This case is, however, unrealistic for the biological problem investigated here.Indeed, the proton expression by hypoxic cells is much reduced and there must be at least some very weak acid buffering, lest all cells (and surrounding tissue) become apoptotic.
Likewise, the steady state p1, h ˚q is unstable already in the case without diffusion and taxis if B h gp1, h ˚q ą 0. This situation may occur at least in a transient manner, e.g. when the cells can still extrude protons while their environment is quite acidic and if the cells are at their carrying capacity and the proton buffering is relatively low.That can lead, e.g., to a choice of the form gpu, hq " u `uh ´γh 2 with γ ď 4{5.

Numerical simulations
In this section we perform numerical simulations of system (5.1), in order to illustrate the solution behavior.The equations are discretized by using the algorithm in [40]; the motility terms were discretized with finite differences (centered for the diffusion, upwind for the drift).The initial conditions are as in [35]: u 0 pxq " # e ´px´x l q 2 , for x l ă x ď 0 e ´x2 l p1 ´x xr q, for 0 ă x ď x r , with x l " ´5, x r " 5.
In a first test we took β " 1, µ " 1, along with the logistic kernel Jpxq " 1{p2 `ex `e´x q (see, e.g., [34]) and the uniform kernel Jpx; ρq " 1 2ρ 1 r´ρ,ρs .The first two columns of Figure 1 show simulation results for α " 2, which is the 'limit value' in (3.4).The solution ceased (in finite time) to exist for sufficiently large α in each of these situations (α " 6.25 and α " 8.2, respectively), u exhibiting strong aggregation near the initial bulk of cells, cf. last two columns in Figure 1.This behavior was also observed for increasing values of µ, with the difference of singularities already occuring for smaller α values.Increasing the values of µ and β leads to patterns, the shape of which depends decisively on the interaction kernel J and also on the values of α and d. Figure 2 shows 1D space-time patterns of the cell density u for β " 20, µ " 100, and several combinations of α and J.The results for the proton concentration h are not shown, as there are only small quantitative differences between the respective cases.Figure 2 suggests that, irrespective of the chosen kernel9 , higher cooperative intraspecific interactions (larger α values) or slower diffusion delay the invasion of cells in the whole region, leading instead to enhanced proliferation.On the long run the cells tend to fill the whole space and remain at their carrying capacity.This behavior endorses the results in Section 4 and is particularly well visible for the logistic kernel, which satisfies all conditions in the proofs of the theoretical results of Sections 3 and 4; the process is much slower when a uniform kernel is used, however it has eventually the same outcome.The last row in Figure 2 exhibits the situation of a cell diffusion which is much slower than that of protons.The effect is a delayed filling of the space with cells (and produced protons) and a later formation of the patterns observed in the upper rows.The asymptotic behavior is similar, only it takes longer for the solution to reach the respective states.
To assess the effect of nonlocality we performed simulations with the source term in the u-equation of (5.1) replaced by µphqu α p1 ´uβ q.The results are shown in Figure 3.The first two columns illustrate the case with the same source term for proton concentration as above, namely gpu, hq " up1 ´hq, for which no patterns seem to develop (we tried several combinations of parameters, including those used for the patterns in Figure 2).In fact, decreasing the value of ρ in the uniform kernel Jpx; ρq eventually leads to the local version of the system.The plots in the leftmost column were produced with d " D H , while those in the middle column used d !D H .The behavior of u and h is the same, with the difference of the second case inferring a slower spread of cells and protons.The last column in Figure 2 already shows the tendency of disappearing patterns when approaching the local case.The last column of Figure 3 shows the case where the source term in the h-equation is replaced by gpu, hq " u `uh ´γh 2 , as proposed in Remark 5.1. 10o patterns for u were observed for the local model, which, together with the simulations performed for intermediary values of ρ, suggests that the patterns are driven by the nonlocality of cell-cell interactions,

Discussion
In this note we investigated a model describing pH-tactic behavior of cells with nonlocal source terms.As such, this work is extending the one in [35], which studied the Fisher-KPP equation with nonlocal intraspecific competition with various powers of the solution.In contrast to [35] we handled here a problem in a bounded domain, and the population dynamics was coupled to that of the proton concentration, which also led to a taxis term.The proof of our results concerning global well-posedness and long time behavior relied, however, to a substantial extent on the methods in [35].We also dealt here with space-dependent tensor coefficients in the motility terms, which involve myopic rather than Fickian diffusion.The dissipative effect of the repellent pH-taxis contributed to reducing some of the difficulties in the analysis -as long as the required conditions on the functions involved in the system are satisfied.
Among the relatively few existing models with nonlocal source terms, the one in [49] is closely related, however it features several differences: the cells perform attractive haptotaxis towards gradients of extracellular matrix (ECM), the nonlocal source terms are contained in both equations, do not involve any powers, and the Fickian diffusion of cells has a constant coefficient.Our model requires less regularity for the interaction kernel and the motility coefficients involve a tensor and are more general.On the other hand, the nonexploding solution behavior is favorized in our case by repellent chemotaxis.We also provided an informal model deduction and an assessment of the long time solution behavior.The analysis done in [41] for a model with standard motility and with nonlocal source terms as in [49], but with one or two species performing chemotaxis towards the same attractant imposes certain requirements on the forcing term of the latter, mainly in order to obtain the asymptotic behavior of the cell-related solution components.Our condition (4.1) imposed for similar purposes on the source term of the tactic signal looks rather differently.The attraction-repulsion chemotaxis models considered in [46] have closer similarities with our setting, as far as the nonlocal intraspecific interactions are concerned.Major differences occur through our system only featuring two equations, in the source terms of the chemical cues, and in the motility terms: the latter involve in our case the space-dependent tensor Dpxq and myopic diffusion, while the nonlocal reaction term in the proton dynamics is more general.We also prove an explicit long time behavior of both solution components and provide a short analysis of space-time patterns (in 1D), along with numerical simulations.
Our preliminary analysis in Section 5 and the simulation results in Section 6 suggest that patterns occur only in the nonlocal model, are not of Turing type, and seem to be driven by the nonlocal source terms and influenced by the chosen kernel and the combination of parameters in the nonlocal term.This is in line with the pattern behavior observed in [35] and with other works concerning reaction-diffusion problems with nonlocal intra-and/or interspecific competition, cf.e.g.[23,25,44,48,53,58].Those works involved more or less similar source terms and no taxis, however the repellent pH-taxis contained in our model does not seem to have a relevant influence on the patterns.
Open problems relate to a thorough study of patterns depending on the interplay between the parameters α, β, µ and the influence of the kernel J.Moreover, the well-posedness, asymptotic and blow-up behavior, along with patterning are largely unknown in the case of a degenerating motility tensor -the less so in combination with myopic diffusion and/or other types of taxis.Indeed, these can lead in the local case to very complex issues even in 1D, as shown e.g. in [56,57].

Figure 3 :
Figure 3: Simulation results for (5.1) with local source term µphqu α p1 ´uβ q replacing the one in the equation for u.Left and middle column: gpu, hq " up1 ´hq with d " D H and d !D H , respectively.Right column: gpu, hq " u `uh ´γh 2 , d " D H .