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

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.


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.

International Journal of Analysis
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.

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 reactiondiffusion equations for the normal cell density (, x), the tumor cell density (, x) (both in cells/cm 3 ), and the concentration of excessive H + (, x) ion concentration (in Mol): where  ≤ 1/2 denotes the strength parameter for the competition between normal cells and neoplastic tissue.Thereby, Ω ⊂ R  ( = 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 constants   and   are given in 1/s and represent the growth rates, and the constants   and   are expressed in cells/cm 3 and provide the carrying capacities of the normal and tumor cells, respectively.The death rate   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 cm 2 /s.The production rate   of H + protons is expressed in Mol ⋅ cm 3 /(cells ⋅ s), the uptake rate   (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   in cm 2 /s.
We also assume that there is no exchange of cells and H + protons through the boundary of the considered domain; thus where n denotes the outer unit normal vector to Ω.
The initial conditions are given by Thereby  0 (x),  0 (x), and  0 (x) are strictly positive functions, which satisfy the no-flux condition: In order to render the system (1) dimensionless we use the following transformations: along with the notations International Journal of Analysis 3 We obtain the system 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.

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 Definition 1.A weak solution of (1) with boundary conditions (2) and initial data (3) is a triple (, , ) of functions in ×  × , such that for all  ∈  1 (Ω) a.e. in [0, ] the following three equations are satisfied: Theorem 2. There exists  > 0, such that the system (1) with initial data (3) and boundary conditions (2) satisfying has a unique solution (, ) ∈ ( × ) ∩ ( × ) and  ∈ .
We set (11) with   ≤ 1 to be defined below.
In order to prove Theorem 2 we construct a sequence and prove its convergence towards the weak solution of the system.Let ( 0 ,  0 ) ∈ ( × ) ∩ ( × ) and  0 ∈  be the weak solution to the homogeneous system while (  ,   ) ∈N 0 ∈ ( × ) ∩ ( × ) and (  ) ∈N 0 ∈  is the weak solution to with the corresponding initial and boundary conditions ( 2) and (3).The existence and uniqueness of the functions (  ,   ,   ) ∈N 0 in the above sequence are ensured by the following.
Lemma 3 (properties of the iteration sequence).Under assumptions (10) there exists  > 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  ∈ N 0 it holds that International Journal of Analysis (ii) the functions   ,   , and   are positive for all  ∈ N 0 .Moreover, the following inequalities hold: (iii) the functions   ,   , and   satisfy for adequate constants (Ω, ) and for all  ∈ N 0 the estimates Remark 4. From ( 19) it follows that for all  ∈ N 0 .
Proof of Lemma 3. We perform mathematical induction with respect to .
Induction Start.The proof of the claims in Lemma 3 for  = 0 is done separately for each of ( 13)-( 15).(a) With the substitution Equation ( 13) becomes the heat equation thus by the theory of linear parabolic differential equations (see, e.g., [22]) and with the assumption  0 ∈  1 (Ω) it follows that there exists a unique solution  0 of ( 13) such that This weak solution also satisfies Further it is known (see, e.g., [23]) that the solution of ( 13) can be written explicitly with respect to the initial condition  0 and the heat kernel and it is therefore positive.(b) Equation ( 14) is linear and has a positive solution: which depends on  0 (, x).It follows immediately that and thus the estimation (23) for  = 0 is obtained.The corresponding statement (19) for  0 is to be justified below.
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 (, x) =   (1 − ( 0 (, x)/  )) in ( 15) depends on time. Let Considering the symmetric bilinear form the dependence of the coefficient (, x) on  leads in its time derivative to a supplementary summand where for shortness we denoted by  the derivative with respect to .
The rest of the proof of Theorem 7.1.5 in [22] can now be adapted to obtain for an arbitrary  > 0 the estimate Now let (recall (33)) Upon integrating with respect to  one can majorize with (Ω, ) 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  0 ∈  1 (Ω) in order to show that there exists a unique weak solution  0 (, x) to ( 15) such that Since it follows from the weak maximum principle that  0 (, x) > 0 and thus also the positivity of  0 (, x).
The proof of the inequalities (21) for  = 0 does not differ from the one for a general  ∈ N given below and is therefore omitted here.
With (a)-(c) we proved all statements of Lemma 3 for  = 0. Induction Hypothesis.Assume the assertions of the lemma hold for an arbitrary  ∈ N 0 .
In order to establish the lower bound for  +1 define an auxiliary function  +1 (, x) :=  +1 (, x) −    −   , for which it holds For every nonnegative  ∈  1 (Ω) the right-hand side is positive.Further,  +1 (0, x) ≥ 0 by construction; thus it follows with the weak maximum principle that  +1 ≥ 0 a.e. which leads to  +1 (, x) ≥    −   .Now ( 17) is a linear, inhomogeneous differential equation, with solution where In order to prove (19) for  + 1 we have to show that 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 ‖ +1 ‖  ∞ ((0,]×Ω) < ∞, which leads to       +1       ∞ ((0,]×Ω) < ∞ (57) and thus (55) is proved.Now we prove the positivity of  +1 and the corresponding inequality in (21).To this aim use the induction hypothesis to observe that Next notice that there exists a positive constant C such that  +1 (, x) ≥ C for a.e.x ∈ Ω,  ∈ [0, ].This leads to the estimate This in turn immediately implies via (52) the positivity of  +1 .
In the next step we prove the estimate (23) for  +1 (, x).
by (23) and the induction hypothesis.
In order to prove the assertions of Lemma 3 for  +1 (, 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 In order to prove the positivity of  +1 we introduce an auxiliary function for  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  ∈ N 0 on an adequate time interval.

Proof (of the Statement (67))
Induction Start.The proof of (67) for  = 0 is identical to the one for  + 1.
Inductive Step.Upon using (66) in ( 18) we get Since for the right-hand side of (68) we have that holds for  <  3 with correspondingly chosen  3 and  such that  > 1.
In virtue of (67), for  ≤ 1/  the right-hand side in ( 18) is positive.Since by hypothesis  +1 0 > 0, the weak maximum principle implies the positivity of  +1 .This ends the proof of all statements in Lemma 3 for an arbitrary  ∈ N 0 and therefore the proof of the lemma itself.Now we are able to pass to the following.

Proof (of Theorem 2)
Existence.In order to prove the existence of a weak solution to (1) and ( 2) we show that the iterative sequence (  ,   ,   ) ∈N 0 is Cauchy.
Due to the completeness of  1 (Ω) and  2 (Ω), this will imply the convergence of the sequence to some limit functions , , and , these being solutions to (1) and (2).
Consider an arbitrary  ∈ N 0 .Since   0 ,  +1 0 ∈  1 (Ω), and   ,  +1 ∈  2 (0, ;  2 (Ω)), it follows that Next, one can apply Theorem 7.1.5 in [22] to the difference  +1 −   to deduce the estimate The right-hand side above can be further estimated and with the embedding constant  3 :=  3 (Ω, ) it follows that where International Journal of Analysis In order to obtain a corresponding estimate for the sequence (  ) ∈N , consider two consecutive terms in (17) written for   and  +1 and substract.This leads to ).Now multiply with ( +1 −   ) and integrate with respect to x to infer Thus Next we estimate the above terms.
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  parts with  + 1 nodes, whereas in two dimensions each of the axes is divided into  parts, thus obtaining a ( + 1) × ( + 1) mesh in two dimensions.Moreover the nodes are reordered in a row-wise manner, leading to a total of ( + 1) 2 nodes.Thereby the subindex  in the discretized equations denotes the spatial node   , where  = 1, 2, . . .,  + 1 and  = 1, 2, . . ., ( + 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: where  denotes the time level and Δ and Δ are the time and space increments, respectively.The discretized equations (96) for each  leads to the following system of equations of the form with H +1 and   H standing for the vectors containing the values of  and Δ   at the ( + 1)st and th time levels at the discretized space points, respectively; A H being the tridiagonal and block tridiagonal matrix for the one-and twodimensional cases, respectively.The updated values H +1 are used to find the values of the normal cell density at the time level  + 1 by solving for each space point .
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, where () =   (1 − ) and   ⊂  is the index set pointing at the direct neighbours of   on the grid.Thus   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 where A K is the ( + 1) × ( + 1) tridiagonal matrix or the block tridiagonal matrix of size (+1) 2 ×(+1) 2 for the oneand two-dimensional cases, respectively, coming from the nonstandard finite difference discretization,  ñ K is the vector coming from the proliferation term with the entries      (1 −    −  +1  ), and K +1 and K  are the vectors containing the  values at the discretized points for the time levels  + 1 and , respectively.
Throughout our simulations we use the biological parameter values from Table 1 which is reproduced from [24].
According to a linear stability analysis performed in [24], the parameter   plays a crucial role in characterizing the aggressivity of the tumor.There,   = 1 was shown to be the crossover value; for   < 1 the tumor is less aggressive, whereas for   > 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 Δ = 0.1 and Δ = 0.01, respectively.
In Figures 1 and 2 we present the simulations with   = 50 (an aggressive tumor) and   = 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).
In the next set of graphs we consider the negative effect of the aggressivity parameter   on the normal cell density.We know from the nondimensionalization addressed in Section 2 that the (nondimensionalized) parameter   is proportional to the death rate   and inversely proportional to the proton reabsorption rate   .Figure 3(a) shows the evolution of the normal cell density with respect to   for the times up to  = 40 at a fixed space point  = 0.4.One can notice that a more aggressive tumor (larger   or equivalently larger   ) 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   and larger   , 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.
Figure 3(b) illustrates the normal cell density at  = 10 depending on the concentration of H + protons for several different values of   .In an organism whose normal cells are more sensitive to pH variations (larger   ) 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.

The Two-Dimensional
Case.We perform 2D simulations in the unit square [0, 1] × [0, 1] using   = 70 and still with the parameter values in Table 1.Analogously to our computations in 1D we consider two different cases:   = 12.5 (an aggressive tumor) and   = 0.5 < 1 (a lessaggressive one).We use Δ = 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 Figure 4.At an earlier time (e.g.,  = 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  = 10, see Figures 4(b) and 4(e)), whereas by  = 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.
Also observe that the proton concentration varies proportionally to the tumor cell density and is inversely proportional to the normal cell density (Figures 5 and 6, resp.).Moreover, the healthy tissue is completely destroyed in the case of an aggressive tumor (Figure 6(c)) by  = 50.

Table 1 :
Parameter values used in the model.