OPEN ACCESS Reproducible Model

The classic Boron and De Weer (1976) paper provided the ﬁrst evidence of active regulation of pH in cells by an energy-dependent acid-base transporter. These authors also developed a quantitative model — comprising passive ﬂuxes of acid-base equivalents across the cell membrane, intracellular reactions, and an active transport mechanism in the cell membrane (modelled as a proton pump) — to help interpret their measurements of intracellular pH under perturbations of both extracellular CO 2 / HCO − 3 and extracellular NH 3 / NH +4 . This Physiome paper seeks to make that model, and the experimental conditions under which it was developed, available in a reproducible and well-documented form, along with a software implementation that makes the model easy to use and understand. We have also taken the opportunity to update some of the units used in the original paper, and to provide a few parameter values that were missing in the original paper. Finally, we provide an historical background to the Boron and De Weer (1976) proposal for active pH regulation and a commentary on subsequent work that has enriched our understanding of this most basic aspect of cellular physiology.


Introduction
In 1976 Boron & De Weer published their landmark paper on "Intracellular pH transients in squid giant axons caused by CO 2 , NH 3 , and metabolic inhibitors" (Boron and De Weer, 1976). The authors used a squid giant axon preparation and a mathematical model of pH buffering and the transport of protons, bicarbonate (HCO − 3 ) and CO 2 to establish the experimental evidence for active regulation of intracellular pH (pH i ) by a transporter in the plasma membrane that -at the expense of energy -either moves acid out of the cell, or base into the cell. Today, we refer to such a transporter generically as an acid-extrusion mechanism. For simplicity, Boron & De Weer modelled it as a proton pump, although the result would have been almost indistinguishable had they modelled it as the uptake of HCO − 3 or carbonate (CO 2− 3 ). The paper reported on the consequences of adding and then removing extracellular CO 2 /HCO − 3 , NH 3 /NH + 4 (where NH + 4 is ammonium), or the metabolic inhibitors, cyanide, azide and dinitrophenol (DNP).
In the present paper, we re-formulate the models from Boron and De Weer (1976), henceforth referred to as 'BDW', and specify the simulation using the Physiome modelling standards CellML (Cuellar et al., 2003) and SED-ML (Bergmann et al., 2017) in order to ensure that the model reproduces the graphs in the original paper and that the model is fully curated. 2 Note that this effort requires the specification of some parameters used in BDW's simulations, but not described in the BDW paper. The curated and annotated model is made available in a form that users can run with OpenCOR 3 to understand the model and to explore the effect of changes in parameter values. A short exposure of the axon to extracellular NH 3 /NH + 4 causes a rapid rise in pH i , followed by a pH i decay that modestly undershoots (lower short arrow) its initial resting value upon removal of extracellular NH 3 /NH + 4 . A longer exposure of squid giant axons to extracellular NH 3 /NH + 4 causes a rapid rise in pH i , followed by a slow and sustained pH i decay. Removing extracellular NH 3 /NH + 4 causes pH i to undershoot substantially its initial resting value (long arrow). Both the plateau-phase acidification (upper short arrow) and the undershoot (long arrow) are indicative of net acid loading during the period of NH 3 /NH + 4 exposure. (B) Cartoon illustrating the processes underlying the initial alkalinisation phase in (A) for both short and long exposures to extracellular NH 3 /NH + 4 . The initial entry of NH 3 leads to the intracellular consumption of H + (and thus to the observed pH i rise) via the reaction NH 3 + H + −→ NH + 4 . (C) Cartoon illustrating the processes underlying the plateau-phase acidification during the long NH 3 /NH + 4 exposure in (A). After NH 3 equilibration across the plasma membrane (pH i peak in panel (A)), the slow entry of NH + 4 -which has always been present but overwhelmed by the influx of NH 3 -leads to the production of H + (and thus to the observed slow pH i decay during the plateau phase) via the reaction NH + 4 −→ NH 3 + H + . The newly formed NH 3 then exits the cell. The pH i undershoots observed upon removal of extracellular NH 3 /NH + 4 , during both short and long NH 3 /NH + 4 exposures, are the result of the accumulation of NH + 4 during exposure to extracellular NH 3 /NH + 4 . BDW used the mathematical model to postulate the above sequence of events, including both the plateau-phase acidification and the pH i undershoot. (A), modified from Boron and De Weer (1976). (B)-(C), modified from Boron (2010).
The total weak acid concentration, [TA], is the sum of [HA] and [A − ]. Note that [TA] is one of the two main unknowns in the BDW model for weak acids.
The dissociation of the cationic weak acid (BH + ) to the uncharged weak base (B) is described by the equilibrium reaction, where the equilibrium constant is An example is the NH + 4 dissociation reaction, The total weak base concentration, [TB], is the sum of [BH + ] and [B]. Note that [TB] is one of the two main unknowns in the BDW model for weak bases.

4/32
The CO 2 /HCO − 3 buffer pair. The formation of HCO − 3 and H + from CO 2 by hydration is given by the equilibrium reaction where the equilibrium constant is Taking logarithms of both sides of Equation 5, and recognising from Henry's law that where s is the solubility coefficient for CO 2 and p CO 2 is the partial pressure of CO 2 , we obtain the familiar Henderson-Hasselbalch equation Here, pH = − log[H + ] 5 and pK CO 2 = − log(K CO 2 ).
In terms of the nomenclature above, one might regard CO 2 as the weak acid HA 6 , and HCO − 3 as its conjugate base A − .
Buffering power (β ). By definition, β is the amount of strong base (e.g., NaOH), or the negative of the amount of strong acid (e.g., HCl), that one must add to 1 L of solution to raise pH by one pH unit: The units of β are mM. For additional details, refer to Roos and Boron (1981); Boron and Boulpaep (2016). Note that BDW defined β as a negative number, as did Koppel and Spiro in their original definition of buffering (Koppel, 1914;Roos and Boron, 1980), rather than as a now-conventional positive quantity, as did Van Slyke in his later work (Van Slyke, 1922). BDW's definition, which they consistently applied, has no effect on the outcome of their simulations. In the present paper, we will follow the definition of Van Slyke -defining β as a positive number -and make appropriate sign changes to the derived equations.

The Boron & De Weer Model for the Permeation by an Uncharged Weak Acid and its Conjugate, Anionic Weak Base
The BDW mathematical model consists of two time-dependent ordinary differential equations (ODEs), one describing the time-course of the concentration of total intracellular buffer ([TA] i = [HA] i + [A − ] i ) and the other the time-course of the intracellular free H + concentration (which is related to pH i ). BDW derived these two equations for the general cases in which any buffer pair HA/A − , or any buffer pair B/BH + , can move passively across the plasma membrane of a prototype cell. Then, they applied these two general equations to their specific experimental conditions, namely exposure of a cell (a squid giant axon) to equilibrated extracellular CO 2 /HCO − 3 or to equilibrated extracellular NH 3 /NH + 4 . Here, following BDW's approach, we begin by deriving the equations for HA/A − . In the next section, we apply the same general formalism to B/BH + . 5 The definition of pH as a pH scale based on powers of 10 was introduced by Sørensen in an attempt to simplify the notation of [H + ] and to avoid having to resort to decimals for tiny amounts of [H + ] (Sørensen, 1909). 6 Note that although CO 2 is often regarded as an acid, the true weak acid is H 2 CO 3 , the product of the reaction of CO 2 with H 2 O.
The constant field equation -also known as the Goldman, Hodgkin, Katz (GHK) (Goldman, 1943;Hodgkin and Huxley, 1952) where P A − (m · s −1 ) is the membrane permeability to the charged conjugate base A − , V m is the membrane potential (intracellular relative to extracellular potential), and is a shorthand for e −V m F /RT . Note that J HA and J A − have units of mol · m −2 · s −1 .
Although HA and A − can interconvert in the cytosol, BDW assumed that the intracellular concentration of total weak acid [TA] i only can change due to the transmembrane fluxes of HA and A − (see Figure 3). Thus, the time rate of change of [TA] i is where ρ (m −1 ) is the area-to-volume ratio for the cell, and converts the transmembrane flux per unit area (in units of mol · m −2 · s −1 ) to a time rate of change per unit cell volume (mol · m −3 · s −1 or mM · s −1 ). Equation 11 is the first of two ODEs of the BDW model for the buffer pair HA/A − . ] i depends on the net rate dQ /dt at which acids are added into the cytosol. BDW assumed that dQ /dt depends on (i) the release of H + by some fraction x of the entering HA (i.e., x J HA ), (ii) the consumption of H + by some fraction y of the entering A − (i.e., y J A − ), and (iii) the additional rate of intracellular H + consumption via metabolism or active acid extrusion (J H + ).
Later, Bevensee and Boron defined the time rate of change per unit volume (e.g., d[TA] i /dt) as a 'pseudoflux' φ, with the area-to-volume ratio folded into the value of φ (Bevensee and Boron, 2013 In their simple system of a squid giant axon exposed to CO 2 (i.e., HA) and HCO − 3 (i.e., A − ), BDW assumed that only three general processes affect dQ /dt : (i) the release of H + by some fraction (x ) of the entering HA (i.e., x J HA ), (ii) the consumption of H + by some fraction (y ) of the entering A − (i.e., y J A − ), and (iii) the "additional" rate of intracellular consumption or active extrusion of H + (J H + ; see Figure 3) above the fixed background rate of H + extrusion necessary to balance the fixed background rate of acid loading (i.e., addition of H + or equivalent acid, or removal of OH − or equivalent base) in the absence of HA/A − . Thus, This is also the fraction of entering A − that combines with H + and becomes HA. Combining the above expression with Equation 2, which BDW defined as α. Conversely, the fraction x of A − that remains A − is This is also the fraction of entering HA that dissociates to form A − and H + . Combining the above expression with Equation 2, In summary, Equation 12 becomes: BDW modelled J H + (mol · m −2 · s −1 ) as the additional proton-extrusion rate above the fixed background rate of pH i ). Note that k /ρ has units of (m · s −1 ), consistent with the membrane permeability terms P HA and P A − in Equation 9 and Equation 10. 8 The BDW authors used the definition of buffering power, in its infinitesimal form, to derive the relation between d[H + ] i /dt and dQ /dt , as shown in the following steps.
Our first goal is to obtain an expression for dpH/dt in terms of dQ /dt . According to the chain rule: By definition (see Equation 8), β = −dQ /dpH (mol · m −3 ), or equivalently Combining Equation 19 and Equation 20 where J HA (from Equation 9), and J A − (from Equation 10) are given by: where Simulation for CO 2 /HCO − 3 experiments. BDW employed Equation 28 and Equation 29 to simulate the experiments in which they exposed a squid giant axon to a solution containing equilibrated CO 2 /HCO − 3 . Their simulation protocol was a step change in (a) extracellular p CO 2 from 0 to 5% CO 2 (37 mmHg or, with s = 0.0321 mM · mmHg −1 , [CO 2 ] o = s .p CO 2 = 1.1877 mM) and (b) extracellular [HCO − 3 ] from 0 to 59.5260 mM (the value that [HCO − 3 ] o has in a solution containing 5% CO 2 at pH o of 7.70) 9 . The step change is applied for 2700 s (45 min) at constant pH o = 7.70. Table 1 and Table 2 report the parameter values used by BDW. Table 1 provides parameter values that are common to both the CO 2 /HCO − 3 and the NH 3 /NH + 4 experiments. Table 2 provides parameter values exclusive to the CO 2 /HCO − 3 experiments only. In the present work, the differential Equation 28 and Equation 29 -when coded in CellML and solved with OpenCOR -produce the plots in Figure 4. The simulation file Boron-CO2.sedml contains the computational setting for running the model. Open the .sedml file in OpenCOR and click Run Simulation. The initial conditions are [TA] i = 0 mM and pH i = 7.40. Note that Figure 4 illustrates the time courses not only of pH i -as presented by BDW -but also of quantities (e.g., various solute concentrations and fluxes) not displayed in the original paper; these values are useful for understanding the processes that contribute to the pH i transient. Moreover, our curated and annotated version of the BDW model also allows one to alter the parameter values from those originally chosen by BDW, thereby extending the ability of the user to investigate the predictive power of the computational model.   3 ] i continues to rise as pH i rises at a constant [CO 2 ] i (the proton pumping rate k is set to 300 s −1 , thus k /ρ = 0.0375 m · s −1 ). Note also that, after the removal of CO 2 /HCO − 3 , pH i rises to a higher value (∼ 8.15) than its starting value (∼ 7.4), indicating the net extrusion of acid from the cell during the CO 2 /HCO − 3 exposure.

11/32 4 The Boron & De Weer Model for the Permeation by an Uncharged Weak Base and its Conjugate, Cationic Weak Acid
Following an approach analogous to the one outlined above for weak acids, BDW derived two timedependent ODEs. Derivation for weak bases. Imagine that a cell is exposed to a solution containing equilibrated B/BH + , and that both B and BH + initially move into the cell -because of the chemical gradient in the case of B, and because of the electrochemical gradient in the case of BH + .
where ρ (m −1 ) is again the area-to-volume ratio for the cell. The equation is an integrated form of Fick's first law of diffusion that describes the net passive flux of B, and describes the net passive influx of BH + according to the GHK equation. In the previous two equations, P B (m · s −1 ) is the membrane permeability to the uncharged weak base B, P BH + (m · s −1 ) is the membrane permeability to the charged conjugate weak acid BH + , and is a shorthand for e V m F /RT . Equation 33 is the first of two ODEs of the BDW model for the buffer pair B/BH + .
The second equation of the BDW model for a weak base -analogous to Equation 27 above -is

12/32
where J H+ is the same as in Equation 18 and and Substituting for α, where  to simulate the experiments in which they exposed a squid giant axon to equilibrated NH 3 /NH + 4 . Their simulation protocol was a step change in extracellular NH 4 Cl from 0 to 9 mM (that is, a step change in [NH + 4 ] o from 0 to 8.86 mM, and in [NH 3 ] i from 0 to 0.14 mM) applied for 1500 s (25 min) at constant pH o = 7.70. 10 Table 1 and Table 3 report the parameter values used by BDW. Note that in the NH 3 /NH + 4 simulations, k is always zero, that is, J H + does not affect these processes.    figure 6B of the BDW paper, P NH + 4 had values of 0, 10 −6 , 10 −5 , and 10 −4 cm −1 . 3 BDW did not report the value of V m , but in their code used −55 mV, which matched the measured mean value. 4 In their code BDW used the value of 7.32 and not 7.30, which they used when constructing the Davenport diagram in their figure 5B.

Discussion
The publication of this retrospective paper provides an opportunity to clarify some concepts in the original paper that have benefitted from subsequent experimental and theoretical advances. We also provide some additional parameters missing from BDW and correct some minor errors. Most importantly, the curated model is now freely available on the Physiome website 11 in standardised form (CellML) that can be run in the open source software OpenCOR. A follow-up paper will be written that recasts the equations in bond graph form to facilitate their incorporation into more complex models where pH regulation is coupled with other cellular processes.

Historical Context
The 1976 Boron & De Weer paper introduced the first models to simulate the time course of pH i . By developing a predictive mathematical model based on first principles, BDW provided a quantitative basis for interpreting their new data on the time-dependent response of pH i to step changes in extracellular CO 2 /HCO − 3 (and HA/A − pairs in general) and NH 3 /NH + 4 (and B/BH + pairs in general). The models also provided a clear, quantitative basis for interpreting BDW's new data on how cells regulate their pH i , which BDW modelled as a pH i -dependent H + -extrusion mechanism. Below, we will introduce a broader concept termed "acid extrusion" (Boron, 1977). The first of the two BDW models elucidates how the transmembrane fluxes of a neutral weak acid and its anionic conjugate weak base affects pH i , with the acid-extrusion becoming increasingly important as pH i falls. The second model simulates how the transmembrane fluxes of a neutral base and its cationic conjugate weak acid affects pH i . Figure 6. Solution of the BDW model during and following a 1500 s period of externally applied NH 4 Cl. In these simulations [NH + 4 ] o = 9 mM. The intracellular fluid becomes alkaline as NH 3 enters (note the J NH 3 time course) and hydrates to form NH + 4 and OH − . Additional passive NH + 4 entry (note J NH + 4 time course) down its electrochemical gradient opposes the effect of the NH 3 entry and slightly reduces the pH i increase. Upon removal of NH 4 Cl, [NH 3 ] i and [NH + 4 ] i decay towards their original values, but pH i drops well below its original value of 7.40.
Of course, BDW were not the first to undertake a quantitative assessment of how acids or bases 15/32 affect, or are affected by, the pH of a solution. Below, we divide the earlier work into two major categories, (a) analyses of how neutral weak acids (and their anion conjugate weak bases) or neutral weak bases (and their cationic conjugate weak acids) affect pH in simple systems, and (b) analyses of how the distribution of HA/A − (or B/BH + ) across a barrier, such as a cell membrane, are affected by pH i and V m .

Development of the Concept of Buffering in Simple and Complex Systems
Buffering power. , Koppel (1914 introduced the first modern definition of the chemical buffering (i.e., "magnitude of moderation", or P ) of H + by a weak-acid/weak-base conjugate pair, and -based on first principles -derived an expression for buffering power. Because these authors defined P in terms of the amount of strong acid that one must add to a solution to produce a pH change, P is a negative number. Roos and Boron (1980) translated the Koppel-Spiro paper from its original German, and provided a historical context.
Initially unaware of the work of Koppel and Spiro, Van Slyke (1922) independently defined buffering power -to which he assigned the Greek letter β . Because he defined β in terms of the amount of strong base that one must add to a solution to produce a pH change (see Equation 8), β is a positive number. Although modest differences exist between the efforts of Koppel and Spiro on the one hand and Van Slyke on the other, they are quite similar. Nevertheless, it is Van Slyke's definition of β that has become the modern convention throughout chemistry and physiological chemistry.
Koppel & Spiro and Van Slyke quantitatively described how -in a one-compartment solution -weak acids, weak bases, ampholytes, and weak-acid/base mixtures can buffer added strong acid or strong base. In their analysis, the system both begins and ends in an equilibrium state. Of course, in their pioneering work, these authors had no reason to contemplate time courses or barriers separating more than one compartment.
In their work, BDW defined β as a negative number, as Koppel and Spiro defined their P .
The "Davenport" diagram. This nomogram (Boron and Boulpaep, 2016) is a powerful tool for graphically computing the effects of respiratory acid-base disorders (caused by changes in [CO 2 ] in a system open to CO 2 ) and metabolic acid-base disorders (caused by the addition or removal of HCO − 3 or a strong acid or base). The underlying assumption for the Davenport diagram is that the system is in equilibrium. The first component of a Davenport diagram (see Figure 7) is a plot of [HCO − 3 ] vs pH i for one or more values of [CO 2 ] -these are the CO 2 isopleths that describe the equilibrium among CO 2 , HCO − 3 , and H + . At any pH on any isopleth, the slope is the open-system CO 2 /HCO − 3 buffering power (β CO 2 ). The second component is a linear plot, on the same axes, of the concentration of all protonated forms of all non-HCO − 3 buffers vs pH. At any pH, the slope is the buffering power of all non-HCO − 3 buffers (β non−HCO − 3 ), and the line is termed the non-HCO − 3 buffer line. Its linearity implies that β non−HCO − 3 is insensitive to changes in pH. The intersection of an isopleth with the non-HCO − 3 buffer line describes the current state of the system, when both CO 2 /HCO − 3 buffer and non-HCO − 3 buffers are simultaneously in equilibrium. Davenport developed a series of rules for using this paradigm to interpret acid-base disorders, and these rules are well founded in physical chemistry.
The Davenport diagram traces its origins to the analyses of blood by several eminent investigators about a century ago. It was Henderson (1921) -as far as we are able to ascertain -who in figure 4 of his paper was the first to plot [HCO − 3 ] vs pH for two H 2 CO 3 (rather than CO 2 ) isopleths, and for two different values of β non−HCO − 3 (i.e., those produced by 10% and 100% HbO 2 ). The power of the Davenport approach is that, knowing the initial conditions and the pH dependence of β non−HCO − 3 , one can compute with fair accuracy (using the nomogram) or great accuracy (using a computer to solve the equations numerically) the result of virtually any acid-base disorder in the pathophysiological range for a system containing CO 2 /HCO − 3 and a mixture of non-HCO − 3 buffers. In their figure 5A (reproduced here as Figure 7), BDW used a Davenport approach to describe the initial steady state of a squid giant axon (point A at pH i ≈ 7.32, [CO 2 ] i = 0.1%), the initial effect of an exposure to increased [CO 2 ] i (point B, an example of intracellular respiratory . This nomogram consists of two kinds of plots. The first kind of plot is represented by the four CO 2 isopleths that slope upwards from left to right. Each isolpleth represents all possible combinations of [HCO − 3 ] i and pH i for a given [CO 2 ] i (described here as the % of the air phase that is CO 2 ). The second kind of plot is represented by the two lines that slope downwards from left to right. The slope of these parallel lines describes the buffering power of non-CO 2 /HCO − 3 buffers. BDW assumed that the experiment starts at point A, at pH i = 7.32 and 0.1% CO 2 . The addition of 5% CO 2 causes the pH i at equilibrium to fall to the point represented by point B. The extrusion of acid during the CO 2 /HCO − 3 exposure causes the system to move along the 5% CO 2 isopleth from point B to point C. Finally, upon removal of CO 2 /HCO − 3 , the system returns to the original CO 2 /HCO − 3 isopleth, but now at point D. The difference between points D and A represents the pH i overshoot. In the BDW paper, β -the slope of the lines in this figure -appeared to be 25 mM/pH unit. In fact, the value of β determined in the NH 3 /NH + 4 experiments (also shown as a Davenport-like diagram in figure 5B of BDW) was 9 mM/pH unit. The reason for this discrepancy was probably that BDW delivered the CO 2 /HCO − 3 solution to the axon using a peristaltic pump and Silastic tubing, which they later realised has a high CO 2 permeability. Thus, the [CO 2 ] reaching the axon was < 5%, accounting for the artificially inflated value for β . In their follow-up paper (Boron and De Weer, 1976), BDW delivered the CO 2 /HCO − 3 solutions from glass syringes and through Saran tubing, which has an extremely low CO 2 permeability. acidosis), the effect of the plateau-phase pH i recovery (point C, an example of intracellular compensatory metabolic alkalosis), and finally the effect of removing extracellular CO 2 (point D, an example of metabolic alkalosis) to account for the pH i overshoot. If one does not know β non−HCO − 3 , the Davenport diagram allows one to compute it from the initial and final pH . The numerical integration of the BDW equations -when P HCO − 3 and the acid extrusion rate are both zeroshould in principle yield, at infinite time, results consistent with the Davenport diagram.
In their paper (their figure 5B), BDW introduced a novel Davenport-like diagram for the NH 3 /NH + 4 buffer system, with [NH + 4 ] on the ordinate (replacing [HCO − 3 ]), NH 3 isopleths (replacing CO 2 isopleths), and a line describing non-NH 3 /NH + 4 buffering power (replacing the line describing non-CO 2 /HCO − 3 buffering power). Like the classical Davenport diagram, this one (or others like it, constructed for other buffer pairs) can be a useful tool for interpreting -in the steady stateproblems in acid-base chemistry.
In the Davenport analysis, the initial and final conditions represent equilibria. Davenport-like diagrams provides no information about the time course of pH between the initial and final states. Nor can Davenport-like diagrams deal with time course, barriers (e.g., cell membranes), permeabilities to substances other than the neutral molecule (e.g., CO 2 , NH 3 ), or active transport. Of course, in their pioneering work, Henderson, Davenport, and other authors contributing to this nomogram had no reason to contemplate these future complexities.

Pre-BDW Analyses of Transmembrane Distributions of Weak Acids and Bases
As summarised by Roos and Boron (1981), about a century ago, several authors -who assumed that CO 2 equilibrates across the plasma membrane but that HCO − 3 is impermeant -used the sum 3 ] i to compute the steady-state pH i of several cell types. Then, beginning in 1940, a series of authors introduced three successively more sophisticated mathematical analyses for the steady-state transmembrane distribution of a neutral weak acid and its anionic conjugate weak base (i.e., TA), and three more for the distribution of a neutral weak base and its cationic, conjugate weak acid (i.e., TB). We will now present these analyses in order of increasing complexity, and according to their sequence in time (see Figure 8). They all have in common the assumption that the system is either in equilibrium or at least in a steady state supported by the input of energy. Figure 8. Timeline of acid-base/pH models prior to BDW. Note that the time-dependent BDW model of weak acid/conjugate weak base (HA/A − ) collapses to Roos (1965Roos ( , 1975) steady-state model for HA/A − . The Roos (1965Roos ( , 1975 Boron and Roos (1976) steady-state model for B/BH + . The Boron and Roos (1976) model in turn reduces to the Orloff and Berliner (1956) model for B/BH + when V m approaches zero. Finally, the Orloff and Berliner (1956) model reduces to the Jacobs-like model -where only one uncharged species B is permeant -as (P B ) approaches zero.
Jacobs (HA, neutral weak acid). Jacobs (1940) presented a general mathematical model that describes the equilibrium transmembrane distribution of total weak acid (TA), assuming that only the neutral species (HA), but not the conjugate weak base (A − ) can cross the membrane: [TA] i [TA] o = 10 pH i −pK + 1 10 pH o −pK + 1 .
In an accompanying document, we show the derivation of the above equation 12 , the linear form of which is: We also present the derivation of this linear version in an accompanying document 13 . A major assumption in the derivations of both Equation 44 and Equation 45 is that HA does not merely move but fully equilibrates across the cell membrane. That is, the system is in equilibrium.
The notion that only the neutral member of the buffer species (i.e., HA) can traverse the membrane is known as nonionic diffusion.
In an accompanying document, we show the derivation of Equation 46, 14 the linear form of which is: We also present a derivation of Equation Orloff & Berliner (B/BH + , V m = 0). Orloff and Berliner (1956) extended the Jacobs-like model to a neutral weak base and its cationic conjugate weak acid by permitting not just B but also BH + to permeate a barrier separating compartments 1 and 2. They avoided the complication of BH + electrodiffusion by assuming a transmembrane voltage of zero (i.e., V m = 0). Because they recognised that the flux of BH + across the barrier would cause pH to drift in opposite directions in the two compartments, they assumed that independent, energy-requiring processes would stabilise pH in the two compartments and establish a steady state described by

19/32
As P BH + approaches zero, Equation 49 reduces to Equation 47 -the Jacobs-like equation for a weak base. The accompanying document 17 also shows the derivation of Equation 49, as well as the mathematics that shows the limit of this expression as P BH + approaches zero.
Note that the flux of BH + -which leads to a flux of B in the opposite direction -tends to push the system off equilibrium. As noted above, the Orloff-Berliner equation requires that the system be in a steady-state, which can be achieved, as they recognised, only by an input of energy to maintain all relevant concentrations constant over time.
Milne et al (HA/A − , V m = 0). Milne and colleagues (1958) developed a steady-state expression similar to that of Orloff & Berliner, but for the transmembrane distribution of a neutral weak acid and its anionic conjugate weak: As P A − approaches zero, Equation 51 reduces to Equation 45 -the Jacobs equation for a weak acid. The accompanying document 19 also shows the derivation of Equation 51, as well as the mathematics that shows the limit of this expression as P A − approaches zero.

Roos (HA/A − , non-zero V m ).
In 1965 and 1975, Roos extended the Irvine model by allowing V m to assume non-zero values (Roos, 1965(Roos, , 1975, and derived the following equation, where has the same meaning as in the derivation of the BDW equations: e −V m F /RT . An accompanying document 20 shows the derivation of Equation 52. This document also shows that, at the limits of certain parameters, Equation 52 simplifies to the expected equation: 1. As V m → 0, the Roos equation simplifies to the equation of Milne et al (which assumes V m = 0). 2. As P A − → 0, the Roos equation simplifies to the Jacobs equation (which assumes P A − = 0). 3. As P HA → 0, the Roos equation simplifies to the Nernst equation (which assumes permeability to only A − ).
The Roos equation was important historically because it allowed one to assess possible errors in pH i values computed from the transmembrane distribution of a neutral weak acid (e.g., CO 2 , DMO) and its monovalent anion conjugate weak base. These errors could in principle arise from membrane permeability to A − (as already considered by Milne et al), as influenced by V m .
conjugate weak acid: where has the same meaning as in the derivation of the BDW equations: e V m F /RT . An accompanying document 21 shows the derivation of Equation 53. This document also shows that, at the limits of certain parameters, Equation 53 simplifies to the expected equation: 1. As V m → 0, the Boron-Roos equation simplifies to the equation of Orloff & Berliner (which assumes V m = 0). 2. As P BH + → 0, the Boron-Roos equation simplifies to the Jacobs-like equation for a neutral weak base (which assumes P BH + = 0). 3. As P B = 0, the Boron-Roos equation simplifies to the Nernst equation (which assumes permeability to only BH + ).
The Boron-Roos equation was important historically because it allowed one to assess possible errors in pH i values computed from the transmembrane distribution of a neutral weak base (e.g., methylamine) and its monovalent anion conjugate weak base. These errors could in principle arise from membrane permeability to BH + , as influenced by V m . In their paper, Boron and Roos used the transmembrane distribution of methylamine/methylammonium to monitor a downward drift in pH i caused by the passive influx of methylammonium. This was the first use of a chemicaldistribution technique to follow pH i changes over time.
We have already noted for the Orloff-Berliner equation that the derivation requires that the system be in a steady-state, which can be achieved only by an input of energy to maintain all relevant concentrations constant over time. The same is true for the equations of Milne et al, Roos, and Boron & Roos.

Comparison of the Pre-BDW Analyses with the BDW Equations
The BDW equations build on the earlier work on buffering and transmembrane distributions of weak acids and bases, presented in the previous two sections. Of course, the work of Koppel and Spiro, and Van Slyke, as well as their predecessors who developed the physico-chemical principles of acid-base chemistry, is at the heart of the BDW model.
An important aspect of the Davenport diagram is its predictive power. For example, given β non−HCO − 3 as well as the initial pH and [CO 2 ], the Davenport approach can predict the effect of an increase in [CO 2 ] on the final equilibrium conditions. However, the Davenport approach makes no statement about mechanism or time course. The BDW approach has all the predictive power of Davenport, but also addresses mechanism and time course.
The six pre-BDW approaches for assessing transmembrane distributions of HA/A − and B/BH + all start with the weak acid/base present and the system in an equilibrium or at least in a steady state. Somehow, the system -the cell and its surrounding fluid -went from a condition with no weak acid/base present to a condition with the weak acid/base present at equilibrium/steady state. The older models make no attempt to describe how and how fast the system achieved the new state, and -unlike the Davenport approach -have no predictive value for relating initial and final conditions. It is worth noting that the investigators who developed these six approaches were interested mainly in using tracer levels of weak acids/bases to compute pH i . The intention was that tracer levels would have minimal effects on the state of the system -hence, the minimal interest in the prediction. Note that, at infinite time (and with no acid extrusion), the BDW equations reduce to those six transmembrane-distribution analyses presented in the previous section.
To some extent, the BDW models represent a merger of the Davenport and the six pre-BDW approaches for assessing transmembrane distributions of HA/A − and B/BH + . Like the Davenport approach, the BDW approach is predictive. However, unlike Davenport's approach, the BDW models provide insight into mechanism and time course, and are applicable even when the system is far from equilibrium/steady state. Like the six transmembrane-distribution models, the BDW models provide insight into how pH i and V m affect [HA] i , [A − ] i , [B] i , and [BH + ] i . Unlike the six transmembrane-distribution models, the BDW models are predictive and provide insight into mechanism and time course.
It is worth noting that BDW derived their equations under the influence of Albert Roos, who had derived the transmembrane-distribution model for HA/A − (Roos, 1975) and who inspired Boron's derivation of the B/BH + model (Boron and Roos, 1976).

Post-BDW Development of Concepts
Fundamental law of pH i regulation. Recall that one of the intermediate steps of the derivation of the BDW equations for the HA/A − system was Equation 21, which we reproduce here: Another intermediate step was the description of dQ /dt in Equation 17, which we reproduce here: Combining the above two equations yields a primitive form of the fundamental law of pH i regulation: In the later literature, Boron and colleagues in effect dissected the dQ /dt term -the net rate at which H + appear in the cytosol (mol · m −3 · s −1 ) -into two concepts that are more general than those considered by BDW in Equation 56: the intracellular acid-loading rate (J L ) and the intracellular acid-extrusion rate (J E ).
In the narrow definition of Equation 56 the only J L term is (1 − α)J HA . In the post-BDW literature, Boron and colleagues defined J L to comprise every process that adds the equivalent of H + to or removes the equivalent of OH − from the cytosol, including H + channels (mediating passive H + influx), a variety of transporters (e.g., Cl-HCO 3 exchangers mediating HCO − 3 efflux), and the metabolic production of acid. Boron et al. (1979) measured and introduced the term "rate of acid introduction". Roos and Boron (1981) later replaced this term when they coined "acid-loading rate".
In the narrow definition of Equation 56, the J E terms are α J A − and J H + . Note that BDW defined J H + as the H + -extrusion rate above baseline. In the post-BDW literature, Boron and colleagues defined J E to comprise every process that removes the equivalent of H + from or adds the equivalent of OH − to the cytosol, including a variety of transporters (e.g., H + pumps, Na-H exchangers) that mediate H + efflux and others (e.g., Na + -driven HCO − 3 transporters, H + /lactate cotransporters) that mediate base efflux. It appears that Boron in 1977 provided the first definition of acid extrusion.
Recasting Equation 56 in terms of J L and J E ,

22/32
The modern version of the fundamental law of pH i regulation is: Here, with the explicit inclusion of ρ, J E and J L have the units (moles · cm −2 · s −1 ).
Equation 58 tells us that pH i is stable when J E = J L , rises when J E > J L , and falls with J E < J L . BDW did not derive the above equation in their 1976 paper. The first statement of the concept of Equation 58 was in sentence form in the review by Roos and Boron (1981), who also provided a graphical example in their figure 12. By 1989, Boron presented a version of Equation 58 that lacks the term ρ. Thus, he implicitly defined J E and J L in terms of (moles · cm −3 · s −1 ). This 1989 paper also included examples of how Equation 58 could help one interpret the time course of pH i in many experimental settings. For example, a bolus introduction of acid into a cell -an "acute" intracellular acid load, coined by Boron by the time of this 1989 review -rapidly lowers pH i , but also increases J E and decreases J L . Because J E now exceeds J L , pH i recovers to its initial value. By 1992, Boron coined the term "fundamental law of pH i regulation" to describe Equation 58.
Note that an acute or bolus intracellular acid load is to be distinguished from a chronic intracellular load (the rate of which is J L ). The amount of an acute acid load can instantly be quantitated in acid equivalents. The amount of a chronic acid load can only be quantitated in acid equivalents by integrating J L over time.
Later, Boron (2004) provided a more detailed description of pH i regulation, as well as a tongue-incheek analogy -which he had been using for years in lectures -between pH i regulation and the temperature regulation of a house, complete with multiple furnaces (acid-extruders), multiple air conditioners (acid loaders), heat capacity (buffering power), a thermostat (pH i sensitivity of the transporters), and weather radar (extracellular sensors for CO 2 , HCO − 3 , and pH). Flux vs pseudoflux. One can define the acid-loading and acid-extrusion rates as strict fluxes, with units of moles/(membrane area × time). Such of use of J E or J L -for example J NBC (the acid-extruding flux mediated by a Na/HCO 3 cotransporter, NBC) or J Cl−HCO 3 (the acid-loading flux mediated by a Cl-HCO 3 exchanger) -is most appropriate in cases in which the cell has a simple geometry (e.g., a squid axon or Xenopus oocyte). However, for cells with complex geometrieswhere surface-to-volume ratios are difficult to define -physiologists often present experimental data in the units "moles/(volume of cell water)". To avoid confusion between the two systems of measurement, Bevensee and Boron (1995) introduce the term "pseudoflux" and the symbol φ (rather than J ). Thus, in their paper, the authors referred to φ E and φ L (rather than J E and J L ), so that the fundamental law of pH i regulation becomes if both B and BH + cross the membrane in the same direction, how will their isodirectional fluxes affect pH on the cis side (the side from which they depart) and the trans side (the side to which they go)?

Post-BDW Models of Acid-Base Fluxes/Chemistry
Following the development of the BDW model, the field of pH regulation has seen several modelling efforts. Here, we summarise several, with emphasis on those models that, in our opinion and to the best of our knowledge, have contributed to advance the field, either by extending the BDW model or by introducing new modelling paradigms.  (Keifer and Roos, 1981). These authors were interested in the transmembrane fluxes of the neutral weak acid DMO and its conjugate weak base. Keifer & Roos assumed that -at the end of the infinitesimal time increment during which HA and A − entered/left the cellthe cytosolic HA, A − and H + re-equilibrated. They used this variation on the BDW approach to estimate P HA and P A − , finding that the plasma membrane permeability to HA is ≈ 10 3 higher than to A − . In the appendix of their paper, Keifer The first of these two equations is identical to that presented by BDW, whereas the second includes the new terms (K new ).
In the appendix of their paper, Keifer Note that as [TA] i approaches zero, or as β approaches infinity, the new terms (K new ) approach zero, and thus the refined BDW model collapses to the original BDW model. Keifer & Roos did not consider net H + efflux (i.e., active acid extrusion).
In order to test whether the Keifer-Roos refinement really improves the predictions of the original BDW model, we implemented the BDW model with and without the Keifer-Roos refinement and employed the two models to simulate how pH i changes when (a) only CO 2 enters the cell (i.e., HCO − 3 permeability is zero) and (b) no proton pumping is present (i.e., J H is zero). With assumptions 'a' and 'b', the BDW model will take the system to equilibrium, where we can compare the predicted final pH i with that produced by the Davenport diagram -our gold standard for the value of pH i at equilibrium. Using the parameter values of Table 1 and Table 2, we find that both the refined BDW model and the Davenport diagram predict final pH i values of 6.972, whereas the original BDW model predicts a slightly lower pH i of 6.964 -lower because the original BDW does not incorporate as much self-buffering (as Keifer & Roos termed it).
Extending BDW to an epithelial cell. As part of their study of the CO 2 permeability of gastric-gland cells,   /dt ). Second, after each step of the integration, they re-equilibrated HA, A − , and H + in the cytosol -a maneuver equivalent to the Keifer-Roos extension of the BDW model. Third, they modelled two separate extracellular fluids (each an infinite reservoir), the equivalent of a luminal solution facing one half of the cell that represented the apical membrane, and a basolateral solution facing the other half of the cell that represented the basolateral membrane. The geometry of the cell was still cylindrical, like the squid axon. In their review,  provide additional information about this modelling, which showed that -to account for the physiological data in the paper by  -the (membrane area) × (CO 2 permeability) product must be > 1000-fold greater for the basolateral than the apical side of the epithelial cell. Thus, this modelling was a critical part of the main conclusion of the paper by , which identified the apical membranes of the gastric chief and parietal cells as the first known membranes with negligible CO 2 permeability.
Model of pH i regulation by the Vaughan-Jones group. Following BDW's approach for modelling the net rate of acid addition into the cytosol, Leem et al. (1999) developed the first comprehensive mathematical model of pH i regulation in cardiac myocytes. They simulated the experimentally observed pH i recovery from an acid load (obtained with the ammonium prepulse technique, introduced by BDW) or a base load (obtained with the acetate-prepulse technique). In their model, Leem and coworkers incorporated acid-base fluxes mediated by four sarcolemma transporters, two acid extruders (i.e., Na-H exchangers and Na/HCO 3 cotransporters) and two acid loaders (i.e., Cl-HCO 3 exchangers and hypothetical Cl-OH exchangers) as well as a time-dependent intracellular buffering by the CO 2 /HCO − 3 system. Their simulations nicely reproduced the pH i changes that they observed in rather complex physiological experiments.
Model of H + diffusion by the Vaughan-Jones group. Richard Vaughan-Jones and his colleagues have extensively investigated the spatial variability of intracellular H + diffusion in a series of experimental studies (e.g., confocal measurements of pH i ) complemented by mathematical modelling (Vaughan-Jones et al., 2002;. These authors developed two-dimensional diffusion models of intracellular H + diffusion from a constant source to estimate the apparent diffusion constant of H + in cardiac ventricular myocytes. They found that the apparent H + mobility is low in cardiac cells because of the presence of mobile and immobile buffers. Later, Swietach and coworkers developed two-dimensional reaction-diffusion models, which they solved numerically in the longitudinal direction only (Swietach et al., , 2005. This approach allowed them to investigate further the spatial variability of intracellular H + diffusion via a proposed H + shuttling among intracellular mobile buffers. More recently, this group developed a reaction-diffusion model of pH regulation in tumor spheroids to study the role of carbonic anhydrase (CA) IX in facilitating CO 2 excretion from tumor cells, thus favoring tumor-cell survival and proliferation (Swietach et al., 2008(Swietach et al., , 2009). Models of the Gros group. Gerolf Gros and his colleagues have employed mathematical modelling of acid-base physiology to estimate the apparent membrane permeability of the red blood cell (RBC) membrane to CO 2 (P CO 2 ) and HCO − 3 (P HCO − 3 ), as well as extra-and intracellular CA activity, from mass-spectrometric data and the 18 O-exchange technique (Wunder et al., 1997). Their compartmental model comprises a system of six ordinary differential equations that account for the 18 O-exchange reactions among HCO − 3 , CO 2 and H 2 O in both the intracellular and extracellular compartments, and for transfer of CO 2 and H 2 O across the intracellular and extracellular compartments.
In 2005, Endeward & Gros used this modelling framework to estimate P CO 2 and P HCO − 3 of the apical membrane of the colonic epithelium in guinea pigs. They found that the apical membrane has a quite low P CO 2 , possibly because of the absence of membrane proteins (e.g., aquaporin 1, AQP1) that act as conduits (i.e., gas channels) for the movement of CO 2 .
In 2009, these same authors extended their previous compartmental model of 18 O-exchange to a one-dimensional reaction-diffusion model of an RBC surrounded by an extracellular unconvected layer that can exchange solutes with a well-stirred bulk solution (Endeward and Gros, 2009). They employed the model to assess the influence of intracellular and extracellular unconvected layers in the estimation of P CO 2 RBC membranes. The combination of physiological experiments and modelling allowed them to estimate the extra-and intracellular unstirred layers, and confirm their earlier conclusions that the P CO 2 of RBC membranes is low in the absence of gas channels.

26/32
Models of the CWRU group. With the goal of assessing the role of CO 2 channels in producing decreases in pH i and transient alkaline transients in extracellular-surface pH (pH S ) as CO 2 enters a Xenopus oocyte (Endeward et al., 2006;Musa-Aziz et al., 2009a, 2010, the group at Case Western Reserve University developed a three-dimensional reaction-diffusion model of CO 2 influx into a spherical cell (Somersalo et al., 2012). The model accounts for a multitude of buffer reactions, as well as solute diffusion within the unconvected intracellular fluid (ICF) and the extracellular unconvected fluid (EUF) that surrounds the cell. Because the electrophysiologists established experimental conditions in which only CO 2 can cross the plasma membrane, the modellers allowed only CO 2 to diffuse between the ICF and EUF. However, all solutes can diffuse within the ICF as well as between the EUF and the bulk extracellular fluid (bECF) that surrounds the EUF. Somersalo and coworkers employed the model to investigate a variety of theoretical conditions. For example, they investigated how changes in the width of the EUF or in the value of P CO 2 affect pH i and pH S transients. The model suggests that, in oocytes, the background P CO 2 (i.e., in the absence of channels) must be low in order for CO 2 channels to have the observed effect on the maximal rate of intracellular acidification, (dpH i /dt ) max , or on the maximal height of the alkaline pH S transient, (∆pH S ) max .
In 2014, Occhipinti and coworkers extended the above theoretical model and employed it to investigate the role that CA II and CA IV play in enhancing transmembrane CO 2 fluxes (Musa-Aziz et al., 2014a,b;Occhipinti et al., 2014). The model was able to mimic the pH i and pH S experiments under a variety of experimental conditions. Moreover, it provided a novel observation that the effects on cytosolic CA II and extracellular CA IV on (dpH i /dt ) max and (∆pH S ) max are supra-additive. More recently, the CWRU group further extended the oocyte model to include the special microenvironment that the pH S electrode creates when pushed against the oocyte membrane to record the pH S transient (Calvetti et al., 2020). This model, solved using finiteelement method, predicts that the special microenvironment between the blunt tip of the pH S electrode and the oocyte membrane greatly amplifies the alkaline (∆pH S ) max as CO 2 enters the cell.

The Future
The advice, "It's tough to make predictions, especially about the future" is attributed to Yogi Berra, amateur philosopher and legend of American baseball. His advice certainly applies here. We imagine that the grand goal of the acid-base segment of world-wide biomedical science would be to model whole-body acid-base regulation on a cell-by-cell basis for a wide range of physiological and pathophysiological challenge. Of course, we all expect modelling to provide new insights that we can test in physiological experiments. However, the models should ultimately provide medical professionals with powerful insights into pathophysiology. These goals represent an enormous challenge in multi-scale modelling -multiscale in both time and space. That lofty goal is probably many decades away, even with appropriate grant support and continued development of computer hardware and software.
The process could begin at the single-cell level. Even here, the challenge is daunting because, as pointed out by Occhipinti and Boron (2015), the movements of acid-base equivalents across the plasma membrane create complex interdependencies among buffer reactions, diffusion processes and transporter mechanisms. Intuition alone cannot explain or predict the consequences of these numerous simultaneous processes on pH. Advanced mathematical models -including details on processes or more complex geometries -performed in conjunction with physiological experiments can provide a useful tool to make predictions and provide mechanistic explanations on observed pH changes. The pH community has begun taking the first steps with cells of simple geometry, like oocytes (i.e., spheres). Even so, current models only simulate passive diffusion of substances like CO 2 and lactic acid. Current reaction-diffusion models do not yet incorporate electrodiffusion, let alone electrodiffusion through particular channels. The next step might be to incorporate specific acid-base transport systems (e.g., Na/HCO 3 cotransporters, Cl-HCO 3 exchangers) with characteristic kinetic properties. Ultimately, one would like to include all traffic (i.e., including non-acid-base traffic) across the cell membrane, and the regulation of this traffic. Model validation would require extensive physiological studies on these simple cell systems to verify that simulations of ion conductances (and ultimately channels) and transporters are reasonable. Moreover, even these relatively simple first steps would require efficient computational methods for solving the governing equations.
From the level of the single cell, the field must move in two opposite directions. In the reductionist direction, we must move from whole cells to nanodomains (i.e., the mesoscopic level) and ultimately to single molecules. At first, we guessed at even the number of different kinds of such proteins. Now we know the identities of specific proteins and their amino-acid sequences, and sometimes even their structures. Even so, cells are often capable of making many protein variants from one gene, and regulating these in response to myriad influences. Even the protein structures are merely an early step in understanding transporter mechanism.
In the opposite, integrative direction, the community must extend the models to more complex geometries, like the cells of a simple epithelium, spindle-shaped cells (e.g., muscle), cylinders (e.g., axons, dendritic processes), and ultimately cells of very complex geometry (e.g., neurons, astrocytes). In real life, even a simple epithelium is far more complex than a cuboid. The apical membrane (nearest the tight junctions) of a cell of the proximal tubule (PT) has microvilli. The basolateral membrane (nearest the blood) has infoldings. These complexities create nanoenvironments on both the cytosolic and extracellular sides of the membrane. The tight junctions that separate the PT cells allow the paracellular diffusion of water and certain solutes. The lateral interspaces between adjacent cells also creates nanoenvironments. All of these nanoenvironments are likely to be extremely important for the physiology of the epithelium. Further complexities include the myriad cellular organelles. These not only affect diffusion of water and solutes through the cytosol (by creating tortuosity), but also can engage in their own collection of transport processes that can affect pH i and contribute to buffering. Several groups have developed either steady-state or time-dependent compartmental models of acid-base transport in PTs and, more generally, solute reabsorption in various segments of the nephron (Thomas and Dagher, 1994;Krahn and Weinstein, 1996;Thomas et al., 2006;Weinstein et al., 2007;Weinstein and Sontag, 2009;Layton and Edwards, 2014;Edwards and Layton, 2017). Bransen and coworkers have used a finite-element approach to solve the partial differential equations that describe their reaction-diffusion model of Na-H exchange in the microvilli of the PT (Brasen et al., 2014).
Beyond cells, we must create whole tissues (e.g., a cylindrical PT) from many cells, and create whole organs (e.g., a kidney) from these tissues. Models of a whole kidney, for example, would have to consider the interstitial fluid and capillaries, the complex exchanges of substances among them, and the changes in composition that occur as fluids flow along the lumens of the tubule and capillaries. Comparable modelling must extend to every tissue and organ, and consider the complex intercommunications among them via the circulatory and neuro-endocrine systems as well as metabolism. In the case of acid-base physiology, one must model the specific roles that certain organs -the brain, the pulmonary system, and the kidneys -play in whole-body pH regulation.
All along the way, community members must cooperate because no one group of investigators can possibly accomplish the ultimate goal single-handedly. The models must be open source, modular and sharable -and the community must share. Finally, in order for modelling approaches to be consistent, the community must establish standard approaches for gathering physiological data, and agree on standard values for physiological parameters. Here, an organisation like the International Union of Physiological Sciences (IUPS) 22 can play a critical role.