From 7c80ca04d28b4c00f0e1b3e2cafa32a9d5b270c2 Mon Sep 17 00:00:00 2001
From: Peter Schubert <Peter.Schubert@hhu.de>
Date: Tue, 7 Jun 2022 17:17:00 +0200
Subject: [PATCH] OptResults added and small adjustments
---
xbanalysis/_version.py | 2 +-
xbanalysis/model/xba_parameter.py | 3 +-
xbanalysis/model/xba_reaction.py | 40 +-
xbanalysis/model/xba_species.py | 23 +-
xbanalysis/problems/gba_problem.py | 823 ++++++++++++++++++++++-------
xbanalysis/utils/__init__.py | 5 +-
xbanalysis/utils/hessians.py | 125 +++++
xbanalysis/utils/opt_results.py | 66 +++
8 files changed, 845 insertions(+), 242 deletions(-)
create mode 100644 xbanalysis/utils/hessians.py
create mode 100644 xbanalysis/utils/opt_results.py
diff --git a/xbanalysis/_version.py b/xbanalysis/_version.py
index 60596f0..2a1d8ea 100644
--- a/xbanalysis/_version.py
+++ b/xbanalysis/_version.py
@@ -1,5 +1,5 @@
"""
Definition of version string.
"""
-__version__ = "0.1"
+__version__ = "0.2.0"
program_name = 'xbanalysis'
diff --git a/xbanalysis/model/xba_parameter.py b/xbanalysis/model/xba_parameter.py
index 656dd5a..812301c 100644
--- a/xbanalysis/model/xba_parameter.py
+++ b/xbanalysis/model/xba_parameter.py
@@ -3,6 +3,7 @@
Peter Schubert, HHU Duesseldorf, May 2022
"""
+
class XbaParameter:
def __init__(self, s_parameter):
@@ -11,4 +12,4 @@ class XbaParameter:
self.sboterm = s_parameter.get('sboterm', '')
self.value = s_parameter['value']
self.constant = s_parameter['constant']
- self.units = s_parameter['units']
\ No newline at end of file
+ self.units = s_parameter['units']
diff --git a/xbanalysis/model/xba_reaction.py b/xbanalysis/model/xba_reaction.py
index 1303bf6..985af26 100644
--- a/xbanalysis/model/xba_reaction.py
+++ b/xbanalysis/model/xba_reaction.py
@@ -31,14 +31,14 @@ class XbaReaction:
def set_enzyme(self, species):
for sid in self.modifiers:
- if species[sid].sboterm == 'SBO:0000250' or species[sid].sboterm == 'SBO:0000252':
+ if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
self.enzyme = sid
break
def add_species_usage(self, species):
# for metabolic reactions add to the species consumed/produced the metabolic reaction id
# !!! we could also enter the Reaction Object references instead !!!!!!!!!!
- if (self.sboterm == 'SBO:0000167') or (self.sboterm == 'SBO:0000179') or (self.sboterm == 'SBO:0000182'):
+ if self.sboterm in ('SBO:0000167', 'SBO:0000179', 'SBO:0000182'):
for sid, val in self.products.items():
species[sid].add_m_reaction(self.id, val)
for sid, val in self.reactants.items():
@@ -46,35 +46,31 @@ class XbaReaction:
# protein degradation reaction, assume that first protein in reactants is getting degraded
# a protein degradation reaction is still a metabolic reaction
if self.sboterm == 'SBO:0000179':
- for sid in self.reactants:
- if (species[sid].sboterm == 'SBO:0000252') or (species[sid].sboterm == 'SBO:0000250'):
- species[sid].set_p_degrad_rid(self.id)
+ for sid, val in self.reactants.items():
+ if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+ species[sid].add_p_degrad_rid(self.id, -val)
break
- # protein interconversion reactions
+ # protein transitions
if self.sboterm == 'SBO:0000182':
- for sid in self.reactants:
- if (species[sid].sboterm == 'SBO:0000252') or (species[sid].sboterm == 'SBO:0000250'):
- species[sid].set_p_conversion_rid(self.id)
+ for sid, val in self.reactants.items():
+ if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+ species[sid].add_p_transition_rid(self.id, -val)
break
- for sid in self.products:
- if (species[sid].sboterm == 'SBO:0000252') or (species[sid].sboterm == 'SBO:0000250'):
- species[sid].set_p_conversion_rid(self.id)
+ for sid, val in self.products.items():
+ if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+ species[sid].add_p_transition_rid(self.id, val)
break
# for protein synthesis reactions identify the enzyme being produced from reaction.products
# to species being consumed add the enzyme product with corresponding stoichiometry
# to species being produced add the enzyme product with corresponding stoichiometry, unless
# the species is itself the enzyme product, in which case we add the ps reaction producing it.
- if (self.sboterm == 'SBO:0000589') or (self.sboterm == 'SBO:0000205') :
- enz_sid = ''
- for sid in self.products:
- if (species[sid].sboterm == 'SBO:0000252') or (species[sid].sboterm == 'SBO:0000250'):
- enz_sid = sid
- self.ps_product = enz_sid
- species[enz_sid].set_ps_rid(self.id)
+ if self.sboterm in ('SBO:0000589', 'SBO:0000205'):
+ for sid, val in self.products.items():
+ if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+ self.ps_product = sid
break
for sid, val in self.products.items():
- if sid != enz_sid:
- species[sid].add_ps_reaction(enz_sid, val)
+ species[sid].add_ps_reaction(self.id, val)
for sid, val in self.reactants.items():
- species[sid].add_ps_reaction(enz_sid, -val)
+ species[sid].add_ps_reaction(self.id, -val)
diff --git a/xbanalysis/model/xba_species.py b/xbanalysis/model/xba_species.py
index e18fdc1..cf94b57 100644
--- a/xbanalysis/model/xba_species.py
+++ b/xbanalysis/model/xba_species.py
@@ -15,14 +15,18 @@ class XbaSpecies:
self.compartment = s_species['compartment']
self.initial_conc = s_species['initialConcentration']
self.constant = s_species['constant']
+ self.boundary = s_species['boundaryCondition']
self.units = s_species['substanceUnits']
attrs = sbmlxdf.misc.extract_xml_attrs(s_species['xmlAnnotation'], token='molecule')
- self.weight = float(attrs['weight_Da'])
+ self.mw = float(attrs['weight_Da'])
+ self.ps_rid = None
self.m_reactions = {}
self.ps_reactions = {}
- self.ps_rid = None
- self.p_degrad_rid = None
- self.p_convert_rid = None
+ self.p_degrad_rids = {}
+ self.p_trans_rids = {}
+
+ def set_ps_rid(self, rid, stoic):
+ self.ps_rid = {'rid': rid, 'stoic': stoic}
def add_m_reaction(self, rid, stoic):
self.m_reactions[rid] = stoic
@@ -30,11 +34,8 @@ class XbaSpecies:
def add_ps_reaction(self, rid, stoic):
self.ps_reactions[rid] = stoic
- def set_ps_rid(self, rid):
- self.ps_rid = rid
-
- def set_p_degrad_rid(self, rid):
- self.p_degrad_rid = rid
+ def add_p_degrad_rid(self, rid, stoic):
+ self.p_degrad_rids[rid] = stoic
- def set_p_conversion_rid(self, rid):
- self.p_convert_rid = rid
+ def add_p_transition_rid(self, rid, stoic):
+ self.p_trans_rids[rid] = stoic
diff --git a/xbanalysis/problems/gba_problem.py b/xbanalysis/problems/gba_problem.py
index 1f242fb..3f73fc7 100644
--- a/xbanalysis/problems/gba_problem.py
+++ b/xbanalysis/problems/gba_problem.py
@@ -1,62 +1,93 @@
-"""Implementation of XBAmodel class.
+"""Implementation of GBA Problem class.
Peter Schubert, HHU Duesseldorf, May 2022
"""
-import re
import numpy as np
-import sympy
+import re
import types
+import sympy
+
+from xbanalysis.utils.hessians import (hess_from_grad_grad, hess_from_v_hesss, hess_from_jac_jac,
+ hesss_from_matrix_hesss, hesss_from_v_hess,
+ hesss_from_grad_jac)
+
+
+# Decorator to cache last result
+def cache_last(func):
+ def wrapper(self, x):
+ fname = wrapper.__name__
+ if fname not in GbaProblem.last_values or np.array_equal(GbaProblem.last_values[fname][0], x) is False:
+ GbaProblem.last_values[fname] = [x, func(self, x)]
+ return GbaProblem.last_values[fname][1]
+ wrapper.__name__ = func.__name__
+ return wrapper
class GbaProblem:
+ last_values = {}
+
def __init__(self, xba_model):
- self.xba_model = xba_model
+ """Initialize GbaProblem from a xba_model.
+ :param xba_model: xba model with all required GBA parameters
+ :type xba_model: XbaModel
+ """
+ self.xba_model = xba_model
+ # self.cache = True
+ self.last_gr0_x = np.zeros(1)
+ self.last_gr0 = None
# metabolites and enzymes in the model
metab_sids = [s.id for s in xba_model.species.values() if s.sboterm == 'SBO:0000247']
- enz_sids = [s.id for s in xba_model.species.values() if (s.sboterm == 'SBO:0000250' or
- s.sboterm == 'SBO:0000252')]
+ enz_sids = [s.id for s in xba_model.species.values() if s.sboterm in ('SBO:0000250', 'SBO:0000252')]
# optimization variables
self.var_metab_sids = [sid for sid in metab_sids if xba_model.species[sid].constant is False]
- self.var_enz_sids = [sid for sid in enz_sids if xba_model.species[sid].constant is False]
+ self.var_enz_sids = enz_sids
self.var_sids = self.var_metab_sids + self.var_enz_sids
- self.var_n_metabs = len(self.var_metab_sids)
+ self.n_vars = len(self.var_sids)
self.var_sid2idx = {sid: idx for idx, sid in enumerate(self.var_sids)}
+ self.var_mask_metab = [True if sid in self.var_metab_sids else False for sid in self.var_sids]
+ self.var_mask_enz = [True if sid in self.var_enz_sids else False for sid in self.var_sids]
self.var_initial = np.array([xba_model.species[sid].initial_conc for sid in self.var_sids])
- # self.var_lower_bounds = np.zeros(len(self.var_sids)) * 1
- # self.var_upper_bounds = np.ones(len(self.var_sids)) * 0.1
+ self.var_vols = np.array([xba_model.compartments[xba_model.species[sid].compartment].size
+ for sid in self.var_sids])
# constant concentrations
self.const_metab_sids = [sid for sid in metab_sids if xba_model.species[sid].constant is True]
self.const_sid2idx = {sid: idx for idx, sid in enumerate(self.const_metab_sids)}
- self.const_enz_sids = [sid for sid in enz_sids if xba_model.species[sid].constant is True]
-
- self.all_enz_sids = self.var_enz_sids + self.const_enz_sids
- # reactions - yet to implement, accept child SBO terms, e.g. here of 00167
+ # TODO: accept child SBO terms (i.e. SBO hierarchy parsing)
self.mr_rids = [r.id for r in xba_model.reactions.values()
- if ((r.sboterm == 'SBO:0000167') or
- (r.sboterm == 'SBO:0000179') or
- (r.sboterm == 'SBO:0000182'))]
+ if r.sboterm in ('SBO:0000167', 'SBO:0000179', 'SBO:0000182')]
self.ps_rids = [r.id for r in xba_model.reactions.values()
- if (r.sboterm == 'SBO:0000589') or (r.sboterm == 'SBO:0000205')]
- spont_rids = [rid for rid in self.mr_rids if xba_model.reactions[rid].enzyme == '']
-
- # calculate stoichiometric for mass balance equality constraints
- s_mat = self.get_stoic_matrix(self.var_metab_sids, self.mr_rids)
- var_sids_vol = np.array([self.xba_model.compartments[self.xba_model.species[sid].compartment].size
- for sid in self.var_metab_sids])
- s_scaled_mat = s_mat/var_sids_vol.reshape(-1, 1)
-
- all_enz_ps_rids = [xba_model.species[sid].ps_rid for sid in self.all_enz_sids]
- p_mat = self.get_stoic_matrix(self.var_metab_sids, all_enz_ps_rids)
- all_enz_vol = np.array([self.xba_model.compartments[self.xba_model.species[sid].compartment].size
- for sid in self.all_enz_sids])
- vol_scale_p_mat = np.outer(1.0/var_sids_vol, all_enz_vol)
- p_scaled_mat = p_mat * vol_scale_p_mat
+ if r.sboterm in ('SBO:0000589', 'SBO:0000205')]
+ # spont_rids = [rid for rid in self.mr_rids if xba_model.reactions[rid].enzyme == '']
+ # calculate stoichiometric sub-matrices and scale them accoring to volume
+ # Note: reordering of P submatrix
+ # - vector of enzyme concentrations is extracted from optimization vector x
+ # - sequenze of enzyme concentrations is according to var_enz_sids
+ # - there is a one-to-one relationship between enzymes and protein synthesis reactions
+ # - columns of the P submatrix have to be ordered along var_enz_sids for correct P.dot(ce)
+ # TODO: sparce matrix storage and calculations
+
+ self.var_enz_ps_rids = [list(xba_model.species[sid].ps_reactions.keys())[0]
+ for sid in self.var_enz_sids]
+ self.ps_stoic = np.array([list(xba_model.species[sid].ps_reactions.values())[0]
+ for sid in self.var_enz_sids])
+
+ metab_x_mrs = self.get_stoic_matrix(self.var_metab_sids, self.mr_rids)
+ metab_x_psrs = self.get_stoic_matrix(self.var_metab_sids, self.var_enz_ps_rids)
+ self.enz_x_mrs = self.get_stoic_matrix(self.var_enz_sids, self.mr_rids)
+ enz_vols = self.var_vols[self.var_mask_enz]
+ metab_vols = self.var_vols[self.var_mask_metab]
+
+ self.pstma_scale = enz_vols / self.ps_stoic
+ self.metab_x_mras = metab_x_mrs / metab_vols.reshape((-1, 1))
+ self.metab_x_psras = metab_x_psrs / metab_vols.reshape((-1, 1))
+ self.enz_x_mras = self.enz_x_mrs / enz_vols.reshape((-1, 1))
+ self.pstm_scale = enz_vols / self.ps_stoic
# model parameters
self.model_params = {pid: p.value for pid, p in xba_model.parameters.items()}
@@ -64,28 +95,24 @@ class GbaProblem:
self.model_params['macro'] = self.model_params.get('macro', 0.0) > 0.0
self.model_params['ext_conc'] = np.array([xba_model.species[sid].initial_conc
for sid in self.const_metab_sids])
- self.model_params['const_enz_conc'] = np.array([xba_model.species[sid].initial_conc
- for sid in self.const_enz_sids])
- self.model_params['S_scaled_mat'] = s_scaled_mat
- self.model_params['P_scaled_mat'] = p_scaled_mat
- # self.model_params['P_scaled_csr'] = scipy.sparse.csr_matrix(p_mat * vol_scale_p_mat)
self.model_params['np'] = np
- self.model_params['weights'] = np.array([xba_model.species[sid].weight for sid in self.var_sids])
- self.model_params['const_enz_weights'] = np.array([xba_model.species[sid].weight
- for sid in self.const_enz_sids])
- self.model_params['weights'] = np.array([xba_model.species[sid].weight for sid in self.var_sids])
+ self.model_params['mws'] = np.array([xba_model.species[sid].mw for sid in self.var_sids])
# create functions
- self.mrs, self.mrs_jac = self.get_reactions(self.mr_rids)
- self.spontrs, self.spontrs_jac = self.get_reactions(spont_rids)
- self.pss, self.pss_jac = self.get_reactions(self.ps_rids)
- self.inv_ps_rates, self.inv_ps_rates_jac = self.get_inv_ps_rates()
- # self.gammas, self.gammas_jac = self.get_gammas()
- self.p_degrads, self.p_degrads_jac = self.get_p_degrad_rates()
- self.enz_conv_mb, self.enz_conv_mb_jac = self.get_enz_interconv_mb()
+ self.mras, self.mras_jac, self.mras_hesss = self.get_reactions_apt(self.mr_rids)
+ self.pstmas, self.pstmas_jac, self.pstmas_hesss = self.get_psm_times_apt()
+ self.enz_trans, self.enz_trans_jac, self.enz_trans_hesss = self.get_enz_transitions()
def get_stoic_matrix(self, sids, rids):
- # matrix for free metabolites and list of reactions
+ """retrieve stoichiometric sub-matrix [sids x rids].
+
+ :param sids: species ids - rows of sub-matrix
+ :type sids: list of strings
+ :param rids: reaction ids - columns of sub-matrix
+ :type rids: list of strings
+ :returns: stoic_mat, stoichmetric sub-matrix
+ :rtype: np.ndarray
+ """
stoic_mat = np.zeros((len(sids), len(rids)))
map_metab_row = {sid: idx for idx, sid in enumerate(sids)}
for col, rid in enumerate(rids):
@@ -101,94 +128,147 @@ class GbaProblem:
self.model_params[key] = val
def _to_vectors(self, rs_str):
- # replace free and constant metabolites by vector elements
+ """replace free and constant metabolites by vector elements.
+
+ :param rs_str: reaction string (extracted from the xba model).
+ :type rs_str: list of strings
+ :returns: rs_v_str
+ :rtype: list of reaction stings with some parameters replaced by vector elements
+ """
rs_v_str = []
for r_v_str in rs_str:
for idx, sid in enumerate(self.var_sids):
r_v_str = re.sub(r'\b' + sid + r'\b', f'x[{idx}]', r_v_str)
for idx, sid in enumerate(self.const_metab_sids):
r_v_str = re.sub(r'\b' + sid + r'\b', f'ext_conc[{idx}]', r_v_str)
- for idx, sid in enumerate(self.const_enz_sids):
- r_v_str = re.sub(r'\b' + sid + r'\b', f'const_enz_conc[{idx}]', r_v_str)
rs_v_str.append(r_v_str)
return rs_v_str
def _gradients(self, rs_v_str):
- # generate derivatives
+ """generate code object for Jacobian and Hessians (lower triangle) of reactions.
+
+ :param rs_v_str: reaction strings with concentrations replaced by vector elements.
+ :type rs_v_str: list of string
+ :returns: jac
+ :rtype: code object to generate Jacobian
+ :returns: hess
+ :rtype: code object to generate Hessians (lower triangle)
+ """
x = sympy.IndexedBase('x')
ext_conc = sympy.IndexedBase('ext_conc')
- const_enz_conc = sympy.IndexedBase('const_enz_conc')
grads_str = []
+ hesss_str = []
for i, r_v_str in enumerate(rs_v_str):
- func_sym = sympy.sympify(r_v_str, locals={'x': x,
- 'ext_conc': ext_conc,
- 'const_enz_conc': const_enz_conc})
- partial_der = [str(func_sym.diff(x[idx])) for idx in range(len(self.var_sids))]
- grads_str.append('[' + ', '.join(partial_der) + ']')
+ func_sym = sympy.sympify(r_v_str, locals={'x': x, 'ext_conc': ext_conc})
+ partials_x = []
+ hess = []
+ for idx1 in range(self.n_vars):
+ partials_x.append(str(func_sym.diff(x[idx1])))
+ partials_xx = []
+ for idx2 in range(self.n_vars):
+ if idx2 <= idx1:
+ partials_xx.append(str(func_sym.diff(x[idx1], x[idx2])))
+ else:
+ partials_xx.append('0')
+ hess.append('[' + ', '.join(partials_xx) + ']')
+ grads_str.append('[' + ', '.join(partials_x) + ']')
+ hesss_str.append('[' + ', '.join(hess) + ']')
# progress_bar(i, len(rs_v_str))
jac = compile('def jac(x): return np.array([' + ', '.join(grads_str) + '])',
'<string>', 'exec')
- return jac
-
- def get_reactions(self, rids):
+ hesss = compile('def hess(x): return np.array([' + ', '.join(hesss_str) + '])',
+ '<string>', 'exec')
+ return jac, hesss
+
+ def get_reactions_apt(self, rids):
+ """create function objects for reactions rates, Jacobian and Hessian
+
+ extracted from kinetic laws for SBML model
+
+ Note: reactions in the SBML model are in amounts per time (not concentration per time)
+ I.e. for mass balances equations, we have to devide by volumes where metabolite is located
+
+ :param rids: reaction ids for reactions in amount per time to be processed.
+ :type rids: list of strings
+ :returns: ras
+ :rtype: function object, returning an 1D array of function results
+ :returns: ras_jac
+ :rtype: function object, returning a 2D array with the Jacobian matrix ∂/∂xi
+ :returns: ras_hesss
+ :rtype: function object, returning a 3D array (or sparce/lower triangele) of the Hessian
+ matrices for each individual function
+ """
kinetics_str = [self.xba_model.reactions[rid].expanded_kl for rid in rids]
v_str = self._to_vectors(kinetics_str)
- func_code = compile('def reactions(x): return np.array([' + ', '.join(v_str) + '])',
- '<string>', 'exec')
- rs = types.FunctionType(func_code.co_consts[0], self.model_params)
- rs_jac = types.FunctionType(self._gradients(v_str).co_consts[0], self.model_params)
- return rs, rs_jac
-
- def get_inv_ps_rates(self):
- # inverse of max protein synthesis rates as function of variable enzymes and constant enzymes
- ps_inv_rs_str = []
- for sid in self.all_enz_sids:
- sp = self.xba_model.species[sid]
- ps_rid = sp.ps_rid
- ps_inv_rs_str.append(f'{sp.compartment}/({self.xba_model.reactions[ps_rid].expanded_kl})')
- v_str = self._to_vectors(ps_inv_rs_str)
- func_code = compile('def ps_inv(x): return np.array([' + ', '.join(v_str) + '])',
+ func_code = compile('def rs_apt(x): return np.array([' + ', '.join(v_str) + '])',
'<string>', 'exec')
- inv_ps_rates = types.FunctionType(func_code.co_consts[0], self.model_params)
- inv_ps_rates_jac = types.FunctionType(self._gradients(v_str).co_consts[0], self.model_params)
- return inv_ps_rates, inv_ps_rates_jac
-
- def get_p_degrad_rates(self):
- # protein degradation rates
- degrads_str = []
- for sid in self.all_enz_sids:
- sp = self.xba_model.species[sid]
- if sp.p_degrad_rid is not None:
- degrads_str.append(f'({self.xba_model.reactions[sp.p_degrad_rid].expanded_kl})/' +
- f'{sp.compartment}')
- else:
- degrads_str.append('0')
- v_str = self._to_vectors(degrads_str)
- func_code = compile('def p_degrads(x): return np.array([' + ', '.join(v_str) + '])',
+ jac_code, hesss_code = self._gradients(v_str)
+ ras = types.FunctionType(func_code.co_consts[0], self.model_params)
+ ras_jac = types.FunctionType(jac_code.co_consts[0], self.model_params)
+ ras_hesss = types.FunctionType(hesss_code.co_consts[0], self.model_params)
+ return ras, ras_jac, ras_hesss
+
+ def get_psm_times_apt(self):
+ """Create function objects for reaction times, Jacobian and Hessian sorted along var_enz_sids
+
+ Maximum protein synthesis rates from model converted to minimum protein synthesis times
+
+ Note: reactions in the SBML model are in amounts per time (not concentration per time)
+
+ :returns: pstmas, returning a 1D array of function results
+ :rtype: function object
+ :returns: pstmas_jac, returning a 2D array with the Jacobian matrix ∂/∂xi
+ :rtype: function object
+ :returns: pstmas_hesss , returning a 3D array (or sparce/lower triangele) of the Hessian
+ matrices for each individual function
+ :rtype: function object
+ """
+ pstmas_str = []
+ for idx, sid in enumerate(self.var_enz_sids):
+ ps_rid = self.var_enz_ps_rids[idx]
+ pstmas_str.append(f'1.0/({self.xba_model.reactions[ps_rid].expanded_kl})')
+ v_str = self._to_vectors(pstmas_str)
+ func_code = compile('def pstms_apt(x): return np.array([' + ', '.join(v_str) + '])',
'<string>', 'exec')
- degrads = types.FunctionType(func_code.co_consts[0], self.model_params)
- degrads_jac = types.FunctionType(self._gradients(v_str).co_consts[0], self.model_params)
- return degrads, degrads_jac
-
- def get_enz_interconv_mb(self):
+ jac_code, hesss_code = self._gradients(v_str)
+ pstmas = types.FunctionType(func_code.co_consts[0], self.model_params)
+ pstmas_jac = types.FunctionType(jac_code.co_consts[0], self.model_params)
+ pstmas_hesss = types.FunctionType(hesss_code.co_consts[0], self.model_params)
+ return pstmas, pstmas_jac, pstmas_hesss
+
+ def get_enz_transitions(self):
+ """Create function objects for protein interconversion, Jacobian and Hessian
+
+ Assumed, that there are only zero of few such conversion in the model
+
+ Note: reactions in the SBML model are in amounts per time (not concentration per time)
+
+ :returns: enz_trans, returning a 1D array of function results
+ :rtype: function object
+ :returns: enz_trans_jac, returning a 2D array with the Jacobian matrix ∂/∂xi
+ :rtype: function object
+ :returns: enz_trans_hesss, returning a 3D array of the Hessian
+ matrices for each individual function
+ :rtype: function object
+ """
# mass balance for interconvering enzymes
- # identify enzymes in ExMr submatrix (enzymes x metabolic reactions)
+ # identify enzymes in enz x mra submatrix (enzymes x metabolic reactions)
# identify sets of enzymes that interconvert (consider unique sets)
enz_in_mrs_sids = []
- for enz_sid in self.all_enz_sids:
+ for enz_sid in self.var_enz_sids:
sp = self.xba_model.species[enz_sid]
partner_enz_sids = {enz_sid}
- for mr_rid, enz_stoic in sp.m_reactions.items():
- rid = self.xba_model.reactions[mr_rid]
+ for p_trans_rid, enz_stoic in sp.p_trans_rids.items():
+ p_trans_r = self.xba_model.reactions[p_trans_rid]
if enz_stoic < 0:
- for other_sid, other_stoic in rid.products.items():
+ for other_sid, other_stoic in p_trans_r.products.items():
sp_other = self.xba_model.species[other_sid]
- if (sp_other.sboterm == 'SBO:0000250') or (sp_other.sboterm == 'SBO:0000252'):
+ if sp_other.sboterm in ('SBO:0000250', 'SBO:0000252'):
partner_enz_sids.add(other_sid)
if enz_stoic > 0:
- for other_sid, other_stoic in rid.reactants.items():
+ for other_sid, other_stoic in p_trans_r.reactants.items():
sp_other = self.xba_model.species[other_sid]
- if (sp_other.sboterm == 'SBO:0000250') or (sp_other.sboterm == 'SBO:0000252'):
+ if sp_other.sboterm in ('SBO:0000250', 'SBO:0000252'):
partner_enz_sids.add(other_sid)
# look for enzyme interconversion
if len(partner_enz_sids) > 1:
@@ -202,128 +282,463 @@ class GbaProblem:
enz_sid = partner_enz_sids.pop()
sp = self.xba_model.species[enz_sid]
parts_str = []
- for mr_rid, enz_stoic in sp.m_reactions.items():
- parts_str.append((f'({enz_stoic} * {self.xba_model.reactions[mr_rid].expanded_kl})/' +
+ for p_trans_rid, enz_stoic in sp.p_trans_rids.items():
+ parts_str.append((f'({enz_stoic} * {self.xba_model.reactions[p_trans_rid].expanded_kl})/' +
f'{sp.compartment}'))
conversions_str.append(' + '.join(parts_str))
v_str = self._to_vectors(conversions_str)
- func_code = compile('def enz_conv_mb(x): return np.array([' + ', '.join(v_str) + '])',
- '<string>', 'exec')
- enz_conv_mb = types.FunctionType(func_code.co_consts[0], self.model_params)
- enz_conv_mb_jac = types.FunctionType(self._gradients(v_str).co_consts[0], self.model_params)
- return enz_conv_mb, enz_conv_mb_jac
-
- def get_gammas(self):
- # ratios between protein degradation and max protein synthesis rates
- gammas_str = []
- for sid in self.all_enz_sids:
- sp = self.xba_model.species[sid]
- if sp.p_degrad_rid is not None:
- degrad_rid = sp.p_degrad_rid
- ps_rid = sp.ps_rid
- gammas_str.append(f'({self.xba_model.reactions[degrad_rid].expanded_kl})/' +
- f'({self.xba_model.reactions[ps_rid].expanded_kl})')
- else:
- gammas_str.append('0')
- v_str = self._to_vectors(gammas_str)
- func_code = compile('def gammas(x): return np.array([' + ', '.join(v_str) + '])',
+ if len(v_str) == 0:
+ v_str = '0'
+ func_code = compile('def enz_trans(x): return np.array([' + ', '.join(v_str) + '])',
'<string>', 'exec')
- gammas = types.FunctionType(func_code.co_consts[0], self.model_params)
- gammas_jac = types.FunctionType(self._gradients(v_str).co_consts[0],
- self.model_params)
- return gammas, gammas_jac
+ jac_code, hesss_code = self._gradients(v_str)
+ enz_trans = types.FunctionType(func_code.co_consts[0], self.model_params)
+ enz_trans_jac = types.FunctionType(jac_code.co_consts[0], self.model_params)
+ enz_trans_hesss = types.FunctionType(hesss_code.co_consts[0], self.model_params)
+ return enz_trans, enz_trans_jac, enz_trans_hesss
+
+# @cache_last
+ def get_growth_rate_0(self, x):
+ """Calculate growth rate at x, assuming no protein degradation
+
+ considering stoichiometric coefficiens for enzyme turnover
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: growth rate at x (without protein degradation)
+ :rtype: float
+ """
+ # if self.cache:
+ # # fname = inspect.stack()[0][3]
+ # if np.array_equal(self.last_gr0_x, x):
+ # return self.last_gr0
+ # gr0 = 1/tau
+ # tau = ce_1 * pstms_1 + ce_2 * pstms_2 + ... + ce_r * pstms_r = sum(ce_i * pstms_i)
+ ce_var = x[self.var_mask_enz]
+ pstms = self.pstma_scale * (self.pstmas(x))
+ gr0 = 1.0 / pstms.dot(ce_var)
+ # if self.cache:
+ # self.last_gr0_x = x.copy()
+ # self.last_gr0 = gr0
+ return gr0
+
+ def get_growth_rate_0_grad(self, x):
+ """Calculate growth rate gradient at x, assuming no protein degradation
+
+ gr0 = 1/tau(x)
+ tau = ce.dot(pstms(x))
+ gr0_grad = ∂(1/tau)/∂x = -(∂tau/∂x)/tau**2 = -gr0**2 * tau_grad
+ tau_grad = ce_jac.dot(pstms(x)) + ce.dot(pstms_jac(x))
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: growth rate gradient at x (without protein degradation)
+ :rtype: numpy.ndarray (1D)
+ """
+ gr0 = self.get_growth_rate_0(x)
+ ce = x[self.var_mask_enz]
+ ce_jac = np.hstack((np.zeros((len(ce), len(x) - len(ce))), np.eye(len(ce))))
+ enz_to_pstms = self.pstma_scale * self.pstmas(x)
+ enz_to_pstms_jac = self.pstma_scale.reshape((-1, 1)) * self.pstmas_jac(x)
+ tau_grad = ce.dot(enz_to_pstms_jac) + enz_to_pstms.dot(ce_jac)
+ gr0_grad = - gr0 ** 2 * tau_grad
+ return np.nan_to_num(gr0_grad, nan=0.0)
+
+ def get_growth_rate_0_hess(self, x):
+ """Calculate growth rate Hessian at x, assuming no protein degradation
+
+ gr0 = 1/tau(x)
+ tau = ce.dot(pstms(x))
+
+ gr0_grad = -gr0**2 * tau_grad
+ tau_grad = ce.dot(pstms_jac(x)) + pstms(x).dot(ce_jac)
+
+ gr0_hess = -2 * gr0 * gr0_grad * tau_grad - gr0**2 * tau_hess
+ gr0_hess = -2 gr0 * -gr0**2 * tau_grad^2 - gr0**2 * tau_hess
+ gr0_hess = gr0**2 * (2 * gr0 * tau_grad^2 - tau_hess)
+ tau_hess = ∂(tau_grad)/∂x = ce.dot(pstms_hesss(x)) + ce_jac.dot(pstms_jac(x))
+ + pstms_jac(x).dot(ce_jac) + pstms.dot(ce_hess)
+ ce_hess = 0
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: growth rate Hessian at x (without protein degradation)
+ :rtype: numpy.ndarray (2D) with shape [x,x]
+ """
+ gr0 = self.get_growth_rate_0(x)
+ gr0_grad = self.get_growth_rate_0_grad(x)
+
+ ce = x[self.var_mask_enz]
+ ce_jac = np.hstack((np.zeros((len(ce), len(x) - len(ce))), np.eye(len(ce))))
+ enz_to_pstms = self.pstma_scale * self.pstmas(x)
+ enz_to_pstms_jac = self.pstma_scale.reshape((-1, 1)) * self.pstmas_jac(x)
+ enz_to_pstms_hesss = (self.pstma_scale ** 2).reshape((-1, 1, 1)) * self.pstmas_hesss(x)
+ tau_grad = ce.dot(enz_to_pstms_jac) + enz_to_pstms.dot(ce_jac)
+
+ hess_a = hess_from_grad_grad(gr0_grad, tau_grad)
+ hess_b1 = hess_from_v_hesss(ce, enz_to_pstms_hesss)
+ hess_b2 = hess_from_jac_jac(ce_jac, enz_to_pstms_jac)
+ hess_b3 = hess_from_jac_jac(enz_to_pstms_jac, ce_jac)
+ # gr0_hess = -2 * gr0 * gr0_grad * tau_grad - gr0**2 * tau_hess
+ # gr0_hess = -2 gr0 * -gr0**2 * tau_grad^2 - gr0**2 * tau_hess
+ # gr0_hess = gr0**2 * (2 * gr0 * tau_grad^2 - tau_hess)
+ # tau_hess = ∂(tau_grad)/∂x = ce.dot(pstms_hesss(x)) + ce_jac.dot(pstms_jac(x))
+ # + pstms_jac(x).dot(ce_jac) + pstms.dot(ce_hess)
+
+ gr0_hess = - 2 * gr0 * hess_a - gr0 ** 2 * (hess_b1 + hess_b2 + hess_b3)
+ return np.nan_to_num(gr0_hess, nan=0.0)
def get_growth_rate(self, x):
- # gr = gr0 * (1 - sum(gammas))
- # gr0 = 1/tau
- # tau = sum(ce/psme) = sum(ce * 1/psme)
- ce_all = np.append(x[self.var_n_metabs:], self.model_params['const_enz_conc'])
- tau = self.inv_ps_rates(x).dot(ce_all)
- gr0 = 1.0/tau
- # gr = gr0 * (1.0 - sum(gammas(x)))
- gr = gr0 * (1.0 - sum(self.p_degrads(x) * self.inv_ps_rates(x)))
+ """Calculate growth rate at x, considering protein degradation
+
+ considering stoichiometric coefficiens unequal unity for enzyme turnover
+ stoichiometric sub-matrices do not require volume scaling (as Vol_enz gets factored out)
+
+ gr = gr0 * (1 + enz_to_mras(x).dot(enz_to_pstmas(x)))
+ enz_to_mras = enz_x_mrs.dot(mras(x))
+ enz_to_pstmas = 1/ps_stoic * pstmas(x)
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: growth rate at x (considering protein degradation)
+ :rtype: float
+ """
+ gr0 = self.get_growth_rate_0(x)
+ if np.count_nonzero(self.enz_x_mrs) == 0:
+ gr = gr0
+ else:
+ enz_to_mras = self.enz_x_mrs.dot(self.mras(x))
+ enz_to_pstmas = 1.0 / self.ps_stoic * self.pstmas(x)
+ gr = gr0 * (1.0 + enz_to_mras.dot(enz_to_pstmas))
return gr
def get_growth_rate_grad(self, x):
- # gr0_grad = d(1/tau)dx = (-dtau/dx)/tau**2 = - dtau/dx * gr0**2
- # gr_grad = gr0_grad * (1 - sum(gammas)) - gr0 * sum(gammas_jac)
- ce_all = np.append(x[self.var_n_metabs:], self.model_params['const_enz_conc'])
- tau = self.inv_ps_rates(x).dot(ce_all)
- gr0 = 1.0/tau
- dtau_dx = ce_all.dot(self.inv_ps_rates_jac(x))
- dtau_dx[self.var_n_metabs:] += self.inv_ps_rates(x)[:len(self.var_enz_sids)]
- gr0_grad = -dtau_dx * gr0**2
- deg_term = (1.0 - sum(self.p_degrads(x) * self.inv_ps_rates(x)))
- deg_term_jac = (self.p_degrads(x).reshape(-1, 1) * self.inv_ps_rates_jac(x) +
- self.inv_ps_rates(x).reshape(-1, 1)
- * self.p_degrads_jac(x)).sum(axis=0)
- gr_grad = gr0_grad * deg_term - gr0 * deg_term_jac
+ """Calculate growth rate gradient at x, considering protein degradation
+
+ gr = gr0 * (1 + enz_to_mras(x).dot(enz_to_pstmas(x)))
+ enz_to_mras = enz_x_mrs.dot(mras(x))
+ enz_to_pstmas = 1/ps_stoic * pstmas(x)
+ gr_grad = gr0_grad * (1 + enz_to_mras(x).dot(enz_to_pstmas(x)))
+ + gr0 * (enz_to_mras(x).dot(enz_to_pstmas_jac(x)) + enz_to_pstmas(x).dot(enz_to_mras_jac))
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: growth rate gradient at x (considering protein degradation)
+ :rtype: numpy.ndarray (1D)
+ """
+ gr0_grad = self.get_growth_rate_0_grad(x)
+ if np.count_nonzero(self.enz_x_mrs) == 0:
+ gr_grad = gr0_grad
+ else:
+ gr0 = self.get_growth_rate_0(x)
+ enz_to_mras = self.enz_x_mrs.dot(self.mras(x))
+ enz_to_mras_jac = self.enz_x_mrs.dot(self.mras_jac(x))
+ enz_to_pstmas = 1.0 / self.ps_stoic * self.pstmas(x)
+ enz_to_pstmas_jac = 1.0 / self.ps_stoic.reshape((-1, 1)) * self.pstmas_jac(x)
+ gr_grad = (gr0_grad * (1.0 + enz_to_mras.dot(enz_to_pstmas))
+ + gr0 * (enz_to_mras.dot(enz_to_pstmas_jac) + enz_to_pstmas.dot(enz_to_mras_jac)))
return np.nan_to_num(gr_grad, nan=0.0)
+ def get_growth_rate_hess(self, x):
+ """Calculate growth rate Hessian at x, considering protein degradation
+
+ gr = gr0 * (1 + enz_to_mras(x).dot(enz_to_pstmas(x)))
+ enz_to_mras = enz_x_mrs.dot(mras(x))
+ enz_to_pstmas = 1/ps_stoic * pstmas(x)
+
+ gr_grad = gr0_grad * (1 + enz_to_mras(x).dot(enz_to_pstmas(x)))
+ + gr0 * (enz_to_mras(x).dot(enz_to_pstmas_jac(x)) + enz_to_pstmas(x).dot(enz_to_mras_jac(x)))
+
+ gr_hess = gr0_hess * (1 + enz_to_mras(x).dot(enz_to_pstmas(x)))
+ + gr0_grad * (enz_to_mras.dot(enz_to_pstmas_jac) + enz_to_pstmas.dot(enz_to_mras_jac))
+ + gr0_grad * (enz_to_mras.dot(enz_to_pstmas_jac) + enz_to_pstmas.dot(enz_to_mras_jac))
+ + gr0 * (enz_to_mras_jac.dot(enz_to_pstmas_jac + enz_to_mras.dot(enz_to_pstmas_hesss)
+ + enz_to_pstmas_jac.dot(enz_to_mras_jac) + enz_to_pstmas.dot(enz_to_mras_hesss) )
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: growth rate Hessian at x (considering protein degradation)
+ :rtype: numpy.ndarray (2D) with shape [x,x]
+ """
+ gr0_hess = self.get_growth_rate_0_hess(x)
+ if np.count_nonzero(self.enz_x_mrs) == 0:
+ gr_hess = gr0_hess
+ else:
+ gr0 = self.get_growth_rate_0(x)
+ gr0_grad = self.get_growth_rate_0_grad(x)
+ enz_to_mras = self.enz_x_mrs.dot(self.mras(x))
+ enz_to_mras_jac = self.enz_x_mrs.dot(self.mras_jac(x))
+ enz_to_mras_hesss = hesss_from_matrix_hesss(self.enz_x_mrs, self.mras_hesss(x))
+ enz_to_pstmas = 1.0 / self.ps_stoic * self.pstmas(x)
+ enz_to_pstmas_jac = 1.0 / self.ps_stoic.reshape((-1, 1)) * self.pstmas_jac(x)
+ enz_to_pstmas_hesss = 1.0 / (self.ps_stoic ** 2).reshape((-1, 1, 1)) * self.pstmas_hesss(x)
+
+ hess_a = gr0_hess * (1.0 + enz_to_mras.dot(enz_to_pstmas))
+ tau_grad = enz_to_mras.dot(enz_to_pstmas_jac) + enz_to_pstmas.dot(enz_to_mras_jac)
+ hess_b = hess_from_grad_grad(gr0_grad, tau_grad)
+ hess_c1 = hess_from_jac_jac(enz_to_mras_jac, enz_to_pstmas_jac)
+ hess_c2 = hess_from_v_hesss(enz_to_mras, enz_to_pstmas_hesss)
+ hess_c3 = hess_from_jac_jac(enz_to_pstmas_jac, enz_to_mras_jac)
+ hess_c4 = hess_from_v_hesss(enz_to_pstmas, enz_to_mras_hesss)
+ # gr_hess = gr0_hess * (1 + enz_to_mras(x).dot(enz_to_pstmas(x)))
+ # + gr0_grad * (enz_to_mras.dot(enz_to_pstmas_jac) + enz_to_pstmas.dot(enz_to_mras_jac))
+ # + gr0_grad * (enz_to_mras.dot(enz_to_pstmas_jac) + enz_to_pstmas.dot(enz_to_mras_jac))
+ # + gr0 * (enz_to_mras_jac.dot(enz_to_pstmas_jac + enz_to_mras.dot(enz_to_pstmas_hesss)
+ # + enz_to_pstmas_jac.dot(enz_to_mras_jac) + enz_to_pstmas.dot(enz_to_mras_hesss))
+
+ gr_hess = hess_a + 2 * hess_b + gr0 * (hess_c1 + hess_c2 + hess_c3 + hess_c4)
+ return np.nan_to_num(gr_hess, nan=0.0)
+
def get_alphas(self, x):
- # alpha_e = gr * c_e/psm_e + gamma_e
- ce_all = np.append(x[self.var_n_metabs:], self.model_params['const_enz_conc'])
- alpha0 = self.get_growth_rate(x) * self.inv_ps_rates(x)*ce_all
- deg_term = self.p_degrads(x) * self.inv_ps_rates(x)
- alphas = alpha0 + deg_term
+ """Calculate ribosomal allocations at x, considering protein degradation
+
+ alpha_e = gr * ce_e * pstm_e + p_degr_e * pstm_e
+ alpha = gr * ce.dot(pstms(x)) + p_degrs(x).dot(pstms(x))
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: ribosomal allocations at x (considering protein degradation)
+ :rtype: numpy.ndarray (1D)
+ """
+ ce = x[self.var_mask_enz]
+ gr = self.get_growth_rate(x)
+
+ enz_to_pstms = self.pstma_scale * self.pstmas(x)
+ alphas0 = gr * ce * enz_to_pstms
+ if np.count_nonzero(self.enz_x_mrs) == 0:
+ alphas = alphas0
+ else:
+ enz_to_mras = self.enz_x_mrs.dot(self.mras(x))
+ enz_to_pstmas = 1.0 / self.ps_stoic * self.pstmas(x)
+ alphas = alphas0 - enz_to_mras * enz_to_pstmas
return alphas
def heq_mass_balance(self, x):
- # considering syntesis costs for protein degradation
- # mass balance for species involved in metabolic reactions:
- # dMb/dt: S.dot(mrs) + deg_term + gr*(P.dot(ce) - cm) = 0
- # deg_term: P.dot(pdegs)
- # assuming protein degradation functions (pdegs) have soichiomtry of -1 for enzyme
- cm_var = x[:self.var_n_metabs]
- ce_all = np.append(x[self.var_n_metabs:], self.model_params['const_enz_conc'])
- degrad_costs = self.model_params['P_scaled_mat'].dot(self.p_degrads(x))
- metabolism = self.model_params['S_scaled_mat'].dot(self.mrs(x))
- pem = self.model_params['P_scaled_mat'].dot(ce_all) - cm_var
- dilution = self.get_growth_rate(x) * pem
- metabs_mb = metabolism + dilution + degrad_costs
- return metabs_mb
+ """Calculate mass balance equations at x for metabolites
+
+ incuding protein degradation , which consume metabolites
+ Note: protein degradation is also part of metabolic reactions (mras), which return
+ metabolites.
+
+ metab_mb = dcm/dt = metabolic_balance + growth_dilution + sythesis_enzyme_balance
+ metabolic_balance = metab_x_mras.dot(mras(x))
+ growth_dilution = λ * (metab_x_psras.dot(ce_scaled) - cm)
+ ce_scaled = pstm_scale * ce # convert concentration to amounts, considering stoichiometry
+ sythesis_enzyme_balance = - metab_x_psras.dot(pstm_scale * enzyme_balance)
+ enzyme_balance = enz_x_mras.dot(mras(x))
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: mass balance calculated at x
+ :rtype: numpy.ndarray (1D)
+ """
+ gr = self.get_growth_rate(x)
+ cm = x[self.var_mask_metab]
+ ce_scaled = self.pstm_scale * x[self.var_mask_enz]
+
+ metabolic_mb = self.metab_x_mras.dot(self.mras(x))
+ ps_drain = self.metab_x_psras.dot(ce_scaled)
+
+ if np.count_nonzero(self.enz_x_mras) == 0:
+ metab_mb = metabolic_mb + gr * (ps_drain - cm)
+ else:
+ enz_mb = self.enz_x_mras.dot(self.mras(x))
+ synth_enz_mb = - self.metab_x_psras.dot(self.pstm_scale * enz_mb)
+ metab_mb = metabolic_mb + gr * (ps_drain - cm) + synth_enz_mb
+ return metab_mb
def heq_mass_balance_jac(self, x):
- cm_var = x[:self.var_n_metabs]
- ce_all = np.append(x[self.var_n_metabs:], self.model_params['const_enz_conc'])
-
- metabolism_jac = self.model_params['S_scaled_mat'].dot(self.mrs_jac(x))
- degrad_costs_jac = self.model_params['P_scaled_mat'].dot(self.p_degrads_jac(x))
-
- # jacobian of dilution term: d/dx gr * (P.dot(ce) - cm)
- # jac = gr_jac * (P.dot(ce) - cm) + gr * (jac(P.dot(ce)) - jac(cm))
- pem = self.model_params['P_scaled_mat'].dot(ce_all) - cm_var
- # construct jacobian of dilution term (P.dot(ce) - cm):
- # d/dx (P.dot(ce) - cm) = d(P.dot(ce) - cm)/dx = P.dot(dce/dx) - dcm/dx
- cm_jac = -np.eye(self.var_n_metabs)
- pe_jac = self.model_params['P_scaled_mat'][:, :len(self.var_enz_sids)]
- pem_jac = np.hstack([cm_jac, pe_jac])
- dilution_jac = np.outer(pem, self.get_growth_rate_grad(x)) + self.get_growth_rate(x) * pem_jac
- metabs_mb_jac = metabolism_jac + dilution_jac + degrad_costs_jac
- return np.nan_to_num(metabs_mb_jac, nan=0.0)
-
- def heq_enz_conv_mb(self, x):
- return self.enz_conv_mb(x)
-
- def heq_enz_conv_mb_jac(self, x):
- return self.enz_conv_mb_jac(x)
+ """Calculate Jacobian of mass balance equations at x for metabolites
+
+ incuding protein degradation , which consume metabolites
+ Note: protein degradation is also part of metabolic reactions (mras), which return
+ metabolites.
+
+ metab_mb = dcm/dt = metabolic_mb + growth_dilution + sythesis_enzyme_balance
+ metabolic_mb = metab_x_mras.dot(mras(x))
+ growth_dilution = λ * (metab_x_psras.dot*(ce_scaled) - cm)
+ ce_scaled = pstm_scale * ce # convert concentration to amounts, considering stoichiometry
+ sythesis_enzyme_balance = - metab_x_psras.dot(pstm_scale * enzyme_balance)
+ enzyme_balance = enz_x_mras.dot(mras(x))
+
+ ∂metab_mb/∂x = metabolic_mb_jac + dilution_jac + sythesis_enzyme_balance_jac
+ metabolic_mb_jac = metab_x_mras.dot(mras_jac(x))
+ dilution_jac = gr_grad * (metab_x_psras.dot(ce_scaled) - cm)
+ + gr * (metab_x_psras.dot(ce_scaled_jac) - cm_jac)
+ ce_scaled = pstm_scale * ce # convert concentration to amounts, considering stoichiometry
+ ce_jac = np.hstack( (np.zeros((len(ce), len(cm))), np.eye(len(ce)) ))
+ cm_jac = np.hstack( (np.eye(len(cm)), np.zeros((len(cm), len(ce))) ))
+ ce_scaled_jac = pstm_scale.reshape((-1,1)) * ce_jac
+ sythesis_enzyme_balance_jac = - metab_x_psras.dot(pstm_scale * enzyme_balance_jac)
+ enzyme_balance_jac = enz_x_mras.dot(mras_jac(x))
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: Jacobian of mass balance equations at x
+ :rtype: numpy.ndarray (2D)
+ """
+ gr = self.get_growth_rate(x)
+ gr_grad = self.get_growth_rate_grad(x)
+
+ cm = x[self.var_mask_metab]
+ ce = x[self.var_mask_enz]
+ ce_scaled = self.pstm_scale * ce
+ cm_jac = np.hstack((np.eye(len(cm)), np.zeros((len(cm), len(x) - len(cm)))))
+ ce_jac = np.hstack((np.zeros((len(ce), len(x) - len(ce))), np.eye(len(ce))))
+ ce_jac_scaled = self.pstm_scale.reshape((-1, 1)) * ce_jac
+
+ metabolic_jac = self.metab_x_mras.dot(self.mras_jac(x))
+ ps_drain = self.metab_x_psras.dot(ce_scaled)
+ dilute_1_jac = np.outer((ps_drain - cm), gr_grad)
+ dilute_2_jac = gr * (self.metab_x_psras.dot(ce_jac_scaled) - cm_jac)
+
+ if np.count_nonzero(self.enz_x_mras) == 0:
+ metab_mb_jac = metabolic_jac + dilute_1_jac + dilute_2_jac
+ else:
+ enz_mb_jac = self.enz_x_mras.dot(self.mras_jac(x))
+ synth_enz_mb_jac = -self.metab_x_psras.dot(self.pstm_scale.reshape((-1, 1)) * enz_mb_jac)
+ metab_mb_jac = metabolic_jac + dilute_1_jac + dilute_2_jac + synth_enz_mb_jac
+ return np.nan_to_num(metab_mb_jac, nan=0.0)
+
+ def heq_mass_balance_hess(self, x):
+ """Calculate Hessians of mass balance equations at x
+
+ incuding protein degradation , which consume metabolites
+ Here: assuming protein degradation has stoichometry of -1 for relevant enzyme
+ Note: protein degradation is also part of metabolic reactions (mras), which return
+ metabolites.
+
+ metab_mb = dcm/dt = metabolic_mb + growth_dilution + sythesis_enzyme_balance
+ metabolic_mb = metab_x_mras.dot(mras(x))
+ growth_dilution = λ * (metab_x_psras.dot*(ce_scaled) - cm)
+ ce_scaled = pstm_scale * ce # convert concentration to amounts, considering stoichiometry
+ sythesis_enzyme_balance = - metab_x_psras.dot(pstm_scale * enzyme_balance)
+ enzyme_balance = enz_x_mras.dot(mras(x))
+
+ ∂metab_mb/∂x = metabolic_mb_jac + dilution_jac + sythesis_enzyme_balance_jac
+ metabolic_mb_jac = metab_x_mras.dot(mras_jac(x))
+ dilution_jac = gr_grad * (metab_x_psras.dot(ce_scaled) - cm)
+ + gr * (metab_x_psras.dot(ce_scaled_jac) - cm_jac)
+ ce_scaled = pstm_scale * ce # convert concentration to amounts, considering stoichiometry
+ ce_jac = np.hstack( (np.zeros((len(ce), len(cm))), np.eye(len(ce)) ))
+ cm_jac = np.hstack( (np.eye(len(cm)), np.zeros((len(cm), len(ce))) ))
+ ce_scaled_jac = pstm_scale.reshape((-1,1)) * ce_jac
+ sythesis_enzyme_balance_jac = - metab_x_psras.dot(pstm_scale * enzyme_balance_jac)
+ enzyme_balance_jac = enz_x_mras.dot(mras_jac(x))
+
+ array of n Hessians
+ metab_mb_hesss = metabolic_mb_hess + dilution_hess + sythesis_enzyme_balance_hess
+ metabolic_mb_hesss = metab_x_mras.dot(mras_hesss(x))
+ dilution_hessd = gr_hess * (metab_x_psras.dot(ce_scaled) - cm)
+ + gr_grad * (metab_x_psras.dot(ce_scaled_jac) - cm_jac)
+ + gr_grad * (metab_x_psras.dot(ce_scaled_jac) - cm_jac)
+ + gr * (metab_x_psras.dot(ce_scaled_hess) - cm_hess)
+ dilution_hesss = gr_hess * (metab_x_psras.dot(ce_scaled) - cm)
+ + 2 * gr_grad * (metab_x_psras.dot(ce_scaled_jac) - cm_jac)
+
+ ce_scaled = pstm_scale * ce # convert concentration to amounts, considering stoichiometry
+ ce_jac = np.hstack( (np.zeros((len(ce), len(cm))), np.eye(len(ce)) ))
+ cm_jac = np.hstack( (np.eye(len(cm)), np.zeros((len(cm), len(ce))) ))
+ ce_scaled_jac = pstm_scale.reshape((-1,1)) * ce_jac
+ cm_hess = np.zeros((x,x)) # as cm_jac only contains constants
+ ce_hess = np.zeros((x,x))
+
+ sythesis_enzyme_balance_hesss = - metab_x_psras.dot(pstm_scale * enzyme_balance_hesss)
+
+ enzyme_balance_jac = enz_x_mras.dot(mras_hesss(x))
+
+ sythesis_enzyme_balance_jac = - metab_x_psras.dot(pstm_scale * enzyme_balance_jac)
+ enzyme_balance_jac = enz_x_mras.dot(mras_jac(x))
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: Hessians of mass balance equations at x
+ :rtype: numpy.ndarray (3D) with shape [n,x,x]
+ """
+ gr_grad = self.get_growth_rate_grad(x)
+ gr_hess = self.get_growth_rate_hess(x)
+
+ cm = x[self.var_mask_metab]
+ ce = x[self.var_mask_enz]
+ ce_scaled = self.pstm_scale * ce
+ cm_jac = np.hstack((np.eye(len(cm)), np.zeros((len(cm), len(x) - len(cm)))))
+ ce_jac = np.hstack((np.zeros((len(ce), len(x) - len(ce))), np.eye(len(ce))))
+ ce_jac_scaled = self.pstm_scale.reshape((-1, 1)) * ce_jac
+
+ metabolic_hesss = hesss_from_matrix_hesss(self.metab_x_mras, self.mras_hesss(x))
+ ps_drain = self.metab_x_psras.dot(ce_scaled)
+ dilute_1_hesss = hesss_from_v_hess(ps_drain - cm, gr_hess)
+ dilute_2_hesss = hesss_from_grad_jac(gr_grad, self.metab_x_psras.dot(ce_jac_scaled) - cm_jac)
+
+ if np.count_nonzero(self.enz_x_mras) == 0:
+ metab_mb_hesss = metabolic_hesss + dilute_1_hesss + 2 * dilute_2_hesss
+ else:
+ enz_mb_hesss = hesss_from_matrix_hesss(self.enz_x_mras, self.mras_hesss(x))
+
+ synth_enz_mb_hesss = -hesss_from_matrix_hesss(self.metab_x_psras,
+ (self.pstm_scale ** 2).reshape((-1, 1, 1)) * enz_mb_hesss)
+ metab_mb_hesss = metabolic_hesss + dilute_1_hesss + 2 * dilute_2_hesss + synth_enz_mb_hesss
+ return np.nan_to_num(metab_mb_hesss, nan=0.0)
+
+ def heq_enz_trans(self, x):
+ return self.enz_trans(x)
+
+ def heq_enz_trans_jac(self, x):
+ return self.enz_trans_jac(x)
+
+ def heq_enz_trans_hess(self, x):
+ return self.enz_trans_hesss(x)
def heq_density(self, x):
- ce_var = x[self.var_n_metabs:]
+ """Calculate density (g/l) constraint at x
+
+ Could be extended having several density constraints
+
+ macromolecular density is total protein mass density
+ alternalively, dry mass density also adds total metabolite mass density
+ formulated in a way to
+
+ using ipopt, the constraint bounds could be set to rho (only add up masses)
+ density = sum(species_concentration * molecular_weight)
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: mass balance calculated at x
+ :rtype: numpy.ndarray (1D)
+ """
+ ce = x[self.var_mask_enz]
if self.model_params.get('macro', True) is True:
- density = self.model_params['weights'][self.var_n_metabs:].dot(ce_var)
+ density = self.model_params['mws'][self.var_mask_enz].dot(ce)
else:
- density = self.model_params['weights'].dot(x)
- density += self.model_params['const_enz_weights'].dot(
- self.model_params['const_enz_conc'])
- return self.model_params['rho'] - density
+ density = self.model_params['mws'].dot(x)
+ # !!! when using Ipopt, the rho could be specified as a bound!
+ return np.array(self.model_params['rho'] - density)
+ # noinspection PyUnusedLocal
def heq_density_jac(self, x):
+ """Calculate gradient of density constraint at x
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: mass balance jacobian at x
+ :rtype: numpy.ndarray (1D)
+ """
if self.model_params.get('macro', True) is True:
- density_grad = np.zeros(len(self.model_params['weights']))
- density_grad[self.var_n_metabs:] = \
- -self.model_params['weights'][self.var_n_metabs:]
+ density_grad = np.zeros(len(x))
+ density_grad[self.var_mask_enz] = -self.model_params['mws'][self.var_mask_enz]
else:
- density_grad = -self.model_params['weights']
+ density_grad = -self.model_params['mws']
return density_grad
+
+ # noinspection PyUnusedLocal
+ def heq_density_hess(self, x):
+ """Calculate hessian of density constraint at x
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: mass balance Hessian at x
+ :rtype: numpy.ndarray (2D), shape [x,x]
+ """
+ hess = np.zeros((self.n_vars, self.n_vars))
+ return np.tril(hess)
diff --git a/xbanalysis/utils/__init__.py b/xbanalysis/utils/__init__.py
index d75f80f..7a3e1e8 100644
--- a/xbanalysis/utils/__init__.py
+++ b/xbanalysis/utils/__init__.py
@@ -1,6 +1,5 @@
"""Subpackage with utility functions """
-#rom .xba_model import XBAmodel
-
-__all__ = []
+from .opt_results import OptResults
+__all__ = ['OptResults']
diff --git a/xbanalysis/utils/hessians.py b/xbanalysis/utils/hessians.py
new file mode 100644
index 0000000..9c63f41
--- /dev/null
+++ b/xbanalysis/utils/hessians.py
@@ -0,0 +1,125 @@
+"""Implementation of Hessians calculations.
+
+Peter Schubert, HHU Duesseldorf, May 2022
+"""
+
+import numpy as np
+
+
+def hess_from_grad_grad(grad_l, grad_r):
+ """Calculate Hessian from two gradiends wrt x
+
+ :param grad_l: left scalar gradient wrt x variables
+ :param grad_l: numpy.ndarray (1D) with length of x
+ :param grad_r: rigth scalar gradient wrt x variables
+ :param grad_r: numpy.ndarray (1D) with length of x
+ :returns: result_hess Hessian wrt to x variables, lower triangle matrix
+ :rtype: numpy.ndarray (2D) with shape [x,x]
+ """
+ assert grad_l.shape[0] == grad_r.shape[0]
+ size = grad_l.shape[0]
+ result_hess = np.zeros((size, size))
+ for i in range(size):
+ for j in range(i+1):
+ result_hess[i, j] = grad_l[i] * grad_r[j]
+ return result_hess
+
+
+def hess_from_v_hesss(v, hesss):
+ """Calculate Hessian from a vector and Hessians
+
+ :param v: vector of length n
+ :param v: numpy.ndarray (1D)
+ :param hesss: Hessians of n vectors wrt to x variables
+ :param hesss: numpy.ndarray (3D) with shape [n,x,x]
+ :returns: result_hess Hessian wrt to x variables, lower triangle matrix
+ :rtype: numpy.ndarray (2D) with shape [x,x]
+ """
+ assert hesss.shape[1] == hesss.shape[2]
+ size = hesss.shape[1]
+ result_hess = np.zeros((size, size))
+ for i in range(size):
+ for j in range(i + 1):
+ result_hess[i, j] = v.dot(hesss[:, i, j])
+ return result_hess
+
+
+def hess_from_jac_jac(jac_l, jac_r):
+ """Calculate Hessian from two Jacobians
+
+ :param jac_l: left Jacobian for vector of size n wrt to x variables
+ :param jac_l: numpy.ndarray (2D) with shape [n,x]
+ :param jac_r: right Jacobian for vector of size n wrt to x variables
+ :param jac_r: numpy.ndarray (2D) with shape [n,x]
+ :returns: result_hess Hessian wrt to x variables, lower triangle matrix
+ :rtype: numpy.ndarray (2D) with shape [x,x]
+ """
+ assert jac_l.shape == jac_r.shape
+ size = jac_l.shape[1]
+ result_hess = np.zeros((size, size))
+ for i in range(size):
+ for j in range(i + 1):
+ result_hess[i, j] = jac_l[:, i].dot(jac_r[:, j])
+ return result_hess
+
+
+def hesss_from_matrix_hesss(matrix, hesss):
+ """Calculate dot product of a stoichiometric matrix with Hessians
+
+ :param matrix: stoichiometric sub-matrix with shape (m, n)
+ :param matrix: numpy.ndarray (2D) with shape [m,n]
+ :param hesss: Hessians of n vectors wrt to x variables
+ :param hesss: numpy.ndarray (3D) with shape [n,x,x]
+ :returns: result_hess Hessians wrt to x variables, lower triangle matrix
+ :rtype: numpy.ndarray (3D) with shape [n,x,x]
+ """
+ assert matrix.shape[1] == hesss.shape[0], "expect as many Hessians as there are matrix columns"
+ assert hesss.shape[1] == hesss.shape[2]
+ size = hesss.shape[1]
+ n_rows, n_cols = matrix.shape
+ result_hesss = np.zeros((n_rows, size, size))
+ for i in range(n_rows):
+ for j in range(n_cols):
+ if matrix[i, j] != 0:
+ result_hesss[i] += matrix[i, j] * hesss[j]
+ return result_hesss
+
+
+def hesss_from_v_hess(v, hess):
+ """Calculate Hessians from vector and a Hessian
+
+ :param v: vector of length m
+ :param v: numpy.ndarray (1D)
+ :param hess: Hessian of scalar wrt to x variables
+ :param hess: numpy.ndarray (2D) with shape [x,x]
+ :returns: result_hesss Hessians wrt to x variables, lower triangle matrix
+ :rtype: numpy.ndarray (3D) with shape [n,x,x]
+ """
+ assert hess.shape[0] == hess.shape[1]
+ size = hess.shape[1]
+ n_rows = v.shape[0]
+
+ result_hesss = np.zeros((n_rows, size, size))
+ for i in range(n_rows):
+ result_hesss[i] = v[i] * hess
+ return result_hesss
+
+
+def hesss_from_grad_jac(grad, jac):
+ """Calculate Hessians from gradient and a Jacobian
+
+ :param grad: gradient of scalar wrt to x variables
+ :param grad: numpy.ndarray (1D) with shape [x]
+ :param jac: Jacobian of a vector size m wrt to x variables
+ :param jac: numpy.ndarray (2D) with shape [m,x]
+ :returns: result_hesss Hessians wrt to x variables, lower triangle matrix
+ :rtype: numpy.ndarray (3D) with shape [m,x,x]
+ """
+ assert grad.shape[0] == jac.shape[1]
+ size = jac.shape[1]
+ n_rows = jac.shape[0]
+
+ result_hesss = np.zeros((n_rows, size, size))
+ for i in range(n_rows):
+ result_hesss[i] = hess_from_grad_grad(jac[i], grad)
+ return result_hesss
diff --git a/xbanalysis/utils/opt_results.py b/xbanalysis/utils/opt_results.py
new file mode 100644
index 0000000..df73f03
--- /dev/null
+++ b/xbanalysis/utils/opt_results.py
@@ -0,0 +1,66 @@
+"""Implementation of OptResults class.
+
+Peter Schubert, HHU Duesseldorf, June 2022
+"""
+
+import numpy as np
+
+
+class OptResults:
+
+ def __init__(self, problem, n_points):
+ """Initialize OptResults.
+
+ :param problem: GbaProblem from where to extract relevant parameters
+ :type problem: GbaProblem
+ :param n_points: number of points to record
+ :type n_points: integer
+ """
+ self.problem = problem
+ self._n_points = n_points
+ self._idx = 0
+ self._n_m = len(self.problem.var_metab_sids)
+ self._n_e = len(self.problem.var_enz_sids)
+
+ self.ress = {}
+ self.ress['metabolites'] = self.problem.var_metab_sids
+ self.ress['enzymes'] = self.problem.var_enz_sids
+ self.ress['mr_reactions'] = self.problem.mr_rids
+ self.ress['enz_sid2idx'] = {sid: idx for idx, sid in enumerate(self.problem.var_enz_sids)}
+ self.ress['var_mask_metab'] = self.problem.var_mask_metab
+ self.ress['var_mask_enz'] = self.problem.var_mask_enz
+
+ self.ress['grs_per_h'] = np.zeros(n_points)
+ self.ress['status'] = np.zeros(n_points)
+ self.ress['x_opt'] = np.zeros((n_points, self.problem.n_vars))
+ self.ress['cm_M'] = np.zeros((n_points, self._n_m))
+ self.ress['rho_m'] = np.zeros((n_points, self._n_m))
+ self.ress['ce_M'] = np.zeros((n_points, self._n_e))
+ self.ress['rho_e'] = np.zeros((n_points, self._n_e))
+ self.ress['alphas'] = np.zeros((n_points, self._n_e))
+ self.ress['mr_fluxes'] = np.zeros((n_points, len(self.problem.mr_rids)))
+ self.ress['psm_fluxes'] = np.zeros((n_points, self._n_e))
+
+ def add(self, info):
+ if self._idx < self._n_points:
+ i = self._idx
+ cx = info['x']
+ self.ress['x_opt'][i] = cx
+ self.ress['status'][i] = info['status']
+ self.ress['grs_per_h'][i] = -3600 * info['obj_val']
+ self.ress['cm_M'][i] = cx[self.problem.var_mask_metab]
+ self.ress['ce_M'][i] = cx[self.problem.var_mask_enz]
+ rho = self.problem.model_params['mws'] * cx
+ self.ress['rho_m'][i] = rho[self.problem.var_mask_metab]
+ self.ress['rho_e'][i] = rho[self.problem.var_mask_enz]
+
+ # !!! what volume to consider, possibly the volume of the ribosome ??
+ scale_factor = 3600 * 1e3 / self.problem.model_params['rho']
+ self.ress['mr_fluxes'][i] = scale_factor * self.problem.mras(cx) / self.problem.var_vols[-1]
+ self.ress['psm_fluxes'][i] = scale_factor / (self.problem.pstm_scale * self.problem.pstmas(cx))
+
+ self.ress['alphas'][i] = self.problem.get_alphas(cx)
+ self._idx += 1
+
+ def get(self):
+ return self.ress
--
GitLab