diff --git a/setup.py b/setup.py
index 2c04d6734c30e6928996e9ff4f15d7dc5485d051..fa9a43be13f533c5d79c8105852341927632c2ec 100755
--- a/setup.py
+++ b/setup.py
@@ -38,6 +38,7 @@ setup(
install_requires=['pandas>=1.3.0',
'numpy >= 1.20.0',
'sympy >= 1.9.0',
+ 'scipy >= 1.7.0',
'sbmlxdf>=0.2.5'],
python_requires=">=3.7",
keywords=['modeling', 'standardization', 'SBML'],
diff --git a/xbanalysis/_version.py b/xbanalysis/_version.py
index 120eb33ded3939e49674a86660c986ede1a23ee9..88646c942b503d2d5b6382a0b91033a5ed23f740 100644
--- a/xbanalysis/_version.py
+++ b/xbanalysis/_version.py
@@ -1,5 +1,5 @@
"""
Definition of version string.
"""
-__version__ = "0.3.0"
+__version__ = "0.3.1"
program_name = 'xbanalysis'
diff --git a/xbanalysis/model/xba_model.py b/xbanalysis/model/xba_model.py
index 6e127885b53eb4da529d343648d05644401e0b7a..c72928d636f154d2abdd3d34722af143524b9f38 100644
--- a/xbanalysis/model/xba_model.py
+++ b/xbanalysis/model/xba_model.py
@@ -35,7 +35,7 @@ class XbaModel:
for fid, row in model_dict['funcDefs'].iterrows()}
self.species = {sid: XbaSpecies(row)
for sid, row in model_dict['species'].iterrows()}
- self.reactions = {rid: XbaReaction(row, self.functions)
+ self.reactions = {rid: XbaReaction(row, self.functions, self.compartments)
for rid, row in model_dict['reactions'].iterrows()}
for r in self.reactions.values():
diff --git a/xbanalysis/model/xba_reaction.py b/xbanalysis/model/xba_reaction.py
index 985af2609255305bc7a13ce10bc766b3d4a6894c..94b3f966bbf8549e9eb88796e34052ac5d880844 100644
--- a/xbanalysis/model/xba_reaction.py
+++ b/xbanalysis/model/xba_reaction.py
@@ -10,7 +10,7 @@ from xbanalysis.utils.utils import get_species_stoic, get_local_params, expand_k
class XbaReaction:
- def __init__(self, s_reaction, functions):
+ def __init__(self, s_reaction, functions, compartments):
self.id = s_reaction.name
self.name = s_reaction.get('name', '')
self.sboterm = s_reaction['sboterm']
@@ -26,6 +26,9 @@ class XbaReaction:
for pid, val in self.local_params.items():
expanded_kl = re.sub(r'\b' + pid + r'\b', str(val), expanded_kl)
self.expanded_kl = expand_kineticlaw(expanded_kl, functions)
+ for cid in compartments.keys():
+ if re.search(r'\b'+cid+r'\b', self.expanded_kl) is not None:
+ self.compartment = cid
self.ps_product = ''
self.enzyme = ''
@@ -37,40 +40,33 @@ class XbaReaction:
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 !!!!!!!!!!
+ # !!! we could also enter the Reaction Object references instead
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():
species[sid].add_m_reaction(self.id, -val)
- # protein degradation reaction, assume that first protein in reactants is getting degraded
- # a protein degradation reaction is still a metabolic reaction
+ # protein degradation reactions
if self.sboterm == 'SBO:0000179':
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 transitions
if self.sboterm == 'SBO:0000182':
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, 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.
+ # protein synthesis reactions
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():
- species[sid].add_ps_reaction(self.id, val)
+ species[sid].set_ps_rid(self.id)
+ else:
+ species[sid].add_ps_reaction(self.id, val)
for sid, val in self.reactants.items():
species[sid].add_ps_reaction(self.id, -val)
diff --git a/xbanalysis/model/xba_species.py b/xbanalysis/model/xba_species.py
index cf94b5735abe0cd44c100999829998579551f527..436537012955b0e481588c15d52e226482d3a1ea 100644
--- a/xbanalysis/model/xba_species.py
+++ b/xbanalysis/model/xba_species.py
@@ -25,8 +25,8 @@ class XbaSpecies:
self.p_degrad_rids = {}
self.p_trans_rids = {}
- def set_ps_rid(self, rid, stoic):
- self.ps_rid = {'rid': rid, 'stoic': stoic}
+ def set_ps_rid(self, rid):
+ self.ps_rid = rid
def add_m_reaction(self, rid, stoic):
self.m_reactions[rid] = stoic
diff --git a/xbanalysis/problems/gba_problem.py b/xbanalysis/problems/gba_problem.py
index eaeb62b36f17d1342f30a1032cfd481b06685337..0eccbb6e93ab4464aa8b99b1404d0f7479abaf55 100644
--- a/xbanalysis/problems/gba_problem.py
+++ b/xbanalysis/problems/gba_problem.py
@@ -7,25 +7,32 @@ import numpy as np
import re
import types
import sympy
+import scipy
-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)
+from xbanalysis.utils.utils import progress_bar
# Decorator to cache last result
def cache_last(func):
def wrapper(self, *args):
fname = wrapper.__name__
- if fname not in GbaProblem.last_values or np.array_equal(GbaProblem.last_values[fname][0], args[0]) is False:
+ if (fname not in GbaProblem.last_values or
+ np.array_equal(GbaProblem.last_values[fname][0], args[0]) is False or
+ GbaProblem.cache is False):
GbaProblem.last_values[fname] = [args[0], func(self, args[0])]
return GbaProblem.last_values[fname][1]
wrapper.__name__ = func.__name__
return wrapper
+# TODO: cache last results from rate functions, Jacobians and Hessians
+# TODO: call get_reaction_rates from solver ?
+# TODO: split get_reaction_rates, get_reaction_jacs, get_reaction_hesses
+# TODO: number of mb equations, enzyme transitions density constraints
+
class GbaProblem:
+ cache = True
last_values = {}
def __init__(self, xba_model):
@@ -36,8 +43,7 @@ class GbaProblem:
"""
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 in ('SBO:0000250', 'SBO:0000252')]
@@ -58,36 +64,40 @@ class GbaProblem:
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)}
- # TODO: accept child SBO terms (i.e. SBO hierarchy parsing)
+ # protein synthesis reactions get ordered according to the order of their enzyme products
+ # to facilitate matrix multiplication without further remapping
+ # - here we expect a one-to-one relationship between enzyme and corresponding protein synthesis reaction
self.mr_rids = [r.id for r in xba_model.reactions.values()
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 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])
+ self.ps_rids = [xba_model.species[sid].ps_rid for sid in self.var_enz_sids]
+ self.mr_vols = np.array([xba_model.compartments[xba_model.reactions[rid].compartment].size
+ for rid in self.mr_rids])
+ self.ps_vols = np.array([xba_model.compartments[xba_model.reactions[rid].compartment].size
+ for rid in self.ps_rids])
+ # calculate stoichiometric sub-matrices and scale them according to volume
+ # TODO: sparce matrix storage and calculations
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.metab_x_psrs = self.get_stoic_matrix(self.var_metab_sids, self.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]
+ self.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
+ self.metab_x_psras = self.metab_x_psrs / metab_vols.reshape((-1, 1))
+ self.enz_x_mras = self.enz_x_mrs / self.enz_vols.reshape((-1, 1))
+
+ # identify reactions ids for enzyme degradation and corresponding submatrix
+ enz_degrad_idxs = self.get_export_idxs(self.enz_x_mrs)
+ enz_degrad_rids = np.array(self.mr_rids)[enz_degrad_idxs]
+ self.enz_x_degrads = self.enz_x_mrs[:, enz_degrad_idxs]
+
+ # identify reactions ids for enzyme transitions and corresponding submatrix
+ # for submatrix perform LU decomposition and remove empty rows from U matrix
+ enz_trans_idxs = self.get_trans_idxs(self.enz_x_mras)
+ enz_trans_rids = np.array(self.mr_rids)[enz_trans_idxs]
+ u = scipy.linalg.lu((self.enz_x_mras[:, enz_trans_idxs]))[2]
+ self.subenz_x_trans = u[abs(u).sum(axis=1) > 0]
# model parameters
self.model_params = {pid: p.value for pid, p in xba_model.parameters.items()}
@@ -99,9 +109,10 @@ class GbaProblem:
self.model_params['mws'] = np.array([xba_model.species[sid].mw for sid in self.var_sids])
# create functions
- 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()
+ self.mras, self.mras_jac, self.mras_hesss = self.get_reactions(self.mr_rids, inv=False)
+ self.pstmas, self.pstmas_jac, self.pstmas_hesss = self.get_reactions(self.ps_rids, inv=True)
+ self.degradas, self.degradas_jac, self.degradas_hesss = self.get_reactions(enz_degrad_rids, inv=False)
+ self.enz_transas, self.enz_transas_jac, self.enz_transas_hesss = self.get_reactions(enz_trans_rids, inv=False)
def get_stoic_matrix(self, sids, rids):
"""retrieve stoichiometric sub-matrix [sids x rids].
@@ -128,12 +139,14 @@ class GbaProblem:
self.model_params[key] = val
def _to_vectors(self, rs_str):
- """replace free and constant metabolites by vector elements.
+ """replace variables and constant metabolites by vector elements.
+
+ Variables converted to x[], external metabolites to ext_conc[].
- :param rs_str: reaction string (extracted from the xba model).
+ :param rs_str: reaction strings (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
+ :rtype: list of reaction stings with parameters replaced by vector elements
"""
rs_v_str = []
for r_v_str in rs_str:
@@ -163,25 +176,26 @@ class GbaProblem:
partials_x = []
hess = []
for idx1 in range(self.n_vars):
- partials_x.append(str(func_sym.diff(x[idx1])))
+ partial_xsp = func_sym.diff(x[idx1])
+ partials_x.append(str(partial_xsp))
partials_xx = []
for idx2 in range(self.n_vars):
if idx2 <= idx1:
- partials_xx.append(str(func_sym.diff(x[idx1], x[idx2])))
+ partials_xx.append(str(partial_xsp.diff(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))
+ progress_bar(i, len(rs_v_str))
jac = compile('def jac(x): return np.array([' + ', '.join(grads_str) + '])',
'<string>', 'exec')
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
+ def get_reactions(self, rids, inv=False):
+ """create function objects for reactions rates/times, Jacobian and Hessian
extracted from kinetic laws for SBML model
@@ -190,6 +204,8 @@ class GbaProblem:
:param rids: reaction ids for reactions in amount per time to be processed.
:type rids: list of strings
+ :param inv: flag indicating if inverse of reaction required (i.e. reaction times)
+ :type inv: boolean (default: False)
:returns: ras
:rtype: function object, returning an 1D array of function results
:returns: ras_jac
@@ -198,9 +214,12 @@ class GbaProblem:
: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]
+ if inv is False:
+ kinetics_str = [self.xba_model.reactions[rid].expanded_kl for rid in rids]
+ else:
+ kinetics_str = [f'1.0/({self.xba_model.reactions[rid].expanded_kl})' for rid in rids]
v_str = self._to_vectors(kinetics_str)
- func_code = compile('def rs_apt(x): return np.array([' + ', '.join(v_str) + '])',
+ func_code = compile('def reactions(x): return np.array([' + ', '.join(v_str) + '])',
'<string>', 'exec')
jac_code, hesss_code = self._gradients(v_str)
ras = types.FunctionType(func_code.co_consts[0], self.model_params)
@@ -208,129 +227,74 @@ class GbaProblem:
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
+ @staticmethod
+ def get_trans_idxs(smat):
+ """Identify metabolic reactions in S matrix with species transitions.
- Maximum protein synthesis rates from model converted to minimum protein synthesis times
+ I.e. columns containing both positive and negative entries
- 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
+ :param smat:
+ :type smat:
+ :return: bolean vector identifying relevant columns
+ :rtype: 1D ndarray with dtype bool
"""
- 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')
- 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
+ mask_trans_idxs = np.zeros(smat.shape[1])
+ for col_idx in range(smat.shape[1]):
+ if ((np.count_nonzero(smat[:, col_idx] < 0) > 0) and
+ np.count_nonzero(smat[:, col_idx] > 0) > 0):
+ mask_trans_idxs[col_idx] = 1.0
+ return mask_trans_idxs > 0
+
+ @staticmethod
+ def get_export_idxs(smat):
+ """Identify metabolic reactions in S matrix with consumption only.
+
+ I.e. columns containing only negative entries
+
+ :param smat:
+ :type smat:
+ :return: bolean vector identifying relevant columns
+ :rtype: 1D ndarray with dtype bool
"""
- # mass balance for interconvering enzymes
- # 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.var_enz_sids:
- sp = self.xba_model.species[enz_sid]
- partner_enz_sids = {enz_sid}
- 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 p_trans_r.products.items():
- sp_other = self.xba_model.species[other_sid]
- 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 p_trans_r.reactants.items():
- sp_other = self.xba_model.species[other_sid]
- 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:
- # consider unique interconversion sets
- if partner_enz_sids not in enz_in_mrs_sids:
- enz_in_mrs_sids.append(partner_enz_sids)
-
- # for each set interconverting enzymes create a mass balance equation
- conversions_str = []
- for partner_enz_sids in enz_in_mrs_sids:
- enz_sid = partner_enz_sids.pop()
- sp = self.xba_model.species[enz_sid]
- parts_str = []
- 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)
- if len(v_str) == 0:
- v_str = '0'
- func_code = compile('def enz_trans(x): return np.array([' + ', '.join(v_str) + '])',
- '<string>', 'exec')
- 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
+ mask_remove_idxs = np.zeros(smat.shape[1])
+ for col_idx in range(smat.shape[1]):
+ if ((np.count_nonzero(smat[:, col_idx] < 0) > 0) and
+ np.count_nonzero(smat[:, col_idx] > 0) == 0):
+ mask_remove_idxs[col_idx] = 1.0
+ return mask_remove_idxs > 0
@cache_last
def get_growth_rate_0(self, x):
"""Calculate growth rate at x, assuming no protein degradation
- considering stoichiometric coefficiens for enzyme turnover
+ Equations:
+ gr0(x) = 1/tau(x)
+ tau(x) = pstms(x).dot(ce)
+ pstms(x) = enz_vols * ptsmas(x)
: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))
+ pstms = self.enz_vols * 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
@cache_last
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))
+ Equations:
+ gr0(x) = 1/tau(x)
+ tau(x) = pstms(x).dot(ce)
+ pstms(x) = enz_vols * ptsmas(x)
+
+ gr0_grad(x) = -tau_grad(x)/tau(x)^2 = -gr0(x)^2 * tau_grad(x)
+ tau_grad(x) = ce.dot(pstms_jac(x)) + pstms(x).dot(ce_jac)
+ pstms_jac(x) = enz_vols.reshape(-1,1) * pstmas(x)
:param x: optimization variables
:type x: numpy.ndarray
@@ -340,28 +304,29 @@ class GbaProblem:
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
+ pstms = self.enz_vols * self.pstmas(x)
+ pstms_jac = self.enz_vols.reshape((-1, 1)) * self.pstmas_jac(x)
+
+ tau_grad = ce.dot(pstms_jac) + pstms.dot(ce_jac)
+ gr0_grad = - gr0**2 * tau_grad
return np.nan_to_num(gr0_grad, nan=0.0)
@cache_last
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))
+ Equations:
+ gr0(x) = 1/tau(x)
+ tau(x) = pstms(x).dot(ce)
+ pstms(x) = enz_vols * ptsmas(x)
- gr0_grad = -gr0**2 * tau_grad
- tau_grad = ce.dot(pstms_jac(x)) + pstms(x).dot(ce_jac)
+ gr0_grad(x) = -tau_grad(x)/tau(x)^2 = -gr0(x)^2 * tau_grad(x)
+ tau_grad(x) = ce.dot(pstms_jac(x)) + pstms(x).dot(ce_jac)
+ pstms_jac(x) = enz_vols.reshape(-1,1) * pstmas(x)
- 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
+ gr0_hess(x) = -2*gr0(x) * gr0_grad(x) * tau_grad(x) - gr0(x)^2 * tau_hess(x)
+ tau_hess = ce_jac.dot(pstms_jac(x)) + ce.dot(pstms_hesss(x)) + pstms_jac(x).dot(ce_jac)
+ ce_hess = 0
:param x: optimization variables
:type x: numpy.ndarray
@@ -370,26 +335,19 @@ class GbaProblem:
"""
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)
+ pstms = self.enz_vols * self.pstmas(x)
+ pstms_jac = self.enz_vols.reshape((-1, 1)) * self.pstmas_jac(x)
+ pstms_hesss = self.enz_vols.reshape((-1, 1, 1)) * self.pstmas_hesss(x)
+
+ tau_grad = ce.dot(pstms_jac) + pstms.dot(ce_jac)
+ hess_a = np.outer(gr0_grad, tau_grad)
+ hess_b1 = (ce.reshape((-1, 1, 1))*pstms_hesss).sum(axis=0)
+ hess_b2 = ce_jac.T.dot(pstms_jac)
+ hess_b3 = pstms_jac.T.dot(ce_jac)
+ gr0_hess = - 2 * gr0 * hess_a - gr0**2 * (hess_b1 + hess_b2 + hess_b3)
+ return np.nan_to_num(np.tril(gr0_hess), nan=0.0)
@cache_last
def get_growth_rate(self, x):
@@ -397,10 +355,11 @@ class GbaProblem:
considering stoichiometric coefficiens unequal unity for enzyme turnover
stoichiometric sub-matrices do not require volume scaling (as Vol_enz gets factored out)
+ With dilution (negative stoichiometric coefficients in enz_x_mrs), growth rate reduces
- 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)
+ Equations:
+ gr(x) = gr0(x) * (1 + enz_to_mras(x).dot(pstmas(x)))
+ enz_to_mras(x) = enz_x_mrs.dot(mras(x))
:param x: optimization variables
:type x: numpy.ndarray
@@ -411,20 +370,23 @@ class GbaProblem:
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))
+ pstmas = self.pstmas(x)
+ enz_to_mras = self.enz_x_degrads.dot(self.degradas(x))
+
+ gr = gr0 * (1.0 + enz_to_mras.dot(pstmas))
return gr
@cache_last
def get_growth_rate_grad(self, x):
"""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))
+ Equations:
+ gr(x) = gr0(x) * (1 + enz_to_mras(x).dot(pstmas(x)))
+ enz_to_mras(x) = enz_x_mrs.dot(mras(x))
+
+ gr_grad(x) = gr0_grad(x) * (1 + enz_to_mras(x).dot(pstmas(x)))
+ + gr0(x) * (pstmas(x).dot(enz_to_mras_jac(x)) + enz_to_mras(x).dot(pstmas_jac(x)))
+ enz_to_mras_jac(x) = enz_to_mrs.dot(mras_jac(x))
:param x: optimization variables
:type x: numpy.ndarray
@@ -436,30 +398,35 @@ class GbaProblem:
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)))
+ pstmas = self.pstmas(x)
+ pstmas_jac = self.pstmas_jac(x)
+ enz_to_mras = self.enz_x_degrads.dot(self.degradas(x))
+ enz_to_mras_jac = self.enz_x_degrads.dot(self.degradas_jac(x))
+
+ gr_grad = (gr0_grad * (1.0 + enz_to_mras.dot(pstmas))
+ + gr0 * (pstmas.dot(enz_to_mras_jac) + enz_to_mras.dot(pstmas_jac)))
return np.nan_to_num(gr_grad, nan=0.0)
@cache_last
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)
+ Equations:
+ gr(x) = gr0(x) * (1 + enz_to_mras(x).dot(pstmas(x)))
+ enz_to_mras(x) = enz_x_mrs.dot(mras(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_grad(x) = gr0_grad(x) * (1 + enz_to_mras(x).dot(pstmas(x)))
+ + gr0(x) * (pstmas(x).dot(enz_to_mras_jac(x)) + enz_to_mras(x).dot(pstmas_jac(x)))
+ enz_to_mras_jac(x) = enz_to_mrs.dot(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) )
+ gr_hess(x) = gr0_hess(x) * (1 + enz_to_mras(x).dot(pstmas(x)))
+ + (pstmas(x).dot(enz_to_mras_jac(x)) + enz_to_mras(x).dot(pstmas_jac(x))) * gr0_grad(x)
+ + gr0_grad(x) * (pstmas(x).dot(enz_to_mras_jac(x)) + enz_to_mras(x).dot(pstmas_jac(x)))
+ + gr0(x) * (pstmas_jac(x).dot(enz_to_mras_jac(x))
+ + enz_to_mras_hess(x).dot(pstmas(x))
+ + enz_to_mras_jac(x).dot(pstmas_jac(x))
+ + pstmas_hess(x).dot(enz_to_mras(x))
+ )
:param x: optimization variables
:type x: numpy.ndarray
@@ -472,28 +439,24 @@ class GbaProblem:
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)
+ pstmas = self.pstmas(x)
+ pstmas_jac = self.pstmas_jac(x)
+ pstmas_hesss = self.pstmas_hesss(x)
+ enz_to_mras = self.enz_x_degrads.dot(self.degradas(x))
+ enz_to_mras_jac = self.enz_x_degrads.dot(self.degradas_jac(x))
+ enz_to_mras_hesss = np.array([(row.reshape((-1, 1, 1)) * self.degradas_hesss(x)).sum(axis=0)
+ for row in self.enz_x_degrads])
+
+ hess_a = gr0_hess * (1.0 + enz_to_mras.dot(pstmas))
+ grad = pstmas.dot(enz_to_mras_jac) + enz_to_mras.dot(pstmas_jac)
+ hess_b1 = np.outer(gr0_grad, grad)
+ hess_b2 = np.outer(grad, gr0_grad)
+ hess_c1 = pstmas_jac.T.dot(enz_to_mras_jac)
+ hess_c2 = (pstmas.reshape((-1, 1, 1)) * enz_to_mras_hesss).sum(axis=0)
+ hess_c3 = enz_to_mras_jac.T.dot(pstmas_jac)
+ hess_c4 = (enz_to_mras.reshape((-1, 1, 1)) * pstmas_hesss).sum(axis=0)
+ gr_hess = hess_a + hess_b1 + hess_b2 + gr0 * (hess_c1 + hess_c2 + hess_c3 + hess_c4)
+ return np.nan_to_num(np.tril(gr_hess), nan=0.0)
@cache_last
def get_alphas(self, x):
@@ -509,50 +472,46 @@ class GbaProblem:
"""
ce = x[self.var_mask_enz]
gr = self.get_growth_rate(x)
+ pstms = self.enz_vols * self.pstmas(x)
- enz_to_pstms = self.pstma_scale * self.pstmas(x)
- alphas0 = gr * ce * enz_to_pstms
+ alphas0 = gr * ce * 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
+ enz_to_mras = self.enz_x_degrads.dot(self.degradas(x))
+ alphas = alphas0 - enz_to_mras * self.pstmas(x)
return alphas
@cache_last
def heq_mass_balance(self, x):
"""Calculate mass balance equations at x for metabolites
- incuding protein degradation , which consume 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)
+ growth_dilution = λ * (metab_x_psrs.dot(ce) - cm)
+ ce_scaled = enz_vols * ce # convert concentration to amounts, considering stoichiometry
+ sythesis_enzyme_balance = - metab_x_psras.dot(enz_vols * 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)
+ :rtype: numpy.ndarray (1D) with shape (len(var_mask_metab), 1)
"""
gr = self.get_growth_rate(x)
cm = x[self.var_mask_metab]
- ce_scaled = self.pstm_scale * x[self.var_mask_enz]
+ ce = 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:
+ dilution = self.metab_x_psrs.dot(ce) - cm
+ metab_mb = metabolic_mb + gr * dilution
+ if np.count_nonzero(self.enz_x_mras) > 0:
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
+ metab_mb -= self.metab_x_psrs.dot(enz_mb)
return metab_mb
@cache_last
@@ -566,19 +525,19 @@ class GbaProblem:
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)
+ ce_scaled = enz_vols * ce # convert concentration to amounts, considering stoichiometry
+ sythesis_enzyme_balance = - metab_x_psras.dot(enz_vols * 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_scaled = enz_vols * 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)
+ ce_scaled_jac = enz_vols.reshape((-1,1)) * ce_jac
+ sythesis_enzyme_balance_jac = - metab_x_psras.dot(enz_vols * enzyme_balance_jac)
enzyme_balance_jac = enz_x_mras.dot(mras_jac(x))
:param x: optimization variables
@@ -588,25 +547,21 @@ class GbaProblem:
"""
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
+ mras_jac = self.mras_jac(x)
- 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)
+ metabolic_jac = self.metab_x_mras.dot(mras_jac)
+ dilution = self.metab_x_psrs.dot(ce) - cm
+ dilute_1_jac = np.outer(dilution, gr_grad)
+ dilute_2_jac = gr * (self.metab_x_psrs.dot(ce_jac) - 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
+ metab_mb_jac = metabolic_jac + dilute_1_jac + dilute_2_jac
+ if np.count_nonzero(self.enz_x_mras) > 0:
+ enz_mb_jac = self.enz_x_mras.dot(mras_jac)
+ metab_mb_jac -= self.metab_x_psrs.dot(enz_mb_jac)
return np.nan_to_num(metab_mb_jac, nan=0.0)
@cache_last
@@ -618,24 +573,6 @@ class GbaProblem:
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))
@@ -646,18 +583,18 @@ class GbaProblem:
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_scaled = enz_vols * 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
+ ce_scaled_jac = enz_vols.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)
+ sythesis_enzyme_balance_hesss = - metab_x_psras.dot(enz_vols * 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)
+ sythesis_enzyme_balance_jac = - metab_x_psras.dot(enz_vols * enzyme_balance_jac)
enzyme_balance_jac = enz_x_mras.dot(mras_jac(x))
:param x: optimization variables
@@ -667,50 +604,71 @@ class GbaProblem:
"""
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)
+ mras_hesss = self.mras_hesss(x)
+ metabolic_hesss = np.array([(row.reshape((-1, 1, 1)) * mras_hesss).sum(axis=0)
+ for row in self.metab_x_mras])
+ dilution = self.metab_x_psrs.dot(ce) - cm
+ dilution_jac = self.metab_x_psrs.dot(ce_jac) - cm_jac
+
+ dilution_1a_hesss = dilution.reshape((-1, 1, 1)) * gr_hess
+ dilution_1b_hesss = np.array([np.outer(row, gr_grad) for row in dilution_jac])
+ dilution_2a_hesss = np.array([np.outer(gr_grad, row) for row in dilution_jac])
+ # dilution_2b_hesss = np.zeros_like(metabolic_hesss)
+
+ metab_mb_hesss = metabolic_hesss + dilution_1a_hesss + dilution_1b_hesss + dilution_2a_hesss
+ if np.count_nonzero(self.enz_x_mras) > 0:
+ enz_mb_hesss = np.array([(row.reshape((-1, 1, 1)) * mras_hesss).sum(axis=0)
+ for row in self.enz_x_mras])
+ metab_mb_hesss -= np.array([(row.reshape((-1, 1, 1)) * enz_mb_hesss).sum(axis=0)
+ for row in self.metab_x_psrs])
+ return np.nan_to_num(np.tril(metab_mb_hesss), nan=0.0)
@cache_last
def heq_enz_trans(self, x):
- return self.enz_trans(x)
+ """Calculate enzyme transition equations at x
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: Jacobian of mass balance equations at x
+ :rtype: numpy.ndarray (1D) with shape (len(subenz_x_trans.shape[1]), 1)
+ """
+ return self.subenz_x_trans.dot(self.enz_transas(x))
@cache_last
def heq_enz_trans_jac(self, x):
- return self.enz_trans_jac(x)
+ """Calculate Jacobian for enzyme transition equations at x
+
+ :param x: optimization variables
+ :type x: numpy.ndarray
+ :returns: Jacobian of mass balance equations at x
+ :rtype: numpy.ndarray (2D)
+ """
+ return self.subenz_x_trans.dot(self.enz_transas_jac(x))
@cache_last
def heq_enz_trans_hess(self, x):
- return self.enz_trans_hesss(x)
+ """Calculate Hessians of enzyme transition equations at 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]
+ """
+ return np.array([(row.reshape((-1, 1, 1)) * self.enz_transas_hesss(x)).sum(axis=0)
+ for row in self.subenz_x_trans])
@cache_last
def heq_density(self, x):
"""Calculate density (g/l) constraint at x
- Could be extended having several density constraints
+ Could be extended having several density constraints (not yet implemented)
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)
@@ -718,20 +676,23 @@ class GbaProblem:
:param x: optimization variables
:type x: numpy.ndarray
:returns: mass balance calculated at x
- :rtype: numpy.ndarray (1D)
+ :rtype: numpy.ndarray (1D) with shape (1,) for now
"""
ce = x[self.var_mask_enz]
if self.model_params.get('macro', True) is True:
density = self.model_params['mws'][self.var_mask_enz].dot(ce)
else:
density = self.model_params['mws'].dot(x)
- # TODO check if rho should be specified as a bound!
return np.array([self.model_params['rho'] - density])
@cache_last
def heq_density_jac(self, x):
"""Calculate gradient of density constraint at x
+ for single density, this returns a gradient vector
+ for multiple densities (not implemented), this
+ would return a Jacobian
+
:param x: optimization variables
:type x: numpy.ndarray
:returns: mass balance jacobian at x
diff --git a/xbanalysis/utils/hessians.py b/xbanalysis/utils/hessians.py
index 9c63f41d805f768374338c4e0c0a4ea2a6c2df55..b3e7f97c0107cc7f6b21685415f170c4f48d303b 100644
--- a/xbanalysis/utils/hessians.py
+++ b/xbanalysis/utils/hessians.py
@@ -9,9 +9,11 @@ import numpy as np
def hess_from_grad_grad(grad_l, grad_r):
"""Calculate Hessian from two gradiends wrt x
+ actually: np.tril(np.outer(grad_l, grad_r))
+
: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: right 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]
@@ -28,6 +30,8 @@ def hess_from_grad_grad(grad_l, grad_r):
def hess_from_v_hesss(v, hesss):
"""Calculate Hessian from a vector and Hessians
+ actually: (v.reshape((-1,1,1))*hesss).sum(axis=0)
+
:param v: vector of length n
:param v: numpy.ndarray (1D)
:param hesss: Hessians of n vectors wrt to x variables
@@ -47,6 +51,8 @@ def hess_from_v_hesss(v, hesss):
def hess_from_jac_jac(jac_l, jac_r):
"""Calculate Hessian from two Jacobians
+ actually: # replaced by jac_l.T.dot(jac_r)
+
: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
@@ -66,6 +72,8 @@ def hess_from_jac_jac(jac_l, jac_r):
def hesss_from_matrix_hesss(matrix, hesss):
"""Calculate dot product of a stoichiometric matrix with Hessians
+ actually: np.array([(row.reshape((-1,1,1)) * hesss).sum(axis=0) for row in matrix])
+
: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
@@ -76,18 +84,18 @@ def hesss_from_matrix_hesss(matrix, hesss):
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
+ n_rows = matrix.shape[0]
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]
+ for i, v in enumerate(matrix):
+ result_hesss[i] = (v.reshape((-1, 1, 1)) * hesss).sum(axis=0)
return result_hesss
def hesss_from_v_hess(v, hess):
"""Calculate Hessians from vector and a Hessian
+ actually: v.reshape((-1,1,1)) * hess
+
:param v: vector of length m
:param v: numpy.ndarray (1D)
:param hess: Hessian of scalar wrt to x variables
@@ -108,6 +116,8 @@ def hesss_from_v_hess(v, hess):
def hesss_from_grad_jac(grad, jac):
"""Calculate Hessians from gradient and a Jacobian
+ actually: np.array([np.outer(grad, row) for row in jac])
+
: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
@@ -121,5 +131,27 @@ def hesss_from_grad_jac(grad, jac):
result_hesss = np.zeros((n_rows, size, size))
for i in range(n_rows):
- result_hesss[i] = hess_from_grad_grad(jac[i], grad)
+ result_hesss[i] = np.outer(grad, jac[i])
+ return result_hesss
+
+
+def hesss_from_jac_grad(jac, grad):
+ """Calculate Hessians from Jacobian and gradient
+
+ actually: np.array([np.outer(row, grad) for row in jac])
+
+ :param jac: Jacobian of a vector size m wrt to x variables
+ :param jac: numpy.ndarray (2D) with shape [m,x]
+ :param grad: gradient of scalar wrt to x variables
+ :param grad: numpy.ndarray (1D) with shape [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] = np.outer(jac[i], grad)
return result_hesss
diff --git a/xbanalysis/utils/opt_results.py b/xbanalysis/utils/opt_results.py
index ec6b9afce3d641274ffe4a67f39aedaf8352a28e..5a03ca3247d347639f5488c890d9ca20a64ddd04 100644
--- a/xbanalysis/utils/opt_results.py
+++ b/xbanalysis/utils/opt_results.py
@@ -22,24 +22,22 @@ class OptResults:
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))
+ self.ress = {'metabolites': self.problem.var_metab_sids,
+ 'enzymes': self.problem.var_enz_sids,
+ 'mr_reactions': self.problem.mr_rids,
+ 'enz_sid2idx': {sid: idx for idx, sid in enumerate(self.problem.var_enz_sids)},
+ 'var_mask_metab': self.problem.var_mask_metab,
+ 'var_mask_enz': self.problem.var_mask_enz,
+ 'grs_per_h': np.zeros(n_points),
+ 'status': np.zeros(n_points),
+ 'x_opt': np.zeros((n_points, self.problem.n_vars)),
+ 'cm_M': np.zeros((n_points, self._n_m)),
+ 'rho_m': np.zeros((n_points, self._n_m)),
+ 'ce_M': np.zeros((n_points, self._n_e)),
+ 'rho_e': np.zeros((n_points, self._n_e)),
+ 'alphas': np.zeros((n_points, self._n_e)),
+ 'mr_fluxes': np.zeros((n_points, len(self.problem.mr_rids))),
+ 'psm_fluxes': np.zeros((n_points, self._n_e))}
def add(self, info):
if self._idx < self._n_points:
@@ -54,10 +52,9 @@ class OptResults:
self.ress['rho_m'][i] = rho[self.problem.var_mask_metab]
self.ress['rho_e'][i] = rho[self.problem.var_mask_enz]
- # TODO extract reaction volume from kinetics to scale fluxes here
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['mr_fluxes'][i] = scale_factor * self.problem.mras(cx) / self.problem.mr_vols
+ self.ress['psm_fluxes'][i] = scale_factor / self.problem.pstmas(cx) / self.problem.ps_vols
self.ress['alphas'][i] = self.problem.get_alphas(cx)
self._idx += 1