The control of acidity in tumor cells: a biophysical model.

Acidosis of the tumor microenvironment leads to cancer invasion, progression and resistance to therapies. We present a biophysical model that describes how tumor cells regulate intracellular and extracellular acidity while they grow in a microenvironment characterized by increasing acidity and hypoxia. The model takes into account the dynamic interplay between glucose and [Formula: see text] consumption with lactate and [Formula: see text] production and connects these processes to [Formula: see text] and [Formula: see text] fluxes inside and outside cells. We have validated the model with independent experimental data and used it to investigate how and to which extent tumor cells can survive in adverse micro-environments characterized by acidity and hypoxia. The simulations show a dominance of the [Formula: see text] exchanges in well-oxygenated regions, and of [Formula: see text] exchanges in the inner hypoxic regions where tumor cells are known to acquire malignant phenotypes. The model also includes the activity of the enzyme Carbonic Anhydrase 9 (CA9), a known marker of tumor aggressiveness, and the simulations demonstrate that CA9 acts as a nonlinear [Formula: see text] equalizer at any [Formula: see text] level in cells that grow in acidic extracellular environments.

Introduction

Acid homeostasis in animal tissues is achieved by active dynamic processes. In physiological conditions, the pH of tissues is maintained between 7.35 and 7.45 in spite of constant metabolic acid production by cells. At the microscopic level, cells must finely regulate their own internal pH to around 7.2 to avoid death[1-3]. Cellular acid homeostasis is carried out by active transport of acid/base equivalents across the cell membranes into the extracellular spaces.

Dysregulation of pH is a well-known hallmark of solid tumors[1-3]. The tissue of solid tumors is characterized by the presence of an irregular network of blood vessels, causing a spatially heterogeneous delivery of nutrients such as glucose and oxygen to tumor cells[1-4]. As the consequence, the inner regions of solid cancers that are distant from blood vessels become hypoxic and acidic. Cancer cells adapt to such adverse environments through a series of molecular changes that involve an increased expression of nutrient and ion transporters and enzymes (reviewed in[1,3,5]). For example, hypoxia activates the Hypoxia Inducible Factor-1 (HIF-1) that up-regulates the transcription of glucose transporters and of enzymes involved in glucose metabolism. Because of hypoxia, glucose is converted mainly to lactic acid through the glycolytic pathway to produce energy under the form of ATP, and the increased production of lactate reduces the pH of the extracellular spaces. A drop in intracellular pH, in turn, increases the activity of lactate and of various ion transporters that collectively contribute to recover intracellular acid homeostasis[1,3,5]. Hypoxia also causes the increased expression of some membrane-bound enzymes such as Carbonic Anhydrase (CA) that, on the cell surface, catalyzes the hydration of carbon dioxide () to protons () and bicarbonate () ions. While the ions contribute to the acidity of the extracellular milieu, ions can be transported back into the cells and increase the buffering potential of the intracellular environment[1,3,5], further contributing to maintain the intracellular pH at normal values.

It has recently been pointed out[1,3] that changes in the control of intracellular and extracellular acidity in the tissue of solid tumors are associated with many phenotypic changes of cancer cells with important implications in tumorigenesis, cancer progression, cancer diffusion, escape from immune surveillance and resistance to therapies. For example, microscopic examination of the tumor/normal tissue interface shows that peritumoral acidity drives tumor invasion in the surrounding normal tissue, with the regions of highest tumor invasion corresponding to those of lowest pH. In these regions the environmental pH reaches values that are toxic for normal but not for tumor cells[2].

Biophysical models can help to disentangle the intricate relationships between regulatory biochemical networks and give support to the interpretation of experimental evidence which is rapidly accumulating in this field. In this paper we describe a comprehensive biophysical model of the control of acidity in tumor cells. We study the action of key molecular actors in acid homeostasis of cancer cells, and investigate to which extent hypoxia and environmental acidosis influence their behavior. We focus on the dynamic interplay between lactate, proton, bicarbonate transporters and CA enzyme, and their regulation by oxygen and both extracellular and intracellular pH. The model includes the bicarbonate buffer that acts both in the extracellular and intracellular milieux and it incorporates results from our previous modeling efforts concerning tumor cell metabolism[6-8]. In particular, our previous models provide values for the rates of glucose and oxygen uptake, lactate and production and lactate/ transport across cell membranes through specific transporters that have already been validated with experimental data. Finally, we fix the model parameters by combining information from a number of experiments carried out with different tumor cell systems.

Results

Preliminary considerations, model assumptions and parameters

We start from the rather detailed model of tumor cell metabolism and growth that we developed in our previous research[6-8] which successfully reproduces the observed behavior of tumor cells in both liquid (e.g. blood tumors) and solid tumors. In particular, for the current work we have excerpted from that model the part that describes the rates of glucose conversion to lactic acid and oxygen consumption. We remark that the model in[6-8] has been set up with the minimal set of chemical and biochemical pathways that drive the dynamics of metabolism and that are common to most, if not all, tumor cells.

Unlike the metabolic model in[6-8], here we must follow the dynamics of , and , both inside and outside a tumor cell. The inputs of the model are the rates of lactate and production (Fig. 1) that depend on how cells take up nutrients, such as glucose, and convert them to ATP through the glycolytic and the oxidative phosphorylation pathways. Lactic acid dissociates immediately to lactate and ions, and both ions are transported through the cell membrane by means of the bi-directional monocarboxylate transporters MCT[6-8]. We remark that this part of the model impacts the rate of change of both intracellular and extracellular pH (from now on and , respectively), and oxygen is assumed to diffuse freely through the cell membrane and its consumption rate is used to determine the rate of production.

Figure 1

Layout of the model of acidity control in tumor cells. A cell takes up from the environment nutrients and oxygen which are then converted by cell metabolism to lactic acid, and ATP. Lactic acid dissociates to lactate and ions, whereas reversibly hydrates to and . These chemical species diffuse through cell membranes () or are actively transported outside and eventually inside the cell by means of specific protein transporters. We consider monocarboxylate (MCT), sodium-hydrogen exchanger (NHE) and generic bicarbonate () transporters. We also model the activity of the membrane-bound enzyme Carbonic Anhydrase 9 (CA9). Chemical reactions are indicated by solid lines and the regulatory pathways by dashed lines. Proton concentrations inside and outside the cell are used to compute the intracellular () and the extracellular () pH. Detailed information on each pathway is given in the main text.

Intracellular ions are transported outside the cell by means of unidirectional sodium-hydrogen exchangers NHE[1]. Different transporters on the other hand are known to drive the flux of bicarbonate ions through the cell membrane. Some of them import or export by exchanging anions and the transport may depend or not on the presence of cations[1]. Experimental works, however, have shown that the efficiency of transport in different cell systems is quite similar, and that the import of is fundamental in tumor cells where it is dominated by the activity of the -dependent exchanger[9,10]. Therefore, we consider the import activity of a generic transporter ( in Fig. 1) which, as a first approximation, assumes the average biochemical characteristics of the -dependent exchanger. We finally model the activity of the membrane-bound Carbonic Anhydrase 9 (CA9) enzyme that catalyzes, on the cell surface, the hydration of . This is an important path since CA9 has been found to be expressed by many solid tumors of different histotypes, and its activity has been correlated to tumor progression and growth[11-13].

It should be noted in Fig. 1 that we do not take into account a possible effect of the extracellular pH on CA9 activity. Previous work has shown that CA9 in cell membrane extracts is sensitive to low pH and is completely inhibited at pH 6.0[14], its pH sensitivity being much steeper than that of other CA isoforms[15]. This observation, however, is at odd with findings obtained using high-resolution techniques with purified enzyme: they showed that the catalytic domain of human CA9, but not of other isoforms, is stable and active still at very low non-physiologic pH but inactive at [16]. Because of these discrepancies, and since we do not want to focus on some specific cell system but rather to keep the model as general as possible, we decided to leave off the possible pH sensitivity of CA9 from the present model. Our modelling strategy is flexible enough to incorporate additional specific details when available, provided they are based on firm experimental conclusions.

We model the kinetics of ion transporters, and of CA9 activity as well, with the Michaelis-Menten/Hill formalism that is described by the following general equation:where is the molar concentration of a given chemical species inside () or outside () the cell, and are the Michaelis-Menten parameters and h is the Hill exponent ().

We assume that:In this work the variables take the following units for length, mass and time, respectively: , pg and s. Molar concentrations (M) have always been converted to mass units by taking into account the volume of the cell (, cell volume is computed by approximating a cell to a sphere of given radius ) or of the environment () and the molecular mass (MW) of chemical species.

can freely diffuse through the cell membrane;

diffusion is driven by the concentration gradient across the membrane and its only important component is the one directed normally with respect to the cell membrane;

the diffusion kinetics of charged ions through the cell membrane are much slower than the kinetics of the other processes in which they are involved, and thus the diffusion of charged ions is negligible;

the mixing of all chemical species in the cell and in the external environment is instantaneous;

within the short characteristic times of the considered chemical reactions the cell volume is constant.

The model defined by the set of differential equations 11 has several parameters. We extensively searched the scientific literature to find their values, and when these values were not directly available they were obtained by fit of specific equations to reported experimental data. Experimental evidence was also used to model regulatory functions given by Eqs. 3, 5, 7 and 10 that tune the activity of transporters and CA9 enzyme as the function of local pH, ATP and/or oxygen availability. The full strategy is detailed in the Supplementary Material and all parameter values are listed in Table 1.

Table 1

Values of model parameters

Parameter	Value	Unit	Reference
MW_H	1	gmol^-1	–
MW_CO_2	44	gmol^-1	–
MW_O_2	32	gmol^-1	–
MW_HCO_3	61	gmol^-1	–
MWAcLa	90.1	gmol^-1	–
PM,CO_2b	3.2×10^4	μms^-1	[17]
gAcL	3.8×10^-4	pgs^-1	[6]
qO_2	3.5×10^-5	pgs^-1	[6]
k_1	0.144	s^-1	[18]
k_2	1.9×10^5	M^-1s^-1	[18]
V_maxAcL	9.58×10^-5	pgs^-1μm^-2	[8]
K_mAcL	0.405×10^-3	pgμm^-3	[8]
a2c_H_slope	1.5	–	[8]
a2c_H_thr	7	–	[8]
c2a_H_slope	1.5	–	[8]
c2a_H_thr	7	–	[8]
V_maxNHE	5.15×10^-7	pgs^-1μm^-2	Fit of data in[9]
K_mNHE	0.196×10^-6	M	Fit of data in[9]
h	2.67	–	Fit of data in[9]
λ_NHE	0.076	–	Fit of data in[9]
pH_0,NHE	7.1	–	Fit of data in[9]
V_maxTHCO3	2.02×10^-5	pgs^-1μm^-2	Fit of data in[19]
K_mTHCO3	7.38×10^-3	M	Fit of data in[19]
λ_THCO3	1.63	–	Fit of data in[9]
pHe_0,THCO3	6.85	–	Fit of data in[9]
γ_THCO3	4.2	–	Fit of data in[10]
pHi_0,THCO3	6.90	–	Fit of data in[10]
V_maxCA9	9.47×10^-2	pgs^-1μm^-2	[20]
K_mCA9	7.2×10^-3	M	[21]
δ_CA9	7.3	–	Fit of data in[22]

AcL=lactic acid/lactate.

Parameter values have been determined and fixed as described in the Supplementary Material section

Once determined, parameter values were fixed and no further tuned to adapt model outputs to data. This means that the model has no free parameters and is strictly predictive. As explained in the next section, for validation purposes we first used it to predict how the intracellular pH () varies when cells are grown into environments with increasing acidity.

Model validation with independent experimental data

Model validation was performed with independent experimental data, i.e. data that were not used to set parameter values. To this end we used the data in the paper by Song et al.[23]. In this paper Song et al. investigated the dependence of on in SCK cells (human choloangiocarcinoma cell line) in standard in vitro cultures. To the best of our knowledge no data concerning the direct expression of specific ion transporters and CA9 in these cells are available. However, the of SCK cells was measured in experiments where cells were also treated with Amiloride and DIDS inhibitors. Amiloride inhibits channels and thus inhibits sodium-hydrogen exchangers, whereas DIDS inhibits all bicarbonate-dependent transport mechanisms (see Song et al.[23] and references cited therein). Thus, the expression of proton and bicarbonate transporters was functionally demonstrated in SCK cells. We do not know if SCK cells express CA9 but, as we shall see below (see Fig. 5), CA9 activity becomes negligible for regulation when the extracellular volume becomes higher than cell volumes, i.e. when the extracellular volume exceeds (the volume of 1 cell of radius is ). The experiments were carried out with cells kept under standard culture conditions where the extracellular volume is much higher, and thus it is irrelevant whether SCK cells express CA9 or not. Data obtained with SCK cells can therefore be used to validate the core model as far as the regulation of due to the activity of ion transporters is concerned.

Figure 5

regulation by CA9 for decreasing size of the extracellular volume. Cell radius was set to the average size of 6.5 . The inset shows values and the main panel to ratio for varying values (i.e. ratio of extracellular to cell volumes) when CA9 activity is turned on or off. In these simulations the extracellular environment is physically closed, i.e. the extracellular volume is unstirred and the diffusion of chemical species toward an “external reservoir” is not allowed.

The radius of SCK cells is not reported nor, to the best of our knowledge, it has been measured previously. This is important because our model equations take into account both cell volume (see Eqns. 1–11) and the cell surface (see e.g. diffusion, Eq. 1) that are computed from cell radius under the assumption that cell geometry can be approximated by a sphere. Thus we run simulations for different cell radii whose values were taken within a reasonable range for animal cells.

Figure 2 shows the model prediction for intracellular pH vs. cell size, under standard culture conditions. At equilibrium there is a difference of in pH between small and large cells ( and , respectively, i.e. a volume ratio of ) but levels reach values that have actually been observed in tumor cells[23]. With the initial conditions discussed above, the simulations approach equilibrium quite fast and this indicates that the numerical solution of model equations is stable.

Figure 2

Plot of as the function of time for cells with the indicated cell radii. We take values that are in the observed range for human tumor cells[24]. The model equations have been solved with the parameter values listed in Table 1. After an initial transient, reach an equilibrium at physiological values and this shows that the model (and its numerical solution) is stable and provides quantitative results in good agreement with actual experimental observations. We also plot for comparison. The extracellular pH does not vary because these runs were carried out for a limited time span and for cells growing in a large volume (1 mL) filled with fresh medium at physiological pH to mimic standard experimental conditions.

The model predictions for values in SCK cells grown in media with increasing acidity are shown in Fig. 3. We ran simulations with varying cell radius within a range of values which is reasonable for tumor cells, i.e. between 4.5 and 9 [24], and computed at equilibrium. As shown in Fig. 2 the numerical solutions approach equilibrium with slower kinetics for increasing cell radii. We chose a conservative criterion to define the equilibrium condition and we halted the simulations when was reached. In these simulations, the volume of the environment was set to , i.e. large enough to assure nearly constant values throughout the simulation runs. Figure 3 shows that model predictions are in excellent agreement with the experimental data.

Figure 3

Plot of for SCK cells grown in media with different values. Experimental data have been redrawn from figure 2 in[23] (closed circles). The lines show values at equilibrium as predicted by our model for the indicated cell radii. It is important to note that these are not fits because our model does not have free parameters. Equilibrium was reached at . The volume of the environment was set at .

Contribution of NHE and THCO3 transporters to in normoxic or hypoxic environments

We have used the model to study the biochemical mechanisms that allow tumor cells to survive to adverse environments. We have investigated the role of NHE and THCO3 transporters in the control of intracellular acidity by tumor cells exposed to normoxic or hypoxic environments. We ran several simulations by alternatively switching off the activity of NHE and THCO3 transporters, i.e. by setting the respective parameters to 0. The results are shown in Fig. 4 where we plot the values at equilibrium (see the previous section) as the function of environmental pH for cells grown under standard oxygen level or at 0.1 fraction thereof.

Figure 4

Contribution of NHE and THCO3 transporters to in normoxic (left panel) or hypoxic (right panel) environments. Simulations were run with the following parameters: cells radius and environmental volume . The intracellular pH was calculated at equilibrium (see also the legend to Fig. 3) as the function of the indicated values. The activity of NHE and THCO3 transporters was switched off by setting the respective parameters to 0. Environmental oxygen levels were tuned by setting the SensO2 parameter to 1 or to 0.1 (see the “Methods” section and the Supplementary Material for details). In both panels, dashed lines have been drawn to show the value at which , a value largely compatible with cell life (see also the experimental data in Fig. 3 for a comparison).

The simulations clearly show that under normoxic condition the contribution of the THCO3 transporter to is negligible. Under this condition is maintained to physiological levels thanks to the activity of NHE transporter that export ions outside the cells. On the contrary, THCO3 activity dominates in hypoxic environments.

Role of Carbonic Anhydrase 9

As previously noted by Swietach et al.[11] regulation is not affected by CA9 expression in isolated tumor cells, but its role becomes important when cells are grown as three-dimensional aggregates (tumor spheroids). When expressed by cells grown as tumor spheroids CA9 induces a near uniform intracellular pH throughout the structure[11], an observation that was explained by diffusion-reaction modeling as follows: CA9 coordinates spatially by facilitating diffusion in the unstirred extracellular space of the spheroid[11]. This intriguing conclusion, supported by experimental evidence, suggests that CA9 activity becomes important for the control of by tumor cells at critical sizes of the extracellular volume. We tested this hypothesis with our model, and the results are shown in Fig. 5.

The role of CA9 in regulation starts to become important at the extracellular to cell volume ratio and reaches a maximum at . It is important to note that we simulated cells that grow in a closed environment. This means that at small extracellular volumes the acidity of the environment becomes too high and runs out of control (see also Fig. 4). However, the results in Fig. 5 show that when and CA9 is active the extracellular pH at equilibrium is around 5.5 and , well within the physiological range.

Simulations in Fig. 5 do not take into account the oxygen levels in the tumor environment. As discussed above (see the “Methods” section) CA9 expression is regulated by hypoxia[22] and thus it is interesting to investigate how is regulated by cells growing in small environments, i.e. when the CA9 role is not negligible, and when levels are lower and lower. Figure 6 shows that when , CA9 acts as a nonlinear equalizer at any levels.

Figure 6

Role of CA9 on regulation for cells grown in a small environment with decreasing oxygen levels. In these simulations the extracellular volume was set to and cell radius to so that . Left panel: plot of at equilibrium as the function of for the indicated fractions of environmental . Right panel: same simulations as those shown in the left panel, but here we plot , i.e. the difference in when CA9 is turned on or off. This plot clearly shows the nonlinear character of CA9 activity in the regulation of .

The model as a tool for exploratory data analysis

Germ-line mutations that inactivate the von Hippel-Lindau (vhl) gene cause the VHL syndrome, a rare inherited disorder characterized mainly, but not only, by renal cancers[25,26]. The VHL protein drives ubiquitination and finally degradation of the hypoxia-inducible factor alpha (HIF) which in turn regulates a number of intracellular pathways that collectively confer resistance to hypoxia to cancer cells[25,26]. However, experimental findings suggest that both HIF-dependent and HIF-independent mechanisms are essential for VHL-mediated tumor suppressor effects[25,26].

Stable transfection of 786-O renal cancer cells with a full-length human vhl gene significantly decreased proton and bicarbonate fluxes with respect to vhl-null cells in spite of increased or unaltered expression of ion transporters[27]. In particular, experiments showed that the rate of change () upon alkali or acid load was reduced to in cells with respect to cells. A number of control experiments were carried out to test possible effects of VHL proteins in these cells, but the effects of VHL protein on ion fluxes remained unexplained[27]. Here we modify our model to provide a possible interpretation of these experimental observations.

In the experiments with and cells, proton fluxes were measured in cells exposed to -deprived solutions, during recovery from -induced cell acidification or subjected to hypertonic shock[27]. Simulations of -induced cell acidification and hypertonic shock would require major revision of the model to include a number of chemical, biochemical and morphological details such as, e.g., dissociation kinetics and intra- and extracellular flows of all involved ionic species, cell volume dynamics during osmotic shock and a detailed description of how cell shrinkage and swelling activate ions transport. In addition, quantitative information which is required to set the values of specific model parameters is not fully available, further hampering the development of specific detailed models. We therefore focus on cell treatments with -deprived solutions.

The rationale behind cell treatment with -deprived solutions was the discovery that VHL expression in 786-O renal cancer cells increased mRNA and protein levels of AE2 anion exchanger by 3.5 fold, although the apparent cell surface expression of AE2 was similar in and cells as evaluated by immunostaining[27]. The AE2 transporter exchange with , and when the cells are exposed to -deprived solutions can only exit from the cells thus forcing import[27]. In other words, the treatment makes an otherwise bidirectional transport unidirectional. Our simplified model takes into account only a generic unidirectional transporter that shuttle from the environment into the cell and that is described by Eq. 6 (see the “Methods” section). To model the AE2 exchanger we introduce one more equation to describe also the rate of efflux (see Eq. 8 in the “Methods” section). We then perturbed the system at equilibrium by suddenly switching to 0 the rate of efflux to model cells placed in -deprived baths or switching it to normal values to model cells re-placed under standard environmental conditions (see Fig. 7).

Figure 7

Rate of intracellular pH change in simulated and cells upon removal (left panel) and restoration (right panel) of environmental anions. Values have been normalized with respect to the maximum (left panel) or to the absolute value of the minimum (right panel) rate of change calculated for cells. Simulations were run until reached equilibrium. removal or restoration was modeled by suddenly (arrows) switching to 0 or to normal values, respectively, the rate of efflux (Eq. 8, see the “Methods” section) as described in the text. The initial rate of change in is reduced to of that of cells as observed in actual experiments.

The expression of many genes is altered in cells and VHL protein is known to affect several physiologic pathways[28]. Quantitative data are not fully available and therefore it is impossible with the present knowledge to reproduce the whole complex phenotype of these cells in silico. However, we note that among the physiologic pathways altered in 786-O cells expressing VHL proteins glycolysis and respiration are prominent[28]. Glycolysis was observed to be approximately one half of that measured for cells, a finding that was paralleled by a corresponding two fold downmodulation of glucose transporters, whereas respiration was found to be increased by a factor of two[28]. VHL expression was also observed to dramatically reduce (i.e. a -fold change) lactate transport in other cell systems[29]. Our model can easily take into account the phenotype of cells as far as these pathways are concerned. We multiplied specific rates by appropriate factors: the rate of proton production ( in Eq. 11) was divided by 2 to model the reduced lactate/ production by glycolysis; the rate of production ( in Eq. 11) was multiplied by 2 to model the increased respiration rate; finally the maximum rate of lactate transport through MCT transporters ( in Eq. 2) was divided by 80 to model the observed reduction of lactate transport.

The simulations show that the initial rate of intracellular pH change () is reduced to in cells with respect to cells (Fig. 7) in agreement with experimental observations[27]. As shown in Fig. 8, a reduced glycolytic rate is mainly responsible for this effect. This shows that the present model, although simplified, can still be adapted to simulate different cell phenotypes and used to suggest novel interpretation of otherwise paradoxical[27] and yet unexplained experimental observations.

Figure 8

Reduced glycolytic rates in cells might explain the effect shown in Fig. 7. VHL protein expression was observed to downmodulate the expression of glucose transporters and to reduce the glycolytic rate, and hence lactate/ production, in renal cancer cells[28]. Here we plot the rate of change in simulated cells as the function of proton production rate through glycolysis (rate in Eq. 11, see the “Methods” section and the Supplementary material). The standard value of is calculated (see Eq. 11) as where the value of gAcl is given in Table 1. When , in cells is equal to that of cells.

Discussion

We have developed a biophysical model to explore the complex molecular mechanisms that allow tumor cells to regulate both intracellular and extracellular acidity, but we are not alone, other modeling efforts have tried to capture the essential features of the biochemical pathways that lead to acid homeostasis in tumor cells (see e.g.[30-33]). We have taken the remarkable models described in[32] and[33] as our starting point, because of their direct applicability to the analysis of experimental data. The former provides a fully tractable quantitative description of the interplay between and transporters with -ATPase and , and ion fluxes, while the latter investigates the interaction of MCT transporters and CA9. We go a few steps further and model the network of important paths that connect together cell metabolism and hypoxia with transport of and ions and CA9 activity (see Fig. 1). The coupling of ion transport mechanisms with metabolism and hypoxia is essential if we want to understand how tumor cells grow and shape their microenvironment, an interplay that is of fundamental importance for the adaptation and evolution of cancer cells within a solid tumor. As mentioned at the beginning of the Results section, we have developed a computer program that successfully reproduces the growth and the behavior of tumor cells in both liquid and solid cancers[6-8]. It is a lattice-free model that contains a rather detailed description of tumor cell metabolism and of the cell cycle, as well as many other biochemical and biophysical features (e.g. cell mechanics, cell division, etc.)[6-8]. This has already allowed us to characterize new biophysical properties of tumors and of their microenvironment[34-37], but the program still contains an excindingly simplified description of how cells control their intracellular pH. The program has an incremental structure, and we add new parts as soon as they are independently validated. The present model is one of these parts, and once integrated in our previous software it will further increase its descriptive and predictive potential. We hope in this way to understand key biological features such as cell adaptation and evolution in tumor microenvironments and explore important aspects such as tumor cell resistance to therapies. Here we show that the present model can nonetheless be used as a tool for exploratory data analysis and for quantitative purposes.

We remark that with the model described here we are able to give a quantitative assessment of the importance of specific molecular mechanisms. For instance, simulations show that efflux from tumor cells dominates the control of intracellular acidity in normoxic environments, whereas import in hypoxic tumor areas (in our simulation where the fraction of oxygen decrease to 0.1 of standard values). Experiments have shown that in in vivo tumor micro-environments oxygen reaches 10% of its normal value at a distance of , i.e. cell diameters, from blood vessels[38]. Thus, within this short distance the control of is attained by tumor cells through a switch from export to import pathways. This observation gives further support to recent work that has shown that inhibition of fluxes inhibits the growth of experimental tumors by increasing intracellular acidity and cell death[39]. When we recall that the hypoxic regions are those where tumor cells show higher resistance to therapies, such as e.g. radiotherapy, then we see that approaches that aim at inhibiting fluxes would target the very cells that colonize the inner tumor regions and that would otherwise be resistant to therapies, and improve cancer control.

Finally, the model singles out the important role of CA9. The simulations show that CA9 acts as a nonlinear equalizer at any level in cells that grow in acidic extracellular environments. This result is in agreement with the experimental observations by Swietach and colleagues[11], collected with tumor spheroids. They observed near-uniform values throughout the spheroid structure due to CA9 activity in spheroids grown up to diameter. It has long been recognized that tumor spheroids of this size show steep gradients of oxygen with fractions that go as far down as 0 at the center of the spheroid[40]. Our simulations show that this is due to the concerted action of CA9 and of hypoxia that up-regulates CA9 expression. These two mechanisms collectively help cells to keep their intracellular pH under control because of increased production followed by import through THCO3 transporters.

Conclusion

While tumor cell adaptation and survival to extreme microenvironments are key concepts in oncogenesis[1-3], we remark that acid homeostasis is central to cellular adaptation in a much wider context. Active transport of acid/base equivalents across cell membranes into the extracellular spaces may cause transient and rapid changes of microenvironmental and cellular pHs like those observed for other ions involved in cell signalling. Indeed, pH transients have been shown to be important in intra- and inter-cellular communication in the nervous system and are known to affect a number of essential functions, like e.g. neuronal excitability and synaptic transmission[41]. This in turn implies that animal cells could sense and adapt to pH changes. The underlying molecular mechanisms are still not well understood, but the role of G-protein coupled receptors in proton sensing is increasingly investigated also in relation to pathological conditions, besides cancer, that result in an increased extracellular acidity, such as infarction and inflammation[42]. We conclude that our model can be used as an essential building block of more comprehensive in silico research on solid tumors[43], but it may also help understanding how other cells can sense and dynamically adapt to pH changes.

Methods

Bicarbonate buffer and initial conditions

Central to the whole scheme of reactions shown in Fig. 1 is the hydration of . It is well known that at physiologic temperature (i.e. ) carbonic acid dissociates very quickly and represents less than 0.5% of the total carbon dioxide and bicarbonate ion[44]. Thus, the hydration of can be approximated by the following chemical reaction:The values of the two rate constants and have been determined in cells under standard culture conditions in two independent experiments with good agreement[11,18]. We take the values in[18]: and .

We compare model outputs with experimental data obtained with cell cultures in vitro, in a standard atmosphere at and 5% at 1 atm pressure. To compute the initial density of dissolved in water under these conditions we use Henry’s law where c is the molar concentration of the gas in water, P the pressure and k(T) is a function of temperaturewith , and (see ref.[45] for further details); we find that the initial density of in cell medium under standard culture conditions is:Finally, given the concentration we find the density of ions from the Henderson-Hasselbach equation:where .

Where not otherwise specified, we fixed the standard intracellular and extracellular pH at 7.4, which determines the initial value of the molar concentration of ions inside and outside the cells.

diffusion through the cell membrane

Given the assumptions above, the component of normal to the cell membrane is described by the Fick’s first law:where is the flux from 1 to 2 in units of concentration over time and surface area , is the permeability of the carbon dioxide and is the concentration of in the i-th volume. Since we model cells grown in an incubator at constant pressure, the concentration can reach values far from equilibrium only inside cells because of the oxygen consumption by cell metabolism and of the equivalent production. This means that in the present model there is only a net outward flux of carbon dioxide from cells to the environment. Thus, the net flux of due to diffusion is:

MCT transporters

The MCTs are a family of bidirectional and lactate co-transporters expressed at the cell membrane and their activity has been shown to depend on the pH values on both sides of the cell membrane (see refs.[6-8] and references therein). We model their activity with parameter values extrapolated from experimental observations[6-8] and we use the following equations and parameters to describe the rate of transport of inside and outside the cell:where , and where the ratio of molecular weights is used to rescale the equations from concentrations to masses.

In Eq. 2, and depend, respectively, on extracellular and intracellular pH, and phenomenologically describe the dependency of MCT transport activity on acidity (for a complete analysis see[6-8]):

NHE transporters

Sodium-hydrogen exchangers (NHE) are membrane transport proteins that exploit the influx of to export ions. The sodium concentration gradient is maintained by the ATP-dependent pump[19,46] so that the activity of NHE indirectly depends on ATP availability. This implies that as long as ATP is available the flux of due to NHE is essentially unidirectional. It has also been reported that NHE activity is inhibited by hypoxia[10,19] and that, in the long-term, hypoxia inhibits the expression of NHE proteins. Energy and oxygen tune NHE activity and as in the previous model of tumor cell metabolism and growth[6-8], here we take into account these regulatory circuits by means of the two variables SensATP and that assume real values in the interval [0, 1].

Experimental observations indicate that NHE activity is described by a Hill equation[9,47,48] and hence the unidirectional flux of from the cell to the environment due to NHE transport is modeled by the equation:where and is a phenomenological function that tunes the activity of NHE transport as a function of extracellular pH:Indeed, it has been observed that extracellular acidity enhances transport through NHE[9,19,49]. In the Supplementary Material we discuss how we fix parameter values and define the function fPHe on the basis of experimental observations.

Transport of bicarbonate ions

As discussed above, we model the activity of a generic bicarbonate ion importer (THCO3). The -dependent exchanger appears to dominate fluxes in tumor cells[9,10], and therefore we take this transporter as a reference to set the values of parameters and fix general biochemical characteristics. This is an important part of the model, because it has been shown that tumor cells do actively import ions to buffer their internal pH[9,10], and that this is a common property of different cancer cells. Experimental studies have demonstrated that import is regulated by both intracellular and extracellular pH but not by hypoxia and that the transport follows a simple Michaelis-Menten kinetics. In the scientific literature there are no indications, as far as we can tell, that transport depends on ATP availability. However, just as observed for proton export by NHE transporters, transport proceeds by parallel fluxes of ions, like and , along their electrochemical gradients that are actively maintained by cells through energy-consuming paths. Thus, it is likely that even transport is controlled by ATP availability, albeit indirectly. On the basis of these considerations we model import as follows:where and the two functions and phenomenologically describe how import is affected by extracellular and intracellular pH, respectively. These functions have been fit to actual experimental data (see the Supplementary Material) and are modeled by the following equations:In in silico experiments with and cells we make transport bidirectional by considering efflux from cells as follows:

Activity of Carbonic Anhydrase 9

The enzyme CA9 is expressed by cells of many different solid tumors, and in general its expression correlates with cancer aggressiveness and poor therapeutic outcome[11-13]. It is a membrane-tethered enzyme and it is mainly found at the external surface of cells where it catalyses the hydration of [11-13]. Importantly, its expression is regulated by hypoxia and indeed CA9 is a marker of hypoxia[22]. Again, experimental observations show that CA9 activity follows a Michaelis-Menten kinetics. Thus:where where and is a phenomenological functions that describe how hypoxia tunes CA9 expression:This is a function of the fraction of available oxygen which, in our model, is defined by SensO2, and it describes the fold change in CA9 expression as observed in actual experiments (see the Supplementary Material).

The full model and its numerical integration

The full model is represented by the following set of differential equations:where and are, respectively, the rates of and production that are proportional to the rate of lactate production gAcL and oxygen consumption as defined in our previous work[6-8], and all the other rates, and regulatory functions, are given in equations 1–10. The multiplicative factor that appears in the right-hand side of equations 11 above comes from the conversion of standard molar concentration units to the units used here where masses are expressed in pg and volumes in .

In in silico experiments with and cells, where transport is bidirectional, the differential equations in the set 11 that describe kinetics were modified as follows:The system of differential equations 11 is nonlinear and stiff because it incorporates processes with different kinetics, from the fast kinetics of hydration and diffusion to the relatively slow kinetics of ion transport and enzyme activity. The system cannot be solved analytically and appropriate numerical approaches are required. We previously investigated this aspect within the context of complex large-scale biophysical models[50] and found that the implicit Euler method is well-suited for the numerical integration of models of this kind. We solved the discretized system of differential equation 11 using the implicit Euler algorithm followed by the Newton-Raphson method to solve numerically the resulting system of nonlinear equations. The code has been implemented in C++ using the computational framework provided by the GNU Scientific Library[51]. We used the standard Newton-Raphson solver gsl_multiroot_fsolver_dnewton and the gsl_multiroot_test_residual library to test the convergence of the algorithm (threshold ) within a maximum number of iterations fixed at .

Supplementary information.