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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {O}_2$$\end{document}O2 consumption with lactate and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CO}_2$$\end{document}CO2 production and connects these processes to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}^+$$\end{document}H+ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {HCO}_3^-$$\end{document}HCO3- 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {H}^+$$\end{document}H+ exchanges in well-oxygenated regions, and of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {HCO}_3^-$$\end{document}HCO3- 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {pH}_i$$\end{document}pHi equalizer at any \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {O}_2$$\end{document}O2 level in cells that grow in acidic extracellular environments.

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][2][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][2][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][2][3][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 ( CO 2 ) to protons ( H + ) and bicarbonate ( HCO − 3 ) ions. While the H + ions contribute to the acidity of the extracellular milieu, HCO − 3 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 .
Scientific RepoRtS | (2020) 10:13613 | https://doi.org/10.1038/s41598-020-70396-1 www.nature.com/scientificreports/ 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][7][8] . In particular, our previous models provide values for the rates of glucose and oxygen uptake, lactate and CO 2 production and lactate/H + 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][7][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][7][8] , here we must follow the dynamics of CO 2 , HCO − 3 and H + , both inside and outside a tumor cell. The inputs of the model are the rates of lactate and CO 2 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 H + ions, and both ions are transported through the cell membrane by means of the bi-directional monocarboxylate transporters MCT [6][7][8] . We remark that this part of the model impacts the rate of change of both intracellular and extracellular pH (from now on pH i and pH e , respectively), and oxygen is assumed to diffuse freely through the cell membrane and its consumption rate is used to determine the rate of CO 2 production.
Intracellular H + ions are transported outside the cell by means of unidirectional sodium-hydrogen exchangers NHE 1 . Different HCO − 3 transporters on the other hand are known to drive the flux of bicarbonate ions through the cell membrane. Some of them import or export HCO − 3 by exchanging Cl − anions and the transport may depend or not on the presence of Na + cations 1 . Experimental works, however, have shown that the efficiency of HCO − 3 transport in different cell systems is quite similar, and that the import of HCO − 3 is fundamental in tumor cells where it is dominated by the activity of the Na + -dependent Cl − /HCO − 3 exchanger 9,10 . Therefore, we consider the import activity of a generic HCO − 3 transporter ( THCO 3 in Fig. 1) which, as a first approximation, assumes the average biochemical characteristics of the Na + -dependent Cl − /HCO − 3 exchanger. We finally model the activity 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, CO 2 and ATP. Lactic acid dissociates to lactate and H + ions, whereas CO 2 reversibly hydrates to HCO − 3 and H + . These chemical species diffuse through cell membranes ( CO 2 ) 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 ( THCO 3 ) 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 ( pH i ) and the extracellular ( pH e ) pH. Detailed information on each pathway is given in the main text. www.nature.com/scientificreports/ of the membrane-bound Carbonic Anhydrase 9 (CA9) enzyme that catalyzes, on the cell surface, the hydration of CO 2 . 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][12][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 nonphysiologic pH but inactive at pH > 8.0 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 [X] C,c is the molar concentration of a given chemical species inside ( [X] C ) or outside ( [X] c ) the cell, V max and K m are the Michaelis-Menten parameters and h is the Hill exponent ( h > 0).
We assume that: • CO 2 can freely diffuse through the cell membrane; • CO 2 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.
In this work the variables take the following units for length, mass and time, respectively: µm , pg and s. Molar concentrations (M) have always been converted to mass units by taking into account the volume of the cell ( V C , cell volume is computed by approximating a cell to a sphere of given radius r C ) or of the environment ( V c ) and the molecular mass (MW) of chemical species.
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.
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 ( pH i ) 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 pH i on pH e 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 pH i of SCK cells was measured in experiments where cells were also treated with Amiloride and DIDS inhibitors. Amiloride inhibits Na + 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 pH i regulation when the extracellular volume becomes higher than 10 4 cell volumes, i.e. when the extracellular volume exceeds ∼ 0.02 µl (the volume of 1 cell of radius ∼ 7 µm is ∼ 2 pl ). 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 pH i due to the activity of ion transporters is concerned.
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. CO 2 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 ≈ 0.1 in pH between small and large cells ( r C = 5.5 and 8.0 µm , respectively, i.e. a volume ratio of ≃ 3 ) but pH i levels reach values that have actually been observed in tumor cells 23 . With the www.nature.com/scientificreports/ initial conditions discussed above, the simulations approach equilibrium quite fast and this indicates that the numerical solution of model equations is stable. The model predictions for pH i 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 µm 24 , and computed pH i 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 �pH i /�t < 10 −5 was reached. In these simulations, the volume of the environment was set to V c = 10 12 µm 3 = 1 ml , i.e. large enough to assure nearly constant pH e values throughout the simulation runs. Figure 3 shows that model predictions are in excellent agreement with the experimental data.

Contribution of NHE and THCO3 transporters to pH
i 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 ν max parameters to 0. The results are shown in Fig. 4 where we plot the pH i 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. The simulations clearly show that under normoxic condition the contribution of the THCO3 transporter to pH i is negligible. Under this condition pH i is maintained to physiological levels thanks to the activity of NHE transporter that export H + ions outside the cells. On the contrary, THCO3 activity dominates in hypoxic environments.  www.nature.com/scientificreports/ Role of Carbonic Anhydrase 9. As previously noted by Swietach et al. 11 pH i 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 pH i spatially by facilitating CO 2 diffusion in the unstirred extracellular space of the spheroid 11 . This intriguing conclusion, supported by experimental evidence, suggests that CA9  Table 1. After an initial transient, pH i 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 pH e 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.   www.nature.com/scientificreports/ activity becomes important for the control of pH i 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 pH i regulation starts to become important at the extracellular to cell volume ratio V c /V C ≈ 10 4 and reaches a maximum at V c /V C ≈ 100 . 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 pH i runs out of control (see also Fig. 4). However, the results in Fig. 5 show that when V c /V C ≈ 100 and CA9 is active the extracellular pH at equilibrium is around 5.5 and pH i ≈ 6.6 , 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 pH i is regulated by cells growing in small environments, i.e. when the CA9 role is not negligible, and when O 2 levels are lower and lower. Figure 6 shows that when pH e ≥ 5.8 , CA9 acts as a nonlinear pH i equalizer at any O 2 levels. 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 pH i change ( dpH i /dt ) upon alkali or acid load was reduced to ∼ 25−45% in VHL + cells with respect to VHL − 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 VHL + and VHL − cells, proton fluxes were measured in cells exposed to Cl − -deprived solutions, during recovery from NH + 4 -induced cell acidification or subjected to hypertonic shock 27 . Simulations of NH + 4 -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., NH 4 Cl 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 Cl − -deprived solutions. www.nature.com/scientificreports/ The rationale behind cell treatment with Cl − -deprived solutions was the discovery that VHL expression in 786-O renal cancer cells increased mRNA and protein levels of Cl − /HCO − 3 AE2 anion exchanger by 3.5 fold, although the apparent cell surface expression of AE2 was similar in VHL + and VHL − cells as evaluated by immunostaining 27 . The AE2 transporter exchange Cl − with HCO − 3 , and when the cells are exposed to Cl − -deprived solutions Cl − can only exit from the cells thus forcing HCO − 3 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 HCO − 3 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 HCO − 3 efflux (see Eq. 8 in the "Methods" section). We then perturbed the system at equilibrium by suddenly switching to 0 the rate of HCO − 3 efflux to model cells placed in Cl − -deprived baths or switching it to normal values to model cells re-placed under standard environmental conditions (see Fig. 7).
The expression of many genes is altered in VHL + 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 VHL − 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 ∼ 80 − 100-fold change) lactate transport in other cell systems 29 . Our model can easily take into account the phenotype of VHL + cells as far as these pathways are concerned. We multiplied specific rates by appropriate factors: the rate of proton production ( gH + in Eq. 11) was divided by 2 to model the reduced lactate/H + production by glycolysis; the rate of CO 2 production ( gCO 2 in Eq. 11) was multiplied by 2 to model the increased respiration rate; finally the maximum rate of lactate transport through MCT transporters ( ν maxMCT 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 ( dpH i /dt ) is reduced to ∼ 40% in VHL + cells with respect to VHL − 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.

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][31][32][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 H + and HCO − 3 transporters with Na + /K + -ATPase and Na + , K + and Cl − ion fluxes, while 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][7][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][7][8] . This has already allowed us to characterize new biophysical properties of tumors and of their microenvironment [34][35][36][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 H + efflux from tumor cells dominates the control of intracellular acidity in normoxic environments, whereas HCO − 3 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 ≈ 150 µm , i.e. ≈ 10 cell diameters, from blood vessels 38 . Thus, within this short distance the control of pH i is attained by tumor cells through a switch from H + export to HCO − 3 import pathways. This observation gives further support to recent work that has shown that inhibition of HCO − 3 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 HCO − 3 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 pH i equalizer at any O 2 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 pH i values throughout the spheroid structure due to CA9 activity in spheroids grown up to ≈ 500 µm 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  Fig. 7. VHL protein expression was observed to downmodulate the expression of glucose transporters and to reduce the glycolytic rate, and hence lactate/H + production, in renal cancer cells 28 . Here we plot the rate of pH i change in simulated VHL + cells as the function of proton production rate through glycolysis (rate gH + in Eq. 11, see the "Methods" section and the Supplementary material). The standard value of gH + is calculated (see Eq. 11) as gAcL · (MW H /MW AcL ) where the value of gAcl is given in Table 1. When gH + = 1 , dpH i /dt in VHL + cells is equal to that of VHL − cells.
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 CO 2 . It is well known that at physiologic temperature (i.e. ∼ 37 • C ) carbonic acid dissociates very quickly and represents less than 0.5% of the total carbon dioxide and bicarbonate ion 44 . Thus, the hydration of CO 2 can be approximated by the following chemical reaction: The values of the two rate constants k 1 and k 2 have been determined in cells under standard culture conditions in two independent experiments with good agreement 11,18 . We take the values in 18 : k 1 ≃ 0.144 s −1 and k 2 ≃ 1.9 · 10 5 M −1 s −1 .
We compare model outputs with experimental data obtained with cell cultures in vitro, in a standard atmosphere at 37 • C and 5% CO 2 at 1 atm pressure. To compute the initial density of CO 2 dissolved in water under these conditions we use Henry's law c = k(T)P where c is the molar concentration of the gas in water, P the pressure and k(T) is a function of temperature with T = 298.15 K , k = 3.3 · 10 −4 mol m −3 Pa −1 and − sol k R = 2400 K (see ref. 45 for further details); we find that the initial density of CO 2 in cell medium under standard culture conditions is: Finally, given the CO 2 concentration we find the density of HCO − 3 ions from the Henderson-Hasselbach equation: where pKa = − log 10 (k 1 /k 2 ) ≃ 6.12.
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 H + ions inside and outside the cells. where J 1→2 is the flux from 1 to 2 in units of concentration over time and surface area S C , P M,CO 2 is the permeability of the carbon dioxide and C i is the concentration of CO 2 in the i-th volume. Since we model cells grown in an incubator at constant CO 2 pressure, the CO 2 concentration can reach values far from equilibrium only inside cells because of the oxygen consumption by cell metabolism and of the equivalent CO 2 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 CO 2 due to diffusion is: MCT transporters. The MCTs are a family of bidirectional H + 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 www.nature.com/scientificreports/ (see refs. [6][7][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 H + inside and outside the cell: where ν maxMCT = V maxAcL · MW H MW AcL · S C , K mMCT = K mAcL · MW H MW AcL and where the ratio of molecular weights is used to rescale the equations from concentrations to masses.
In Eq. 2, a2c H and c2a H depend, respectively, on extracellular and intracellular pH, and phenomenologically describe the dependency of MCT transport activity on acidity (for a complete analysis see [6][7][8] ): nHe transporters. Sodium-hydrogen exchangers (NHE) are membrane transport proteins that exploit the influx of Na + to export H + ions. The sodium concentration gradient is maintained by the ATP-dependent Na + /K + 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 H + 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 SensO 2 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 H + from the cell to the environment due to NHE transport is modeled by the equation: where ν maxNHE = V maxNHE · S C and fPHe NHE 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 H + 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 Na + -dependent Cl − /HCO − 3 exchanger appears to dominate HCO − 3 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 HCO − 3 ions to buffer their internal pH 9,10 , and that this is a common property of different cancer cells. Experimental studies have demonstrated that HCO − 3 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 HCO − 3 transport depends on ATP availability. However, just as observed for proton export by NHE transporters, HCO − 3 transport proceeds by parallel fluxes of ions, like Na + and Cl − , along their electrochemical gradients that are actively maintained by cells through energy-consuming paths. Thus, it is likely that even HCO − 3 transport is controlled by ATP availability, albeit indirectly. On the basis of these considerations we model HCO − 3 import as follows: where ν maxTHCO3 = V maxTHCO3 · S C and the two functions fpHe THCO3 and fpHi THCO3 phenomenologically describe how HCO − 3 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 VHL + and VHL − cells we make HCO − 3 transport bidirectional by considering HCO − 3 efflux from cells as follows: (2) www.nature.com/scientificreports/ 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][12][13] . It is a membrane-tethered enzyme and it is mainly found at the external surface of cells where it catalyses the hydration of CO 2 [11][12][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 ν maxCA9 = V maxCA9 · S C and h CA9 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 gH + = gAcL · MW H /MW AcL and gCO 2 = qO 2 · MW CO 2 /MW O 2 are, respectively, the rates of H + and CO 2 production that are proportional to the rate of lactate production gAcL and oxygen consumption qO 2 as defined in our previous work [6][7][8] , and all the other rates, and regulatory functions, are given in equations 1-10. The multiplicative factor 10 3 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 µm 3 .
In in silico experiments with VHL + and VHL − cells, where HCO − 3 transport is bidirectional, the differential equations in the set 11 that describe HCO − 3 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 CO 2 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++ (8) ν in→out THCO3 = SensATP · fpHe THCO3 · fpHi THCO3 · ν maxTHCO3 · m HCO − 3 ,C V C · MW HCO3 · K mHCO3 + m HCO − 3 ,C (9) ν CA9 = h CA9 · ν maxCA9 · m CO 2 ,c V c · MW CO2 · K mCA9 + m CO 2 ,c (10) h CA9 = 3 + 2 · tanh (−δ CA9 · SensO2) (11) dm CO 2 ,C dt = gCO 2 − k 1 · m CO 2 ,C + k 2 · m H + ,C · m HCO − 3 ,C · 10 3 · MW CO 2 V C · MW H · MW HCO3 − J C→c · S C · MW CO 2 dm H + ,C dt = gH + + k 1 · m CO 2 ,C · MW H MW CO 2 − k 2 · m H + ,C · m HCO − 3 ,C · www.nature.com/scientificreports/ 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 ǫ < 10 −6 ) within a maximum number of iterations fixed at N max = 1000.