Commit 0da1bfbf authored by eladnoor's avatar eladnoor

implementing phases

parent a265ad7a
......@@ -29,9 +29,10 @@ from typing import Dict, List, Tuple
import numpy as np
from component_contribution import GibbsEnergyPredictor
from equilibrator_cache.reaction import Reaction
from . import Compound, ureg
from .phased_reaction import PhasedReaction
from .phased_reaction import AQUEOUS_PHASE, PHASED_COMPOUND_DICT, PhasedReaction
from .settings import FARADAY, R, default_I, default_pH, default_pMg, default_T
......@@ -74,12 +75,20 @@ class ComponentContribution(object):
U as the covariance of a Gaussian used for sampling
(where dG0_cc is the mean of that Gaussian).
"""
return self.predictor.standard_dg_prime_multi(
reactions,
rxn_aquoues, dg_prime_nonaq = zip(
[reaction.split_aqueous_phase(
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature)
for reaction in reactions])
standard_dg_prime, dg_sigma = self.predictor.standard_dg_prime_multi(
rxn_aquoues,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
return standard_dg_prime + dg_prime_nonaq, dg_sigma
def standard_dg_prime(
self,
......@@ -92,13 +101,21 @@ class ComponentContribution(object):
calculate the confidence interval, use the range -1.96 to 1.96 times
this value
"""
return self.predictor.standard_dg_prime(
reaction,
rxn_aquoues, dg_prime_nonaq = reaction.split_aqueous_phase(
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
standard_dg_prime, dg_sigma = self.predictor.standard_dg_prime(
rxn_aquoues,
p_h=self.p_h,
ionic_strength=self.ionic_strength,
temperature=self.temperature
)
return standard_dg_prime + dg_prime_nonaq, dg_sigma
def dg_prime(
self,
reaction: PhasedReaction,
......@@ -118,9 +135,7 @@ class ComponentContribution(object):
(which is the value shown on eQuilibrator)
"""
dg_prime, dg_uncertainty = self.standard_dg_prime(
reaction
)
dg_prime, dg_uncertainty = self.standard_dg_prime(reaction)
dg_prime += self.RT * reaction.dg_correction(compound_to_conc)
return dg_prime, dg_uncertainty
......
......@@ -26,32 +26,64 @@
import logging
from typing import Dict
from collections import namedtuple
from typing import Dict, Tuple
import numpy as np
from equilibrator_cache.reaction import Reaction
from equilibrator_cache.thermodynamic_constants import legendre_transform
from . import Compound, ccache
from . import Q_, Compound, R, ccache, ureg
# list of compounds that are non aqueous (we keep their formation energies
# as they appear in component-contribution, and only override their possible
# phases to the one indicated)
AQUEOUS_PHASE = 'aq'
GAS_PHASE = 'g'
SOLID_PHASE = 's'
LIQUID_PHASE = 'l'
PhaseInfo = namedtuple("phase_info", "name "
"standard_abundance "
"physiolical_abundance "
"dimensionality")
PHASE_INFO_DICT = {
AQUEOUS_PHASE: PhaseInfo("aqueous",Q_("1 M"),Q_("1 mM"),"[concentration]"),
GAS_PHASE: PhaseInfo("gas", Q_("1 bar"), Q_("1 mbar"), "[pressure]"),
SOLID_PHASE: PhaseInfo("solid", None, None, "[dimensionless]"),
LIQUID_PHASE: PhaseInfo("liquid", None, None, "[dimensionless]")
}
# dictionary of compounds that can be non aqueous. The values are tuples of
# possible phases, starting with the default phase
NON_AQUEOUS_COMPOUND_DICT = {
"MNXM2": "liquid", # H2O
"MNXM89": "solid", # sulfur
"MNXM2": (LIQUID_PHASE, ), # H2O
"MNXM89": (SOLID_PHASE, ), # sulfur
"MNXM4": (AQUEOUS_PHASE, GAS_PHASE), # O2
"MNXM195": (AQUEOUS_PHASE, GAS_PHASE), # H2
"MNXM724": (AQUEOUS_PHASE, GAS_PHASE), # N2
"MNXM13": (AQUEOUS_PHASE, GAS_PHASE), # CO2
"MNXM10881": (AQUEOUS_PHASE, GAS_PHASE), # CO
}
# list of compounds that can also be in gas phase (besides the default
# aqueous phase). The values are the formation energy and number of protons
# in gas phase (charge must always be 0).
GAS_COMPOUND_DICT = {
"MNXM4": (0, 0), # O2
"MNXM195": (0, 2), # H2
"MNXM724": (0, 0), # N2
"MNXM13": (-394.36, 0), # CO2
"MNXM15": (None, 3), # NH3
"MNXM10881": (-137.17, 0), # CO
MicroSpecie = namedtuple("microspecie", "standard_dgf num_protons charge")
# list of compounds that can also be in non aqueous phases. The values are the
# formation energy, number of protons and charge. Values are from [Alberty's
# Thermodynamics of Biochemical Reactions, 2003]. If the standard_dgf is
# None, use the value from component_contribution
PHASED_COMPOUND_DICT = {
("MNXM4", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("0 kJ/mol"), num_protons=0, charge=0), # O2
("MNXM195", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("0 kJ/mol"), num_protons=2, charge=0), # H2
("MNXM724", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("0 kJ/mol"), num_protons=0, charge=0), # N2
("MNXM13", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("-394.36 kJ/mol"), num_protons=0, charge=0), # CO2
("MNXM10881", GAS_PHASE): MicroSpecie(
standard_dgf=Q_("-137.17 kJ/mol"), num_protons=0, charge=0), # CO
("MNXM89", SOLID_PHASE): MicroSpecie(
standard_dgf=Q_("0 kJ/mol"), num_protons=0, charge=0), # S8
}
......@@ -73,12 +105,114 @@ class PhasedReaction(Reaction):
self.reaction_id = 'COCO:R%05d' % PhasedReaction.REACTION_COUNTER
PhasedReaction.REACTION_COUNTER += 1
self.phase_dict = dict({compound: 'aqueous' for compound in
sparse.keys()})
# TODO - complete the work on the phases, i.e. make sure the default
# phase is used in gas and liquid compounds. Also, when calculating
# the dG'0 we should be aware of the phase and use the appropriate
# formation energy in case it's not aqueous
self.compound_to_phase = {cpd: PhasedReaction.get_default_phase(cpd)
for cpd in sparse.keys()}
@staticmethod
def get_default_phase(compound: Compound):
if compound.mnx_id in NON_AQUEOUS_COMPOUND_DICT:
return NON_AQUEOUS_COMPOUND_DICT[compound.mnx_id][0]
else:
return AQUEOUS_PHASE
def change_phase(self, compound: Compound, phase: str) -> None:
if compound.mnx_id in NON_AQUEOUS_COMPOUND_DICT:
assert phase in NON_AQUEOUS_COMPOUND_DICT[compound.mnx_id]
else:
assert phase == AQUEOUS_PHASE
self.compound_to_phase[compound] = phase
@ureg.check(None, None, "[concentration]", "[temperature]")
def split_aqueous_phase(
self,
p_h: float,
ionic_strength: float,
temperature: float
) -> Tuple[Reaction, float]:
"""
:param p_h:
:param ionic_strength:
:param temperature:
:return: a tuple (reaction_aqueous, dg_prime_nonaq)
"""
dg_prime_nonaq = 0
aquoues_sparse = {}
for compound, phase in self.compound_to_phase.items():
if (compound.mnx_id, phase) in PHASED_COMPOUND_DICT:
# for compound-phase pairs in this dictionary use the stored
# value and transform directly using the Legendre transform
ms = PHASED_COMPOUND_DICT[(compound.mnx_id, phase)]
dgf_prime = ms.standard_dgf + legendre_transform(
p_h=p_h, ionic_strength=ionic_strength,
temperature=temperature, num_protons=ms.num_protons,
charge=ms.charge
) * R * temperature
dg_prime_nonaq += self.get_coeff(compound) * dgf_prime
else:
# for compounds in the aqueous phase, keep them in the
# reaction object so we can later use component_contribution
# on them
aquoues_sparse[compound] = self.get_coeff(compound)
rxn_aquoues = Reaction(aquoues_sparse, arrow=self.arrow, rid=self.rid)
return rxn_aquoues, dg_prime_nonaq
def dg_correction(
self,
compound_to_conc: Dict[Compound, float]
) -> float:
"""
Calculate the concentration adjustment in the dG' of reaction.
Arguments:
compound_to_conc - a dictionary mapping Compound objects
to concentration in M (default is 1M)
Returns:
ddg_over_rt - correction for delta G in units of RT
"""
# here we check that each concentration is in suitable units,
# depending on the phase of that compound
ddg_over_rt = Q_(0.0)
for compound, conc in compound_to_conc.items():
if conc is None:
continue
if compound == Reaction.PROTON_COMPOUND:
continue
if compound in self.compound_to_phase:
p_info = PHASE_INFO_DICT[self.compound_to_phase[compound]]
if p_info.standard_abundance is None:
raise ValueError(f"compounds in {p_info.name} phase "
f"cannot have a relative abundance")
if not conc.check(p_info.dimensionality):
raise ValueError(f"abundance of {p_info.name} compounds "
f"must be in {p_info.dimensionality}")
ddg_over_rt += self.get_coeff(compound) * np.log(
compound_to_conc[compound] / p_info.standard_abundance
)
return ddg_over_rt
def physiological_dg_correction(self) -> float:
"""
Calculate the concentration adjustment in the dG' of reaction
assuming all reactants are in the default physiological
concentrations (i.e. 1 mM)
Returns:
ddg_over_rt - correction for delta G in units of RT
"""
compound_to_conc = {}
for compound, phase in self.compound_to_phase.items():
p_info = PHASE_INFO_DICT[phase]
compound_to_conc[compound] = p_info.physiolical_abundance
return self.dg_correction(compound_to_conc)
def balance_by_oxidation(self):
"""
......@@ -127,7 +261,6 @@ class PhasedReaction(Reaction):
"""
return PhasedReaction({kegg_id: -1}).balance_by_oxidation()
if __name__ == '__main__':
import sys
......
......@@ -26,8 +26,8 @@
from equilibrator_cache.thermodynamic_constants import (
FARADAY, POSSIBLE_REACTION_ARROWS, Q_, R, Rlog10, default_I, default_pH,
default_T, physiological_concentration, standard_concentration, ureg)
FARADAY, POSSIBLE_REACTION_ARROWS, Q_, R, default_I, default_pH, default_T,
physiological_concentration, standard_concentration, ureg)
default_pMg = Q_("14.0")
......
......@@ -27,43 +27,70 @@
import warnings
import numpy
import pytest
from equilibrator_api import Q_, ComponentContribution, Reaction, ccache
from equilibrator_api.phased_reaction import GAS_PHASE
from . import approx_unit
@pytest.fixture(scope="module")
def comp_contribution() -> ComponentContribution:
return ComponentContribution(p_h = Q_("7"), ionic_strength = Q_("0.1M"))
return ComponentContribution(p_h = Q_("7"), ionic_strength = Q_("0.25M"))
@pytest.fixture(scope="module")
def atp_hydrolysis() -> Reaction:
formula = "KEGG:C00002 + KEGG:C00001 = KEGG:C00008 + KEGG:C00009"
return Reaction.parse_formula(formula)
@pytest.fixture(scope="module")
def fermentation_gas() -> Reaction:
formula = "KEGG:C00031 = 2 KEGG:C00469 + 2 KEGG:C00011"
rxn = Reaction.parse_formula(formula)
rxn.change_phase(ccache.get_compound("KEGG:C00011"), GAS_PHASE)
return rxn
def test_gibbs_energy_atp_hydrolysis(comp_contribution):
def test_atp_hydrolysis_physiological_dg(comp_contribution, atp_hydrolysis):
warnings.simplefilter('ignore', ResourceWarning)
assert atp_hydrolysis.is_balanced()
assert float(atp_hydrolysis.physiological_dg_correction()) == \
pytest.approx(numpy.log(1e-3), rel=1e-3)
# ATP + H2O = ADP + Pi
formula = "KEGG:C00002 + KEGG:C00001 = KEGG:C00008 + KEGG:C00009"
reaction = Reaction.parse_formula(formula)
compound_to_conc = {ccache.get_compound('KEGG:C00001'): None,
ccache.get_compound('KEGG:C00002'): Q_("1uM"),
ccache.get_compound('KEGG:C00008'): Q_("1uM"),
ccache.get_compound('KEGG:C00009'): Q_("1uM")}
assert reaction.is_balanced()
assert float(atp_hydrolysis.dg_correction(compound_to_conc)) == \
pytest.approx(numpy.log(1e-6), rel=1e-3)
def test_fermentation_gas(comp_contribution, fermentation_gas):
assert fermentation_gas.is_balanced()
assert float(fermentation_gas.physiological_dg_correction()) == \
pytest.approx(numpy.log(1e-9), rel=1e-3)
def test_atp_hydrolysis_dg(comp_contribution, atp_hydrolysis):
warnings.simplefilter('ignore', ResourceWarning)
compound_to_conc = {ccache.get_compound('KEGG:C00002'): Q_("1mM"),
compound_to_conc = {ccache.get_compound('KEGG:C00001'): None,
ccache.get_compound('KEGG:C00002'): Q_("1mM"),
ccache.get_compound('KEGG:C00008'): Q_("10mM"),
ccache.get_compound('KEGG:C00009'): Q_("10mM")}
standard_dg_prime, dg_uncertainty = comp_contribution.standard_dg_prime(
reaction)
atp_hydrolysis)
approx_unit(standard_dg_prime, Q_("-26.4 kJ/mol"), abs=0.1)
approx_unit(standard_dg_prime, Q_("-25.8 kJ/mol"), abs=0.1)
approx_unit(dg_uncertainty, Q_("0.3 kJ/mol"), abs=0.1)
dg_prime, _ = comp_contribution.dg_prime(reaction, compound_to_conc)
approx_unit(dg_prime, Q_("-32.1 kJ/mol"), abs=0.1)
dg_prime, _ = comp_contribution.dg_prime(atp_hydrolysis, compound_to_conc)
approx_unit(dg_prime, Q_("-31.5 kJ/mol"), abs=0.1)
physiological_dg_prime, _ = comp_contribution.physiological_dg_prime(
reaction)
approx_unit(physiological_dg_prime, Q_("-43.5 kJ/mol"), abs=0.1)
atp_hydrolysis)
approx_unit(physiological_dg_prime, Q_("-42.9 kJ/mol"), abs=0.1)
def test_gibbs_energy_pyruvate_decarboxylase(comp_contribution):
......@@ -78,7 +105,7 @@ def test_gibbs_energy_pyruvate_decarboxylase(comp_contribution):
standard_dg_prime, dg_uncertainty = comp_contribution.standard_dg_prime(
reaction)
approx_unit(standard_dg_prime, Q_("-18.4 kJ/mol"), abs=0.1)
approx_unit(standard_dg_prime, Q_("-18.0 kJ/mol"), abs=0.1)
approx_unit(dg_uncertainty, Q_("3.3 kJ/mol"), abs=0.1)
......@@ -94,7 +121,7 @@ def test_reduction_potential(comp_contribution):
standard_e_prime, e_uncertainty = comp_contribution.standard_e_prime(
reaction)
approx_unit(standard_e_prime, Q_("-175.0 mV"), abs=1.0)
approx_unit(standard_e_prime, Q_("-177.0 mV"), abs=1.0)
approx_unit(e_uncertainty, Q_("3.3 mV"), abs=1.0)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment