A nonlocal sample dependence SDE-PDE system modeling proton dynamics in a tumor

A nonlocal stochastic model for intra- and extracellular proton dynamics in a tumor is proposed. 
The intracellular dynamics is governed by an SDE coupled to a reaction-diffusion 
equation for the extracellular proton concentration on the macroscale. In a more general context 
the existence and uniqueness of solutions for local and nonlocal 
SDE-PDE systems are established allowing, in particular, to analyze the proton dynamics model both, 
in its local version and the case with nonlocal path dependence. 
Numerical simulations are performed 
to illustrate the behavior of solutions, providing some insights into the effects of randomness on tumor acidity.


1.
Introduction. This work is motivated by modeling the interactions between extracellular and intracellular proton dynamics in the context of tumor growth. Hypoxia is a characteristic of invasive tumors, resulting from an imbalance between oxygen supply and its consumption at the cellular level. The increased glycolysis metabolism in cancer cells leads to acidification of the peritumoral region, hence conferring an advantage against normal cells, which -unlike neoplastic tissue-are known to have a reduced capability of surviving at low pH values, see e.g. [18,22]. The reversed pH gradient between tumors and normal tissue promotes tumor invasion and proliferation [25]. Starting with the model by Gatenby and Gawlinski [6] involving reaction-diffusion equations for the dynamics of extracellular protons in interaction with tumor and normal cell densities, several classes of models extending that setting have been proposed and analyzed, see e.g., [17,22] or [8,9,24,19] for more recent, multiscale approaches coupling the dynamics of intra-and extracellular protons with the evolution of tumor cells and normal tissue. Focusing on where D ⊂ R n , n = 1, 2, 3, is a bounded domain with smooth boundary ∂D and W t is a standard scalar Wiener process.
The model variables are: X, the extracellular proton concentration, and Y , the intracellular proton concentration. Both are normalized w.r.t. the maximum concentration, i.e., X and Y take values within the unit interval [0, 1]. Thus, X t is a deterministic quantity satisfying a reaction-diffusion equation, while Y t is a stochastic process evolving according to an Itô SDE. The evolution of X t is influenced by diffusion, the source term r(X t , E(Y t )) modeling the proton extrusion through the cell membrane into the extracellular environment, and a decay term with a rate α characterizing the loss of extracellular protons by processes other than membrane based transport into the tumor cells (e.g., uptake by vasculature, normal cells, buffering etc.). The dependence of r on E(Y t ) instead of Y t highlights that the effects of H + coming from the intracellular regions of a set of tumors are averaged when considering the proton concentration in the extracellular space. In the intracellular space, however, Y t is seen as a genuine stochastic process influencing as such the proton extrusion. The function ϕ models production of Y by glycolysis (which in cancer cells is much amplified when compared to normal cells and hence non-negligible). The decay term βY t describes the loss of protons by intracellular buffering (e.g., by organelles), and the diffusion coefficient g(Y ) = γY (1 − Y ) quantifies the (stochastic) variability in the production/decay of intracellular H + . It accounts for the Y t values ranging between 0 (complete alkalinization) and 1 (maximum acidification), both bounds being lethal for the cell and hence leading to no variability. Also, a maximum threshold is achieved in the middle of this interval, suggesting that larger spreads of the Y t distribution are not allowed.
For a concrete choice of the membrane-based transport terms we use the setting in [24], which in turn was motivated by the choices in [8,28] using the quantitative information in [1], where rates of H + flux due to NDCBE, NHE, and AE transporters were measured 1 . Hence, the above function r will take the form where a i , i = 1, . . . , 4, are positive constants.
Thus, for the interaction functions we have The space dependent function ζ is positive and bounded by one, and the random variable η can only have realizations within the unit interval, as well.
In the following section we establish the well-posedness for a more general class of models including the above mathematical description of proton dynamics as a particular case. Namely, we consider coupled, nonlocal SDE-RPDE systems of the 2236 PETER E. KLOEDEN, STEFANIE SONNER AND CHRISTINA SURULESCU with notations corresponding to the previous ones. We will first analyze local SDE-RPDE systems in Section 3 and then apply these results to show the well-posedness of system (1).
3. Well-posedness of a local SDE-PDE system.
Existence and uniqueness of the solution. Let (Ω, F, P) be a complete probability space with a normal filtration (F t ) t≥0 , let (W t ) t≥0 be a standard scalar Wiener process and dW t denote the corresponding Itô differential. In this section we prove the well-posedness of stochastic systems of the form where D ⊂ R n , n = 1, 2, 3, is a bounded domain with smooth boundary ∂D. Moreover, ∆ = ∆ x denotes the Laplace operator with respect to the spatial variable x ∈ D, and ∂ ν the outward normal derivative on the boundary. The initial data ζ, η are F 0 -adapted random variables in L 2 (D) such that E ζ 2 L 2 (D) < ∞ and E η 2 L 2 (D) < ∞. We denote by A the operator −∆ in D with homogeneous Neumann boundary conditions and by e −At , t ≥ 0, the analytic semigroup in L 2 (D) generated by A.
We remark that the stochastic integral equation in Definition 3.1 is equivalent to the identity Assumptions.
Theorem 3.2. We assume (A 1 ) is satisfied. Then, for every T > 0, p ≥ 1 and for some constant c ≥ 0 depending on T.
Proof. We formulate the proof for p > 1; the case p = 1 follows similarly. For a constant λ > 0 that will be chosen below and p > 1 we denote by Ξ p,λ the space of F t -adapted, continuous processes in H such that Then, Ξ p,λ is a Banach space equipped with the norm show that the mapping Φ : X → X , is well-defined, Lipschitz-continuous, and has a unique fixed point in X .
From now on, the letter C will always denote a non-negative constant, independent of T , that may vary in each occurrence and from line to line.
Step 1. Φ is a well-defined, bounded operator and satisfies the estimate for some constant c T > 0 depending on T .
Here and in the sequel, we use the notation Z t = (X t , Y t ). By assumption (A 1 ) we obtain where we used the estimate e −A(t−s) L(H) ≤ Ce −C(t−s) and L(H) denotes, as usual, the space of linear and continuous functions on H. Now taking the supremum and expectation value in the above inequality it follows that Similarly, we derive the estimate Theorem I.7.2, p.40, in [16] and hypothesis (A 1 ) further imply that Summing up we obtain for some constant C ≥ 0, which implies estimate (4). It remains to prove continuity of the image function t → Φ(Z) t . Let 0 < s < t < T. For the second equation we obtain and consequently, by a version of the Burkholder-Gundy-Davis inequality (see Theorem 7.1, p.39, in [16]) it follows that Kolmogorov's continuity criterion (Theorem 2.1 and the subsequent remark, p.6, in To prove the continuity of Φ 1 (X, Y ) we first recall the following estimates Hence, for the first component we obtain The above estimates imply and by (A 1 ), for the second term it follows that

PETER E. KLOEDEN, STEFANIE SONNER AND CHRISTINA SURULESCU
The last integral can be estimated by using (5), (A 1 ) and Hölder's inequality Taking the 2p-th power and expectation value and summing up we obtain for 0 < s < t ≤ T. As before, the continuity of t → Φ 1 (X, Y ) t now follows from Kolmogorov's continuity criterion.
Step 2. For λ sufficiently large Φ is a contraction in X . We assume Z = (X, Y ) and Z = ( X, Y ) are processes in X . Then, Φ(X, Y ) − Φ( X, Y ) satisfies the system Assumption (A 1 ) implies for the first equation , and taking the supremum and expectation value we obtain for some constant C T ≥ 0. Similarly, we derive a bound for the first integral in the second equation To estimate the remaining term we use Theorem I.7.2, p.40, in [11] and hypothesis (A 1 ), Adding the inequalities deduced above we obtain for some constant C ≥ 0. Consequently, if λ > 0 is sufficiently large Φ is a contraction in X = Ξ p,λ × Ξ p,λ .
Step 3. Existence and uniqueness. For λ large enough Φ is a contraction in X = Ξ p,λ × Ξ p,λ . By Banach's fixed point theorem it possesses a unique fixed point (X, Y ), which is the unique mild solution of the initial value problem (2).
Let (X, Y ) = Z = Φ(Z) be the mild solution of (2). By Hölder's inequality and (A 1 ) it follows that To estimate the stochastic integral we use again Theorem I.7.2, p.40, in [16] and (A 1 ), Adding the relevant inequalities we obtain for some constant C ≥ 0, which by Gronwall's lemma implies (3).

3.1.
Invariance. The solutions of mathematical models in biology often describe quantities that are non-negative and bounded by a certain maximum value. The following conditions ensure that solutions emanating from initial data in a given admissible range remain within this range for t > 0.
, and the stochastic perturbation fulfillŝ respectively. An admissible stochastic perturbation is, e.g., the function ensuring invariance of the set [0, m * 1 ] × [0, m * 2 ]. Theorem 3.3. In addition to the hypotheses of Theorem 3.2 we assume that (A 2 ) holds and the initial data (ζ, η) is deterministic and such that 0 ≤ ζ ≤ m * 1 , 0 ≤ η ≤ m * 2 in D. Then, the solutions are almost surely non-negative and uniformly bounded by m * 1 and m * 2 , respectively. Proof. Let (ζ, η) be given initial data satisfying the stated assumptions. We first assume thatf = (f 1 ,f 2 ) andĝ satisfy the conditions in (A 2 ) for all (x 1 , x 2 ) ∈ R 2 , t ∈ [0, T ]. More precisely, we denote byf = (f 1 ,f 2 ) andg modified functions that coincide on [0, m * 1 ] × [0, m * 2 ] withf andĝ and satisfy (A 2 ) in R 2 × [0, T ]. Moreover, we denote the solutions of the corresponding modified system by X and Y . Due to the continuity of solutions shown in Step 1 the first equation can be considered pathwise. By deterministic comparison principles for scalar parabolic equations it follows that X remains pathwise non-negative and bounded by m * 1 . On the other hand, for every fixed x ∈ D the SDE in (2) fulfills the hypothesis of the stochastic invariance criterion (see [20] or [3]), which implies that Y takes values within the interval [0, m * 2 ] with probability 1. Finally, the solutionsX andỸ satisfy the original system withf = (f 1 ,f 2 ) and g, and by the uniqueness of solutions we conclude that the set [0, m * 1 ] × [0, m * 2 ] is invariant for the original problem (2). Consequently, solutions corresponding to initial data within the given range are non-negative, uniformly bounded and exist globally.

3.2.
Boundedness properties of the solutions. We will need the following properties of solutions of (2) to show the well-posedness of stochastic mean-field models.
for some constant c 1 ≥ 0 depending on T .
Proof. Let Z = (X, Y ) be the mild solution of (2). Itô's formula implies that Y t 2 satisfies the SDE taking the expectation value in the above equality we obtain where we used hypothesis (A 1 ). On the other hand, multiplying the PDE in (2) by X and integrating over Ω leads to 1 2 We integrate the equation from 0 to t, disregard the negative term, and use (A 1 ) to conclude Taking expectation values we obtain Adding both estimates leads to and the proposition follows from Gronwall's Lemma.

Proposition 2.
Under the hypotheses of Theorem 3.2, for every p ∈ N the solutions of system (2) satisfy the estimate for 0 < s < t ≤ T and some constant c 2 ≥ 0 depending on T .
Proof. Let Z = (X, Y ) be the mild solution of (2) and t, s ∈ [0, T ] be such that t > s. For the first component we obtain Similar to the proof of Theorem 3.2, Step 1, we obtain and using (A 1 ) the second term can be estimated by Taking the 2p-th power and expectation value and using Hölder's inequality we obtain Furthermore, using (5), (A 1 ) and Hölder's inequality it follows that Taking the 2p-th power and expectation value it follows by Hölder's inequality that The solution of the SDE satisfies The first integral can be estimated similarly to J 2 , and the stochastic integral by applying Theorem I.7.1, p.39, in [16], Summing up all previous estimates and using Proposition 1 it follows that for some constant C T ≥ 0, hence obtaining the stated inequality.

4.1.
Well-posedness of the nonlocal model. To prove the well-posedness of (1) we use the results for local SDE-PDE systems in Section 3 and ideas applied in [13] to scalar SDEs with nonlocal sample dependence.
Assumptions. (H 2 ) For all r > 0 there exists a constant c r such that for all t ∈ [0, T ], x ∈ R 2 and y,ỹ ∈ R 2 such that |y| ≤ r.
for some constant c ≥ 0 depending on T . Moreover, if the initial data belongs to L 2p (Ω; H) with p > 1, then the solution (X, Y ) has continuous sample paths.
Proof. We construct approximate solutions on equidistant partitions of the interval [0, T ]. This will provide a Cauchy sequence in C([0, T ]; L 2 (Ω; H)) that converges to the solution of the original model (1). For n ∈ N such that 2 n > T we set h n := T 2 n and t n k := kh n for k = 0, . . . , 2 n . We denote by Z n = (X n , Y n ) the approximate solution obtained by successively solving the local systems dX n We use finite differences and an explicit Euler scheme to solve the PDE and the Euler-Mayurama method for the SDE. The parameters used in the simulations are summarized in Table 1, and the initial data is chosen as Observe that for this choice of interaction functions and parameter values the conditions ensuring invariance are satisfied, i.e., for any initial data such that 0 ≤ ζ ≤ 1, 0 ≤ η ≤ 1 the proton concentrations X and Y remain non-negative and bounded by 1. The shape of the initial conditions is motivated by the assumption that the tumor is located at the left end of the spatial interval and the cancer cells emit protons, rendering the tumor microenvironment acidic. The acidity is supposed to decrease towards the tumor edge (hence with advancing space). The concentration of intracellular protons should be lower than the one of their extracellular counterparts, in order to allow the survival and proliferation of tumor cells (thus the maintenance of proton dynamics). Figure 1 shows the time evolution of extra-and intracellular protons. We plotted five sample paths for solutions of the local model (9). Notice the inter-path differences in both proton populations, suggesting a corresponding tumor-to-tumor variability w.r.t. acidity, although the same type of cancer is considered. Figure 2 shows the behavior of solutions of the deterministic model and of the nonlocal system (8). The expectation values involved in the PDE for extracellular protons were computed by averaging over 20 ( Figure 2b) and 1000 sample paths (Figure 2c). Figure 2a shows the solutions of the deterministic model, i.e., the model with γ = 0, while each of Figures 2b and 2c shows the solution X of the mean-field SDE-PDE system and one sample path for the intracellular protons Y . As expected, when averaging over a large number of tumors the differences between the nonlocal model and the pure deterministic one are small. However, in the case of a much reduced number of tumors (corresponding to a less frequent cancer type) the randomness (inter-tumor variations) seems to play a significant role in acidification. This can be relevant for the sensitivity against therapies and for tumor aggressiveness, as it is well known by now that the low extracellular pH (pH e ) and the gradients between intracellular pH (pH i ) and pH e significantly influence the response of tumors to various treatments like radiotherapy and chemotherapy [5,23,4]. Our findings seem to endorse the necessity of paying particular attention to individualized treatment, especially for rare tumors, where it is difficult to rely on clinical experience acquired with a rather small number of patients.
The next step will be to include the dynamics of tumor cells and normal tissue and to assess the behavior of the two cell populations under the influence of the stochastic proton dynamics studied here. This is ongoing work.