pH Buffering and the Bicarbonate System

To solve this problem, express the above sum for total carbon dioxide concentration [$\Sigma$CO$_2$] in terms of the unknown [CO$_2$] and the known pH:

$$ [\Sigma\textrm{CO}_2] =  [\textrm{CO}_2] + K_{eq} [\textrm{CO}_2] + \frac{K_2[\textrm{H}_2\textrm{CO}_3]}{[\textrm{H}^+]} $$

or

$$ [\Sigma\textrm{CO}_2] =  [\textrm{CO}_2] + K_{eq} [\textrm{CO}_2] + \frac{K_2K_{eq}[\textrm{CO}_2]}{[\textrm{H}^+]} = [\textrm{CO}_2]\left( 1 + K_{eq} + \frac{K_2K_{eq}}{[\textrm{H}^+]} \right) .$$

We can solve this problem in MATLAB using the following commands:

% Input constants

K1 = 10^(-6.8);

K2 = 10^(-3.60);

Keq = 2.7e-3;

pH = 7.4;

H = 10^(-pH); % H+ activity

SCO2 = 30e-3; % total CO2 (M)

% Computed aqueous CO2

CO2 = SCO2 / (1 + Keq + K2*Keq./H ) 

Your solution will require simultaneously solving three conservation relationships. The first is for total phosphate

$$ P_0 = [\textrm{H}_2\textrm{PO}_3^-] + [\textrm{H}\textrm{PO}_3^{2-}] . \; \; \; (1)$$

The second is for total carbon dioxide

$$ [\Sigma\textrm{CO}_2] = [\Sigma\textrm{CO}_2]_0 + \Delta\textrm{CO}_2 = [\textrm{CO}_2] + [\textrm{H}_2\textrm{CO}_3] + [\textrm{HCO}_3^-],  \; \; \;  (2)$$

where $\Delta\textrm{CO}_2 $ is the additional CO$_2$ injected into the system and $[\Sigma\textrm{CO}_2]_0$ is the initial total CO$_2$ (before injecting additional).

The third conservation relationship is for total exchangeable hydrogen ion, which is the sum of free plus bound hydrogen ions:

$$ \textrm{H} = [\textrm{H}^+] + [\textrm{H}_2\textrm{PO}_3^-] +  [\textrm{H}_2\textrm{CO}_3], \; \; \; (3)$$

Where the value of H must be obtained based on the amount of exchangeable hydrogen ion created stoichiometrically by the carbonic anhydrase reaction. 

The value of H in the third equation is obtained based on the initial concentrations plus the turnover of the carbonic anhydrase reaction:

$$ \textrm{H} = [\textrm{H}^+]_0 + [\textrm{H}_2\textrm{PO}_3^-]_0 +  [\textrm{H}_2\textrm{CO}_3]_0 + \Delta\textrm{H}_{CA} = [\textrm{H}^+] + [\textrm{H}_2\textrm{PO}_3^-] +  [\textrm{H}_2\textrm{CO}_3] , $$

where $\Delta\textrm{H}_{CA}$ is obtained from the initial and final bicarbonate concentrations:

$$ \Delta\textrm{H}_{CA} = [\textrm{H}\textrm{CO}_3^-] + [\textrm{H}_2\textrm{CO}_3] -  \left( [\textrm{H}\textrm{CO}_3^-]_0 + [\textrm{H}_2\textrm{CO}_3]_0 \right). $$

(In other words, what the above equation does is count how many times the CA reaction has turned over. There is a hydrogen ion generated every time that reaction turns over.)

Substitution yields

$$ [\textrm{H}^+]_0 + [\textrm{H}_2\textrm{PO}_3^-]_0 - [\textrm{H}\textrm{CO}_3^-]_0  = [\textrm{H}^+] + [\textrm{H}_2\textrm{PO}_3^-] -  [\textrm{H}\textrm{CO}_3^-] $$

for the final conservation relation. Expressing the three conservation relations in terms of the unknown final concentrations [H$^+$], [CO$_2$], and [HPO$_3^{2-}$], we get

$$ \textrm{P}_0 = [\textrm{HPO}_3^{2-}]\left( 1 + \frac{[\textrm{H}^+]}{K_1} \right) , $$

$$ [\Sigma\textrm{CO}_2] + \Delta\textrm{CO}_2 = [\textrm{CO}_2]\left(1 + K_{eq} + \frac{K_2 K_{eq}}{[\textrm{H}^+]} \right) , \; \textrm{and}$$

$$ [\textrm{H}^+]_0 + [\textrm{H}_2\textrm{PO}_3^-]_0 - [\textrm{H}\textrm{CO}_3^-]_0  = [\textrm{H}^+] + \frac{P_0}{\left(1 + \frac{[\textrm{H}^+]}{K_1} \right)} - [\textrm{CO}_2]\frac{K_2 K_{eq}}{[\textrm{H}^+]}.$$

Here we have three equations and three unknowns. In MATLAB we can solve systems like this using the function fsolve. The first step is to write a function that computes the three equations as a function of the three variables:

function f = conserve3(x,delta_CO2,Po,SCO2o,Ho,H2PO3o,HCO3o)

% Chemical constants.

K1 = 10^(-6.8);

K2 = 10^(-3.60);

Keq = 2.7e-3;

% The three unknowns

HPO3 = x(1);

CO2 = x(2);

H = x(3);

% Conservation equations (should equal 0)

f(1,:) = Po - HPO3*(1 + H/K1);

f(2,:) = SCO2o + delta_CO2 - CO2*(1 + Keq + K2*Keq/H);

f(3,:) = Ho + H2PO3o - HCO3o - H - Po*(H/K1)/(1+H/K1) + CO2*K2*Keq/H;

Notice that the output of the function, 'f', is a vector that is equal to 0 when Equation (4) is satisfied. In addition to the vector of unknowns, x, the function accepts as inputs the injected CO2, the total phosphate, initial total CO2, initial hydrogen ion activity, initial protonated phosphate, and initial protonated bicarb: 'delta_CO2', 'Po', 'SCO2o', 'Ho', 'H2PO3o', 'HCO3o'. 

To call the 'fsolve' function and compute the equilibrium concentration, we use the following syntax

delta_CO2 = +10e-3;

x = fsolve(@conserve3,x,[],delta_CO2,Po,SCO2o,Ho,H2PO3o,HCO3o);

pH = -log10(x(3))

This calculation yields a pH of approximately 6.6, illustrating that injection of CO$_2$ causes a drop in pH. (What happens to the predicted drop in pH if the total phosphate concentration is increased? You should find that increasing the phosphate concentration results in smaller changes in pH for a given amount of CO$_2$ injected.)

To make a plot of pH versus CO$_2$ injected we can use a simple loop

for i = 1:25

  delta_CO2(i) = (i-1)*1e-3;

  x = fsolve(@conserve3,x,[],delta_CO2(i),Po,SCO2o,Ho,H2PO3o,HCO3o);

  pH(i) = -log10(x(3));

end

plot(delta_CO2*1e3,pH); xlabel('\Delta CO_2 (mM)'); ylabel('pH');