Skip to content
Snippets Groups Projects
Commit 798da8ba authored by Peter Schubert's avatar Peter Schubert
Browse files

start with scaled Flux GBA, step 1: model normalization

parent bd0ccaeb
Branches
No related tags found
No related merge requests found
......@@ -3,5 +3,6 @@
from .gba_problem import GbaProblem
from .rba_problem import RbaProblem
from .fba_problem import FbaProblem
from .gba_sflux_problem import GbaSfluxProblem
__all__ = ['GbaProblem', 'FbaProblem', 'RbaProblem']
__all__ = ['GbaProblem', 'FbaProblem', 'RbaProblem', 'GbaSfluxProblem']
"""Implementation of GBA scaled flux Problem class.
Based on Hugo Dourado scaled flux GBA optimization
Ref: Presentation
Growth Mechanics: the Economy, Control and Optimality of Self-Replication
Metabolic Pathway Analysis 2021, Hugo Dourado
First step: convert XBA into Hugo RGBA format
Second step: implement optimization (inline with Hugo formalism)
Peter Schubert, HHU Duesseldorf, September 2022
"""
import numpy as np
import pandas as pd
import re
import sympy as sp
import sbmlxdf
from xbanalysis.model.xba_unit_def import XbaUnit
from xbanalysis.utils.utils import expand_kineticlaw
class GbaSfluxProblem:
def __init__(self, xba_model):
"""Initialize GbaSfluxProblem from a xba_model.
:param xba_model: xba model with all required GBA parameters
:type xba_model: XbaModel
"""
self.xba_model = xba_model
# while XBA model contains enzymes for all metabolic reactions with related synthesis reactions,
# in scaled flux GBA optimization one a single protein (here we use the ribosome) is used
metab_sids = [s.id for s in xba_model.species.values() if s.sboterm.is_in_branch('SBO:0000247')]
self.sids_e = [sid for sid in metab_sids if xba_model.species[sid].constant is True]
sids_s = [sid for sid in metab_sids if xba_model.species[sid].constant is False]
sid_p = [s.id for s in xba_model.species.values() if s.sboterm.is_in_branch('SBO:0000250')][0]
self.sids_sp = sids_s + [sid_p]
self.sids_esp = self.sids_e + self.sids_sp
rids_be = [r.id for r in xba_model.reactions.values() if r.sboterm.is_in_branch('SBO:0000167')]
rid_r = xba_model.species[sid_p].ps_rid
self.rids_ber = rids_be + [rid_r]
self.mw_esp = np.array([xba_model.species[sid].mw for sid in self.sids_esp])
self.mw_ber = np.array([xba_model.species[xba_model.reactions[rid].enzyme].mw
for rid in self.rids_ber])
# mw_ber is used for scaling of kcats, weight for protein synthesis needs to be rescaled
s = xba_model.species[sid_p]
if hasattr(s, 'rba_process_costs'):
p_len = s.rba_process_costs['translation']
else:
r = xba_model.reactions[rid_r]
p_len = min(r.reactants.values())
self.mw_ber[-1] *= p_len
# retrieve S-matrix: internal metabos + total protein and metabolic reaction + ribosome synthesis
sm_espxber = xba_model.get_stoic_matrix(self.sids_esp, self.rids_ber)
# calulate mass turnover in gram reactants/products per mol reactions (g/mol)
sm_espxber_react = np.where(sm_espxber > 0, 0, sm_espxber)
sm_espxber_prod = np.where(sm_espxber < 0, 0, sm_espxber)
mass_turnover_react = -sm_espxber_react.T.dot(self.mw_esp)
mass_turnover_prod = sm_espxber_prod.T.dot(self.mw_esp)
self.mass_turnover = np.max((mass_turnover_react, mass_turnover_prod), axis=0)
# mass normalized stoichiometric matrix
self.smmn_espxber = (sm_espxber * self.mw_esp.reshape(-1, 1)) / self.mass_turnover
self.rho = xba_model.parameters['rho'].value
self.cemn = np.array([xba_model.species[sid].initial_conc * xba_model.species[sid].mw
for sid in self.sids_e])
self.cimn = np.array([xba_model.species[sid].initial_conc * xba_model.species[sid].mw
for sid in sids_s])
# identify unit ids for kinetic parameters
units_per_s = [XbaUnit({'kind': 'second', 'exp': -1.0, 'scale': 0, 'mult': 1.0})]
units_mol_per_l = [XbaUnit({'kind': 'litre', 'exp': -1.0, 'scale': 0, 'mult': 1.0}),
XbaUnit({'kind': 'mole', 'exp': 1, 'scale': 0, 'mult': 1.0})]
self.kcat_udid = self.get_unit_def_id(units_per_s)
self.km_udid = self.get_unit_def_id(units_mol_per_l)
self.kcats_global = {pid: p.value for pid, p in xba_model.parameters.items()
if p.units == self.kcat_udid}
self.kms_global = {pid: p.value for pid, p in xba_model.parameters.items()
if p.units == self.km_udid}
def get_unit_def_id(self, query_units):
for udid, ud in self.xba_model.unit_defs.items():
if ud.is_equivalent(query_units):
return udid
return None
def get_kcats(self, rid):
r = self.xba_model.reactions[rid]
kcats_local = {pid: val for pid, [val, units] in r.local_params.items() if units == self.kcat_udid}
# analyse kinetic law to identify kcat_f and kcat_r
expanded_kl = sbmlxdf.misc.mathml2numpy(r.kinetic_law)
kl = expand_kineticlaw(expanded_kl, self.xba_model.functions)
used_kcats = {kcat_id: val for kcat_id, val in (self.kcats_global | kcats_local).items()
if re.search(r'\b' + kcat_id + r'\b', kl)}
kcat_f = 0.0
kcat_r = 0.0
for kcat_id, val in used_kcats.items():
# replace other parameters with value of 1.0 (exclude selected kcat and numbers)
kl_kcat = re.sub(r'\b(?!(' + kcat_id + r'\b))[a-zA-Z_]\w*\b', str(1.0), kl)
# set two different kcat values to check which increases/reduces the rate
kl_kcat_v2 = re.sub(r'\b' + kcat_id + r'\b', str(2.0), kl_kcat)
kl_kcat_v10 = re.sub(r'\b' + kcat_id + r'\b', str(10.0), kl_kcat)
if sp.sympify(kl_kcat_v2) > sp.sympify(kl_kcat_v10):
kcat_r = val
else:
kcat_f = val
return kcat_f, kcat_r
def get_kms(self, rid):
"""
:param rid: reaction id
:type rid: string
:return:
:rtype:
"""
r = self.xba_model.reactions[rid]
kms_local = {pid: val for pid, [val, units] in r.local_params.items() if units == self.km_udid}
# analyse kinetic law to identify Km and Ki values
expanded_kl = sbmlxdf.misc.mathml2numpy(r.kinetic_law)
kl = expand_kineticlaw(expanded_kl, self.xba_model.functions)
used_kms = {km_id: val for km_id, val in (self.kms_global | kms_local).items()
if re.search(r'\b' + km_id + r'\b', kl)}
kms = {}
kis = {}
# check correlation between metabolite and KM value, scaling both equally should not change value of kl
for km_id, val in used_kms.items():
for sid in r.reactants:
# in first step: replace parameters with value of 1.0 (excluding selected km and sid)
# km_id, sid are excluded from matching, e.g. \b(?!(' + km_id + r'\b))
# parameter names starting with letter or underscore [a-zA-Z_] are matched
# numbers are not matched, as they do not start with [a-zA-Z_]
# in second phase: protected km_id and sid are replaced by values
kl_km_m = re.sub(r'\b(?!(' + km_id + r'\b))(?!(' + sid + r'\b))[a-zA-Z_]\w*\b', str(1.0), kl)
kl_km_m_v2 = re.sub(r'\b[a-zA-Z_]\w*\b', str(2.0), kl_km_m)
kl_km_m_v10 = re.sub(r'\b[a-zA-Z_]\w*\b', str(10.0), kl_km_m)
if bool(abs(sp.sympify(kl_km_m_v2) - sp.sympify(kl_km_m_v10)) <= 1e-8):
kms[sid] = val
for sid in r.products:
kl_km_m = re.sub(r'\b(?!(' + km_id + r'\b))(?!(' + sid + r'\b))[a-zA-Z_]\w*\b', str(1.0), kl)
kl_km_m_v2 = re.sub(r'\b[a-zA-Z_]\w*\b', str(2.0), kl_km_m)
kl_km_m_v10 = re.sub(r'\b[a-zA-Z_]\w*\b', str(10.0), kl_km_m)
if abs(sp.sympify(kl_km_m_v2) - sp.sympify(kl_km_m_v10)) <= 1e-8:
kis[sid] = val
return kms, kis
def mass_normalized_model(self):
"""Retrieve mass normalized model
:return: gba model mass normalized as per Hugo Dourado
:rtype: dict of pandas DataFrames
"""
protein_mn = self.rho - sum(self.cimn)
kcats = np.array([self.get_kcats(rid) for rid in self.rids_ber]).T
skcats = kcats * 3600.0 * self.mass_turnover / self.mw_ber
km = np.zeros_like(self.smmn_espxber)
ki = np.zeros_like(self.smmn_espxber)
rid2idx = {mrid: idx for idx, mrid in enumerate(self.rids_ber)}
sid2idx = {msid: idx for idx, msid in enumerate(self.sids_esp)}
for rid in self.rids_ber:
kms, kis = self.get_kms(rid)
for sid, val in kms.items():
km[sid2idx[sid], rid2idx[rid]] = val
for sid, val in kis.items():
ki[sid2idx[sid], rid2idx[rid]] = val
# scale Michaelis constants with molecluar weights
skm = km * self.mw_esp.reshape(-1, 1)
ski = ki * self.mw_esp.reshape(-1, 1)
ska = np.zeros_like(self.smmn_espxber)
mn_model = {
'N': pd.DataFrame(self.smmn_espxber, index=self.sids_esp, columns=self.rids_ber),
'KM': pd.DataFrame(skm, index=self.sids_esp, columns=self.rids_ber),
'KI': pd.DataFrame(ski, index=self.sids_esp, columns=self.rids_ber),
'KA': pd.DataFrame(ska, index=self.sids_esp, columns=self.rids_ber),
'kcat': pd.DataFrame(skcats, index=['kcat_f', 'kcat_b'], columns=self.rids_ber),
'conditions': pd.DataFrame(np.hstack((self.rho, self.cemn)),
index=['rho'] + self.sids_e, columns=[1]),
'lower_c': pd.DataFrame(np.zeros(len(self.sids_sp)),
index=self.sids_sp, columns=['lower']),
'initial': pd.DataFrame(np.hstack((self.cimn, protein_mn)),
index=self.sids_sp, columns=['initial']),
}
return mn_model
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment