# Mathematical Analysis and Numerical Simulations for a System Modeling Acid-Mediated Tumor Cell Invasion

Authors: Christian Märkl; Gülnihal Meral; Christina Surulescu
Journal: International Journal of Analysis
 (2013)

---

## Abstract

This work is concerned with the mathematical analysis of a model proposed
by Gatenby and Gawlinski (1996) in order to support the hypothesis that tumor-induced
alteration of microenvironmental pH may provide a simple but comprehensive mechanism
to explain cancer invasion. We give an intuitive proof for the existence of a solution under
general initial conditions upon using an iterative approach. Numerical simulations are also
performed, which endorse the predictions of the model when compared with experimentally
observed qualitative facts.


---

## Body

## 1. Introduction

Despite major progress in medicine and science there still are incurable diseases which can threaten human lives. Cancer is among the most severe ones and it manifests itself as an uncontrolled growth of cells which are produced by the organism subsequently to mutations. Cancer cells migrate through the surrounding tissue and degrade it on their way toward blood vessels and distal organs where they initiate and develop further tumors, a process known as metastasis.In the last decades various classes of models have been proposed aiming to provide a quantitative description of tumor growth. They range from the microscopic level of intracellular signaling pathways conditioning the growth of neoplastic tissue by stimulation or inhibition of apoptosis (e.g., by the influence of tumor necrosis factors [1]) or tumor cell motility, for example, by restructuring the cytoskeleton or by producing matrix degrading enzymes [2], through the level of cell-cell or cell-tissue interactions and up to the macroscopic level characterizing the behavior of the entire cell population. Multiscale settings like those in [3, 4] involve several of these scales and offer a systemic approach to the modeling process.When ignoring the setups relying on mechanical force balance and/or on the theory of mixtures, in the study of tumor invasion and metastasis one can distinguish between the so-called kinetic approach and the direct modeling at the macroscopic level. In the former a mesoscopic model is considered, consisting of an integro-partial differential equation for the evolution of the cell density, possibly coupled with integro-differential and/or reaction-diffusion equations for the fibre density of the extracellular matrix (ECM) and the chemotactic signal (see, e.g., [5] and the references therein [3, 4]). Then with an appropriate scaling the macroscopic limit is deduced, usually leading to a Keller-Segel-type model or some hyperbolic systems; see, for example, [6]. The macroscopic approach involves the largest class of existing models and directly accounts for processes at the level of cell populations, leading to systems of reaction-diffusion (transport) equations like, for example, in [7, 8].The role of tumor microenvironment in determining cancer malignancy has been put in evidence in several references; see, for example, [9, 10]. For instance hypoxia and acidity are factors that can trigger the progression from benign to malignant growth. In order to survive in the unfavourable environment they create, cancer cells upregulate certain proton extrusion mechanisms [11], the consequence of which is that the extracellular tumor environment has an acidic pH, which boosts apoptosis of normal cells and thus allows the neoplastic tissue to extend in the space becoming available. Hence the pH level directly influences the metastatic potential of tumor cells [12, 13]. These facts led Gatenby and Gawlinski [8, 14] to propose a model for the acid-mediated tumor invasion, which describes the interaction between the density of normal cells, tumor cells, and the concentration of H+ protons produced by the latter via reaction-diffusion equations. Starting from this model, travelling waves have been used to explain the aggressive action of cancer cells on their surroundings [15]. Further settings issued from Gatenby and Gawlinski’s model involve nutrient dynamics influenced by both vascular and avascular growth of multicellular tumor spheroids [16, 17] assuming rotational symmetry and investigating existence and qualitative properties of the solutions.In this work we reconsider the model in [8], whereby also explicitly allowing for crowding effects (due to competition with cancer cells) in the growth of normal cells. (This can be also done for the growth term in the tumor cell equation; however it does not change the analysis nor the qualitative behavior of the solutions to the system. Moreover, its biological motivation is not strong enough, since the competition between the two cell types does not really affect the growth, but rather the invasion of the neoplastic tissue: the acidity increase in the peritumoral environment is a byproduct of the enhanced glycolysis of cancer cells and not produced with the purpose of killing the normal cells and eluding concurrence.) For this setting we perform mathematical analysis and numerical simulations in order to verify the model predictions with respect to experimentally observed qualitative facts. In order to prove the existence of a unique (weak) solution we propose an intuitive method relying on an iterative procedure which has also been applied in [18, 19] in a different context (one of the supplementary difficulties here is the diffusion coefficient being nonconstant, but depending on the solution itself) and allows to avoid the use of operator semigroups.


## 2. Problem Setting

The model by Gatenby and Gawlinski [8] describes the evolution of normal and tumor cell density, respectively, in a domain where these cell types interact on the basis of pH value modifications. The mathematical description of these processes is ensured by the following system of reaction-diffusion equations for the normal cell density N(t,x), the tumor cell density K(t,x) (both in cells/cm3), and the concentration of excessive H+(t,x) ion concentration (in Mol):
(1)∂N∂t=wNN(1-NKN-θKKK)-dNHNin(0,T)×Ω,∂K∂t=wKK(1-KKK)+∇·(DK(1-NKN)∇K)in(0,T)×Ω,∂H∂t=wHK-dHH+DHΔHin(0,T)×Ω,
where θ≤1/2 denotes the strength parameter for the competition between normal cells and neoplastic tissue. Thereby, Ω⊂ℝn(n=1,2,3) is a regular-enough and bounded domain and only microscopically small processes are considered at the interface between tumor and healthy tissue. Observe that the diffusion coefficient of the cancer cells depends on the normal cell density: when the healthy tissue is at its carrying capacity the neoplastic tissue cannot diffuse; thus the tumor is confined. It can only spread if the surrounding normal tissue is diminished from its carrying capacity and this is assumed to happen due to lowering the pH level upon secretion of H+ protons by cancer cells.The constantswN and wK are given in 1/s and represent the growth rates, and the constants KN and KK are expressed in cells/cm3 and provide the carrying capacities of the normal and tumor cells, respectively. The death rate dN of the normal cells is measured in 1/(Mol·s), and the diffusion coefficient of tumor cells in the absence of normal cells is given in cm2/s. The production rate wH of H+ protons is expressed in Mol·cm3/(cells·s), the uptake rate dH (due for instance to proton buffering (see, e.g., [20] and the references therein) and/or various ion exchangers between intracellular and extracellular domains (see e.g. [21])) is measured in 1/s, and their diffusion coefficient of DH in cm2/s.We also assume that there is no exchange of cells andH+ protons through the boundary of the considered domain; thus
(2)∂K∂n=∂H∂n=0in(0,T)×∂Ω,
where n denotes the outer unit normal vector to ∂Ω.The initial conditions are given by(3)K(0,x)=K0(x),N(0,x)=N0(x),H(0,x)=H0(x)inΩ.
Thereby K0(x),N0(x),  and H0(x) are strictly positive functions, which satisfy the no-flux condition:
(4)∂K0∂n=∂H0∂n=0on∂Ω.
In order to render the system (1) dimensionless we use the following transformations:
(5)N~=NKN,K~=KKK,H~=H·dHwHKK,t~=wN·t,x~=wNDH·x
along with the notations
(6)δN=dNwHKKdHwN,ρK=wKwN,ΔK=DKDH,δH=dHwN.We obtain the system(7)∂N∂t=N(1-N-θK)-δNHN,∂K∂t=ρKK(1-K)+∇·(ΔK(1-N)∇K),∂H∂t=δHK-δHH+ΔH,
where for simplicity the tilde notations have been ignored. The stability analysis of this system (with θ=0) has been performed in [8], leading to biologically significant predictions.


## 3. Existence and Uniqueness of Solutions

In this section we provide a natural proof for the existence and uniqueness of a weak solution to the system (1) with initial data (3) and boundary conditions (2). We make use of an iterative procedure instead of the classical approach via semigroup theory; this is more intuitive and allows for a separate treatment of the three equations in each step.Consider the function spaces(8)X∶=L∞(0,T;H1(Ω)),Y∶={u∈L2(0,T;H2(Ω)):ut∈L2(0,T;L2(Ω))},Z∶=L∞(0,T;L2(Ω)).Definition 1.
A weak solution of (1) with boundary conditions (2) and initial data (3) is a triple (H,N,K) of functions in X×Y×Z, such that for all ϕ∈H1(Ω) a.e. in [0,T] the following three equations are satisfied:
(9)∫Ω‍wHKϕdx=∫Ω‍Htϕdx+∫Ω‍DH∇H∇ϕdx+∫Ω‍dHHϕdx,∫Ω‍wNN(1-NKN-θKKK)ϕdx=∫Ω‍Ntϕdx+∫Ω‍NdNHϕdx,∫Ω‍wKK(1-KKK)ϕdx=∫Ω‍Ktϕdx+∫Ω‍DK(1-NKN)∇K∇ϕdx.Theorem 2.
There existsT>0, such that the system (1) with initial data (3) and boundary conditions (2) satisfying
(10)H0∈H1(Ω)∩C(Ω),N0∈L∞(Ω)∩H1(Ω),K0∈H1(Ω),H0≥CH>0,0<N0≤KN2,0<K0≤KK
has a unique solution (H,K)∈(X×X)∩(Y×Y) and N∈Z.We set(11)T∶=∏i=16‍Ti
with Ti≤1 to be defined below.In order to prove Theorem2 we construct a sequence
(12)(Hm,Km)m∈ℕ0∈(X×X)∩(Y×Y),(Nm)m∈ℕ0∈Z
and prove its convergence towards the weak solution of the system.Let(H0,K0)∈(X×X)∩(Y×Y) and N0∈Z be the weak solution to the homogeneous system
(13)∂H0∂t-DHΔH0+dHH0=0,(14)∂N0∂t+dNH0N0=0,(15)∂K0∂t-∇·(DK(1-N0KN)∇K0)=0,
while (Hm,Km)m∈ℕ0∈(X×X)∩(Y×Y) and (Nm)m∈ℕ0∈Z is the weak solution to
(16)∂Hm+1∂t-DHΔHm+1+dHHm+1=wHKm,(17)∂Nm+1∂t+dNHm+1Nm+1=wNNm(1-NmKN-θKmKK),(18)∂Km+1∂t-∇·(DK(1-Nm+1KN)∇Km+1)=wKKm(1-KmKK)
with the corresponding initial and boundary conditions (2) and (3).The existence and uniqueness of the functions(Hm,Km,Nm)m∈ℕ0 in the above sequence are ensured by the following.Lemma 3 (properties of the iteration sequence).
Under assumptions (10) there exists T>0 such that(i)
there exists a unique weak solution to the systems (13)–(15) and (16)–(18) with conditions (3) and (2), and for every m∈ℕ0 it holds that
(19)Nm,Ntm∈L∞((0,T]×Ω),(20)Hm,Km∈L2(0,T;H2(Ω))∩L∞(0,T;H1(Ω)),Htm,Ktm∈L2(0,T;L2(Ω));(ii)
the functionsHm, Nm, and Km are positive for all m∈ℕ0. Moreover, the following inequalities hold:
(21)Hm(t,x)≥CHe-dHt,Nm(t,x)≤KN2,Km(t,x)≤KKfora.e.x∈Ω,t∈[0,T];(iii)
the functionsHm, Nm, and Km satisfy for adequate constants C(Ω,T) and for all m∈ℕ0 the estimates
(22)∥Hm∥X+∥Hm∥L2(0,T;H2(Ω))≤C(Ω,T)(∥K0∥H1(Ω)+∥H0∥H1(Ω)),(23)∥Nm∥X2≤C(Ω,T)∥N0∥H1(Ω)2,(24)∥Km∥X+∥Km∥L2(0,T;H2(Ω))≤2C(Ω,T)∥K0∥H1(Ω).Remark 4.
From (19) it follows that
(25)Nm∈L∞(0,T;L2(Ω))
for all m∈ℕ0.Proof of Lemma3.
We perform mathematical induction with respect tom.
Induction Start. The proof of the claims in Lemma 3 for m=0 is done separately for each of (13)–(15).
(a) With the substitution(26)H~0(t,x)=H0(t,x)edHt.
Equation (13) becomes the heat equation
(27)H~t0-DHΔH~0=0;
thus by the theory of linear parabolic differential equations (see, e.g., [22]) and with the assumption H0∈H1(Ω) it follows that there exists a unique solution H0 of (13) such that
(28)H0∈L2(0,T;H2(Ω))∩L∞(0,T;H1(Ω)),Ht0∈L2(0,T;L2(Ω)).
This weak solution also satisfies
(29)∥H0∥X+∥H0∥L2(0,T;H2(Ω))≤C(Ω,T)∥H0∥H1(Ω).
Further it is known (see, e.g., [23]) that the solution of (13) can be written explicitly with respect to the initial condition H0 and the heat kernel and it is therefore positive.
(b) Equation (14) is linear and has a positive solution:
(30)N0(t,x)=N0e-∫0t‍dNH0(s,x)ds>0
which depends on H0(t,x). It follows immediately that
(31)∥N0∥L∞(0,T;H1(Ω))2=∥N0e-∫0t‍dNH0dt∥L∞(0,T;H1(Ω))2≤∥N0∥H1(Ω)2,
and thus the estimation (23) for m=0 is obtained.
The corresponding statement (19) for N0 is to be justified below.
(c) In order to prove the claims of Lemma3 for K0 we show first that
(32)N0∈L∞((0,T]×Ω),(33)Nt0∈L∞((0,T]×Ω).
The former follows from
(34)∥N0∥L∞((0,T]×Ω)=(30)∥N0e-∫0t‍dNH0dt∥L∞((0,T]×Ω)≤∥N0∥L∞(Ω)<∞.
For t≥δ>0 it is
(35)∥Nt0∥L∞((0,T]×Ω)=(30)dN∥N0·e-∫0t‍dNH0dt·H0∥L∞((0,T]×Ω)<∞.
For t→0 we can consider (27). For its solution it holds (see, e.g., [22]) that
(36)lim(t,x)→(0,x0)H~0(t,x)=H0(x0)foreveryx0∈Ω.
Therefore,
(37)lim(t,x)→(0,x0)H0(t,x)=lim(t,x)→(0,x0)H~0(t,x)e-dHt=H0(x0)
and finally (33) follows, thus also (19) for m=0.
The following proof of (20) and (24) upon starting from (15) relies on Theorem 7.1.5 in Evans [22]. However, that result cannot be directly applied to the present case, since the diffusion coefficient a(t,x)=DK(1-(N0(t,x)/KN)) in (15) depends on time.
Let(38)km(t)∶=∑i=1m‍dmi(t)wi
with functions wi=wi(x) such that
(39){wi}i=1∞isanorthogonalbasisofH1(Ω),{wi}i=1∞isanorthonormalbasisofL2(Ω).
Considering the symmetric bilinear form
(40)A[km,km]∶=∫Ω‍a(t,x)(∇km)2dx,
the dependence of the coefficient a(t,x) on t leads in its time derivative
(41)ddtA[km,km]=∫Ω‍a′(t,x)(∇km)2dx+2∫Ω‍a(t,x)(∇km)′∇kmdx
to a supplementary summand
(42)∫Ω‍a′(t,x)(∇km)2dx=-∫Ω‍DKKN(N0)′(t,x)(∇km)2dx,
where for shortness we denoted by ′ the derivative with respect to t.
The rest of the proof of Theorem7.1.5 in [22] can now be adapted to obtain for an arbitrary ζ>0 the estimate
(43)∥km′∥L2(Ω)2+ddt(12A[km,km])≤Cζ(∥km∥H1(Ω)2+∥f∥L2(Ω)2)+2ζ∥km′∥L2(Ω)2+12∫Ω‍DKKN(N0)′(∇km)2dx.
Now let (recall (33))
(44)MN0∶=DKKN∥Nt0∥L∞((0,T]×Ω).
Upon integrating with respect to t one can majorize
(45)∫0T‍∫Ω‍DKKN(N0)′(∇km)2dxdt≤MN0∫0T‍∥∇km∥L2(Ω)2dt≤MN0∥km∥L2(0,T;H1(Ω))2≤γ(Ω,T)<∞,
with γ(Ω,T) an adequate constant. The rest of the proof can be done as in Theorem  7.1.5 in [22], upon taking into account (32) and K0∈H1(Ω) in order to show that there exists a unique weak solution K0(t,x) to (15) such that
(46)K0∈L2(0,T;H2(Ω))∩L∞(0,T;H1(Ω)),Kt0∈L2(0,T;L2(Ω)),∥K0∥X+∥K0∥L2(0,T;H2(Ω))≤C(Ω,T)∥K0∥H1(Ω).
Since
(47)K00(x)>0,
it follows from the weak maximum principle that K0(t,x)>0 and thus also the positivity of K0(t,x).
The proof of the inequalities (21) for m=0 does not differ from the one for a general m∈ℕ given below and is therefore omitted here.
With (a)–(c) we proved all statements of Lemma3 for m=0.
Induction Hypothesis. Assume the assertions of the lemma hold for an arbitrary m∈ℕ0.
Inductive Step. The proof for m+1 is to be done separately for each of (16)–(18). Since for a corresponding embedding constant c1∶=c1(Ω,T)(48)∫0T‍∥Km∥L2(Ω)2dt≤c1∫0T‍∥Km∥H1(Ω)2dt≤(24)ind.hyp.4c1C2(Ω,T)T∥K0∥H1(Ω)2<∞
and thus
(49)Km∈L2(0,T;L2(Ω)),
the existence of a unique weak solution to (16), (2), and (3) follows from the theory of linear parabolic differential equations. The solution Hm+1(t,x) satisfies
(50)Hm+1∈L2(0,T;H2(Ω))∩L∞(0,T;H1(Ω)),Htm+1∈L2(0,T;L2(Ω)),∥Hm+1∥X+∥Hm+1∥L2(0,T;H2(Ω))≤C1(Ω,T)(2wHC(Ω,T)c1T∥K0∥H1(Ω)+∥H0∥H1(Ω))≤𝒞(Ω,T)(∥K0∥H1(Ω)+∥H0∥H1(Ω)),
with 𝒞(Ω,T)=max{C1(Ω,T),C1(Ω,T)2wHC(Ω,T)c1T}.
In order to establish the lower bound forHm+1 define an auxiliary function ψm+1(t,x)∶=Hm+1(t,x)-CHe-dHt, for which it holds
(51)〈ψtm+1(t),ϕ〉+DH∫Ω‍∇ψm+1∇ϕdx+dH∫Ω‍ψm+1ϕdx=〈wHKm,ϕ〉.
For every nonnegative ϕ∈H1(Ω) the right-hand side is positive. Further, ψm+1(0,x)≥0 by construction; thus it follows with the weak maximum principle that ψm+1≥0 a.e. which leads to Hm+1(t,x)≥CHe-dHt.
Now (17) is a linear, inhomogeneous differential equation, with solution
(52)Nm+1(t,x)=e-α(t,x)(N0(x)+∫0t‍β(s,x)eα(s,x)ds),
where
(53)α(t,x)=∫0t‍dNHm+1(v,x)dv,β(s,x)=wNNm(s,x)(1-Nm(s,x)KN-θKm(s,x)KK).
In order to prove (19) for m+1 we have to show that
(54)Nm+1∈L∞((0,T]×Ω),(55)Ntm+1∈L∞((0,T]×Ω).
Obviously, the first assertion (54) holds, due to the induction hypothesis.
Next, estimate(56)∥Ntm+1∥L∞((0,T]×Ω)≤wN∥Nm∥L∞((0,T]×Ω)∥1-NmKN-θKmKK∥L∞((0,T]×Ω)+dN∥Nm+1∥L∞((0,T]×Ω)∥Hm+1∥L∞((0,T]×Ω)<∞,
due to (19).
Using again the induction hypothesis, the regularity of the initial data, and the properties of the solutions to the heat equations it follows immediately that∥Hm+1∥L∞((0,T]×Ω)<∞, which leads to
(57)∥Ntm+1∥L∞((0,T]×Ω)<∞
and thus (55) is proved.
Now we prove the positivity ofNm+1 and the corresponding inequality in (21). To this aim use the induction hypothesis to observe that
(58)Nm+1(t,x)≤KN2e-α(t,x)+∫0t‍β(s,x)e-(α(t,x)-α(s,x))ds≤KN2e-α(t,x)+wNKN2∫0t‍e-(α(t,x)-α(s,x))ds.
Next notice that there exists a positive constant C~H such that Hm+1(t,x)≥C~H for a.e. x∈Ω, t∈[0,T]. This leads to the estimate
(59)Nm+1(t,x)≤KN2e-dNC~Ht+wNKN21dNC~H(1-e-dNC~Ht)≤KN2((1-wNdNC~H)e-dNC~Ht+wNdNC~H)≤KN2.
This in turn immediately implies via (52) the positivity of Nm+1.
In the next step we prove the estimate (23) for Nm+1(t,x).
Due to (52) we get
(60)∥Nm+1(t)∥H1(Ω)2=∥e-α(t)N0+e-α(t)∫0t‍β(s)eα(s)ds∥H1(Ω)2≤2∥N0∥H1(Ω)2+2wN2∥∫0t‍Nm(s)(1-Nm(s)KN-θKm(s)KK)ds∥H1(Ω)2≤2∥N0∥H1(Ω)2+2wN2∥∫0t‍(Nm(s)-(Nm(s))2KN2)ds∥H1(Ω)2≤2∥N0∥H1(Ω)2+4wN2(∥∫0t‍Nm(s)ds∥H1(Ω)2+1KN2∥∫0t‍(Nm(s))2ds∥H1(Ω)2)≤∥N0∥H1(Ω)2[2+4wN2C(Ω,T)T2]≤𝒞(Ω,T)∥N0∥H1(Ω)2
by (23) and the induction hypothesis.
In order to prove the assertions of Lemma3 for Km+1(t,x) one can apply Theorem 7.1.5 in [22], with (54), (55), and the same justification as for the induction start at (c).
With an adequate embedding constantc2∶=c2(Ω,T)(61)∫0T‍∥Km(1-KmKK)∥L2(Ω)2dt≤∫0T‍(∥Km∥L2(Ω)+∥(Km)2KK∥L2(Ω))2dt≤2∫0T‍∥Km∥L2(Ω)2dt+2∫0T‍∥(Km)2KK∥L2(Ω)2dt≤2c12∫0T‍∥Km∥H1(Ω)2dt+2c24KK2∫0T‍∥Km∥H1(Ω)4dt≤8c12C2(Ω,T)∥K0∥H1(Ω)2T1T2+32c24KK2C4(Ω,T)∥K0∥H1(Ω)4T1T2<∞,
by (24) and the induction hypothesis; therefore Km(1-(Km/KK))∈L2(0,T;L2(Ω)) and K0(x)∈H1(Ω). By applying Theorem 7.1.5 in [22] it follows that (18) has a unique weak solution Km+1(t,x) with
(62)Km+1∈L2(0,T;H2(Ω))∩L∞(0,T;H1(Ω)),Ktm+1∈L2(0,T;L2(Ω)).
Now choose T1 such that max{T1C2(Ω,T),T1C4(Ω,T)}≤1 and
(63)T2∶=min{12,116wK2c12∥K0∥,KK264wK2c24∥K0∥3}.
Then
(64)∫0T‍∥wKKm(1-KmKK)∥L2(Ω)2dt≤∥K0∥H1(Ω)
and thus the estimate
(65)∥Km+1∥X+∥Km+1∥L2(0,T;H2(Ω))≤2C(Ω,T)∥K0∥H1(Ω)
holds.
In order to prove the positivity ofKm+1 we introduce an auxiliary function
(66)ξm+1(t,x)∶=-Atexp(1-ηt)+KK-Km+1(t,x),
for A positive and large enough and η a positive constant to be correspondingly chosen (see below). With the aid of this function we show that for all m∈ℕ0(67)Km+1≤KK,
on an adequate time interval.
Proof (of the Statement (67))
Induction Start. The proof of (67) for m=0 is identical to the one for m+1.
Induction Hypothesis. Assume assertion (67) holds for an arbitrary m∈ℕ0.
Inductive Step. Upon using (66) in (18) we get
(68)∂ξm+1∂t-∇·(DK(1-Nm+1KN)∇ξm+1)=A(ηt-1)exp(1-ηt)-wKKm(1-KmKK).
Since
(69)Km(1-KmKK)≤KK,
for the right-hand side of (68) we have that
(70)A(ηt-1)exp(1-ηt)-wKKm(1-KmKK)≥0
holds for t<T3 with correspondingly chosen T3 and η such that ηt>1.
Since by constructionξm+1(0,x)≥0, we can apply the weak maximum principle for t≤T3 to show that
(71)A(ηt-1)exp(1-ηt)+KK-Km+1(t,x)=ξm+1≥0,
from which it also follows that
(72)Km+1≤KK.
This completes the proof of the statement (67).
In virtue of (67), for T≤1/wK the right-hand side in (18) is positive. Since by hypothesis K0m+1>0, the weak maximum principle implies the positivity of Km+1. This ends the proof of all statements in Lemma 3 for an arbitrary m∈ℕ0 and therefore the proof of the lemma itself.Now we are able to pass to the following.Proof (of Theorem2)Existence. In order to prove the existence of a weak solution to (1) and (2) we show that the iterative sequence (Nm,Km,Hm)m∈ℕ0 is Cauchy.Due to the completeness ofH1(Ω) and L2(Ω), this will imply the convergence of the sequence to some limit functions N, K, and H, these being solutions to (1) and (2).Consider an arbitrarym∈ℕ0. Since H0m,H0m+1∈H1(Ω), and Km,Km+1∈L2(0,T;L2(Ω)), it follows that
(73)H0m+1-H0m∈H1(Ω),Km+1-Km∈L2(0,T;L2(Ω)).
Next, one can apply Theorem  7.1.5 in [22] to the difference Hm+1-Hm to deduce the estimate
(74)∥Hm+1-Hm∥X2≤C(Ω,T)∫0T‍∥wHKm-wHKm-1∥L2(Ω)2dt.
The right-hand side above can be further estimated and with the embedding constant c3∶=c3(Ω,T) it follows that
(75)∥Hm+1-Hm∥X2≤C(Ω,T)wH2c32∫0T‍∥Km-Km-1∥H1(Ω)2dt≤C(Ω,T)wH2c32T4∥Km-Km-1∥X2≤12∥Km-Km-1∥X2,
where
(76)T4=min{14,14C(Ω,T)wH2c32}.
In order to obtain a corresponding estimate for the sequence (Nm)m∈ℕ, consider two consecutive terms in (17) written for Nm and Nm+1 and substract. This leads to
(77)∂∂t(Nm+1-Nm)+dN(Hm+1Nm+1-HmNm)=wN(Km-1KKNm(1-NmKN-KmKK)-Nm-1(1-Nm-1KN-Km-1KK)).
Denote h(Nm,Nm-1)∶=wN(Nm(1-(Nm/KN)-(Km/KK))-Nm-1(1-(Nm-1/KN)-(Km-1/KK))).Now multiply with(Nm+1-Nm) and integrate with respect to x to infer
(78)12∫Ω‍∂∂t(Nm+1-Nm)2dx+dN∫Ω‍(Nm+1-Nm)2Hm+1dx=∫Ω‍(h(Nm,Nm-1)-dNNm(Hm+1-Hm))×(Nm+1-Nm)dx.
Thus
(79)ddt∥Nm+1-Nm∥L2(Ω)2≤2wN∫Ω‍|(Nm(1-NmKN-KmKK)-Nm-1(1-Nm-1KN-Km-1KK))×(Nm+1-Nm)KmKK|dx+2dN∫Ω‍|Nm(Hm+1-Hm)(Nm+1-Nm)|dx≤[2wN∥Nm(1-NmKN-KmKK)-Nm-1(1-Nm-1KN-Km-1KK)∥L2(Ω)+2dN∥Nm(Hm+1-Hm)∥L2(Ω)NmKN]×∥Nm+1-Nm∥L2(Ω).
Next we estimate the above terms.Let (recall (19))
(80)Mmax∶=max{∥Nm-1∥L∞((0,T],×Ω)MNm∶=∥Nm∥L∞((0,T]×Ω),NNm-1∶=∥Nm-1∥L∞((0,T]×Ω)}.
With the embedding constant c4∶=c4(Ω,T) we obtain for the first term on the right-hand side of (79)
(81)2wN∥Nm-Nm-1-(Nm)2KN+(Nm-1)2KN-NmKmKK+Nm-1Km-1KK(Nm-1)2KN∥L2(Ω)≤2wN∥Nm-Nm-1∥L2(Ω)+4wNMmaxKN∥Nm-Nm-1∥L2(Ω)+2wNKNKK∥Km-Km-1∥L2(Ω)+2wN∥Nm-Nm-1∥L2(Ω)≤CN~∥Nm-Nm-1∥L2(Ω)+CK~∥Km-Km-1∥L2(Ω)
with CN~∶=4wN(1+(Mmax/KN)) and CK~∶=2wNKN/KK.Now for the second term on the right-hand side of (79)
(82)2dN∥Nm(Hm+1-Hm)∥L2(Ω)≤dNc4KN2∥Hm+1-Hm∥H1(Ω)=CH~∥Hm+1-Hm∥H1(Ω),
with CH~∶=dNKNc4. The two estimates above thus lead to
(83)ddt∥Nm+1-Nm∥L2(Ω)2≤12(CN~∥Nm-Nm-1∥L2(Ω)+CK~∥Km-Km-1∥L2(Ω)+CH~∥Hm+1-Hm∥H1(Ω))2+12∥Nm+1-Nm∥L2(Ω)2≤CN~2∥Nm-Nm-1∥L2(Ω)2+CK~2∥Km-Km-1∥L2(Ω)2+CH~2∥Hm+1-Hm∥H1(Ω)2+12∥Nm+1-Nm∥L2(Ω)2.
Applying Gronwall’s inequality we deduce
(84)∥Nm+1-Nm∥L2(Ω)2≤et/2∫0t‍(CN~2∥Nm-Nm-1∥L2(Ω)2+CK~2∥Km-Km-1∥L2(Ω)2+CH~2∥Hm+1-Hm∥H1(Ω)2)ds,
and finally with D(Ω,T)=eT/2max{CN~2,CK~2,CH~2} we get
(85)∥Nm+1-Nm∥L∞(0,T;L2(Ω))2≤D(Ω,T)(∥Nm-Nm-1∥L∞(0,T;L2(Ω))2+∥Km-Km-1∥L∞(0,T;L2(Ω))2+∥Hm+1-Hm∥X2∥Nm-Nm-1∥L∞(0,T;L2(Ω))2)T5≤14(∥Nm-Nm-1∥L∞(0,T;L2(Ω))2+∥Hm+1-Hm∥X2+∥Km-Km-1∥L∞(0,T;L2(Ω))2)≤14(54∥Nm-Nm-1∥L∞(0,T;L2(Ω))2+54∥Km-Km-1∥X2).T5 is chosen such that
(86)D(Ω,T)T5≤14.
Now since K0m,K0m+1∈H1(Ω) and Km(1-(Km/KK)),Km+1(1-(Km+1/KK))∈L2(0,T;L2(Ω)), we get
(87)K0m+1-K0m∈H1(Ω),[wKKm+1(1-Km+1KK)-wKKm(1-KmKK)]∈L2(0,T;L2(Ω)).
Theorem  7.1.5 in [22] can be applied to the difference Km+1-Km, leading to
(88)∥Km+1-Km∥X2≤C(Ω,T)∫0T‍∥wKKm(1-KmKK)-wKKm-1(1-Km-1KK)∥L2(Ω)2dt.
The right hand side of this inequality can further be majorized and with the embedding constants c5∶=c5(Ω,T) and c6∶=c6(Ω,T) it follows that
(89)∫0T‍∥wKKm(1-KmKK)-wKKm-1(1-Km-1KK)∥L2(Ω)2dt≤2∫0T‍wK2KK2∥(Km)2-(Km-1)2∥L2(Ω)2+wK2∥Km-Km-1∥L2(Ω)2dt≤2∫0T‍c54wK2KK2∥Km-Km-1∥H1(Ω)2∥Km+Km-1∥H1(Ω)2+wK2c62∥Km-Km-1∥H1(Ω)2dt≤(∫0T‍(4c54wK2KK2[∥Km∥H1(Ω)2+∥Km-1∥H1(Ω)2]+2wK2c62wK2KK2)dt)∥Km-Km-1∥X2≤(32C2(Ω,T)c54wK2KK2T6∥K0∥H1(Ω)2+2wK2c62T6)×∥Km-Km-1∥X2≤14∥Km-Km-1∥X2,
where
(90)T6∶=min{18,18κ,18λ},k∶=32C2(Ω,T)c54wK2KK2∥K0∥H1(Ω)2,λ∶=2wK2c62.
Thus, putting all together,
(91)∥Km+1-Km∥X2+∥Nm+1-Nm∥L∞(0,T;L2(Ω))2+∥Hm+1-Hm∥X2≤14(3∥Km-Km-1∥X2+∥Nm-Nm-1∥L∞(0,T;L2(Ω))2).
Therefore, (Hm,Nm,Km) is a Cauchy sequence in X×L∞(0,T;L2(Ω))×X, from which the existence of a weak solution follows.Uniqueness. Let (K1,N1,H1) and (K2,N2,H2) be two solutions to (1)–(3). Due to the previous estimates
(92)∥K1-K2∥X2≤14∥K1-K2∥X2,(93)∥H1-H2∥X2≤14∥K1-K2∥X2(94)∥N1-N2∥L∞(0,T;L2(Ω))2≤14∥N1-N2∥L∞(0,T;L2(Ω))2+14∥H1-H2∥X2+14∥K1-K2∥X2;
thus K1=K2 (92) and with (93) it follows that H1=H2. Finally, (94) implies that N1=N2. This completes the proof of the uniqueness.Regularity of the Solution. From (20) it follows that (Km,Hm) is uniformly bounded with respect to m in Y×Y; therefore Y is compactly embedded in L2(0,T;H1(Ω)). This implies that for m→∞ we have (K,H)∈Y×Y.Theorem 5.
(The local solution from Theorem2 exists globally).The proof follows upon sequentially extending the time interval on which the solution exists: the previously deduced estimates allow for a bootstrap of the local existence proof in a subsequent step on the time interval[T,2T], then on [2T,T], and so forth. Eventually the existence of a unique solution in shown on [0,𝕋] for any bounded 𝕋.


## 4. Numerical Simulations

In this section we perform the numerical simulation of the system (7). The boundary conditions for K and H are the no-flux boundary conditions given by (2). We assume that initially the normal cells are at half of their carrying capacity, while the tumor cells can be close to theirs and thus prone to invade the surrounding tissue. Since the pH level is lowered by the cancer cells the concentration of protons is taken proportional to the density of the latter. Thus using these assumptions we choose for the initial conditions of the system (7)
(95)K(0,x)=exp(-|x|2ϵ),ϵ>0,N(0,x)=0.5*(1-exp(-|x|2ϵ)),ϵ>0,H(0,x)=ζK(0,x),ζ∈[0,1),
where x∈[0,1] and x∈[0,1]×[0,1] in one- and two-dimensional cases, respectively. In our computations we have taken both ζ and the strength parameter θ to be 0.5.For the discretization of the model the finite difference method is employed. Thereby in the one-dimensional case the interval[0,1] is divided into m parts with m+1 nodes, whereas in two dimensions each of the axes is divided into m parts, thus obtaining a (m+1)×(m+1) mesh in two dimensions. Moreover the nodes are reordered in a row-wise manner, leading to a total of (m+1)2 nodes. Thereby the subindex i in the discretized equations denotes the spatial node xi, where i=1,2,…,m+1 and i=1,2,…,(m+1)2 in one and two dimensions, respectively. We use forward differences for the time derivatives in the system. The central difference is used for the diffusion term in the equation characterizing the ion concentration:
(96)Hin+1-HinΔt=δHKin-δHHin+1+(Hi-1n+1-2Hin+1+Hi+1n+1Δx2),
where n denotes the time level and Δt and Δx are the time and space increments, respectively. The discretized equations (96) for each i leads to the following system of equations of the form
(97)AHHn+1=Hn+ϑHn
with Hn+1 and ϑHn standing for the vectors containing the values of H and ΔtδHK at the (n+1)st and nth time levels at the discretized space points, respectively; AH being the tridiagonal and block tridiagonal matrix for the one- and two-dimensional cases, respectively. The updated values Hn+1 are used to find the values of the normal cell density at the time level n+1 by solving
(98)Nin+1=11+ΔtδNHin+1[Nin+ΔtNin(1-θKin-Nin)]
for each space point i.In order to write the term characterizing the dispersion of the neoplastic tissue into the healthy tissue we make use of the nonstandard finite difference scheme [24], that is,
(99)∇(D(N)∇K)|xi=12(Δx)2∑k∈Ni‍(D(Nkn+1)+D(Nin+1))×(Kkn+1-Kin+1),
where D(N)=δk(1-N) and Ni⊂I is the index set pointing at the direct neighbours of xi on the grid. Thus Ni has two elements for the one-dimensional case and four elements for the two-dimensional case and the matrix-vector form of the scheme reads
(100)AKKn+1=Kn+ϑKn~,
where AK is the (m+1)×(m+1) tridiagonal matrix or the block tridiagonal matrix of size (m+1)2×(m+1)2 for the one- and two-dimensional cases, respectively, coming from the nonstandard finite difference discretization, ϑKn~ is the vector coming from the proliferation term with the entries ρkKin(1-Kin-Nin+1), and Kn+1 and Kn are the vectors containing the K values at the discretized points for the time levels n+1 and n, respectively.Throughout our simulations we use the biological parameter values from Table1 which is reproduced from [24].Table 1
Parameter values used in the model.
Parameters
Range
K
N
5
×
1
0
7
/
c
m
3
K
K
5
×
1
0
7
/
c
m
3
w
N
1
×
1
0
-
6
/
s
w
K
1
×
1
0
-
6
/
s
D
K
2
×
1
0
-
10
/
c
m
2
/
s
D
H
5
×
1
0
-
6
/
c
m
2
/
s
w
H
2.2
×
1
0
-
17
M
·
c
m
3
/
s
d
H
1.1
×
1
0
-
4
/
s
d
N
0
→
10
/
M
·
sAccording to a linear stability analysis performed in [24], the parameter δN plays a crucial role in characterizing the aggressivity of the tumor. There, δN=1 was shown to be the crossover value; for δN<1 the tumor is less aggressive, whereas for δN>1 it becomes highly aggressive. This will be an important factor in our computations, too.


### 4.1. The One-Dimensional Case

We consider the space interval[0,1] and choose for the time and space increments Δt=0.1 and Δx=0.01, respectively.In Figures1 and 2 we present the simulations with δN=50 (an aggressive tumor) and δN=0.5 (a less-aggressive one), respectively. The rest of parameters are taken from Table 1. As time progresses, the difference between these two cases is more visible. In the aggressive case at later times the cancer cells invade a larger region and destroy the healthy tissue much more than in the case of a less-aggressive tumor (Figure 2).Variations ofH+ protons (starred), normal tissue (dashed), and neoplastic tissue (solid) for an aggressive tumor.
(a)
(b)
(c)
(d)Variations ofH+ protons (starred), normal tissue (dashed), and neoplastic tissue (solid) for a less-aggressive tumor comparing to the one in Figure 1.
(a)
(b)
(c)
(d)In the next set of graphs we consider the negative effect of the aggressivity parameterδN on the normal cell density. We know from the nondimensionalization addressed in Section 2 that the (nondimensionalized) parameter δN is proportional to the death rate dN and inversely proportional to the proton reabsorption rate dH.Figure3(a) shows the evolution of the normal cell density with respect to δN for the times up to t=40 at a fixed space point x=0.4. One can notice that a more aggressive tumor (larger δN or equivalently larger dN) leads to a faster decay in the density of the normal cells, as expected. On the other hand, Figure 3(a) also supports the intuitive fact that in an organism which can poorly buffer the issuing excessive protons (smaller dH and larger δN, resp.) the pH value will decay faster as well, thus triggering the decay of normal cell density, which in such an acid environment can no longer be sustained at a physiologically convenient level.(a) Evolution of the normal cell density for several different values ofδN. (b) Normal cell density with respect to the H+ proton concentration for several different values of δN.
(a)
(b)Figure3(b) illustrates the normal cell density at t=10 depending on the concentration of H+ protons for several different values of δN. In an organism whose normal cells are more sensitive to pH variations (larger dN) the density of these cells will decay faster for the same concentration of H+ protons. It can also be seen that a smaller reabsorption rate of excessive protons leads to a faster decay of normal cell density.


### 4.2. The Two-Dimensional Case

We perform 2D simulations in the unit square[0,1]×[0,1] using δH=70 and still with the parameter values in Table 1. Analogously to our computations in 1D we consider two different cases: δN=12.5 (an aggressive tumor) and δN=0.5<1 (a less-aggressive one). We use Δt=0.01 as the time increment and for the spatial discretization we take 11 nodes on each axis, leading to 121 nodes in the computational domain.The evolution of cancer cells is plotted in Figure4. At an earlier time (e.g., t=1, see Figures 4(a) and 4(d)) the difference between the less- and the higher-aggressive tumors is not relevant. However, as time progresses the difference starts to be visible (e.g., at t=10, see Figures 4(b) and 4(e)), whereas by t=50 the more-aggressive tumor (Figure 4(c)) invaded almost the whole domain and the less-aggressive one was not able to penetrate that far.Variations of cancer cells for the cases of an aggressive ((a)–(c)) and a less-aggressive ((d)–(f)) tumor.
(a)
(b)
(c)
(d)
(e)
(f)Also observe that the proton concentration varies proportionally to the tumor cell density and is inversely proportional to the normal cell density (Figures5 and 6, resp.). Moreover, the healthy tissue is completely destroyed in the case of an aggressive tumor (Figure 6(c)) by t=50.Variations of proton concentration for the cases of an aggressive ((a)–(c)) and a less-aggressive ((d)–(f)) tumor.
(a)
(b)
(c)
(d)
(e)
(f)Variations of healthy tissue for the cases of an aggressive ((a)–(c)) and a less-aggressive ((d)–(f)) tumor.
(a)
(b)
(c)
(d)
(e)
(f)


## 4.1. The One-Dimensional Case

We consider the space interval[0,1] and choose for the time and space increments Δt=0.1 and Δx=0.01, respectively.In Figures1 and 2 we present the simulations with δN=50 (an aggressive tumor) and δN=0.5 (a less-aggressive one), respectively. The rest of parameters are taken from Table 1. As time progresses, the difference between these two cases is more visible. In the aggressive case at later times the cancer cells invade a larger region and destroy the healthy tissue much more than in the case of a less-aggressive tumor (Figure 2).Variations ofH+ protons (starred), normal tissue (dashed), and neoplastic tissue (solid) for an aggressive tumor.
(a)
(b)
(c)
(d)Variations ofH+ protons (starred), normal tissue (dashed), and neoplastic tissue (solid) for a less-aggressive tumor comparing to the one in Figure 1.
(a)
(b)
(c)
(d)In the next set of graphs we consider the negative effect of the aggressivity parameterδN on the normal cell density. We know from the nondimensionalization addressed in Section 2 that the (nondimensionalized) parameter δN is proportional to the death rate dN and inversely proportional to the proton reabsorption rate dH.Figure3(a) shows the evolution of the normal cell density with respect to δN for the times up to t=40 at a fixed space point x=0.4. One can notice that a more aggressive tumor (larger δN or equivalently larger dN) leads to a faster decay in the density of the normal cells, as expected. On the other hand, Figure 3(a) also supports the intuitive fact that in an organism which can poorly buffer the issuing excessive protons (smaller dH and larger δN, resp.) the pH value will decay faster as well, thus triggering the decay of normal cell density, which in such an acid environment can no longer be sustained at a physiologically convenient level.(a) Evolution of the normal cell density for several different values ofδN. (b) Normal cell density with respect to the H+ proton concentration for several different values of δN.
(a)
(b)Figure3(b) illustrates the normal cell density at t=10 depending on the concentration of H+ protons for several different values of δN. In an organism whose normal cells are more sensitive to pH variations (larger dN) the density of these cells will decay faster for the same concentration of H+ protons. It can also be seen that a smaller reabsorption rate of excessive protons leads to a faster decay of normal cell density.


## 4.2. The Two-Dimensional Case

We perform 2D simulations in the unit square[0,1]×[0,1] using δH=70 and still with the parameter values in Table 1. Analogously to our computations in 1D we consider two different cases: δN=12.5 (an aggressive tumor) and δN=0.5<1 (a less-aggressive one). We use Δt=0.01 as the time increment and for the spatial discretization we take 11 nodes on each axis, leading to 121 nodes in the computational domain.The evolution of cancer cells is plotted in Figure4. At an earlier time (e.g., t=1, see Figures 4(a) and 4(d)) the difference between the less- and the higher-aggressive tumors is not relevant. However, as time progresses the difference starts to be visible (e.g., at t=10, see Figures 4(b) and 4(e)), whereas by t=50 the more-aggressive tumor (Figure 4(c)) invaded almost the whole domain and the less-aggressive one was not able to penetrate that far.Variations of cancer cells for the cases of an aggressive ((a)–(c)) and a less-aggressive ((d)–(f)) tumor.
(a)
(b)
(c)
(d)
(e)
(f)Also observe that the proton concentration varies proportionally to the tumor cell density and is inversely proportional to the normal cell density (Figures5 and 6, resp.). Moreover, the healthy tissue is completely destroyed in the case of an aggressive tumor (Figure 6(c)) by t=50.Variations of proton concentration for the cases of an aggressive ((a)–(c)) and a less-aggressive ((d)–(f)) tumor.
(a)
(b)
(c)
(d)
(e)
(f)Variations of healthy tissue for the cases of an aggressive ((a)–(c)) and a less-aggressive ((d)–(f)) tumor.
(a)
(b)
(c)
(d)
(e)
(f)



---

*Source: 878051-2013-03-10.xml*