ProteinLigandProg

Robert A. Alberty

Department of Chemistry
Massachusetts Institute of Technology
Cambridge, MA 02139

Introduction

Protein molecules interact with species in their aqueous environment.  They are always in equilibrium with hydrogen ions, but they may also bind other ligands at specific sites.  This package is concerned with proteins like hemoglobin that bind oxygen molecules at specific sites that are linked with acid dissociations in the protein.  The simple system considered here is
          (Subscript[pK, 1M]) (Subscript[pK, 2M])													(1)
    M = HM = (Subscript[H, 2])M
K  ||
    (Subscript[MO, 2])= (Subscript[HMO, 2])= (Subscript[H, 2])(Subscript[MO, 2])
           (Subscript[pK, 1  MO2])   pK\(TraditionalForm\`\_\(2  MO2\)\)
where K = [M] [(Subscript[O, 2])]/[ (Subscript[MO, 2])] and the hydrogen ions have been omitted in the acid disociations.  This is a hypothetical system, but the pKs and their enthalpies of reaction have been chosen to be the same as those determined by Antonini and coworkers [1] for human hemoglobin.
	There are three ways to make thermodynamic calculations on a system like this, depending on the thermodynamic variables that are specified.  When the temperature and pressure are specified, the criterion for spontaneous change and equilibrium is provided by the Gibbs energy G.  In this approach, which is the one used in chemical thermodynamics, the system is described in terms of species, like the eight species involved in the system considered here.  
	In biochemistry it is convenient to specify in addition the pH, and in this case the system is described in terms of reactants, which are sums of species.  When this is done the above system can be described in terms of reactants (pseudoisomer groups), and system (1) is described by a single biochemical reaction:
(Subscript[H, av])(Subscript[MO, 2]) = (Subscript[H, av])M+ (Subscript[O, 2])   												(2)
where the sub av indicates an average number.  The equilibrium constant (Subscript[K, O2])' for this biochemical reaction is referred to as an apparent equilibrium constant because it is a function of pH as well as temperature.
(Subscript[K, O2])'=[ (Subscript[H, av])M][(Subscript[O, 2])] /[(Subscript[H, av])(Subscript[MO, 2])]   											(3)
When the pH is specified, the Gibbs energy no longer provides the criterion for spontaneous change and equilibrium, and it is necessary to use a Legendre transform to define a transformed Gibbs energy G' that does.  
	In discussing a ligand like oxygen molecules, it is often convenient to specify the equilibrium concentration of molecular oxygen.  In this case the system can simply be described in terms of  (Subscript[H, av])(Subscript[MO, 2  av]).  When temperature, pressure, pH, and [(Subscript[O, 2])] are specified, it is necessary to use a Legendre transform to define a further transformed Gibbs energy G'' to provide the criterion for spontaneous change and equilibrium.  These three types of calculations, in the order in which we will discuss them they are:

I.  At specified T, P, and pH, the transformed Gibbs energy G' provides the criterion for equilibrium
II.  At specified T, P, pH, and [(Subscript[O, 2])], the further transformed Gibbs energy G'' provides the criterion for equilibrium
III.  At specified T and P, the Gibbs energy G provides the criterion for equilibrium

The third type of calculation is the most basic, but in the absence of experimental methods for detecting species directly, the approaches in I and II are used. However, approaches I and II can yield the standard thermodynamic properties of species, as described in this package.   
	The Gibbs energies G', G'', and G are the most important thermodynamic properties of these systems because if they can be determined as functions of their natural variables (the variables indicated in I, II, and III), all the other thermodynamic properties can be obtained by taking partial derivatives.  
	This Introduction cannot provide more information on the ideas and equations behind the calculations presented here, but all the background is provided in ref. 2.  More specifically, the thermodynamics of protein-ligand binding is discussed in references 3-5.  This last article emphasizes the importance of Maxwell relations in providing relationships between thermodynamic properties.  Maxwell relations (second cross partial derivatives) apply to exact differentials.

Maxwell relations are important because they show the relationships between different properties.  Since properties like G, G', and G'' cannot be measured directly, it is important that some Maxwell relations yield experimental properties that are related to properties like G''.  A prime example of this is the average number of oxygen molecules bound \(N\&_\_O2\), which can be determined spectrophotometrically.  Another example is (Subscript[Δ, r]) (Subscript[N, H]), which can be obtained as the difference between the acid titration curves of the protein in the absence of oxygen and saturated with oxygen.  If \(N\&_\_O2\) can be determined at a series of pHs and at least two temperatures, all the other thermodynamic properties can be calculated, including the standard properties of the species involved.

Standard further transformed thermodynamic properties like (Subscript[Δ, f]) G'' ° are very complicated functions of T, pH, and [O\(\_2\)]. Taking partial derivatives yields even more complex functions.  These could never be done by hand, but Mathematica makes it possible to make them.
	Data and programs for calculations on the thermodynamics of metabolism  are freely available in MathSource [6] if you have Mathematica in your personal computer.  More Mathematica programs for the thermodynamics of biochemical reactions are available in the Mathematica Journal by Akers and Goldberg [7].. 

References

(1)  E. Antonini, J. Wymna, M. Brunori, C. Fronticelli, E. Bucci, and A. Rossi-Fanelli, J. Biol. Chem. 240, 1090-1103 (1965).
(2)  R. A. Alberty, Thermodynamics of Biochemical Reactions, Wiley, Hoboken, NJ (2003).
(3)  R. A. Alberty, Biophys. Chem. 62, 141-159 (1996).
(4)  R. A. Alberty, Biophys. Chem. 63, 119-132 (1997).
(5)  R. A. Alberty, Biophys. Chem., in press.
(6)  R.  A. Alberty, BasicBiochemData2, (2002) http://www.mathsource.com/cgi/msitem?0211-622.  
(7)  D. L. Akers and R. N. Goldberg, Math. J. 8, 86-113 (2001). http://www.mathematica-journal.com/issue/v8i1/

Initialization

BeginPackage["ProteinLigandProg`"];

Off[General::spell, General::spelll];

calcappeqconst::usage = "calcappeqconst derives the function of \
temperature and pH that yields an apparent equilibrium constant.";
\(\(reactionprops293::usage = \*""reactionprops293 yields the apparent equilibrium (constantx10^5), the base 10 logarithm of the apparent equilibrium constant, the standard transformed reaction Gibbs energy, the change in binding of hydrogen ions, the standard transformed enthalpy of reaction, and the standard transformed entropy of reaction at 293.25 K and five pHs."";\)\)
\(\(reactionprops313::usage = \*""reactionprops313 yields the apparent equilibrium (constantx10^5), the base 10 logarithm of the apparent equilibrium constant, the standard transformed reaction Gibbs energy, the change in binding of hydrogen ions, the standard transformed enthalpy of reaction, and the standard transformed entropy of reaction at 313.25 K and five pHs."";\)\)
\(\(calcstdfrtrGe::usage = \*""calcstdfrtrGe  derives the function of temperature, pH, and concentration of molecular oxygen that yields the standard further transformed Gibbs energy of formation of (Subscript[H, av])(Subscript[MO, 2  av]). Energies are in J (mol^(-1))."";\)\)
calcstdfrtrGe293::usage = "calcstdfrtrGe293 calculates the standard \
transformed Gibbs energy of formation at 293.15 K and a series of oxygen \
concentrations and pHs.";
calcstdfrtrGe313::usage = "calcstdfrtrGe313 calculates the standard \
transformed Gibbs energy of formation at 313.15 K and a series of oxygen \
concentrations and pHs.";
nO2::usage = "nO2 calculates the number of oxygen molecules bound \
as a function of T, pH, and concentration of molecular oxygen.";
taboxycols293::usage = "taboxycols293 calculates the number of \
oxygen molecules bound at 293.15 K and a series of oxygen concentrations and \
pHs.";
 taboxycols313::usage = "taboxycols313 calculates the number of \
oxygen molecules bound at 313.15 K and a series of oxygen concentrations and \
pHs.";
frtrGbyint::usage = "frtrGbyint calculates the further transformed \
Gibbs energy plus an integration constant.";
\(\(calcnegRTlnPM::usage = \*"" calcnegRTlnPM derives the function of temperature, pH, and concentration of molecular oxygen that yields -(Subscript[RTlnP, M]). Energies are in J (mol^(-1))."";\)\)
nH::usage = "nH calculates the number of hydrogen atoms bound as a \
function of T, pH, and concentration of molecular oxygen.";
nH293::usage = "nH293 calculates the number of hydrogen atoms bound \
at 293.15 K as a function of pH and concentration of molecular \
oxygen.";
nH313::usage = "nH293 calculates the number of hydrogen atoms bound \
at 313.15 K as a function of pH and concentration of molecular \
oxygen.";
frtrenthalpy::usage = "frtrenthalpy gives the function of T, pH, \
and concentrations of molecular oxygen that gives the further transformed \
enthalpy in J/mol.";
frtrenthalpy293::usage = "frtrenthalpy293 calculates the further \
transformed enthalpy at 293.15 K and a series of oxygen concentrations and \
pHs.";
frtrenthalpy313::usage = "frtrenthalpy313 calculates the further \
transformed enthalpy at 313.15 K and a series of oxygen concentrations and \
pHs.";
frtrentropy::usage = "frtrentropy gives the function of T, pH, and \
concentrations of molecular oxygen that gives the further transformed entropy \
in J/K mol.";
frtrentropy293::usage = "frtrentropy293 calculates the further \
transformed entropy at 293.15 K and a series of oxygen concentrations and \
pHs.";
frtrentropy313::usage = "frtrentropy313 calculates the further \
transformed entropy at 313.15 K and a series of oxygen concentrations and \
pHs.";
\(\(calcchemrxprops::usage = \*""calcchemrxprops calculates the functions for calculating K and four pKs and the reaction Gibbs energy, stndard reaction enthalpy, and standard reaction entropy at a specified temperature.  The energies are in J (mol^(-1)), and the entropies are in J (K^(-1))(mol^(-1))."";\)\)
t::usage = "t is the absolute temperature in Kelvins.";
pH::usage = "pH is an independent variable.";
co2::usage = "co2 is the specified concentration of molecular \
oxygen.";
appk::usage = "apparent equilibrium constant*10^5.";
logappk::usage = "base 10 logarithm of the apparent equilibrium \
constant.";
\(\(trge::usage = "standard transformed Gibbs energy of reaction.";\)\)
chgNH::usage = "change in the binding of hydrogen ions in a \
biochemical reaction.";
trenthal::usage = "standard transformed enthalpy of a biochemical \
reaction.";
trentro::usage = "standard transformed entropy of a biochemical \
reaction.";

Begin["ProteinLigandProg`Private`"]

I.  At specified T, P, and pH, the transformed Gibbs energy G' provides the criterion for equilibrium

Ia.  Calculate the apparent equilibrium constant for the reaction (Subscript[H, av])(Subscript[MO, 2]) = (Subscript[H, av])M+ (Subscript[O, 2])

The apparent equilibrium constant for reaction (2) is given by
(Subscript[K, O2])' = (Subscript[KP, M])/(Subscript[P, MO2])													(4)
where K = [M] [(Subscript[O, 2])]/[(Subscript[MO, 2])].  (Subscript[P, M])is the binding polynomial for (Subscript[H, av])M that is given by
(Subscript[P, M])= 1 + (10^(-pH + pK1M)) + (10^(-2 pH + pK1M + pK2M))									(5)
The binding polynomial for (Subscript[H, av])M(Subscript[O, 2]) is given by 
(Subscript[P, MO2])= 1 + (10^(-pH + pK1MO2)) + (10^(-2 pH + pK1MO2 + pK2MO2))								(6)
The following program derives the function of T and pH that gives (Subscript[K, O2])'.  The values of the five properties in system (1) are given as functions of temperature in the program calcappeqconst.

calcappeqconst := 
      Module[{pK1M, pK2M, pK1MO2, pK2MO2, k, pM, 
          pMO2},   
pK1M = 
          7.85 - (37.7*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
pK2M = 
          5.46 - (-6.3*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
        pK1MO2 = 
          6.67 - (37.7*10^3/(8.3145*Log[10]))*(1/293.15 - 1/t); 
        pK2MO2 = 
          6.04 - (-6.3*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
k = 
          10^-5*
            Exp[((60.7*10^3)/8.3145)*(1/293.15 - 
                    1/t)]; 
pM = (1 + 
              10^(-pH + pK1M) + 
              10^(-2*pH + pK1M + 
                    pK2M)); 
pMO2 = (1 + 
              10^(-pH + pK1MO2) + 
              10^(-2*pH + pK1MO2 + pK2MO2)); 
k*
          pM/pMO2]

calcappeqconst

calcappeqconst /. t -> 293.15

Plot[Evaluate[{calcappeqconst /. 
                t -> {283.15, 293.15, 303.15, 313.15}}],{pH, 5, 9},AxesLabel->{"pH","(Subscript[(K), O2])'"}];

The highest curve is for the highest temperatire.

Since (Subscript[Δ, r]) H° = 60.7 kJ (mol^(-1)) for K, raising the temperature makes K bigger.  If K is bigger, at a given pH, the binding is smaller.

Ib.  Construct tables that give all the thermodyamic properties of the reaction (Subscript[H, av])(Subscript[MO, 2]) = (Subscript[H, av])M+ (Subscript[O, 2]) at five pHs and two temperatures

Make a table of (Subscript[K, O2])',(Subscript[logK, O2])',(Subscript[Δ, r]) G' °,(Subscript[Δ, r]) (Subscript[N, H]),(Subscript[Δ, r]) H' °,and (Subscript[Δ, r]) S' ° at pHs 5, 6, 7, 8, and 9 at 293.15 K and 313.15 K.

appk = calcappeqconst*10^5;

logappk = Log[10, calcappeqconst];

trge = -8.3145*t*Log[calcappeqconst];

chgNH = (-1/Log[10])*D[Log[calcappeqconst], pH];

trenthal = -t^2*
          D[(-8.3145*Log[calcappeqconst]), t];

trentro = -D[-8.3145*t*Log[calcappeqconst], t];

reactionprops293 = {appk, logappk, trge, chgNH, trenthal, trentro} /. 
            t -> 293.15 /. pH -> {5, 6, 7, 8, 9};

reactionprops313 = {appk, logappk, trge, chgNH, trenthal, trentro} /. 
            t -> 313.15 /. pH -> {5, 6, 7, 8, 9};

Table 1  Thermodynamic Properties of the Biochemical Reaction (Subscript[H, av])(Subscript[MO, 2]) = (Subscript[H, av])M+ (Subscript[O, 2]) at 293.15 K.  The energies are in J (mol^(-1)) and the entropies are in J (K^(-1))(mol^(-1)).

\(TableForm[reactionprops293, 
      TableHeadings -> {{\*""(Subscript[K, O2])'(x10^5)"", \*""(Subscript[logK, O2])'"", \*""(Subscript[Δ, r])G'°"", \*""(Subscript[Δ, r])(Subscript[N, H])"", \*""(Subscript[Δ, r])H'°"", \*""(Subscript[Δ, r])S'°""}, {"pH 5", "pH 6", "pH 7", "pH 8", "pH 9"}}]\)

Table 2  Thermodynamic Properties of the Biochemical Reaction (Subscript[H, av])(Subscript[MO, 2]) = (Subscript[H, av])M+ (Subscript[O, 2]) at 313.15 K.  The energies are in J (mol^(-1)) and the entropies are in J (K^(-1))(mol^(-1)).

\(TableForm[reactionprops313, 
      TableHeadings -> {{\*""(Subscript[K, O2])'(x10^5)"", \*""(Subscript[logK, O2])'"", \*""(Subscript[Δ, r])G'°"", \*""(Subscript[Δ, r])(Subscript[N, H])"", \*""(Subscript[Δ, r])H'°"", \*""(Subscript[Δ, r])S'°""}, {"pH 5", "pH 6", "pH 7", "pH 8", "pH 9"}}]\)

These calculations are based on the assumption that the standard enthalpies and entropies of species are independent of temperature.  The standard transformed enthalpies and entropies do depend on the temperature because (Subscript[H, av])(Subscript[MO, 2]) and (Subscript[H, av])M are pseudoisomer groups.  Reference 5 shows that for such a system, there are three Maxwell relations and two equations for the transformed enthalpy, which are derived from Maxwell relations, that show how these four thermodynamic porperties are related to each other.  The change of (Subscript[Δ, r]) (Subscript[N, H]) with temperature determines the changes of (Subscript[Δ, r]) H' ° and (Subscript[Δ, r]) S' ° with temperature..   
	The values of (Subscript[K, O2])' can be used to calculate the average number of oxygen molecules bound at a specified concentration of molecular oxygen \(N\&_\_O2\).
\(N\&_\_O2\)=1/(1+(Subscript[K, O2])'/[O\(\_2\)])=( [(Subscript[O, 2])]/(Subscript[K, O2])')/(1+[O\(\_2\)]/(Subscript[K, O2])')							(7)

II.  At specified T, P, pH, and [(Subscript[O, 2])], the further transformed Gibbs energy G'' provides the criterion for equilibrium

IIa.  Calculate the standard further transformed Gibbs energy of formation of (Subscript[H, av])(Subscript[MO, 2  av])

When these independent variables are specified, the system can be represented by (Subscript[H, av])(Subscript[MO, 2  av]), where av indicates an average number.  The standard further transformed Gibbs energy of formation of (Subscript[H, av])(Subscript[MO, 2  av]) is given by
(Subscript[Δ, f]) G'' ° = -RTln(Subscript[P, M])- RTln(1 + [O\(\_2\)]/(Subscript[K, O2])')									(8)
The following program derives the function for (Subscript[Δ, f]) G'' °  of (Subscript[H, av])(Subscript[MO, 2  av]):

calcstdfrtrGe := 
      Module[{pK1M, pK2M, pK1MO2, pK2MO2, k, pM, pMO2, 
          kprimeO2},   
pK1M = 
          7.85 - (37.7*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
pK2M = 
          5.46 - (-6.3*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
        pK1MO2 = 
          6.67 - (37.7*10^3/(8.3145*Log[10]))*(1/293.15 - 1/t); 
        pK2MO2 = 
          6.04 - (-6.3*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
k = 
          10^-5*
            Exp[(60.7*10^3/8.3145)*(1/293.15 - 
                    1/t)]; 
pM = (1 + 
              10^(-pH + pK1M) + 
              10^(-2*pH + pK1M + 
                    pK2M)); 
pMO2 = (1 + 
              10^(-pH + pK1MO2) + 
              10^(-2*pH + pK1MO2 + 
                    pK2MO2)); 
kprimeO2 = 
          k*pM/pMO2; 
-8.3145*t*Log[pM] - 
          8.3145*t*Log[1 + co2/kprimeO2]]

calcstdfrtrGe

calcstdfrtrGe /. t -> 293.15

calcstdfrtrGe /. t -> 293.15 /. pH -> 7

Plot3D[Evaluate[calcstdfrtrGe /. t -> 293.15],{co2, 10^-10,  .0001},{pH, 5, 9},AxesLabel->{"[(Subscript[(O), 2])]","pH",""},\(PlotLabel -> \*""(Subscript[Δ, f])G''°""\)];

The slopes in the oxygen direction provide information on the binding of molecular oxygen, and the slopes in the pH direction provide information on the binding of hydrogen ions.

Table 3  Standard Further Transformed Gibbs Energies of Formation of (Subscript[H, av])(Subscript[MO, 2  av])at 293.15 K as a Function of pH and Concentration of Molecular Oxygen.  The energies are in J (mol^(-1))and the concentrations of molecular oxygen are molar.

calcstdfrtrGe293 = calcstdfrtrGe /. t -> 293.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[calcstdfrtrGe293], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

Table 4  Standard Further Transformed Gibbs Energies of Formation of (Subscript[H, av])(Subscript[MO, 2  av])at 313.15 K as a Function of pH and Concentration of Molecular Oxygen.  The energies are in J (mol^(-1))and the concentrations of molecular oxygen are molar.

calcstdfrtrGe313 = calcstdfrtrGe /. t -> 313.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[calcstdfrtrGe313], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

IIb.  Calculate the average number of oxygen molecules bound in (Subscript[H, av])(Subscript[MO, 2  av])

A Maxwell relation shows that the average number of oxygen molecules bound in (Subscript[H, av])(Subscript[MO, 2  av]) is given by -[O\(\_2\)]/RT times the derivative of (Subscript[Δ, f]) G'' ° with respect to [O\(\_2\)].

nO2 = -(co2/(8.3145*t))*D[calcstdfrtrGe, co2]

nO2 /. t -> 293.15

nO2 /. t -> 293.15 /. pH -> 7

Plot3D[Evaluate[nO2 /. t -> 293.15],{co2, 10^-10,  .0001},{pH, 5, 9},AxesLabel->{"[(Subscript[(O), 2])]","pH",""},\(PlotLabel -> \*""\!\(N\&_\_O2\)""\)];

Table 5  Average Number of Oxygen Molecules Bound in (Subscript[H, av])(Subscript[MO, 2  av])at 293.15 K as a Function of pH and Concentration of Molecular Oxygen.  The concentrations of molecular oxygen are molar.

taboxycols293 = nO2 /. t -> 293.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[taboxycols293], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

If the average number of oxygen molecules bound by the protein is measured at a specified pH and a series of concentrations of molecular oxygen, the value of the apparent equilibrium constant at that pH and temperature can be calculated by plotting 1/\(TraditionalForm\`N\&_\_O2\) versus 1/[(Subscript[O, 2])].  For example, at 293.15 K and pH 5 the the average numbers of oxygen molecules bound at the six oxygen concentrations are given by

Transpose[taboxycols293][[1]]

The reciprocal binding numbers are given by

recipnos = 1/Transpose[taboxycols293][[1]]

The reciprocal concentrations are given by

recipconc = 
      1/{ .0000001,  .000001,  .00001,  .0001,  .001,  .01}

Fit[Transpose[{recipconc, recipnos}], {1, x}, x]

The slope of the double reciprocal plot is equal to the apparent equilibrium constant under these conditions.

calcappeqconst /. t -> 293.15 /. pH -> 5

Table 6  Average Number of Oxygen Molecules Bound in (Subscript[H, av])(Subscript[MO, 2  av])at 313.15 K as a Function of pH and Concentration of Molecular Oxygen.  The concentrations of molecular oxygen are molar.

taboxycols313 = nO2 /. t -> 313.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[taboxycols313], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

calcappeqconst /. t -> 293.15

calcstdfrtrGe /. t -> 293.15

nO2 /. t -> 293.15

The average number of molecules of oxygen bound is given by
\(TraditionalForm\`N\&_\_O2\)=([(Subscript[O, 2])]/K')/(1+[(Subscript[O, 2])]/K')												(9)
A Maxwell relation shows that \(TraditionalForm\`N\&_\_O2\)/[(Subscript[O, 2])] can be integrated with respect to [(Subscript[O, 2])]  to obtain G'' + integration constant

Define x = [(Subscript[O, 2])]/K'.

frtrGbyint = -8.3145*t*Integrate[1/(1 + x), x] /. 
        x -> co2/calcappeqconst

frtrGbyint /. t -> 293.15

Since integration yields the second term in equation 8, and so the first term in (Subscript[P, M]) has to be added.

calcnegRTlnPM := 
      Module[{pK1M, pK2M, 
          pM},   
pK1M = 
          7.85 - (37.7*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
pK2M = 
          5.46 - (-6.3*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
pM = (1 + 
              10^(-pH + pK1M) + 
              10^(-2*pH + pK1M + 
                    pK2M)); 
-8.3145*t*
          Log[pM]]

calcnegRTlnPM

(calcnegRTlnPM + frtrGbyint) /. t -> {293.15, 313.15} /. 
        pH -> {12, 13, 14} /. co2 -> 0

calcstdfrtrGe /. t -> {293.15, 313.15} /. 
        pH -> {12, 13, 14} /. co2 -> 0

This shows that experimental data on \(TraditionalForm\`N\&_\_O2\) can be integrated to obtain the oxygen-dependent term in the standard further transformed Gibbs energy.  Measurements of \(TraditionalForm\`N\&_\_O2\) as a function of pH and [(Subscript[O, 2])] at at least two temperatures can be used to obtain all the thermodynamic properties of the system.

IIc.  Calculate the average number of hydrogen ions bound in (Subscript[H, av])(Subscript[MO, 2  av])

According to a Maxwell equation, the average number of hydrogen ions bound in (Subscript[H, av])(Subscript[MO, 2  av]) is given by (1/RTln(10)) times the partial derivative of G'' with respect to pH.

nH = (1/(8.3145*t*Log[10]))*D[calcstdfrtrGe, pH]

nH /. t -> 293.15

nH /. t -> 293.15 /. pH -> 7

Simplify[nH /. t -> 293.15 /. pH -> 7]

Plot3D[Evaluate[nH /. t -> 293.15],{co2, 10^-10,  .0001},{pH, 5, 9},AxesLabel->{"[(Subscript[(O), 2])]","pH",""},\(PlotLabel -> \*""\!\(N\&_\_H\)""\)];

Plot the titration curve at 293.15 K in the absence of molecular oxygen and at co2 = 0.01 M.

\(\(Plot[{\(nH /. t -> 293.15\) /. 
            co2 -> 1, \(nH /. t -> 293.15\) /. co2 -> 0}, {pH, 
          5, 9}, AxesLabel -> {"pH", \*""\!\(N\_H\)""}];\)\)

Table 7  Average Number of Hydrogen Ions Bound in (Subscript[H, av])(Subscript[MO, 2  av])at 293.15 K as a Function of pH and Concentration of Molecular Oxygen.  The concentrations of molecular oxygen are molar.

nH293 = nH /. t -> 293.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[nH293], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

Table 8  Average Number of Hydrogen Ions Bound in (Subscript[H, av])(Subscript[MO, 2  av])at 313.15 K as a Function of pH and Concentration of Molecular Oxygen.  The concentrations of molecular oxygen are molar.

nH313 = nH /. t -> 313.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[nH313], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

IId.  Calculate the standard further transformed enthalpy of formation of (Subscript[H, av])(Subscript[MO, 2  av])

The standard further transformed enthalpy of formation of is calculated using the Gibbs-Helmholtz equation.

frtrenthalpy = -(t^2)*D[(calcstdfrtrGe/t), t]

frtrenthalpy /. t -> 293.15 /. pH -> 7

Plot3D[Evaluate[frtrenthalpy /. t -> 293.15],{co2, 10^-10,  .0001},{pH, 5, 9},AxesLabel->{"[(Subscript[(O), 2])]","pH",""},\(PlotLabel -> \*""(Subscript[Δ, f])H''°""\)];

Note that the standard further transformed enthalpy of formation is with respect to (Subscript[Δ, f]) H°(M) = 0 at each temperature.

Table 9  Standard Further Transformed Enthalpy of Formation of (Subscript[H, av])(Subscript[MO, 2  av])at 293.15 K as a Function of pH and Concentration of Molecular Oxygen.  The enthalpies are in kJ (mol^(-1)).

frtrenthalpy293 = frtrenthalpy /. t -> 293.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[frtrenthalpy293], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

Table 10  Standard Further Transformed Enthalpy of Formation of (Subscript[H, av])(Subscript[MO, 2  av])at 313.15 K as a Function of pH and Concentration of Molecular Oxygen.  The enthalpies are in kJ (mol^(-1)).

frtrenthalpy313 = frtrenthalpy /. t -> 313.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[frtrenthalpy313], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

Plot3D[Evaluate[frtrenthalpy /. t -> 313.15],{co2, 10^-10,  .0001},{pH, 5, 9},AxesLabel->{"[(Subscript[(O), 2])]","pH",""},\(PlotLabel -> \*""(Subscript[Δ, f])H''°""\)];

IIe.  Calculate the standard further transformed entropy of formation of (Subscript[H, av])(Subscript[MO, 2  av])

The standard further transformed entropy of formation of (Subscript[H, av])(Subscript[MO, 2  av]) is calculated by taking the derivative of (Subscript[Δ, f]) G'' ° with respect to T and multiplying by minus 1.

frtrentropy = -D[calcstdfrtrGe, t]

frtrentropy /. t -> 293.15 /. pH -> 7

Plot3D[Evaluate[frtrentropy /. t -> 293.15],{co2, 0,  .0001},{pH, 5, 9},AxesLabel->{"[(Subscript[(O), 2])]","pH",""},\(PlotLabel -> \*""(Subscript[Δ, f])S''°""\)];

Note that these standard further transformed entropies are with reapect to (Subscript[Δ, f]) S°(M) = 0 at each temperature.

Table 11  Standard Further Transformed Entropy of Formation of (Subscript[H, av])(Subscript[MO, 2  av])at 293.15 K as a Function of pH and Concentration of Molecular Oxygen.  The entropies are in J (K^(-1))(mol^(-1)).

frtrentropy293 = frtrentropy /. t -> 293.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[frtrentropy293], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

Table 12  Standard Further Transformed Entropy of Formation of (Subscript[H, av])(Subscript[MO, 2  av])at 313.15 K as a Function of pH and Concentration of Molecular Oxygen.  The entropies are in J (K^(-1))(mol^(-1)).

frtrentropy313 = frtrentropy /. t -> 313.15 /. 
            co2 -> { .0000001,  .000001,  .00001,  .0001,  .001,  .01} /. pH -> {5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9};

\(TableForm[Transpose[frtrentropy313], 
      TableHeadings -> {{5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 
            9}, {".0000001", ".000001", ".00001", ".0001", ".001", ".01"}}]\)

(frtrenthalpy - t*frtrentropy) /. t -> 293.15 /. 
        co2 -> {0,  .0001} /. pH -> {5, 7, 9}

calcstdfrtrGe /. t -> 293.15 /. co2 -> {0,  .0001} /. 
      pH -> {5, 7, 9}

These agree as expected.

III.  At specified T and P, the Gibbs energy G provides the criterion for equilibrium

IIIa  Calculate the standard reaction properties for M(Subscript[O, 2]) = M + O\(TraditionalForm\`\_2\)

The following program can be used to calculate the thermodynamic properties of the five reactions involved in system (1) at any desired temperature in the range approximately 273.25 K to  313.15 K.  Since these reactions are written in terms of species, the program is used here to calculate properties at 298.15 K since that is the usual temperature for chemical thermodynamic tables.

calcchemrxprops := 
      Module[{k, pK1M, pK2M, pK1MO2, pK2MO2, line1rx, line2rx, line3rx, 
          line4rx, 
          line5rx},  (*This program calculates the functions for calculating K and four pKs and the reaction Gibbs energy,  
          stndard reaction enthalpy,  
          and standard reaction entropy at a specified temperature .   The energies are in joules per mole,  
          and the entropies are in joules per K per mole.*) 
k = 
          10^-5*
            Exp[((60.7*10^3)/8.3145)*(1/293.15 - 
                    1/t)]; 
pK1M = 
          7.85 - (37.7*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
pK2M = 
          5.46 - (-6.3*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
pK1MO2 = 
          6.67 - (37.7*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
pK2MO2 = 
          6.04 - (-6.3*10^3/(8.3145*Log[10]))*(1/293.15 - 
                  1/t); 
line1rx = {k, -8.3145* .29815*Log[k], -t^2*
              D[-8.3145*Log[k], t]/1000, -D[-8.3145*t*Log[k], 
                t]}; 
line2rx = {pK1M, 
            8.3145* .29815*Log[10]*pK1M, -t^2*
              D[8.3145*Log[10]*pK1M, t]/1000, -D[8.3145*t*Log[10]*pK1M, 
                t]}; 
line3rx = {pK2M, 
            8.3145* .29815*Log[10]*pK2M, -t^2*
              D[8.3145*Log[10]*pK2M, t]/1000, -D[8.3145*t*Log[10]*pK2M, 
                t]}; 
line4rx = {pK1MO2, 
            8.3145* .29815*Log[10]*pK1MO2, -t^2*
              D[8.3145*Log[10]*pK1MO2, t]/1000, -D[8.3145*t*Log[10]*pK1MO2, 
                t]}; 
line5rx = {pK2MO2, 
            8.3145* .29815*Log[10]*pK2MO2, -t^2*
              D[8.3145*Log[10]*pK2MO2, t]/1000, -D[8.3145*t*Log[10]*pK2MO2, 
                t]}; 
{line1rx, line2rx, line3rx, 
          line4rx, line5rx}]

Table 12 Thermodynamic Properties of Chemical Reactions at 298.15 K

\(TableForm[calcchemrxprops /. t -> 298.15, 
      TableHeadings -> {{"K(MO2=M+O2)", "pK1M(HM=H+M)", "pK2M(H2M=H+HM)", "pK1MO2(HMO2=H+MO2)", "pK2MO2(H2MO2=H+HMO2)"}, {"K or pK", \*""(Subscript[Δ, r])G°"", \*""(Subscript[Δ, r])H°"", \*""(Subscript[Δ, r])S°""}}]\)

IIIb  Calculate standard thermodynamic properties of species

Eight species are involved in these five chemical reactions, and the standard thermodynamic propertie of H\(TraditionalForm\`\(\^+\)\) and O\(TraditionalForm\`\_2\),aq are available in the NBS Tables.  Since the standard thermodynamic properties of the remaining six species cannot be calculated from the properties of the five reactions, the properties of M are assigned to be zero at each temperature.  This is the same as the convention of the usual chemical thermodynamic properties for (SuperPlus[H]).  The entropy of formation of (Subscript[O, 2])(aq) at 298.15 K, which is not given in the NBS Tables is

(-11.7 - 16.4)/ .29815

in J (K^(-1))(mol^(-1)).

The other properties can be calculated from Table 12 by using LinearSolve.  When the species are arranged in the order HM, H2M,MO2,HMO2,H2MO2, the stoichiometric number matrix is

stoichmat = {{0, 0, -1, 0, 0}, {-1, 0, 0, 0, 0}, {1, -1, 0, 
          0, 0}, {0, 0, 1, -1, 0}, {0, 0, 0, 1, -1}}

TableForm[stoichmat]

When the properties are arranged in the order G, H, S, the thermodynamic properties of the reactions are given by.

thpropsrxs = {{11.11, 44.17, 31.27, 37.43, 34.58}, {72.4, 
            37.7, -6.3, 
            37.7, -6.3}, {205.59, -21.68, -126.02, 
            0.91, -137.13}};

TableForm[Transpose[thpropsrxs]]

The values of the species properties are given by

TableForm[LinearSolve[stoichmat, Transpose[thpropsrxs]]]

Table 13  Standard Thermodynamic Properties of Species at 298.15 K

TableForm[{{"","(Subscript[Δ, r])G°/\!\(TraditionalForm\`\*
StyleBox[\(\(J\)\(\\ \)\),
FontWeight->"Plain",
FontSlant->"Plain",
FontTracking->"Plain",
FontVariations->{"Underline"->False,
"Outline"->False,
"Shadow"->False,
"StrikeThrough"->False,
"Masked"->False,
"CompatibilityType"->0,
"RotationAngle"->0}]\)(mol^(-1))","(Subscript[Δ, r])H°/\!\(TraditionalForm\`\*
StyleBox[\(\(J\)\(\\ \)\),
FontWeight->"Plain",
FontSlant->"Plain",
FontTracking->"Plain",
FontVariations->{"Underline"->False,
"Outline"->False,
"Shadow"->False,
"StrikeThrough"->False,
"Masked"->False,
"CompatibilityType"->0,
"RotationAngle"->0}]\)(mol^(-1))","(Subscript[Δ, r])S°/J (K^(-1))(mol^(-1))"},{M, 0, 0, 0},{HM, -44.17, -37.70, 21.68},{(Subscript[H, 2]) M, 75.44, -31.40, 147.70},{(Subscript[MO, 2]), -11.11, -72.40, -205.59},{(Subscript[HMO, 2]), -48.54, -110.1, -206.50},{(Subscript[H, 2]) (Subscript[MO, 2]), -83.12, -103.80, -169.37},\({\(H\^+\), 0, 0, 0}\),{(Subscript[O, 2  aq]), 16.40, -11.70, -94.25}}]

The values of (Subscript[Δ, f]) G° and (Subscript[Δ, f]) H° of the species in Table 13 are the most efficient way to store information on the system discussed here.  The values of (Subscript[Δ, f]) S° are of interest, but are redundant.  The values of (Subscript[Δ, f]) G° and (Subscript[Δ, f]) H° can be used to calculate all the tables in this package and related tables.  It is important to understand that when other ligands (for example CO) are bound at the same site, the entries for M, HM, and (Subscript[H, 2])M remain applicable.  These calculations have all been based on the assumption that (Subscript[Δ, f]) H° values of species are independent of temperature in the range approximately 273.15 to 313.15 K.  The values of (Subscript[Δ, f]) G° at temperatures other than 298.15 K can be calculated using
 (Subscript[Δ, f]) G°(T)=(T/298.15)(Subscript[Δ, f]) G°(298.15) + (1 - T/298.15)(Subscript[Δ, f]) H°(298.15)						(10)

End[];

EndPackage[];