Research Article | Open Access
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.
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 ) or tumor cell motility, for example, by restructuring the cytoskeleton or by producing matrix degrading enzymes , 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.,  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, . 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 , 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 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 . 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 , 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  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 , the tumor cell density (both in ), and the concentration of excessive ion concentration (in Mol): where denotes the strength parameter for the competition between normal cells and neoplastic tissue. Thereby, 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 protons by cancer cells.
The constants and are given in and represent the growth rates, and the constants and are expressed in and provide the carrying capacities of the normal and tumor cells, respectively. The death rate of the normal cells is measured in , and the diffusion coefficient of tumor cells in the absence of normal cells is given in . The production rate of protons is expressed in , the uptake rate (due for instance to proton buffering (see, e.g.,  and the references therein) and/or various ion exchangers between intracellular and extracellular domains (see e.g. )) is measured in , and their diffusion coefficient of in .
We also assume that there is no exchange of cells and protons through the boundary of the considered domain; thus where denotes the outer unit normal vector to .
The initial conditions are given by Thereby , and 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
We obtain the system where for simplicity the tilde notations have been ignored. The stability analysis of this system (with ) has been performed in , 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
We set with 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.
The existence and uniqueness of the functions in the above sequence are ensured by the following.
Lemma 3 (properties of the iteration sequence). Under assumptions (10) there exists 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 it holds that (ii)the functions , , and are positive for all . Moreover, the following inequalities hold: (iii)the functions , , and satisfy for adequate constants and for all the estimates
Remark 4. From (19) it follows that for all .
Proof of Lemma 3. We perform mathematical induction with respect to .
Induction Start. The proof of the claims in Lemma 3 for 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., ) and with the assumption it follows that there exists a unique solution of (13) such that This weak solution also satisfies Further it is known (see, e.g., ) that the solution of (13) can be written explicitly with respect to the initial condition and the heat kernel and it is therefore positive.
(b) Equation (14) is linear and has a positive solution: which depends on . It follows immediately that and thus the estimation (23) for is obtained.
The corresponding statement (19) for is to be justified below.
(c) In order to prove the claims of Lemma 3 for we show first that The former follows from For it is For we can consider (27). For its solution it holds (see, e.g., ) that Therefore, and finally (33) follows, thus also (19) for .
The following proof of (20) and (24) upon starting from (15) relies on Theorem in Evans . However, that result cannot be directly applied to the present case, since the diffusion coefficient in (15) depends on time.
Let with functions such that Considering the symmetric bilinear form the dependence of the coefficient 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 in  can now be adapted to obtain for an arbitrary 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 in , upon taking into account (32) and in order to show that there exists a unique weak solution to (15) such that Since it follows from the weak maximum principle that and thus also the positivity of .
The proof of the inequalities (21) for does not differ from the one for a general given below and is therefore omitted here.
With (a)–(c) we proved all statements of Lemma 3 for .
Induction Hypothesis. Assume the assertions of the lemma hold for an arbitrary .
Inductive Step. The proof for is to be done separately for each of (16)–(18). Since for a corresponding embedding constant and thus the existence of a unique weak solution to (16), (2), and (3) follows from the theory of linear parabolic differential equations. The solution satisfies with .
In order to establish the lower bound for define an auxiliary function , for which it holds For every nonnegative the right-hand side is positive. Further, by construction; thus it follows with the weak maximum principle that a.e. which leads to .
Now (17) is a linear, inhomogeneous differential equation, with solution where In order to prove (19) for we have to show that Obviously, the first assertion (54) holds, due to the induction hypothesis.
Next, estimate 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 , which leads to and thus (55) is proved.
Now we prove the positivity of and the corresponding inequality in (21). To this aim use the induction hypothesis to observe that Next notice that there exists a positive constant such that for a.e. , . This leads to the estimate This in turn immediately implies via (52) the positivity of .
In the next step we prove the estimate (23) for .
Due to (52) we get by (23) and the induction hypothesis.
In order to prove the assertions of Lemma 3 for one can apply Theorem in , with (54), (55), and the same justification as for the induction start at (c).
With an adequate embedding constant by (24) and the induction hypothesis; therefore and . By applying Theorem in  it follows that (18) has a unique weak solution with Now choose such that and Then and thus the estimate holds.
In order to prove the positivity of 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 on an adequate time interval.
Proof (of the Statement (67))
Induction Start. The proof of (67) for is identical to the one for .
Induction Hypothesis. Assume assertion (67) holds for an arbitrary .
Inductive Step. Upon using (66) in (18) we get Since for the right-hand side of (68) we have that holds for with correspondingly chosen and such that .
Since by construction , we can apply the weak maximum principle for to show that from which it also follows that This completes the proof of the statement (67).
In virtue of (67), for the right-hand side in (18) is positive. Since by hypothesis , the weak maximum principle implies the positivity of . This ends the proof of all statements in Lemma 3 for an arbitrary and therefore the proof of the lemma itself.
Now we are able to pass to the following.
Proof (of Theorem 2)
Consider an arbitrary . Since , and , it follows that Next, one can apply Theorem in  to the difference to deduce the estimate The right-hand side above can be further estimated and with the embedding constant it follows that where In order to obtain a corresponding estimate for the sequence , consider two consecutive terms in (17) written for and and substract. This leads to Denote .
Now multiply with and integrate with respect to to infer Thus Next we estimate the above terms.
Now for the second term on the right-hand side of (79) with . The two estimates above thus lead to Applying Gronwall’s inequality we deduce and finally with we get is chosen such that Now since and , we get Theorem in  can be applied to the difference , leading to The right hand side of this inequality can further be majorized and with the embedding constants and it follows that where Thus, putting all together, Therefore, is a Cauchy sequence in , from which the existence of a weak solution follows.
Regularity of the Solution. From (20) it follows that is uniformly bounded with respect to in ; therefore is compactly embedded in . This implies that for we have .
Theorem 5. (The local solution from Theorem 2 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 , then on , and so forth. Eventually the existence of a unique solution in shown on for any bounded .
4. Numerical Simulations
In this section we perform the numerical simulation of the system (7). The boundary conditions for and 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) where and in one- and two-dimensional cases, respectively. In our computations we have taken both and the strength parameter to be .
For the discretization of the model the finite difference method is employed. Thereby in the one-dimensional case the interval is divided into parts with nodes, whereas in two dimensions each of the axes is divided into parts, thus obtaining a mesh in two dimensions. Moreover the nodes are reordered in a row-wise manner, leading to a total of nodes. Thereby the subindex in the discretized equations denotes the spatial node , where and 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 and standing for the vectors containing the values of and at the st and th time levels at the discretized space points, respectively; being the tridiagonal and block tridiagonal matrix for the one- and two-dimensional cases, respectively. The updated values are used to find the values of the normal cell density at the time level 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 , that is, where 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 is the tridiagonal matrix or the block tridiagonal matrix of size for the one- and two-dimensional cases, respectively, coming from the nonstandard finite difference discretization, is the vector coming from the proliferation term with the entries , and and are the vectors containing the values at the discretized points for the time levels and , respectively.
According to a linear stability analysis performed in , the parameter plays a crucial role in characterizing the aggressivity of the tumor. There, was shown to be the crossover value; for the tumor is less aggressive, whereas for 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 and choose for the time and space increments and , respectively.
In Figures 1 and 2 we present the simulations with (an aggressive tumor) and (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).