In [2]:

    

from scipy.misc import logsumexp

dG0_H2O = -237.2 # kJ/mol
R = 8.31e-3 # kJ/(K*mol)
T = 298.15 # K

# Debye-Huckel
# Approximation of the temperature dependency of ionic strength effects
DH_alpha = 1e-3*(9.20483*T) - 1e-5*(1.284668 * T**2) + 1e-8*(4.95199 * T**3)
DH_beta = 1.6
debye_huckel = lambda I : DH_alpha * I**(0.5) / (1.0 + DH_beta * I**(0.5))


Alberty = {'CO2(g)':    {'nH':0, 'z':0,  'dG0_f': -394.36},
           'CO2(sp)':   {'nH':0, 'z':0,  'dG0_f': -385.97},
           'H2O(l)':    {'nH':2, 'z':0,  'dG0_f': -237.13},
           'H2CO3(sp)': {'nH':2, 'z':0,  'dG0_f': -606.33},
           'HCO3-(sp)': {'nH':1, 'z':-1, 'dG0_f': -586.77},
           'CO3-2(sp)': {'nH':0, 'z':-2, 'dG0_f': -527.81}}

KEGG2sp = {'C00288': ['H2CO3(sp)', 'HCO3-(sp)', 'CO3-2(sp)'],
           'C00001': ['H2O(l)'],
           'C00011': ['CO2(sp)']}

def transform(species, pH, I):
    DH = debye_huckel(I)
    dG0_primes = []
    for d in map(Alberty.get, species):
        dG0_prime = d['dG0_f'] + \
                    d['nH'] * (R*T*log(10)*pH + DH) - \
                    d['z']**2 * DH
        dG0_primes.append(dG0_prime)
    
    dG0_prime_tot = -R*T * logsumexp([-x/(R*T) for x in dG0_primes])
    
    P = map(lambda x: exp((dG0_prime_tot - x)/(R*T)), dG0_primes)
    return dG0_prime_tot, dict(zip(species, P))

def calc_dG0_prime(sparse, pH, I):
    dG0_prime = 0
    for kegg_id, coeff in sparse.iteritems():
        dG0_prime += coeff * transform(KEGG2sp[kegg_id], pH, I)[0]
    return dG0_prime


    

In [15]:

    

# calculate bicarbonate pKas
pKa1 = (Alberty['HCO3-(sp)']['dG0_f'] - Alberty['H2CO3(sp)']['dG0_f']) / (R*T*log(10))
pKa2 = (Alberty['CO3-2(sp)']['dG0_f'] - Alberty['HCO3-(sp)']['dG0_f']) / (R*T*log(10))
print pKa1, pKa2

# calculate bicarbonate transformed energy at pH 7, I=0.25M

print transform(KEGG2sp['C00288'], 7, 0.25)


    

    




3.42860397482 10.334892145
(-546.83995564466272, {'H2CO3(sp)': 0.00013929830105068118, 'CO3-2(sp)': 0.0017061512883912138, 'HCO3-(sp)': 0.9981545504105569})

In [53]:

    

data = []
pH_range = linspace(2, 14, 50)
for pH in pH_range:
    dG0_prime, P_dict = transform(KEGG2sp['C00288'], pH, 0.2)
    data.append(P_dict.values())
    
data = matrix(data)
fig = figure(figsize=(6,6))
plot(pH_range, data)
legend(P_dict.keys(), loc='best', fontsize=14)
xlabel('pH', fontsize=14)
ylabel('fractional abundance', fontsize=14)


    

    
Out[53]:





<matplotlib.text.Text at 0x7fef3c4652d0>






    

In [54]:

    

def draw_pH_response(I, figure):
    sparse = {'C00001': -1, 'C00011': -1, 'C00288': 1}
    data = []
    pH_range = linspace(2, 14, 50)
    for pH in pH_range:
        data.append(calc_dG0_prime(sparse, pH, I))
        
    
    plot(pH_range, data, label='I = %.1fM' % I, figure=figure)
    xlabel('pH', fontsize=20)
    ylabel('$\Delta G\'^\circ$', fontsize=14)
    title('CO2 + H2O = bicarbonate', fontsize=14) 

fig = figure(figsize=(6,6))
for I in [0, 0.1, 1]:
    draw_pH_response(I, fig)
legend(fontsize=14)


    

    
Out[54]:





<matplotlib.legend.Legend at 0x7fef3c4435d0>