Ahmed M Fouad
Department of Physics, Temple University, Philadelphia, PA 19122, USA
*Corresponding author: Ahmed M. Fouad, Department of Physics, Temple University, Philadelphia, PA 19122, USA.
Received: June 25, 2019 Published: September 30, 2019
The acid-mediated tumor invasion hypothesis proposes that altered glucose metabolism exhibited by the vast majority of tumors leads to increased acid (H+ ion) production which subsequently facilitates tumor invasion [1-3]. The reaction-diffusion model  that captures the key elements of the hypothesis shows how the densities of normal cells, tumor cells, and excess H+ ions change with time due to both chemical reactions between these three populations and density-dependent diffusion by which they spread out in three-dimensional space. Moreover, it proposes that each cell has an optimal pH for survival; that is, if the local pH deviates from the optimal value in either an acidic or alkaline direction, the cells begin to die, and that the death rate saturates at some maximum value when the microenvironment is extremely acidic or alkaline. We have previously studied in detail how the death-rate functions of the normal and tumor populations depend upon the H+ ion density . Here, we extend previous work by investigating how the equilibrium densities (at which the time rates of change of the cellular densities are equal to zero) reached by the normal and tumor populations in three-dimensional space are affected by the presence of the H+ ions, and we present detailed analytical and computational techniques to analyze the dynamical stability of these equilibrium densities. For a sample set of biological input parameters and within the acid-mediation hypothesis, our model predicts the transformation to a malignant behavior, as indicated by the presence of unstable sets of equilibrium densities.
KEYWORDS: Dynamical Stability Analysis; Fixed Points: Cancer Biology; Acid-Mediated Tumor invasion; Reaction-Diffusion Systems; Partial Differential Equations
The revolutionary work of Warburg almost a century ago demonstrated that tumor cells utilize anaerobic metabolism to produce energy by converting glucose to ATP even in the presence of abundant oxygen . The major disadvantage of anaerobic metabolism (when compared to the regular oxidation of glucose to CO2 and H2O) is the relatively low efficiency (low ATP output); however, tumor cells compensate for that by increasing the glucose flux. Moreover, FDG-PET relies basically on that glucose flux for tumor imaging purposes [6-10]. In fact, it has been well proven that the glucose consumption of tumors (primary and metastatic) is noticeably larger than normal-tissue cells [6-10]. In addition, through PET imaging, a direct relationship between tumor aggressiveness (and prognosis) and the rate of glucose uptake has been shown [6,7]. The transformation from normal to anaerobic metabolism has been directly linked to the increased acid production exhibited by tumor cells . Initially, it was assumed that the lactic acid produced by anaerobic metabolism would result in an acidic intracellular pH; however, studies with NMR spectroscopy have unanimously shown that the intracellular pH of tumors is generally slightly alkaline relative to normal cells . Other measurements have demonstrated that tumor cells actually excrete protons through amplification of the Na+/H+ antiporter and other membrane transporters. As a result, it is now clear that the extracellular pH of tumors is substantially lower (usually by about 0.5 pH unit) than normal [12-15]. Generally, a persistent pH below about 7.1 results in death of normal cells due to a p53-dependent apoptosis pathway triggered by increased caspase activity [16,17]. However, tumor cells are much more tolerant to acidic pH, presumably due to mutations in p53 or other components of the apoptotic pathway. In fact, tumor cells typically exhibit a maximum proliferation rate in relatively acidic mediums (i.e. pH 6.8) [18-20]. As a consequence, it has been concluded that the invasiveness of tumors are directly attributed to the acidification of their microenvironments, which, in turn, makes it ideal for tumor cells to proliferate and, moreover, corroborates the death of normal cells that compete against tumor cells for glucose and substrate; because of their relatively high intolerance of these acidic microenvironments.
An acid-mediated tumor invasion hypothesis (AMTIH) has been developed [1-3], wherein microenvironmental changes at the tumor-host interface caused by altered metabolism in transformed cells (accelerated uptake of glucose due to increased reliance on glycolysis) lead to an excess acid production, thus facilitating tumor invasion. Concurrent with this increase inglucose consumption is the up-regulation of the Na+/H+ antiport system [21,22]. The resulting net increase in tumor intracellular pH leads to normalization, and therefore to a decrease in extracellular pH by about 0.5 pH unit compared to normal tissue [11-13,15,20]. While some investigators have attributed this increased reliance on glycolysis to the existence of a compromised tumor blood supply leading to poor O2 perfusion of tumor tissue , others have found evidence that it is a phenotypical consequence of transformation .
The reaction-diffusion model of the AMTIH  predicts the degree to which excess hydrogen ions in the tumor interstitium diffuse into the surrounding tissue. It also predicts the extent of the zone of altered interstitial pH found at the tumor-host interface. This microenvironment is harmful to normal tissue because of three reasons. First, normal cell viability declines sharply below extracellular pH of 7.1 [11-13,15,20]. Second, an acidic microenvironment stimulates production of enzymes which degrade the extracellular matrix [25,26]. Third, diminished pH results in loss of intercellular gap junctions reducing the cohesion, collaboration, and communication among normal cells at the tumor edge . This altered microenvironment leads to a progressive loss of normal cells and degradation of extracellular matrix. This subsequently allows tumor cells (which are much more resistant to acidic extracellular pH) to invade the host tissue. Eventually, the acid-mediation hypothesis of tumor invasion is supported by recent reports that high lactate levels are directly correlated with the likelihood of metastases, tumor recurrence, and restricted patient survival .
The main research focus in this paper is to investigate how the equilibrium densities reached by the normal and tumor populations are affected by the existence of the H+ ions. To achieve that, given some sample values of biological input parameters (which are patient-specific), we present a detailed approach (that utilizes both analytical and computational techniques) to solving numerically for the equilibrium points (fixed points) that represent the equilibrium densities at which the normal, tumor, and acid populations could coexist in three-dimensional space. Moreover, through solving numerically for the eigenvalues of the Jacobian matrix of the partial differential equations that represent how the three populations evolve with time, we determine the dynamical stability of these fixed points.
This article is outlined as follows. In Sec. 2 we review the previous work on modeling tumor growth and invasion. The canonical Gatenby-Gawlinski reaction-diffusion model that captures the key components of the AMTIH is presented in Sec. 3. In Sec. 4 we present our analytical and computational results on the dynamical stability analysis of the reaction-diffusion model. The conclusions and outlook are presented in Sec. 5.
Tumor growth and other biological systems are particularly amenable to modeling with cellular automata. This approach has been used increasingly since Wolfram (1986) demonstrated that every automaton with local interactions belongs to one of a small number of universality classes . Some of the earliest work in tumor modeling with cellular automata was carried out by Düchting and Vogelsaenger to investigate the effects of different radio-therapeutic strategies on tumor growth [30,31].
More recently, Qi et al.  developed a cellular automaton model of tumor growth which includes immune system surveillance and accounts for mechanical pressure from within the tumor. This model is of interest because it reproduces the Gompertz growth law. The cellular automaton approach was also employed by Smolle and Stettner , and by Smolle [34,35] in an attempt to infer relationships between tumor morphology and basic biological features; such as cell cohesion, motility, and autocrine and paracrine growth factors. Kansal et al. [36,37] have made use of the Voronoi tessellation of space in a cellular automaton model of brain tumor growth dynamics. This elegant approach allows for the preservation of the discrete nature of cells (or groups of cells); but removes the anisotropy of a regular lattice.
Cellular automata have also been used to model neoangiogenesis, which is another crucial aspect of cancer biology. Chaplain and Anderson [38,39] have used cellular automata to solve the discrete version of a system of non-linear partial differential equations which describe the response in space and time of endothelial cells to tumor angiogenesis factor and fibronectin, including migration, proliferation, anastomosis, and branching. Cellular automata have also been used to model the development of leaf veins, insect trachea, and tumor neovascularization .
Some models that use systems of partial differential equations have been developed to simulate tumor growth and invasion; such as tumor angiogenesis [38,39,42,46], post-surgical response to tumor removal , avascular tumor growth [43,44], tumor cell motility during invasion , and drug resistance and vascular structure in chemotherapy . All these models (either automaton-based or continuum, stochastic or deterministic) simulate one or more aspects of cancer biology. Other models, including chemotherapeutic agent transfer  and diffusive vasculature for nutrient transfer , do so with a spatially-uniform term accounting for delivery or removal, not by including individual vessels distributed throughout space.
Reaction-diffusion model of tumor invasion
Gatenby and Gawlinski propose that the alteration of microenvironmental pH induced by the tumor’s anaerobic metabolism may provide a simple, yet complete mechanism for cancer invasion. In this section, we extend a reaction-diffusion model they developed, which describes the spatial distribution and temporal development of tumor tissue, normal tissue, and excess H+ ion density . A pH gradient across the tumor-host interface is proposed by the model, which is verified by experimental data. A transition from benign to malignant growth, analogous to the adenoma-carcinoma sequence, is predicted by the model through investigating the dynamics of the tumor-host interactions. This transition is controlled through the effects of many biological parameters. Moreover, it is crucial to pinpoint that the applicability of their model is limited to the vicinity of the tumor-host interface.
All realistic models should include as fundamental factors properties common to almost all tumors, regardless of their originating genetic instability and heterogeneity. One fundamental cellular dynamic is that tumor cells evolve away from the original differentiated state toward one that is more primitive (less ordered). Metabolic changes adapt to this evolution by increasing reliance on glycolytic metabolism (even in the presence of abundant oxygen) despite the 19-fold reduction in energy production that comes as a consequence . The Gatenby-Gawlinski model proposes that it is specifically this energy-production inefficiency associated with that altered metabolism that corroborates successful tumor invasion, particularly the transformation-induced reversion of neoplastic tissue to primitive glycolytic metabolic pathways with a consequent increased acid production and the diffusion of the resulting acidic protons into neighboring normal tissue. That very acidic microenvironment is what eventually corroborates tumor invasion which is relatively more tolerant to it than normal cells. The following temporal sequence can be concluded after all. Initially, high H+ ion concentrations in tumors diffuse into neighboring normal tissue, exposing these normal cells to tumor-like interstitial pH. Next, normal cells immediately adjacent to the tumor edge are incapable of surviving in this chronically acidic environment. Eventually, the progressive loss of layers of normal tissue at the tumor-host interface facilitates tumor invasion.
From the temporal sequence mentioned above, we can conclude that the invasion mechanism capitalizes primarily on the low interstitial pH of tumors (resulting from the anaerobic metabolism) and the subsequent, relatively low viability of normal tissue in an acidic environment favorable to tumor tissue. The Gatenby-Gawlinski model mathematically frames the acid-mediation hypothesis as a system of reaction-diffusion equations whose solutions could precisely predict the structure and dynamics of the tumor-host interactions.
The acid-mediated tumor invasion hypothesis described above can be represented by the following system of reaction-diffusion equations. Let N1 and N2 denote the cell densities in cell/cm3 of the normal and tumor populations, respectively. Assuming these populations only compete for available space, their temporal evolution is governed by the equations
where r1,2 and K1,2 are the growth rates in s-1 ( s denotes second) and spatial carrying capacities in cell/cm3 of the respective populations. Moreover, if we propose a diffusion model similar to the Fickian one, where the diffusion parameters are themselves density-dependent [maximum value in vacant space and zero in high concentrations (close-packed cells)], then Eq. (1) becomes
Next, the model assumes that each cell has an optimal pH for survival, and that if the local pH deviates from the optimal value (either in the acidic or alkaline direction), the cells start to die. In addition, it proposes that the death rate saturates at some maximum value when the environment is extremely acidic or alkaline. A simple ad-hoc functional form meeting these criteria is the “inverted Gaussian”
where H is the local concentration of H+ ions in mol/L, d1,2 are the saturated death rates in s-1, H opt are the local H+ ion concentrations in mol/L, corresponding to the optimal pH’s, and H width are the half-widths of the “inverted Gaussian” in mol/L. Substituting the death rates (3) into (2) we obtain
The H+ ions are produced at a rate proportional to the local concentration of tumor andremoved by the combined effects of buffering and vascular evacuation, both of which are proportional to micro-vessel areal density. Thus
where H is the H+ ion concentration in mol/cm3, r3 is the H+ ion production rate in mol/(cell *s), d3 is the H+ ion uptake rate in s-1, H0 is the H+ ion concentration in serum, and D3 is the H+ ion diffusion constant in cm2/s.
Eq. (4) and (5) can be rewritten in terms of the variables
In dimensionless form Eq. (4) and (5) become
where ρ2 = r2 / r1, ρ3 = r3 K / (H0 r1),δ3 = d3 / r3 , ∆3=D3 / DN. The death-rate functions are dimensionless and have the form
For in vitro spheroids, doubling times are between 1 and 4 days. Therefore, this yields r2 = ln 2/2.5 days ≈ 2.0 x 10-6/s. For normal tissue wound healing, 4 days seem reasonable for the doubling time; thus, r1 = ln 2/4.0 days ≈ 2.0 × 106/s . The volume-limited carrying capacities of tumor and normal tissues are assumed to be the same and given by K1 = K ≈ 5 × 108 cells/cm3. For vascular evacuation without buffering d3 = αp, where α ≈ 200/cm is the vessel areal density and p ≈ 1.2 × 104 cm/s is the vessel permeability for lactate, resulting in a removal rate of 2.4 ×10-2/s . Local buffering might increase this by 25%. Thus, the final estimate for this rate is d3 ≈ 3.0 ×10-2/s. If the serum pH0 = 7.4 is also the optimal pH for normal tissue growth, then H1opt = H0 3.98 × 10-11 mol/cm3. An optimal pH of 6.8 for tumor growth gives H2opt=1.58 × 10-10 mol/cm3.
As for acid production rate, if a tumor is sufficiently large, then the temporal and spatial derivatives at its core are small. From Equation (5) we see that r3 ≈ d3 (Hcore− H0)/ K2. Assuming a core pH of 6.4, this yields r ≈ 2.2 × 10-20 mol/(cell*s). This last value is remarkably consistent with the curve fit results of the Martin-Jain data to the old model . The lactic acid and cellular diffusion constants are D3 ≈ 5 × 10-6 cm2/s and DN ≈ 2 × 10-10 cm2/s, respectively.
Using the values above, the numerical estimates of the dimensionless parameters in Eq. (7) and (8) are as follows:
and δ2 = d2 / r1 = 2.0. The last four values are admittedly guesswork; but it is reasonable to choose the saturated death rates greater than the growth rates and the width of the optimal pH region for tumor cells greater than that for normal cells.
Equilibrium densities and dynamical stability analysis
Here, we present a detailed technique we developed to solve for the fixed points of the reaction-diffusion system defined in Eq. (7), with the sample numerical values of the biological input parameters given in Sec. 3, and determine their dynamical stability. The significance of the fixed points is that they determine equilibrium densities of the normal, tumor, and acid populations at which they could coexist in three-dimensional space (if no malignant behavior is already in progress). The dynamical stability of the fixed points can be determined by solving numerically for the eigenvalues of the Jacobian matrix. If all eigenvalues are negative, the fixed point is considered to be stable; however, if one or more of the eigenvalues is (are) positive, the fixed point is considered to be unstable. A stable fixed point represents an intermediate balance between the normal, tumor, and acid populations where none of them is likely to undergo any noticeable growth into the other two; however, an unstable fixed point where tumor cells are very likely to invade normal tissue upon any slight perturbations of the cellular densities could potentially develop into a malignancy. Note that in Equation (7), we present multiple reaction and diffusion factors that could contribute to a change in the cellular densities of the three populations. By mathematical definition, if the fixed point is stable, after slight changes in the equilibrium cellular densities (due to e.g. cellular diffusion) have taken place, the cellular densities tend to converge infinitesimally back to their initial equilibrium values (so the time rates of change of the cellular densities remain at zero after all); however, on the other hand, in an unstable fixed point, after the occurrence of slight perturbations in the equilibrium cellular densities, significant divergences off the equilibrium values are very likely to take place (that is, the cellular populations are not necessarily expected to converge back to their initial values), where malignancies are very likely to develop upon any small perturbations of unstable fixed points; because the presence of the acid specifically would always significantly catalyze the process in the direction where tumor invasion into normal tissue takes place (not vice versa).
Note that, when the microenvironment is extremely acidic, the death rate of tumor cells saturates at a value considerably lower than its counterpart for normal cells, and this is interpreted as tumor cells being more capable of surviving in an extremely acidic microenvironment [2,4].
We find the fixed points and analyze their dynamical stability as follows. If we restrict our attention to solutions η1, η2, and η3 ≡ Λ, which depend on time but are independent of the spatial coordinate ξ, Eq. (7) reduce to three coupled ordinary differential equations
The fixed points or points of equilibrium are points (η1*, η2*, η3*) where the right sides of Equation (9) vanish. Near a fixed point, the differential equations may be linearized in the form
Substituting ηj −ηj* = aj eλt, one finds n exponential solutions, with λ and aj given by the eigenvalues and eigenvectors of the matrix If all n of the eigenvalues λ are negative, η −η* decays exponentially with time, and the fixed point is “stable”. If one or more of the λ values is (are) positive on the other hand, the fixed point is unstable. Note that variables η1 (t), η2 (t), and Λ(t) depend on t but not on the coordinate ξ.
Next, we solve numerically for the fixed points. To find the fixed points, we set the derivatives in the cellular equations [the first and second of Eq. (7)] equal to zero after eliminating Λ using the acid equation [the third in Eq. (7)]. This yields two non-linear equations in two unknowns, namely η1 and η2, which are physically contained within the interval [0,1]. Therefore, we partition the domain 0 ≤ η1≤ 1 ⊗ 0 ≤ η2 ≤ 1 into a rectangular grid with ∆η = ∆η = 0.01 (10000-points grid), and use these 10000 points as numerical inputs for a multi-dimensional Newton’s algorithm. For a specific Newton’s solution for η1 and η2, η3 is found using the third equation in (7). We determine the stability of the resulting solutions by numerically computing the eigenvalues of the Jacobian matrix (we have published the software package we used to perform the dynamical stability analysis presented here ). The fixed- points analysis for the parameters given in the previous section is shown below [4,51] in Table 1. Fixed point 3 is stable (all eigenvalues are negative). Fixed point 4 is unstable (some eigenvalues are positive) and is very likely to develop into a malignancy upon any slight perturbations of the equilibrium cellular densities [note that when η1 is zero, this means that all the normal population has died at the interface (microenvironment) joining the three populations]. Fixed points 1, 2, and 5 are all unstable and unphysical (negative density values). Note that the existence of two different physical solutions (either stable or unstable) means that for the biological input parameters given in Sec. 3, wherever the normal, tumor, and acid populations are found in a dynamical equilibrium in three-dimensional space, it must be in one of these two solutions [note that we are not reporting here trivial fixed points such as (0,0,1) and (1,0,1)].
|(− 0.26500634, 0.01500634, 1.14005922)||(− 1.500000e+04, 1.205189e−01, 1.205189e−01)||Unstable (unphysical)|
|(− 0.23499366, − 0.01500634, 0.85994078)||(− 1.500000e+04, 1.295203e−01, 1.295203e−01)||Unstable (unphysical)|
|(0.00000000, 0.37206032, 4.47256298)||(− 1.498978e+04, − 1.081070e+01, −1.372060e+00)||Stable|
The canonical Gatenby-Gawlinski model of the acid-mediated tumor invasion hypothesis presented in Sec. 3 shows how the densities of normal cells, tumor cells, and excess H+ ions change due to two factors: first, local chemical reactions between the three populations and, second, density-dependent diffusion by which these populations spread out in space. According to the hypothesis, the dominant mechanism for cancer invasion is the tumor-induced alteration of microenvironmental pH. Moreover, the hypothesis predicts an acidic pH gradient extending into the peritumoral normal tissue, where the normal mammalian cells (but not tumor cells) are intolerant of acidic interstitial pH in the range typically found within this gradient.
DISSCUSSION, CONCLUSIONS, AND FUTURE WORK
In Sec. 4 we have developed detailed analytical and computational techniques to analyze the dynamical stability of the Gatenby-Gawlinski reaction-diffusion system by solving for the fixed points of the differential equations [see Eq. (7)] that represent how the normal, tumor, and acid populations evolve with time, then determining the stability of these fixed points by solving for the eigenvalues of the Jacobian matrix. The significance of the fixed points is that they determine possible densities at which the normal, tumor, and acid populations could coexist in dynamical equilibrium in three-dimensional space (if no malignant behavior is already in progress; that is, the time rates of change of the densities of the three populations are equal to zero). As for the dynamical stability of the fixed points, if all eigenvalues are negative, the fixed point is stable; however, if one or more of the eigenvalues is (are) positive, the fixed point is unstable. In a stable fixed point, after slight changes in the equilibrium cellular densities have taken place (due to e.g. cellular diffusion), the cellular densities tend to converge infinitesimally back to their initial equilibrium values (so the time rates of change of the cellular densities remain at zero after all); however, on the other hand, in an unstable fixed point, after the occurrence of slight changes in the equilibrium densities, significant divergences off the initial equilibrium values are likely to take place, where malignancies are likely to emerge; because, on balance, the presence of the acid specifically would always significantly corroborate the process into only one direction; that is, tumor invasion into normal tissue (not vice versa). Whenever a malignant behavior is already in progress, this implies that the time rate of change of the density of tumor cells is positive. A fixed point (either stable or unstable) physically represents an intermediate balance between the three populations where none of them is undergoing any noticeable growth into the other two. If the fixed point is stable, a malignant behavior is unlikely to emerge, even after the occurrence of slight perturbations of the equilibrium cellular densities at the interface (microenvironment) joining the three populations. On the other hand, an unstable fixed point where tumor cells are very likely to invade normal tissue upon any slight perturbations of the equilibrium cellular densities could potentially develop into a malignancy. Note that, in the reaction-diffusion model we present here, from a mathematical standpoint and within the acid-mediation hypothesis, only the dynamical stability of the equilibrium points (not the initial equilibrium density values) controls whether a prospective malignant behavior is likely; because the dynamical stability by itself controls what happens after any slight perturbations of the initial equilibrium density values have taken place. We have solved numerically for the fixed points (and their dynamical stability) of the reaction-diffusion system presented in Eq. (7) with the sample numerical values of the biological input parameters provided in Sec. 3. As demonstrated in Table 1, for physical fixed points (all density values are positive), our numerical results show that both stable (fixed point 3) and unstable (fixed point 4) configurations are possible. For the particular set of the biological input parameters we have operated on in Sec. 3, the existence of two different physical solutions means that wherever the three populations are found in a dynamical equilibrium in three-dimensional space, it must be in one of these two solutions.
As for future work, it would be interesting to study systematically how the fixed-point behavior depends on the biological input parameters. In simple systems this can be done analytically; but in systems as complex as the one introduced in Eq. (7), numerical methods are needed. Possibly, one can use numerical algorithms , designed to compute branches of stable and unstable solutions, and compute the Floquet multipliers that determine stability along these branches. The initial starting points for the computation of periodic orbits are generated automatically at the Hopf bifurcation points. Such techniques can also locate folds, branch points, period-doubling bifurcations, and bifurcations to tori. Along branches of periodic solutions, branch switching is possible at branch points and at period-doubling bifurcations.
Once the bifurcation diagrams in parameter space are understood in this way, dynamical simulations can be performed to determine tumor growth or regression rates.
I am grateful to Dr. Edward Gawlinski and Dr. Theodore Burkhardt for very thoughtful comments and discussions.
Copyright: Fouad AM. ©2019. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation: Fouad AM (2019). Dynamical Stability Analysis of Tumor Growth and Invasion: A Reaction- Diffusion Model. Oncogen 2(4): 20.