Commit d7ebdaba authored by eladnoor's avatar eladnoor

allow calculating dGs using CC for pathways (if they are not provided in the SBtab)

parent 6b1892ec
......@@ -39,7 +39,8 @@ from . import Compound, ccache
from .bounds import Bounds
from .component_contribution import ComponentContribution
from .phased_reaction import PhasedReaction
from .settings import Q_, R, default_T, ureg
from .settings import Q_, ureg, R, default_I, default_pH, default_pMg, \
default_T
from .thermo_models import PathwayThermoModel
......@@ -65,13 +66,11 @@ class Pathway(object):
Designed for checking input prior to converting to a stoichiometric model.
"""
@ureg.check(None, None, None, "[energy]/[substance]",
"[energy]/[substance]", None)
def __init__(self,
reactions: List[PhasedReaction],
fluxes: np.array,
standard_dg_primes: np.array,
dg_uncertainties: np.array,
standard_dg_primes: np.array = None,
dg_uncertainties: np.array = None,
bounds: Bounds = None
) -> object:
"""Initialize.
......@@ -88,17 +87,28 @@ class Pathway(object):
Nr = len(reactions)
self.fluxes = fluxes
self.standard_dg_primes = standard_dg_primes
assert self.fluxes.shape == (Nr, )
assert self.standard_dg_primes.shape == (Nr,)
assert dg_uncertainties.shape == (Nr, Nr)
self.dg_uncertainties = dg_uncertainties
if standard_dg_primes is None:
self.set_aqueous_params()
self.standard_dg_primes, self.dg_uncertainties = \
self.comp_contribution.standard_dg_prime_multi(self.reactions)
# TODO: right now MDF does not work with uncertainties, and we
# need to fix that. In the meantime, we use a 0-matrix to
# indicate there is no uncertainty in the dG'0 estimates.
self.dg_uncertainties *= 0.0
else:
self.standard_dg_primes = standard_dg_primes
self.dg_uncertainties = dg_uncertainties
phys_corr = np.array(
[float(r.physiological_dg_correction()) for r in reactions]) * R \
* default_T
assert self.standard_dg_primes.shape == (Nr,)
assert self.standard_dg_primes.check("[energy]/[substance]")
assert self.dg_uncertainties.shape == (Nr, Nr)
assert self.dg_uncertainties.check("[energy]**2/[substance]**2")
phys_corr = np.array([float(r.physiological_dg_correction()) for
r in reactions]) * R * default_T
self.physiological_dg_prime = self.standard_dg_primes + phys_corr
if bounds is None:
......@@ -137,6 +147,26 @@ class Pathway(object):
"""
self.compound_names = list(map(mapping, self.S.index))
@ureg.check(None, None, None, '[concentration]', '[temperature]')
def set_aqueous_params(
self,
p_h: float = default_pH,
p_mg: float = default_pMg,
ionic_strength: float = default_I,
temperature: float = default_T
) -> None:
self.comp_contribution = ComponentContribution(
p_h=p_h, p_mg=p_mg, ionic_strength=ionic_strength,
temperature=temperature)
def calculate_dgs(self) -> None:
self.standard_dg_primes, self.dg_uncertainties = \
self.comp_contribution.standard_dg_prime_multi(self.reactions)
phys_corr = np.array([float(r.physiological_dg_correction()) for
r in self.reactions]) * R * default_T
self.physiological_dg_prime = self.standard_dg_primes + phys_corr
@property
def bounds(self) -> Tuple[Iterable[float], Iterable[float]]:
return self._bounds.GetBounds(self.S.index)
......@@ -158,8 +188,8 @@ class Pathway(object):
dg_over_rt.shape)
@property
def dg_uncertainty_prime_over_rt(self) -> np.array:
u_over_rt = self.dg_uncertainties / (R * default_T)
def sqrt_dg_uncertainty_over_rt(self) -> np.array:
u_over_rt = self.dg_uncertainties**(0.5) / (R * default_T)
return np.array(list(map(float, u_over_rt.flat))).reshape(
u_over_rt.shape)
......@@ -167,8 +197,10 @@ class Pathway(object):
def from_csv_file(cls,
f: TextIO,
bounds: Bounds = None,
pH: float = None,
ionic_strength: float = None):
p_h: float = default_pH,
p_mg: float = default_pMg,
ionic_strength: float = default_I,
temperature: float = default_T) -> object:
"""Returns a pathway parsed from an input file.
Caller responsible for closing f.
......@@ -191,11 +223,12 @@ class Pathway(object):
reactions.append(rxn)
fluxes.append(flux)
cc = ComponentContribution(pH=pH, ionic_strength=ionic_strength)
standard_dg_primes, _ = zip(*map(cc.standard_dg_prime,
reactions))
standard_dg_primes = list(standard_dg_primes)
return Pathway(reactions, fluxes, standard_dg_primes, bounds=bounds)
pp = Pathway(reactions, fluxes, bounds=bounds)
pp.set_aqueous_params(p_h=p_h, p_mg=p_mg,
ionic_strength=ionic_strength,
temperature=temperature)
pp.calculate_dgs()
return pp
@staticmethod
def get_compounds(reactions: List[PhasedReaction]) -> List[Compound]:
......@@ -288,15 +321,21 @@ class Pathway(object):
# reaction energies from the equilibrium constants
keq_df = param_df[param_df['QuantityType'] == 'equilibrium constant']
keq_df.set_index('Reaction', inplace=True)
keq_df.Value = keq_df.Value.apply(float)
keq = keq_df.loc[reaction_ids, 'Value'].values
standard_dg_primes = -np.log(keq) * R * default_T
if not set(reaction_ids).issubset(keq_df.index):
logging.warning("Not all reactions have Keq values in the "
"Parameter table. Therefore, we use "
"Component-Contribution to calculate all of them.")
pp = Pathway(reactions, v, bounds=bounds)
else:
keq_df.Value = keq_df.Value.apply(float)
keq = keq_df.loc[reaction_ids, 'Value'].values
standard_dg_primes = -np.log(keq) * R * default_T
dg_uncertainties = np.zeros((len(reactions), len(reactions))) * Q_(
"kJ/mol")
dg_uncertainties = np.zeros((len(reactions), len(reactions))) * Q_(
"kJ**2/mol**2")
pp = Pathway(reactions, v, standard_dg_primes, dg_uncertainties,
bounds=bounds)
pp = Pathway(reactions, v, standard_dg_primes, dg_uncertainties,
bounds=bounds)
# override the compound names with the ones in the SBtab
compound_to_name = dict(map(reversed, name_to_compound.items()))
......
......@@ -257,7 +257,7 @@ class PathwayThermoModel(object):
# Define and apply the constraints on the concentrations
inds = np.nonzero(np.diag(self.I_dir))[0].tolist()
A11 = self.I_dir[inds] @ self.pathway.dg_uncertainty_prime_over_rt
A11 = self.I_dir[inds] @ self.pathway.sqrt_dg_uncertainty_over_rt
A12 = self.I_dir[inds] @ self.pathway.S.T
A13 = np.ones((len(inds), 1))
......
......@@ -42,7 +42,7 @@ TEST_DIR = os.path.dirname(
def test_mdf():
warnings.simplefilter('ignore', ResourceWarning)
sbtab_fname = os.path.join(TEST_DIR, 'pathway_ethanol_SBtab.tsv')
sbtab_fname = os.path.join(TEST_DIR, 'pathway.tsv')
pp = Pathway.from_sbtab(sbtab_fname)
pp_mdf_data = pp.calc_mdf()
......@@ -53,7 +53,7 @@ def test_mdf():
"2 KEGG:C00002 + 2 KEGG:C00001 + 2 KEGG:C00469 + 2 KEGG:C00011")
assert net_rxn == reference_net_rxn
approx_unit(pp_mdf_data.mdf, Q_("1.69 kJ/mol"), abs=0.1)
approx_unit(pp_mdf_data.mdf, Q_("2.04 kJ/mol"), abs=0.1)
shadow_prices = pp_mdf_data.reaction_df.set_index(
'reaction_id').shadow_price
......
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