ON A COUPLED SDE-PDE SYSTEM MODELING ACID-MEDIATED TUMOR INVASION

. We propose and analyze a multiscale model for acid-mediated tumor invasion accounting for stochastic eﬀects on the subcellular level. The setting involves a PDE of reaction-diﬀusion-taxis type describing the evolution of the tumor cell density, the movement being directed towards pH gradients in the local microenvironment, which is coupled to a PDE-SDE system charac-terizing the dynamics of extracellular and intracellular proton concentrations, respectively. The global well-posedness of the model is shown and numerical simulations are performed in order to illustrate the solution behavior.

1. Introduction. The analytical and numerical theory for stochastic differential equations (SDEs) have strongly developed since the 1940s, and there are many reference books collecting the main achievements in the field. Several results have also been obtained for systems of SDEs, including couplings of forward and backward SDEs and their connections to systems of quasilinear parabolic PDEs, see e.g., [33,25,6]. Stochastic PDEs (SPDEs) have inferred, too, an unprecedented progress during the last decades, although the mathematical tools are less developed than for SDEs. In contrast, references concerned with coupled PDE-SDE systems are rather scarce. In [17] the authors considered a system of a 1D parabolic PDE coupled to an SDE in a half-local way, letting the variable satisfying the SDE intervene in the PDE only by its expectation. The possibility of a coupling via compatibility conditions at the mutual interface between the different domains on which a PDE and an SDE are respectively stated was mentioned as well, however without further dwelling on it. The emphasis of [17] was on numerics.
Models describing gene expression in the zebrafish hindbrain have been introduced in [42] and involve a system of SPDEs with spatio-temporal noise, coupled to two SDEs. The one-way coupling therein allows to reduce the complexity of the problem; for the numerical simulations the setting was further simplified to avoid the simultaneous noise occurrence in the (S)PDEs and (S)ODEs components of the system.
Except for the simpler model in [20] we are not aware of other works handling a genuine coupling between PDEs and SDEs, in which the SDE holds at each point where the PDE does. To our knowledge models featuring the interplay between reaction-diffusion-taxis PDEs and SDEs have not been addressed so far, neither numerically nor analytically. Motivated by a biomedical problem, we propose here such a model and investigate the existence of solutions to the considered system.
One of the hallmarks of cancer is the upregulation of glycolysis, both in aerobic and hypoxic conditions [13]. Among its most prominent consequences is the extrusion of excess protons, triggering the acidification of the extracellular region, see e.g. [39]. This seems to confer tumor cells several advantages: the intracellular region being at the alkaline side of neutrality promotes advancement through the cell cycle towards mitosis, along with evasion of apoptosis and cytoskeletal remodeling [5,40], hence facilitating migration. The invasion into adjacent tissue is further fostered by extracellular matrix (ECM) degradation and induced apoptosis of stroma cells [40]. Several modeling approaches have been proposed and investigated in order to describe acid-mediated tumor growth and invasion. Most of them are continuum deterministic settings involving ordinary and/or partial differential equations, the latter being used whenever spatial effects are accounted for. The PDEs are mainly of reaction-diffusion type, like the model in [9], which seems to be the first to address this biological issue in such mathematical framework. Quite a few works followed, concerning the mathematical analysis of the PDEs involved in the model and some of its extensions and modifications, see e.g., [7,10,30,27,35]. All these models describe the spatio-temporal interaction between tumor cells and normal tissue under the influence of extracellular proton dynamics, [27] also involving the evolution of matrix degrading enzymes. The PDEs are set on the macroscopic scale of concentrations and cell densities depending on time and space only. Biological evidence indicates, however, that subcellular processes (including intracellular proton production, buffering, and transport into the extracellular space) regulate and are in turn influenced by the macroscale dynamics, see e.g. [22,38,40]. This calls for multiscale formulations in which this dependence is addressed, at least in some of its aspects. Deterministic models connecting subcellular level dynamics with those on the macroscale have been proposed and analyzed e.g., in [29,28,36,37]. Of these, [28,36] specifically refer to acid-mediated tumor invasion described by way of pH-taxis. The latter characterizes the biased motion of cancer cells in the direction of an extracellular pH gradient [2,32].
Being inherent to most biological processes, stochasticity is a relevant feature, in particular on the level of individual cells and also of the subcellular dynamics, see e.g., [8,22]. Recently we considered in [20] a two-scale model with nonlocal sample dependence describing the proton dynamics in a tumor, where the intracellular one is governed by an SDE that is coupled to a reaction-diffusion equation for the macroscopic concentration of extracellular protons. The models in [15,14] have a multiscale character, as well; they couple random ODEs with PDEs of reaction-(cross)diffusion-taxis type and show the relevance of stochasticity in explaining transiently observed phenomena like hypocellular gaps between the tumor and the surrounding normal tissue, further infiltrative growth patterns, or tumor aggressivity depending on cell phenotype switching. In this work we consider a model connecting the subcellular scale (dynamics of intracellular protons, described by an SDE) with the macroscopic one (tumor cell density and extracellular proton concentration, each described by a reaction-diffusion PDE -the one for tumor cells also including pH-taxis).
The paper is organized as follows: In Section 2 we introduce the multiscale model for acid-mediated cancer migration. Section 3 is dedicated to the statement of the well-posedness result and its proof. Numerical simulations are presented in Section 4, and in Section 5 we conclude with some comments concerning the results and further issues to be addressed in future research.

2.1.
Modeling the intracellular proton dynamics on the microscale. We describe the dynamics of intracellular proton concentration by an Itô SDE where W denotes a standard scalar Wiener process and dW the corresponding Itô differential. The production of protons by glycolysis is described by the function ϕ, a constant decay rate δ h models the loss of protons by intracellular buffering, while the exchange of protons between the intra-and extracellular regions is determined by the function ψ. Moreover, the function g determines the volatility of the random perturbation.
2.2. Evolution of cancer cells and extracellular protons on the macroscale. Intracellular protons are expulsed into the extracellular region through several membrane-based proton transporters upregulated as a consequence of enhanced glycolysis [40,34]. 1 On the macrolevel, the extracellular protons are transported by diffusion and buffered or uptaken by vasculature. The invasive capacity of the tumor is influenced by the level of the extracellular proton concentration and local cancer cell density. Moreover, the movement of malignant cells is biased towards gradients of the extracellular proton concentration (pH taxis), and the death and growth of cancer cells also depends on the local cell density and concentration of extracellular protons.
Hence, the evolution of c (density of cancer cells) and p (concentration of extracellular protons) is governed by the system of reaction-diffusion(-taxis) equations where ∆ denotes the Laplace operator with respect to the spatial variable and ∇ the gradient. The nonlinear diffusion of tumor cells is described with the aid of the density-dependent coefficient Φ(c, p), also featuring the dependence on proton concentration available in the peritumoral region and favoring invasion. The convective effect of proton gradients is described by the taxis term, the coefficient 1 e.g., NDCBE (Na + dependent Cl − -HCO − 3 exchanger), NHE (Na + and H + exchanger) and AE (Cl − -HCO − 3 or anion exchanger) are specific transporters present on the cell membrane Ψ(c, p) representing pH-tactic sensitivity. The proliferation of cancer cells is described by a logistic function, and δ c denotes the death rate. The extracellular protons are dissolved in the domain and transported by diffusion with a constant diffusion coefficient d. The function ψ determines the proton exchange between the intra-and extracellular regions and the constant decay rate δ p models the loss of protons through buffering and vasculature.
The micro-macro model for acid mediated tumor invasion (1)-(2) extends the stochastic proton dynamics model [20] by taking into account the dynamics of cancer cells, thereby also involving pH-taxis. It is completed with appropriate boundary and initial data, which will be specified in the following section.
3.1. Basic notation and functional spaces. We will use the notation R + 0 := [0, ∞). Moreover, ∂ t denotes the partial derivative with respect to time t > 0, and ∆, ∇, and ∇· are the Laplace, gradient, and divergence operator with respect to the spatial variable x, respectively. The spatial variable belongs to D, where D is a given smooth spatial domain in R n for some n ∈ N. On the boundary of D, ∂ ν is defined as the derivative with respect to the outward unit normal. Furthermore, let (W (t)) t≥0 be a standard scalar Wiener process defined on the filtered probability space (Ω, F, P). We assume that F = (F t ) t≥0 is the usual completion of the natural filtration of (W (t)) t≥0 . The corresponding Itô differential is denoted by dW , and as usual, E denotes the expectation value.
For a Banach space X we denote by L p ((Ω, F t ); X) the L p -space of measurable functions w.r.t. the σ-algebra F t . Moreover, X k , k ∈ N, represents the product space X × · · · × X (k-times) with the usual norm.
To simplify the notation when dealing with purely PDE (SDE) properties which hold P-a.s. in Ω (a.e. in D), we sometimes drop the dependence upon variable ω (variable x) and write, for example, c instead of c(·, ω, ·) (c(·, ·, x)). Moreover, for a stochastic process u : [0, T ] × Ω → V we write u(t) instead of u(t, ·). Thereby, V denotes any space of functions defined for x ∈ D.

3.2.
Problem setting and main result. Motivated by stochastic micro-macro models for acid mediated cancer invasion we consider the following coupled system of two PDEs (with their initial and boundary conditions) P-a.s.
and an Itô SDE where T > 0 and D is a smooth bounded domain in R n . The diffusion and pH-taxis coefficients are given by the functions Φ and Ψ, and We make the following assumptions: , and there exist some constants K 1 , K 2 > 0 such that Furthermore, it holds that (c) h 0 ∈ L p ((Ω, F 0 ); W β,p (D)) for some p ∈ (2, ∞) and β ∈ (0, 1) such that β > n p . Remark 1. Typical choices for the diffusion and pH-taxis coefficients Φ and Ψ such that conditions (5) are satisfied are e.g., (see e.g., [15,29,36] and Section 4) where the constants κ,α 1 ,α 2 are positive and such thatα 2 > 1. Thereby, the diffusivity of tumor cells is enhanced (up to a certain upper limit) by their interactions with an increasingly acidic environment; this is in line with the actually observed behavior of cells whose motion into the normal tissue is favorized by the acidity degrading the latter, see e.g., [27,35] and references therein. The subsequent analysis also holds, however, if instead of the condition on Φ(c, p) in (5) we require A class of functions satisfying such condition is where the constants α 1 , α 2 are positive and such that α 1 ≤ α 2 . This choice accounts again for a (limited) enhancement of cancer cell diffusivity by acidity, but here this enhancement is primarily due to acidity, involving the interactions between cells and protons only to impose a limit on it (thus we have in this situation a mainly 'acidity-driven' diffusivity). We will use this choice in the numerical simulations in Section 4, as it leads to more interesting solution patterns.
holds P-a.s for all t ∈ [0, T ] and x ∈ D. If (c, p, h) is a local solution for all T > 0, then we call it a global solution of (3)-(4).

Proof of Theorem 3.2.
In this Section we prove the global existence result Theorem 3.2. We use the Banach fixed point theorem in order to obtain local wellposedness, but prior to that we show that each solution necessarily satisfies certain a priori estimates. The latter allows us to extend a local solution to a global one. Finally, we prove that solutions depend continuously on the initial data on each finite time interval, which yields the overall uniqueness of solutions.
In order to apply Banach's fixed point theorem, we need to decouple our system. For this purpose, we first study the PDE system (3) assuming that (c 0 , p 0 ) and h are given, and then study the SDE (4) for given h 0 and (c, p).

Remark 4 (Notation).
In the sequel C i denotes (for all indices i) a non-negative constant. The dependence of such quantity upon the parameters of the problem, i.e.: the space dimension n, the domain D, the constants K 1 , K 2 , K 3 , α, β, p, and the structure of the coefficient functions Φ, Ψ, f 1 , f 2 , f 3 (meaning their norms and the norms of their derivatives on compact sets) is subsequently not indicated in an explicit way.
Step 1 (Well-posedness for (3) with respect to (c, p)). Assume that (c 0 , p 0 ) satisfies the regularity and compatibility assumptions (3a) and (3b) in Assumptions 1, and that h fulfils condition 5 in Definition 3.1 for some T > 0. As observed in Remark 3, (3) is a weakly coupled reaction-diffusion-transport system. The semilinear equation (3b) together with the corresponding boundary and initial conditions can be solved with respect to p. The latter can then be plugged into equation (3a) in order to obtain c. Equations (3b) and (3a) are standard semi-and quasilinear parabolic equations, respectively. Equation (3a) is in divergence form. Moreover, the coefficient functions Φ, Ψ are smooth, the functions f 1 , f 2 are continuous, and both h and (c 0 , p 0 ) are Hölder continuous. In this situation, we may apply standard theory of semi-and quasilinear parabolic PDEs (see [21]) and the regularity result [23, Theorem 1.2], which yield the existence, for each ω ∈ Ω, of a classical solution The solution is unique (also among weak solutions) and satisfies conditions 2 and 3 of Definition 3.1, provided that c(·, ω, ·) is a priori non-negative and bounded and 0 ≤ p(·, ω, ·) ≤ 1 in [0, T ] × D. The latter is checked in the subsequent step.

Remark 5 (Weak solution).
If h is not assumed to belong to some Hölder class and only satisfies 0 ≤ h ≤ 1, the classical theory based on parabolic maximal regularity (see [21] and [1]) yields the existence of a unique weak solution (c, p) to (3). As in the case of a classical solution, uniform a priori bounds are required.
Step 2 (A priori estimates for (c, p) as a solution of (3)). Assume that (c 0 , p 0 ) satisfies the regularity and compatibility conditions (3a) and (3b) from Assumptions 1 with and assume that h ∈ L ∞ (Ω; L ∞ ((0, T ) × D)) for some T > 0 (which ensures that h is measurable in (0, T ) × D) P-a.s.). Let (c, p, h) be the pathwise (w.r.t. ω ∈ Ω) classical or weak solution (in the standard PDE sense) to (3). The existence and uniqueness of such solution was discussed in Step 1. Equation (3b) is a semilinear parabolic PDE. Due to (7b) it follows by the parabolic comparison principle that Moreover, the weak maximum principle implies that Using the assumptions on Φ, Ψ, f 1 and c 0 and estimates (12)-(15), we apply once again the regularity result [23, Theorem 1.2] and thus obtain that In particular, due to (13) and (17) we have that Further, standard inner regularity results for quasilinear parabolic PDEs (see [21]) yield that Remark 6. In case of a Hölder continuous h, in (20) one can take ε equal to 1.
Step 3 (Stability of (c, p) as a solution of (3)). Let (c, p, h) and (c,p,h) be two pathwise (w.r.t. ω ∈ Ω) classical solutions (in the standard PDE sense) to (3) for some T > 0 with initial data satisfying the regularity and compatibility conditions (3a) and (3b) from Assumptions 1 and The difference of (c, p) and (c,p) satisfies the initial-/boundary value problem Multiplying (21a) by (c −c) and integrating by parts over D yields Using the hypotheses on Φ, Ψ, f 1 and (18), we obtain where we used Hölder's and Young's inequality in the last step. Multiplying (21b) by (p −p), integrating by parts over D, and using the assumptions on f 2 , it follows that 1 2 Moreover, if we multiply (21b) by ∆(p −p) and integrate, we obtain by the assumptions on f 2 that 1 2 , where we used again Hölder's and Young's inequality. Consequently, it follows that d dt Combining this estimate with (22) and (23) and taking expectation values yields Step 4 (Adaptivity of (c, p) as a solution of (3)). We begin with the case when is a simple càdlàg process for some T > 0. More precisely, the latter means that there exist a partition 0 = t 0 < t 1 < · · · < t n = T , n ∈ N, of the interval [0, T ], a family of disjoint sets Ω ij ∈ F ti , and functions h ij ∈ L ∞ (D) for where 1 denotes the indicator function. Let now t ∈ (0, T ) be arbitrary. Then, there exists 1 ≤ l < n such that t ∈ [t l , t l+1 ), and since Building disjoint intersections, we can replace {Ω ij } 1≤j≤Ni,1≤i≤l by a family of (3) is also a simple function w.r.t. the σ-algebra F t . Since t ∈ (0, T ) was arbitrary, this implies that (c, p) Assume now that (c 0 , p 0 ) satisfies the regularity and compatibility assumptions (3a) and (3b) from Assumptions 1, and that h satisfies condition 5 from Definition 3.1 for some T > 0. Let (c, p) be the corresponding classical solution (its existence and uniqueness were established in Step 1). Then, there exist: a sequence (c n0 , p n0 ) n∈N of simple functions w.r.t the σ-algebra F 0 such that and a sequence of simple càdlàg processes (h n ) n∈N such that If (c n , p n ) denotes the solution corresponding to (c n0 , p n0 ) and h n , then, by the above arguments, it follows that (c n , p n ) : [0, T ] × Ω → [C 1+α (D)] 2 is adapted. Due to the estimate (24), we have that Consequently, for each t ∈ [0, T ] there is a subsequence such that (without relabelling) Moreover, due to the Sobolev interpolation inequality and estimate (20) (take ε ∈ (α, 1)), we conclude with (26) that Therefore, for each t ∈ [0, T ], (c, p)(t) : Ω → [C 1+α (D)] 2 is measurable w.r.t. the σalgebra F t as it is the limit of a P-a.s. converging sequence of simple functions w.r.t. the σ-algebra F t . This means that the process (c, p) Step 5 (Well-posedness for the SDE (4)). Assume that (c, p) satisfies conditions 1 and 2 from We observe next that for every fixed x ∈ D the SDE (4) satisfies due to assumptions on f 3 and g the hypotheses of the stochastic invariance criterion (see [31] or [4]) and, consequently, it holds a priori that Altogether, the classical result on the existence of solutions for an Itô SDE (see, e.g. Theorem 1, p. 50 in [11]) is applicable in this situation and yields for every x ∈ D a unique solution h; in particular, h(·, ·, x) is a continuous adapted process for all x ∈ D.
Step 6 (A priori estimates for h as solution of the SDE (4)). Assume that h 0 satisfies (3c) from Assumptions 1, and that (c, p) satisfies conditions 2 and 3 from Definition 3.1 and estimate (18) for some T, R > 0. Moreover, let h fulfil the stochastic integral equation (11) for all x ∈ D. In this part of the proof, we exploit the regularity of h. In order to shorten notations, we introduce the difference operator for each pair (x, y) ∈ D × D and a function u defined in D. For the difference D x,y h we obtain Taking the expectation value and p-th power (p > 2, from the assumptions on h 0 ) in (29), it follows from Hölder's inequality, a version of the Burkholder-Gundy-Davis inequality [26,Theorem 7.1,p. 39], the assumptions on f 3 and g, and (18) that Applying Gronwall's lemma yields Further, dividing both sides of (31) by |x − y| βp+n and integrating over D × D, we arrive at the estimate Combining (18), the assumptions on h 0 , and the Sobolev embedding C 1 (D) → W β,p (D) (recall that β ∈ (0, 1)), we conclude from (32) that for all t ∈ [0, T ]. Since 0 ≤ h ≤ 1, (33) yields (recall the definition of the Sobolev-Slobodeckij norm for W β,p (D)) that for all t ∈ [0, T ]. Similarly, the difference D x,y (h(t) − h(s)), for 0 ≤ s < t ≤ T , satisfies the stochastic integral equation Using (33) Combining (34) and (36), we arrive at Step 7 (Adaptivity of h as solution to (4)). Assume that h 0 satisfies (3c) from Assumptions 1, and that (c, p) satisfies conditions 1 and 2 from Definition 3.1 for some T > 0. Let h be the solution of the SDE (4) which was defined pointwise for every x ∈ D (the existence and uniqueness of such solution was discussed in Step 5).
For each x ∈ D and t ∈ (0, T ], the random variable h(t, ·, x) can be approximated using the standard Euler-Maruyama scheme (see, e.g., [19, Chapter 9 §1]): Since strong and weak limits coincide, we conclude that (without relabelling the subsequence) it holds that Finally, as weak convergence preserves measurability in Bochner spaces of functions with values in a separable Banach space (and W β,p (D) is separable since p ∈ (2, ∞)), (40) and (43) yield that h : [0, T ] × Ω → W β,p (D) is an adapted process.
In order to construct a corresponding local solution, we use the Banach fixed point argument. To this end, we define the following metric space for an arbitrary T > 0:  We also make use of the following set: with the constants C 2 (T, R) and C 4 (T, R) as in estimates (12)- (13) and (16)- (17), respectively. The results obtained in the preceding steps imply that the solution operators θ 1 for (3) (h given) and θ 2 for (4) ((c, p) given), are well-defined and continuous. Let us now define Combining (24) and (45), we obtain Observe that the particular choice of the 'starting time' was not essential for what we obtained heretofore. Indeed, all previous results remain valid if we replace the interval [0, T ] by [t 0 , t 0 +T ] for any t 0 > 0 and consider the 'shifted' in time filtration (F t+t0 ) t≥0 instead of the original one. It is also clear from the above steps that if a solution defined for t ∈ [0, t 0 ] for some t 0 is prolonged by a solution defined for [t 0 , t 1 ] for some t 1 > t 0 , both solutions being understood in terms of Definition 3.1, then the resulting solution on [0, t 1 ] is also a solution in terms Definition 3.1. We now use these considerations in order to construct a solution which exists globally. Let T, R > 0 be arbitrary and set Step 10 (Uniqueness of solutions of (3)-(4) and continuous dependence on initial data). Let (c, p, h) and (c,p,h) be solutions to system (3)-(4) with Then, estimates (24) and (45) hold. Adding them and using Gronwall's inequality, we obtain the following Lipschitz-type estimate The uniqueness of solutions to (3)-(4) is a direct consequence of (46). This completes the proof of Theorem 3.2.  Table 1.
Before describing the discretization scheme we provide the initial conditions and explicit forms of the functions and coefficients in (3), (4) used for the numerical simulations. The former are illustrated in Figure 1, both for the 1D and the 2D cases. The proton concentrations and tumor cell density are scaled, thus the units are arbitrary.
The parameters chosen for the 1D and 2D simulations are given in Table 2    Concerning the set of functions and coefficients involved in the PDE-SDE system, we first choose some satisfying the conditions for which the well-posedness proof in Section 3 works: It can be easily checked that these functions satisfy the Assumptions 1 (with the condition on Φ in (5) replaced by (9)).
The results of the 1D simulations are plotted in Figure 2; it shows different time snapshots of the same solution sample path. Figure 3 Table 2. Simulation parameters (1D and 2D) the intracellular proton concentration exhibits oscillations due to the stochastic effects, however only infers a higher variation (increase) at the front of the tumor wave, after which it stabilizes again. This less interesting dynamics is due to the choice of functions in (47), however other choices are possible as well -although for them the proof in Section 3 no longer holds. An alternative proof, however, could possibly allow for more general forms (we will provide more comments on this issue in Section 5).
In the following we propose different functions and coefficients which account in a more pronounced way for the nonlinear couplings in the system and which lead to more realistic tumor patterns. Thereby, we also allow the functions f i (i = 1, 2, 3) to depend on all variables c, p, and h, while the volatility coefficient g in the SDE for intracellular protons can depend (beside h) also on c, to model the direct effect of tumor cell density on the perturbations inferred by h.
The proliferation rate of tumor cells is inversely related to the concentration of intracellular protons, and the carrying capacity is reduced by the peritumoral acidity, which both motivate the choice of f 1 . Moreover, we use a function J to weight the proton transport across the cell membrane (modeled by the last terms in f 2 and f 3 ) in such a way that the latter is enhanced during the proliferative regime of tumor cells. This is in line with biological evidence of increased glycolytic activity during the mitotic phase. The first (higher order polynomial) terms in f 2 and f 3 are chosen in order to ensure some multistable dynamics (with 4 or 5 steady states) for p and h, respectively. Figure 4 shows the qualitative behavior of J as a function of c.
The following choices are made for the coefficients in (48): The involved constants can be found in Table 2.
We use an Euler-Maruyama scheme ( [19], Chapter 10) for discretizing the intracellular proton dynamics equation (4). Setting J n k,j := J(c n k,j ), c n k,j := c(t n , x k , y j ), p n k,j := p(t n , x k , y j ), h n k,j := h(t n , x k , y j ),  f + 3 (c n k,j , p n k,j , h n k,j ) : we get for h the following discretization: h n+1 k,j = h n k,j + τ f 3 (c n k,j , p n k,j , h n k,j , h n+1 k,j ) + √ τ g(h n k,j , c n k,j )ξ = h n k,j + τ f + 3 (c n k,j , p n k,j , h n k,j ) − f − 3 (c n k,j , p n k,j , h n k,j , h n+1 k,j ) + √ τ g(h n k,j , c n k,j )ξ, and therefore, The Euler-Maruyama discretization ensures convergence of order 1 2 for τ < 1, while the nonstandard way of expressing the negative terms implicitly ensures positivity of the approximated solutions. Apart from this, if h n k,j ≤ 0 at any given spatial point (x k , y j ) we set h n+1 k,j = 0 at that point. For the extracellular proton dynamics we use an implicit finite difference scheme, wherein the diffusion term is discretized implicitly. In particular, setting we get for p the following discretization: For the cancer cell population dynamics (equation (3a)) we use an IMEX finite difference scheme, wherein an implicit central finite difference scheme is used to discretize the diffusion terms while the taxis and reaction terms are discretized explicitly. Setting f 1,n k,j := f 1 (c n k,j , p n+1 k,j , h n+1 k,j ), Φ n k,j := Φ(c n k,j , p n+1 k,j ), Ψ n k,j := Ψ(c n k,j , p n+1 k,j ). we get for c the following discretization: Figure 5 shows 3 out of 1000 randomly chosen sample paths, each sample path being a solution of (3) whose coefficient functions take the explicit form given in (48). Figure 6 depicts the expectation value of the solution for (3). It was numerically computed by averaging over 1000 sample solutions.
The plots show the more or less strong aggregation of tumor cells at different time points, in regions with high extracellular and low intracellular proton concentrations. This is mainly due to the pH-taxis and the proton exchange through the cell membranes. The weighting of the latter by the function J leads to stronger oscillations in the dynamics of h, thus triggering the formation of cell aggregates.    Although the solution infers very steep increasing and high densities of cells in the aggregates, it does not seem to infer blow-up. Instead, several such high-density aggregates form and the solution remains bounded. This is, however, merely what simulations suggest; a mathematical proof of this conjecture seems to be currently out of reach. Since the observed peaks occur at different time and space points, the expectation curves in Figure 6 maintain some wiggly shapes, hinting on rather irregular patterns for the cancer cells, where high density regions alternate in a quick succession with hypocellular patches. Figure 7 shows 2D simulations for several different sample paths illustrating various patterns both for the tumor cells and the extracellular protons. Since the 1D simulations have been run independently of the 2D ones the sample paths do not correspond and cannot be identified; therefore we show also different time points of the sample solutions than in the 1D case, as we focus on the mentioned patterns. Notice also in this case the sample-to-sample path variability in the solution behavior at the same time points: this applies to the tumor spatial extent, to the pattern (with different localizations of the cell aggregates and different regions of acidity), and to the effective levels of cell density and extracellular proton concentration. The model predicts very heterogeneous tumors, with alternating hyper-and hypocellular regions, with usually higher acidity near the cell aggregates. While time goes on the tumor spread becomes increasingly infiltrative, exhibiting an irregular shape with islands of cells apparently having no connection to the rest of the tumor.
Computing the expectation by averaging over 500 sample solutions in 2D allows to assess the average patterns of cancer cell density and extracellular proton concentration. Several time snapshots are shown in Figure 8, which confirms the biologically well known fact that the highest acidity level is located at the sites with highest tumor cell density. Although averaged over so many sample solutions,    the tumor and acidity patterns exhibit (especially at later times) fractal shapes, transiently forming 'islands' of high density and concentration, respectively.
Next we consider the situation with so-called nonlocal coupling, meaning that the effect of intracellular protons on the dynamics of their extracellular counterparts and on the cancer cell evolution is described not directly by the stochastic process h, but via its expectation E(h). This approach has also been considered in a simpler framework in [20], where -relying on [18]-we called it nonlocal sample dependence. In this nonlocal framework we consider again the system (3)-(4) whose coefficient functions take the explicit forms given in (48), with the exception of the functions f 1 and f 2 taking the following form: where the constants are as given in (48). Figure 9 shows a sample path of the solution, while Figure 10 illustrates the expectation, numerically computed by averaging over 1000 sample solutions. In the nonlocal case the sample paths have a very similar appearance (from path to path), exhibiting-as expected-much smoother curves for h and c than in the previous case accounting for the full stochasticity of h. They all have in common the steep profile of the tumor front, with a rather large accumulation of cell mass at the front side (due to the pH-taxis), followed by smaller local maxima for which the nonlinear diffusion and the proton dynamics play a larger role. The wave of cancer cells eventually invades the whole available space. As in the previous case, for later times we observe cell accumulation in strongly localized, high density aggregates; this is again mainly due to the pH-taxis, but the shape of the involved functions and coefficients (with the multistability caused by f 2 and f 3 and the dependence on J(c)) also plays a role. 2D simulations in the nonlocal case are presented in Figure 11 showing the expectations (for c and p) computed by averaging over 500 sample paths of the solution. As previously in the 2D plots, in order to save computational time the spatial components are each in a smaller interval than in the 1D case. The cell and acidity patterns form quite fast, exhibiting a transient behavior as triggered by the multistability introduced via (48) and (49). Thus, Figure 11 illustrates an initially compact tumor in a highly acidic environment, advancing gradually and changing thereby the proton concentration from high levels at the tumor core towards lower levels surrounded by higher acidity. The acidity distribution is then inverted, as the protons fill the further extending tumor region and the cells change their phenotype from migrating to proliferating or vice versa 2 . Eventually, the tumor cells form 'islands' of very high densities, which usually are highly acidic; however, smaller, less acidic, and apparently disconnected aggregates are possible as well. The very fast accumulation of extremely dense cell aggregates might suggest blowup of the solution to the nonlocal equations considered here, however -as before-a proof confirming or refuting this conjecture is not available.

5.
Discussion. In this work we proposed a micro-macro model for acid-mediated tumor invasion which couples the subcellular scale with the cell population level. It accounts for stochastic effects on the lower scale and for pH-taxis on the macrolevel. For the highly nonlinear SDE-PDE system we proved in Theorem 3.2 the global well-posedness and non-negativity of solutions. The boundedness of one solution component (tumor cell density) in a more general setting is still open. The numerical simulations suggest blow-up for certain choices of the functions and coefficients involved in (3)-(4). Furthermore, we asked these to satisfy Assumptions 1, which are rather restrictive and not fulfilled e.g., by (48) or (49), for which the more   interesting dynamics has been observed, in line with known biological facts. The mathematical model introduced here has a much larger versatility than the proof in Section 3 affords, therefore we also performed the numerical simulations for those more general functions and put in evidence the decisive role played by pH-taxis and stochasticity in reproducing realistic infiltrative growth patterns of acid-mediated tumor invasion.
Alternative proofs relying on more advanced fixed point theorems (as in the deterministic framework) would require compactness arguments which cannot just be translated from the deterministic to the stochastic case. A common approach to well-posedness of parabolic SPDEs relies on some monotonicity properties of the elliptic operator (see e.g., [24] and references therein), however such properties are often not satisfied when dealing with strongly coupled, highly nonlinear systems. Another well studied approach -still within the context of stochastic evolution equations in Hilbert spaces-relies on the theory of semigroups (see e.g., [3]) and requires, too, rather strong smoothness assumptions about the involved operators; in particular, it is not clear how to apply it for SDE-PDE systems including nonlinear diffusion and taxis effects. In [14] we applied such method to a larger system involving nonlinear diffusion and pH-taxis, however coupling PDEs with an ODE and a random ODE; the proof was quite involved, but ensured global well-posedness under the conditions imposed on the coefficients. When only weaker assumptions can be made about the functions and coefficients of the system to be studied in this or related stochastic settings then -as mentioned-an approach relying on compactness and a priori estimates (as in the deterministic framework) would be desirable, since it allows to approximate the highly complex problem at hand with a less complicated one. A first step in this direction has been done in [43].
The model introduced here was motivated by the problem of tumor growth and spread, which led to many challenging deterministic mathematical models, of which we already mentioned a few in the introduction. Extending such models to account for stochastic effects is justified by the necessity to characterize (at least some of) the uncertainties inherent to each living system. As shown by our simulations, there is a large variability between the tumor cell and acidity patterns obtained for different sample paths of the solution. Assuming that each patient has a single tumor (i.e., no metastases) this translates into a substantial interpatient variability, for which a deterministic counterpart of this model is not able to account for. Our stochastic model, however, is able to reproduce 3 the often observed highly infiltrative, irregular patterns of cancer. Although in average the predicted shape of the tumor is more regular, with advancing time the margins of the neoplastic tissue become more and more fractal, as illustrated in Figure 8. Indeed, the cell aggregates (as far as they occur) are formed at different spatial points and at different times; they also have different sizes, and they seem to be distributed over the entire considered spatial domain.
Another related issue to be considered is the effect of tumor cells and acidity on the normal tissue, as the latter acts both as support and hindrance for the tumor cells invading their surroundings. Moreover, it is degraded by the cancer cells (actually by the acidity they cause) but can also be remodeled by them. Including these effects increases the level of detail, but also the complexity of the models. The so-called haptotaxis (motion directed towards the tissue gradient) leads even in the deterministic framework to serious challenges both from an analytical and a numerical point of view, the more so in a multiscale setting (see e.g., [44,45] and [29,37], respectively). Its mathematical investigation in a stochastic context (PDE-ODE-SDE setting) is still to be addressed.