README.md

@@ -2,6 +2,13 @@
 
 XBA problem generation from SBML code metabolic/kinetic models. Support of Growth Balance Analysis (GBA), FBA and RBA
 
+References to be included !!!
+- sbo terms (sbo.owl)
+- SBML and extensions for fbc, groups
+- pFBA, FBA, RBA, FVA
+- sbmlxdf
+
+
 ## Getting started
 
 To make it easy for you to get started with GitLab, here's a list of recommended next steps.

xbanalysis/model/__init__.py

@@ -6,5 +6,7 @@ from .xba_function import XbaFunction
 from .xba_parameter import XbaParameter
 from .xba_reaction import XbaReaction
 from .xba_species import XbaSpecies
+from .sbo_term import SboTerm
 
-__all__ = ['XbaModel', 'XbaCompartment', 'XbaFunction', 'XbaParameter', 'XbaReaction', 'XbaSpecies']
+__all__ = ['XbaModel', 'XbaCompartment', 'XbaFunction', 'XbaParameter',
+           'XbaReaction', 'XbaSpecies', 'SboTerm']

xbanalysis/model/sbo_term.py

+"""Implementation of SboTerm class.
+
+Peter Schubert, HHU Duesseldorf, July 2022
+"""
+import re
+import os
+from functools import lru_cache
+import xml.etree.ElementTree
+
+
+ns = {'owl': 'http://www.w3.org/2002/07/owl#',
+      'rdf': 'http://www.w3.org/1999/02/22-rdf-syntax-ns#',
+      'rdfs': 'http://www.w3.org/2000/01/rdf-schema#'}
+
+local_sbo_tree = {
+    'SBO_0000000': {'parents': [], 'label': 'systems biology representation'},
+    'SBO_0000167': {'parents': ['SBO_0000375'], 'label': 'biochemical or transport reaction'},
+    'SBO_0000176': {'parents': ['SBO_0000167'], 'label': 'biochemical reaction'},
+    'SBO_0000179': {'parents': ['SBO_0000176'], 'label': 'degradation'},
+    'SBO_0000180': {'parents': ['SBO_0000176'], 'label': 'dissociation'},
+    'SBO_0000181': {'parents': ['SBO_0000176'], 'label': 'conformational transition'},
+    'SBO_0000182': {'parents': ['SBO_0000176'], 'label': 'conversion'},
+    'SBO_0000183': {'parents': ['SBO_0000205'], 'label': 'transcription'},
+    'SBO_0000184': {'parents': ['SBO_0000205'], 'label': 'translation'},
+    'SBO_0000204': {'parents': ['SBO_0000205'], 'label': 'DNA replication'},
+    'SBO_0000205': {'parents': ['SBO_0000375'], 'label': 'composite biochemical process'},
+    'SBO_0000231': {'parents': ['SBO_0000000'], 'label': 'occurring entity representation'},
+    'SBO_0000236': {'parents': ['SBO_0000000'], 'label': 'physical entity representation'},
+    'SBO_0000240': {'parents': ['SBO_0000236'], 'label': 'material entity'},
+    'SBO_0000245': {'parents': ['SBO_0000240'], 'label': 'macromolecule'},
+    'SBO_0000246': {'parents': ['SBO_0000245'], 'label': 'information macromolecule'},
+    'SBO_0000247': {'parents': ['SBO_0000240'], 'label': 'simple chemical'},
+    'SBO_0000250': {'parents': ['SBO_0000246'], 'label': 'ribonucleic acid'},
+    'SBO_0000251': {'parents': ['SBO_0000246'], 'label': 'deoxyribonucleic acid'},
+    'SBO_0000252': {'parents': ['SBO_0000246'], 'label': 'polypeptide chain'},
+    'SBO_0000327': {'parents': ['SBO_0000247'], 'label': 'non-macromolecular ion'},
+    'SBO_0000328': {'parents': ['SBO_0000247'], 'label': 'non-macromolecular radical'},
+    'SBO_0000375': {'parents': ['SBO_0000231'], 'label': 'process'},
+    'SBO_0000377': {'parents': ['SBO_0000176'], 'label': 'isomerisation'},
+    'SBO_0000395': {'parents': ['SBO_0000375'], 'label': 'encapsulating process'},
+    'SBO_0000589': {'parents': ['SBO_0000205'], 'label': 'genetic production'},
+    'SBO_0000627': {'parents': ['SBO_0000631'], 'label': 'exchange reaction'},
+    'SBO_0000628': {'parents': ['SBO_0000631'], 'label': 'demand reaction'},
+    'SBO_0000629': {'parents': ['SBO_0000395'], 'label': 'biomass production'},
+    'SBO_0000630': {'parents': ['SBO_0000395'], 'label': 'ATP maintenance'},
+    'SBO_0000631': {'parents': ['SBO_0000375'], 'label': 'pseudoreaction'},
+    'SBO_0000632': {'parents': ['SBO_0000631'], 'label': 'sink reaction'},
+    'SBO_0000656': {'parents': ['SBO_0000182'], 'label': 'activation'},
+    'SBO_0000672': {'parents': ['SBO_0000176'], 'label': 'spontaneous reaction'},
+}
+data_dir = 'data'
+sbo_owl_fname = 'sbo.owl'
+
+
+class SboTerm:
+    """SboTerm variables
+
+    Implemented as a flyweight pattern and caching for is_in_batch lookups
+    Class builds a SBO hierachie tree, loaded from a 'sbo.owl' file,
+    if file not found, tree is constructed from a local dict.
+    Sbo term tree implemented as dictionary
+
+    Class provides methods to parse the tree
+    - retrieving parents
+    - retrieving children
+    - check if a sbo term belongs is within a specified tree branch
+    - retrieving subtrees (could be used to build default tree)
+    """
+    _instances = {}
+    sbo_tree = None
+
+    def __init__(self, sbo_term):
+        """Initialize SboTerm with sbo_term id
+
+        :param sbo_term: SBO term id
+        :type sbo_term: string in format 'SBO:0000xxx' or int
+        """
+        self.sbo_term = sbo_term
+
+    def __new__(cls, sbo_term):
+        """Create only one SboTerm per sbo_term.
+
+        Flyweight pattern
+
+        Improves performance on caching.
+        Also construct initial SBO Tree
+
+        :param sbo_term: SBO term id
+        :type sbo_term: string in format 'SBO:0000xxx' or int
+        """
+        self = cls._instances.get(sbo_term)
+        if self is None:
+            self = cls._instances[sbo_term] = object.__new__(SboTerm)
+            if cls.sbo_tree is None:
+                cls.load_sbo_tree()
+        return self
+
+    @property
+    def sbo_term(self):
+        return self.__sbo_term
+
+    @sbo_term.setter
+    def sbo_term(self, sbo_term):
+        if type(sbo_term) is str:
+            self.__sbo_term = re.sub(':', '_', sbo_term)
+        else:
+            self.__sbo_term = f'SBO_{sbo_term:07d}'
+
+    def __repr__(self):
+        return self.sbo_term
+
+    def get_parents(self):
+        """Get parent hierarchy for a specific node.
+
+        Supporting multiple parents of a node.
+
+        :return: ordered list of unique sbo_terms of parent hierarchy
+        :rtype: list of strings
+        """
+        return [repr(node) for node in type(self)._get_parents(self, [])]
+
+    def get_children(self):
+        """Get children subtree for a specific node.
+
+        :return: ordered list of unique sbo_terms for children subtree
+        :rtype: list of strings
+        """
+        return [repr(node) for node in type(self)._get_children(self, [])]
+
+    @lru_cache
+    def is_in_branch(self, upper_sbo_term):
+        """Check if sbo term is part subtree starting given value.
+
+        :param upper_sbo_term: SBO term that starts the branch
+        :rtype upper_sbo_term: string or int
+        :return: True if sbo_term part of the given branch
+        :rtype: bool
+        """
+        return repr(SboTerm(upper_sbo_term)) in self.get_parents() + [self.sbo_term]
+
+    def get_subtree(self):
+        """Retrieve a subtry including parents and children of give sbo term.
+
+        :return: all nodes forming the subtree with label and hierarchy information
+        :rtype: dict with nodes forming the subtree
+        """
+        nodes = set(type(self)._get_children(self, [self]))
+        nodes = nodes.union(set(type(self)._get_parents(self, [])))
+        sub_tree = {}
+        for node in nodes:
+            sbo_term = node.__sbo_term
+            sub_tree[sbo_term] = {'parents': type(self).sbo_tree[sbo_term]['parents'],
+                                  'label': type(self).sbo_tree[sbo_term]['label']}
+        return sub_tree
+
+    @classmethod
+    def _get_parents(cls, node, branch):
+        """Class method to recursively identify all parents.
+
+        :param node: sbo term of current node
+        :type node: SboTerm
+        :param branch: list of parents already identified
+        :type branch: list of SboTerm
+        :return: list of parents so far identified
+        :rtype: list of SboTerm
+        """
+        if node.__sbo_term in SboTerm.sbo_tree:
+            for parent in SboTerm.sbo_tree[node.__sbo_term]['parents']:
+                sbo_parent = SboTerm(parent)
+                branch.append(sbo_parent)
+                cls._get_parents(sbo_parent, branch)
+        return branch
+
+    @classmethod
+    def _get_children(cls, node, branch):
+        """Class method to recursively identify all children.
+
+        :param node: sbo term of current node
+        :type node: SboTerm
+        :param branch: list of children already identified
+        :type branch: list of SboTerm
+        :return: list of children so far identified
+        :rtype: list of SboTerm
+        """
+        if node.__sbo_term in SboTerm.sbo_tree:
+            for child in SboTerm.sbo_tree[node.__sbo_term]['children']:
+                sbo_child = SboTerm(child)
+                branch.append(sbo_child)
+                cls._get_children(sbo_child, branch)
+        return branch
+
+    @classmethod
+    def load_sbo_tree(cls):
+        """Class method to construct an SBO Tree from file or dict.
+
+        updated sbo file can be downloaded from https://www.ebi.ac.uk/sbo/main/download
+        and stored in the package under data/sbo.owl
+
+        """
+        here = os.path.dirname(os.path.abspath(__file__))
+        sbo_owl_path = os.path.join(here, data_dir, sbo_owl_fname)
+
+        if os.path.exists(sbo_owl_path):
+            tree = xml.etree.ElementTree.parse(sbo_owl_path)
+            root = tree.getroot()
+            sbo_tree = {}
+            for sbo_class in root.findall('owl:Class', ns):
+                sbo_term_id = sbo_class.attrib['{' + ns['rdf'] + '}' + 'about']
+                subclasses_of = []
+                label = ''
+                comment = ''
+                for child in sbo_class:
+                    if child.tag == '{' + ns['rdfs'] + '}' + 'subClassOf':
+                        if ('{' + ns['rdf'] + '}' + 'resource') in child.attrib:
+                            subclasses_of.append(child.attrib['{' + ns['rdf'] + '}' + 'resource'])
+                        else:
+                            for class_node in child:
+                                if class_node.tag == '{' + ns['owl'] + '}' + 'Class':
+                                    subclasses_of.append(class_node.attrib['{' + ns['rdf'] + '}' + 'about'])
+                    if child.tag == '{' + ns['rdfs'] + '}' + 'label':
+                        label = child.text
+                    if child.tag == '{' + ns['rdfs'] + '}' + 'comment':
+                        comment = child.text
+                sbo_tree[sbo_term_id] = {'parents': subclasses_of, 'label': label, 'comment': comment}
+        else:
+            print('loading local sbo subtree')
+            sbo_tree = local_sbo_tree
+
+        # add children information, so we can parse up/downwards through the tree
+        for sbo_term in sbo_tree:
+            sbo_tree[sbo_term]['children'] = []
+        for sbo_term in sbo_tree:
+            for parent in sbo_tree[sbo_term]['parents']:
+                if sbo_term not in sbo_tree[parent]['children']:
+                    sbo_tree[parent]['children'].append(sbo_term)
+        SboTerm.sbo_tree = sbo_tree

xbanalysis/model/xba_compartment.py

@@ -3,6 +3,11 @@
 Peter Schubert, HHU Duesseldorf, May 2022
 """
 
+import sbmlxdf
+
+from xbanalysis.utils.utils import xml_rba_ns
+from xbanalysis.model.sbo_term import SboTerm
+
 
 class XbaCompartment:
 
@@ -15,3 +20,9 @@ class XbaCompartment:
             self.units = s_compartment['units']
         if 'spatialDimension' in s_compartment:
             self.dimension = s_compartment['spatialDimension']
+        if 'xmlAnnotation' in s_compartment:
+            attrs = sbmlxdf.extract_xml_attrs(s_compartment['xmlAnnotation'], ns=xml_rba_ns, token='density')
+            if 'frac_gcdw' in attrs:
+                self.rba_frac_gcdw = float(attrs['frac_gcdw'])
+        if 'sboterm' in s_compartment:
+            self.sboterm = SboTerm(s_compartment['sboterm'])

xbanalysis/model/xba_model.py

@@ -41,7 +41,7 @@ class XbaModel:
 
         self.species = {sid: XbaSpecies(row)
                         for sid, row in model_dict['species'].iterrows()}
-        self.reactions = {rid: XbaReaction(row, self.functions, self.compartments)
+        self.reactions = {rid: XbaReaction(row, self.species, self.functions, self.compartments)
                           for rid, row in model_dict['reactions'].iterrows()}
 
         for r in self.reactions.values():

xbanalysis/model/xba_parameter.py

@@ -2,6 +2,7 @@
 
 Peter Schubert, HHU Duesseldorf, May 2022
 """
+from xbanalysis.model.sbo_term import SboTerm
 
 
 class XbaParameter:
@@ -13,4 +14,4 @@ class XbaParameter:
         self.constant = s_parameter['constant']
         self.units = s_parameter['units']
         if 'sboterm' in s_parameter:
-            self.sboterm = s_parameter['sboterm']
+            self.sboterm = SboTerm(s_parameter['sboterm'])
\ No newline at end of file

xbanalysis/model/xba_reaction.py

@@ -6,28 +6,31 @@ Peter Schubert, HHU Duesseldorf, May 2022
 import re
 import sbmlxdf
 from xbanalysis.utils.utils import get_species_stoic, get_local_params, expand_kineticlaw
+from xbanalysis.utils.utils import xml_rba_ns
+from xbanalysis.model.sbo_term import SboTerm
 
 
 class XbaReaction:
 
-    def __init__(self, s_reaction, functions, compartments):
+    def __init__(self, s_reaction, species, functions, compartments):
         self.id = s_reaction.name
         self.name = s_reaction.get('name', self.id)
         if 'sboterm' in s_reaction:
-            self.sboterm = s_reaction['sboterm']
+            self.sboterm = SboTerm(s_reaction['sboterm'])
         self.reaction_string = s_reaction['reactionString']
         self.reversible = s_reaction['reversible']
-        self.products = get_species_stoic(s_reaction['products'])
         self.reactants = get_species_stoic(s_reaction['reactants'])
+        self.products = get_species_stoic(s_reaction['products'])
+        self.metabs_only = True
+        for sid in list(self.reactants) + list(self.products):
+            if species[sid].sboterm.is_in_branch('SBO:0000247') is False:
+                self.metabs_only = False
         if 'modifiers' in s_reaction:
             self.modifiers = get_species_stoic(s_reaction['modifiers'])
         else:
             self.modifiers = {}
-        if 'fbcLb' in s_reaction and 'fbcUb' in s_reaction:
-            self.fbc_lb = s_reaction['fbcLb']
-            self.fbc_ub = s_reaction['fbcUb']
-        self.local_params = get_local_params(s_reaction)
         if 'kineticLaw' in s_reaction:
+            self.local_params = get_local_params(s_reaction)
             self.kinetic_law = s_reaction['kineticLaw'] if type(s_reaction['kineticLaw']) == str else ''
             # convert to numpy, replace local parameters with numerical values, inline lambda functions:
             expanded_kl = sbmlxdf.misc.mathml2numpy(self.kinetic_law)
@@ -37,41 +40,56 @@ class XbaReaction:
             for cid in compartments.keys():
                 if re.search(r'\b'+cid+r'\b', self.expanded_kl) is not None:
                     self.compartment = cid
+        # retrieve optional flux balance bounds, ensuring consistent bounds
+        if 'fbcLb' in s_reaction and 'fbcUb' in s_reaction:
+            self.fbc_lb = s_reaction['fbcLb']
+            self.fbc_ub = max(self.fbc_lb, s_reaction['fbcUb'])
+            if self.reversible is False:
+                self.fbc_lb = max(0, self.fbc_lb)
+                self.fbc_ub = max(0, self.fbc_ub)
+        if 'xmlAnnotation' in s_reaction:
+            attrs = sbmlxdf.extract_xml_attrs(s_reaction['xmlAnnotation'], ns=xml_rba_ns, token='kinetics')
+            if 'kapp_f_per_s' in attrs:
+                self.rba_kapp_f = float(attrs['kapp_f_per_s'])
+            if 'kapp_r_per_s' in attrs:
+                self.rba_kapp_r = float(attrs['kapp_r_per_s'])
+            if 'Km_ext_M' in attrs:
+                self.rba_km_ext = float(attrs['Km_ext_M'])
         self.ps_product = ''
-        self.enzyme = ''
+        self.enzyme = None
 
     def set_enzyme(self, species):
         for sid in self.modifiers:
-            if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+            if species[sid].sboterm.is_in_branch('SBO:0000246'):
                 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 in ('SBO:0000167', 'SBO:0000179', 'SBO:0000182'):
+        if self.sboterm.is_in_branch('SBO:0000167'):
             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 reactions
-        if self.sboterm == 'SBO:0000179':
+        if self.sboterm.is_in_branch('SBO:0000179'):
             for sid, val in self.reactants.items():
-                if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+                if species[sid].sboterm.is_in_branch('SBO:0000246'):
                     species[sid].add_p_degrad_rid(self.id, -val)
         # protein transitions
-        if self.sboterm == 'SBO:0000182':
+        if self.sboterm.is_in_branch('SBO:0000182'):
             for sid, val in self.reactants.items():
-                if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+                if species[sid].sboterm.is_in_branch('SBO:0000246'):
                     species[sid].add_p_transition_rid(self.id, -val)
             for sid, val in self.products.items():
-                if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+                if species[sid].sboterm.is_in_branch('SBO:0000246'):
                     species[sid].add_p_transition_rid(self.id, val)
 
         # protein synthesis reactions
-        if self.sboterm in ('SBO:0000589', 'SBO:0000205'):
+        if self.sboterm.is_in_branch('SBO:0000205'):
             for sid, val in self.products.items():
-                if species[sid].sboterm in ('SBO:0000252', 'SBO:0000250'):
+                if species[sid].sboterm.is_in_branch('SBO:0000246'):
                     self.ps_product = sid
                     species[sid].set_ps_rid(self.id)
                 else:

xbanalysis/model/xba_species.py

@@ -4,10 +4,8 @@ Peter Schubert, HHU Duesseldorf, May 2022
 """
 
 import sbmlxdf
-
-
-xml_gba_ns = 'http://www.hhu.de/ccb/gba/ns'
-xml_token = 'molecule'
+from xbanalysis.utils.utils import xml_gba_ns, xml_rba_ns
+from xbanalysis.model.sbo_term import SboTerm
 
 
 class XbaSpecies:
@@ -16,7 +14,7 @@ class XbaSpecies:
         self.id = s_species.name
         self.name = s_species.get('name', self.id)
         # set sboterm to 000247 to support Flux Balance Analysis
-        self.sboterm = s_species.get('sboterm', 'SBO:0000247')
+        self.sboterm = SboTerm(s_species.get('sboterm', 'SBO:0000247'))
         self.compartment = s_species['compartment']
         self.constant = s_species['constant']
         self.boundary = s_species['boundaryCondition']
@@ -25,8 +23,19 @@ class XbaSpecies:
         if 'units' in s_species:
             self.units = s_species['substanceUnits']
         if 'xmlAnnotation' in s_species:
-            attrs = sbmlxdf.extract_xml_attrs(s_species['xmlAnnotation'], ns=xml_gba_ns, token=xml_token)
+            attrs = sbmlxdf.extract_xml_attrs(s_species['xmlAnnotation'], ns=xml_gba_ns, token='molecule')
+            if 'weight_Da' in attrs:
                 self.mw = float(attrs['weight_Da'])
+            attrs = sbmlxdf.extract_xml_attrs(s_species['xmlAnnotation'], ns=xml_rba_ns, token='costs')
+            if 'target_mmol_per_gcdw' in attrs:
+                self.rba_target = float(attrs.pop('target_mmol_per_gcdw'))
+            if 'process_machine' in attrs:
+                self.rba_process_machine = attrs.pop('process_machine')
+            if 'process_capacity_per_s' in attrs:
+                self.rba_process_capacity = float(attrs.pop('process_capacity_per_s'))
+            if len(attrs) > 0:
+                self.rba_process_costs = attrs
+
         self.ps_rid = None
         self.m_reactions = {}
         self.ps_reactions = {}

xbanalysis/problems/__init__.py

 """Subpackage with XBA model classes """
 
 from .gba_problem import GbaProblem
+from .rba_problem import RbaProblem
 from .fba_problem import FbaProblem
 from .lp_problem import LpProblem
 
-__all__ = ['GbaProblem', 'FbaProblem', 'LpProblem']
+__all__ = ['GbaProblem', 'FbaProblem', 'LpProblem', 'RbaProblem']

xbanalysis/problems/fba_problem.py

@@ -35,29 +35,33 @@ class FbaProblem:
         :param xba_model: xba model with all required FBA parameters
         :type xba_model: XbaModel
         """
-        self.is_fba_model = False
         self.xba_model = xba_model
 
-        self.metab_sids = [s.id for s in xba_model.species.values() if s.sboterm == 'SBO:0000247']
-        metabs = set(self.metab_sids)
-        # identify reactions that do not use enzymes in reactants and products
-        self.mrids = [rid for rid, rxn in xba_model.reactions.items()
-                      if len((set(rxn.products) | set(rxn.reactants)) - metabs) == 0]
-        self.mrid2idx = {mrid: idx for idx, mrid in enumerate(self.mrids)}
-        self.reversible = np.array([self.xba_model.reactions[rid].reversible for rid in self.mrids])
-        self.s_mat = xba_model.get_stoic_matrix(self.metab_sids, self.mrids)
-
-        try:
-            self.flux_bounds = np.array([[xba_model.reactions[rid].fbc_lb for rid in self.mrids],
-                                         [xba_model.reactions[rid].fbc_ub for rid in self.mrids]])
-        except AttributeError:
-            print('Error: Model does not contain FBA flux bounds!')
-            return
+        # identify reaction excluding protein synthesis and protein degradation
+        self.mr_rids = [r.id for r in xba_model.reactions.values()
+                        if r.sboterm.is_in_branch('SBO:0000205') is False and
+                        r.sboterm.is_in_branch('SBO:0000179') is False]
+        self.mrid2idx = {mrid: idx for idx, mrid in enumerate(self.mr_rids)}
+        self.reversible = np.array([self.xba_model.reactions[rid].reversible for rid in self.mr_rids])
 
-        if hasattr(xba_model, 'fbc_objectives') is False:
-            print('Error: Model does not contain FBA objectives!')
+        sids = set()
+        for rid in self.mr_rids:
+            for sid in xba_model.reactions[rid].reactants:
+                sids.add(sid)
+            for sid in xba_model.reactions[rid].products:
+                sids.add(sid)
+        self.sids = list(sids)
+
+        self.s_mat = xba_model.get_stoic_matrix(self.sids, self.mr_rids)
+
+        self.is_fba_model = self.check_fba_model(xba_model)
+        if self.is_fba_model is False:
+            print('Error: not a FBA model - missing parameters')
             return
-        self.obj_coefs = np.zeros(len(self.mrids))
+
+        self.flux_bounds = np.array([[xba_model.reactions[rid].fbc_lb for rid in self.mr_rids],
+                                     [xba_model.reactions[rid].fbc_ub for rid in self.mr_rids]])
+        self.obj_coefs = None
         for oid, fbc_obj in xba_model.fbc_objectives.items():
             if fbc_obj.active is True:
                 self.obj_id = fbc_obj.id
@@ -65,21 +69,43 @@ class FbaProblem:
                 self.set_obj_coefs(fbc_obj.coefficiants)
                 break
 
-        if hasattr(self, 'obj_id') is False:
-            print('Error: Model does not contain an active FBA objective!')
-            return
+    def check_fba_model(self, xba_model):
+        """Check if fba related data available in the XBA model
 
-        self.is_fba_model = True
+        :param xba_model: Xba model with data extracted from SBML file
+        :type xba_model: XbaModel
+        :return: indicator if model contains data required for FBA problem formulation
+        :rtype: bool
+        """
+        is_fba_model = True
+
+        r_check = {'flux lower bound': 'fbc_lb', 'flux upper bound': 'fbc_ub'}
+        for check, param in r_check.items():
+            missing = [rid for rid in self.mr_rids if hasattr(xba_model.reactions[rid], param) is False]
+            if len(missing) > 0:
+                print(f'missing {check} in:', missing)
+                is_fba_model = False
+
+        if hasattr(xba_model, 'fbc_objectives') is False:
+            print('missing flux objectives')
+            is_fba_model = False
+        else:
+            active_objs = [oid for oid in xba_model.fbc_objectives if xba_model.fbc_objectives[oid].active]
+            if len(active_objs) == 0:
+                print('missing an active flux objective')
+                is_fba_model = False
+
+        return is_fba_model
 
     def combine_fwd_bwd(self, vector):
         """combine forward and backward flux information.
 
         :param vector: contains data on all reactions + backward information for reversible reactions
-        :type vector: 1D ndarray of shape (len(mrids) + sum(reversible))
+        :type vector: 1D ndarray of shape (len(mr_rids) + sum(reversible))
         :return: combined flux information
-        :rtype: 1D ndarray of shape (len(mrids))
+        :rtype: 1D ndarray of shape (len(mr_rids))
         """
-        n_cols = len(self.mrids)
+        n_cols = len(self.mr_rids)
         combined = vector[:n_cols]
         combined[self.reversible] -= vector[n_cols:]
         return combined
@@ -97,14 +123,14 @@ class FbaProblem:
             return None
 
         lp = LpProblem()
-        # self.reversible = np.array([self.xba_model.reactions[rid].reversible for rid in self.mrids])
+        # self.reversible = np.array([self.xba_model.reactions[rid].reversible for rid in self.mr_rids])
         fwd_bwd_s_mat = np.hstack((self.s_mat, -self.s_mat[:, self.reversible]))
         fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs[self.reversible]))
         fwd_bwd_bnds = np.hstack((np.fmax(self.flux_bounds, 0),
                                   np.vstack((np.fmax(-self.flux_bounds[1, self.reversible], 0),
                                              np.fmax(-self.flux_bounds[0, self.reversible], 0)))))
-        row_bnds = np.zeros((2, fwd_bwd_s_mat.shape[0]))
-        lp.configure(fwd_bwd_obj_coefs, fwd_bwd_bnds, fwd_bwd_s_mat, row_bnds, self.obj_dir)
+        row_bnds = np.zeros((fwd_bwd_s_mat.shape[0], 2))
+        lp.configure(fwd_bwd_obj_coefs, fwd_bwd_s_mat, row_bnds, fwd_bwd_bnds, self.obj_dir)
         return lp
 
     def fba_optimize(self):
@@ -115,7 +141,7 @@ class FbaProblem:
         """
         lp = self.create_fba_lp()
         if lp is None:
-            return {'simplex_status': 1, 'simplex_message': 'not an FBA model'}
+            return {'success': False, 'message': 'not an FBA model'}
         res = lp.solve()
         del lp
         res['x'] = self.combine_fwd_bwd(res['x'])
@@ -133,9 +159,9 @@ class FbaProblem:
         """
         lp = self.create_fba_lp()
         if lp is None:
-            return {'simplex_status': 1, 'simplex_message': 'not an FBA model'}
-        res = lp.solve(short_result=True)
-        if res['simplex_status'] == 0:
+            return {'success': False, 'message': 'not an FBA model'}
+        res = lp.solve(short_result=False)
+        if res['success']:
             # fix FBA objective as constraint
             fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs[self.reversible]))
             if self.obj_dir == 'maximize':
@@ -144,7 +170,7 @@ class FbaProblem:
                 row_bds = np.array([np.nan, res['fun'] / fract_of_optimum])
             lp.add_constraint(fwd_bwd_obj_coefs, row_bds)
             # new LP objective is to minimize sum of fluxes
-            n_cols = len(self.mrids)
+            n_cols = len(self.mr_rids)
             lp.set_obj_coefs(np.ones(n_cols + sum(self.reversible)))
             lp.set_obj_dir('minimize')
             res = lp.solve()
@@ -165,13 +191,13 @@ class FbaProblem:
         """
         lp = self.create_fba_lp()
         if lp is None:
-            return {'simplex_status': 1, 'simplex_message': 'not an FBA model'}
+            return {'success': False, 'message': 'not an FBA model'}
         res = lp.solve(short_result=True)
 
         if rids is None:
-            rids = self.mrids
+            rids = self.mr_rids
         flux_ranges = np.zeros((4, len(rids)))
-        if res['simplex_status'] == 0:
+        if res['success']:
             # fix FBA objective as constraint
             fwd_bwd_obj_coefs = np.hstack((self.obj_coefs, -self.obj_coefs[self.reversible]))
             if self.obj_dir == 'maximize':
@@ -180,7 +206,7 @@ class FbaProblem:
                 row_bds = np.array([np.nan, res['fun'] / fract_of_optimum])
             fba_constr_row = lp.add_constraint(fwd_bwd_obj_coefs, row_bds)
 
-            n_cols = len(self.mrids)
+            n_cols = len(self.mr_rids)
             for idx, rid in enumerate(rids):
                 # select reaction maximize/minimize flux
                 new_obj_coefs = np.zeros(n_cols)
@@ -190,20 +216,20 @@ class FbaProblem:
 
                 lp.set_obj_dir('minimize')
                 res = lp.solve(short_result=True)
-                if res['simplex_status'] == 0:
+                if res['success']:
                     flux_ranges[0, idx] = res['fun']
                     flux_ranges[2, idx] = lp.get_row_value(fba_constr_row)
                 else:
-                    print(f"FVA min failed for {rid}, {res['simplex_message']}")
+                    print(f"FVA min failed for {rid}")
                     flux_ranges[0, idx] = np.nan
 
                 lp.set_obj_dir('maximize')
                 res = lp.solve(short_result=True)
-                if res['simplex_status'] == 0:
+                if res['success']:
                     flux_ranges[1, idx] = res['fun']
                     flux_ranges[3, idx] = lp.get_row_value(fba_constr_row)
                 else:
-                    print(f"FVA max failed for {rid}, {res['simplex_message']}")
+                    print(f"FVA max failed for {rid}")
                     flux_ranges[1, idx] = np.nan
         del lp
 
@@ -228,7 +254,7 @@ class FbaProblem:
         :returns: flux bounds configured in the FBA problem
         :rtype: dict (key: reaction id, value: list with two float values for lower/upper bound)
         """
-        return {mrid: [bds[0], bds[1]] for mrid, bds in zip(self.mrids, self.flux_bounds.T)}
+        return {mrid: [bds[0], bds[1]] for mrid, bds in zip(self.mr_rids, self.flux_bounds.T)}
 
     def set_flux_bounds(self, flux_bounds):
         """selectivly set lower and upper flux bounds for given reactions.
@@ -255,6 +281,6 @@ class FbaProblem:
         :param coefficiants: objective coefficients of selected reactions
         :type coefficiants: dict (key: reaction id, value: float)
         """
-        self.obj_coefs = np.zeros(len(self.mrids))
+        self.obj_coefs = np.zeros(len(self.mr_rids))
         for rid, coef in coefficiants.items():
             self.obj_coefs[self.mrid2idx[rid]] = coef

xbanalysis/problems/gba_problem.py

@@ -43,8 +43,8 @@ class GbaProblem:
         # self.cache = True
 
         # 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')]
+        metab_sids = [s.id for s in xba_model.species.values() if s.sboterm.is_in_branch('SBO:0000247')]
+        enz_sids = [s.id for s in xba_model.species.values() if s.sboterm.is_in_branch('SBO:0000246')]
 
         # optimization variables
         self.var_metab_sids = [sid for sid in metab_sids if xba_model.species[sid].constant is False]
@@ -66,7 +66,7 @@ class GbaProblem:
         #  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')]
+                        if r.sboterm.is_in_branch('SBO:0000167')]
         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])

xbanalysis/problems/lp_problem.py

@@ -69,7 +69,6 @@ class LpProblem:
         self.lp = glpk.glp_create_prob()
         self.ncols = 0
         self.nrows = 0
-        self.simplex_result = None
         self.res = {}
 
     @staticmethod
@@ -92,30 +91,30 @@ class LpProblem:
             else:
                 ub = 0.0
                 var_type = glpk.GLP_FX
-        elif not np.isfinite(lb):
+        elif not np.isfinite(lb) and not np.isfinite(ub):
             lb = 0.0
-            var_type = glpk.GLP_UP
-        elif not np.isfinite(ub):
             ub = 0.0
-            var_type = glpk.GLP_LO
-        else:
+            var_type = glpk.GLP_FR
+        elif np.isfinite(ub):
             lb = 0.0
+            var_type = glpk.GLP_UP
+        else:
             ub = 0.0
-            var_type = glpk.GLP_FR
+            var_type = glpk.GLP_LO
         return var_type, lb, ub
 
-    def configure(self, obj_coefs, col_bnds, constr_coefs, row_bnds, direction='maximize'):
+    def configure(self, obj_coefs, constr_coefs, row_bnds, col_bnds, direction='maximize'):
         """Configure the LP Problem.
 
         Note: glpk related indices start at 1 (not zero)
         :param obj_coefs: coefficients of objective function
         :type obj_coefs: 1D ndarray of shape (1, ncols)
-        :param col_bnds: lower and upper bounds of structural (col) variables
-        :type col_bnds: 2D ndarray of shape(2, ncols) (1st row: lower bounds, 2nd row: upper bounds)
         :param constr_coefs: matrix with constraint coefficients
         :type constr_coefs: 2D nparray with shape(nrows, ncols)
         :param row_bnds: lower and upper bounds of auxhiliary (row) variables
-        :type row_bnds: 2D ndarray of shape(2, nrows) (1st row: lower bounds, 2nd row: upper bounds)
+        :type row_bnds: 2D ndarray of shape(nrows, 2) (1st col: lower bounds, 2nd col: upper bounds)
+        :param col_bnds: lower and upper bounds of structural (col) variables
+        :type col_bnds: 2D ndarray of shape(2, ncols) (1st row: lower bounds, 2nd row: upper bounds)
         :param direction: direction for the optimization ('maximize' or 'minimize')
         :type direction: string
         """
@@ -123,7 +122,7 @@ class LpProblem:
 
         assert len(obj_coefs) == self.ncols
         assert col_bnds.shape == (2, self.ncols)
-        assert row_bnds.shape == (2, self.nrows)
+        assert row_bnds.shape == (self.nrows, 2)
 
         # configure objective function
         glpk.glp_add_cols(self.lp, self.ncols)
@@ -137,7 +136,7 @@ class LpProblem:
 
         # configure constraints
         glpk.glp_add_rows(self.lp, self.nrows)
-        for i, lb_ub in enumerate(row_bnds.T):
+        for i, lb_ub in enumerate(row_bnds):
             var_type, lb, ub = self._get_variable_type(lb_ub)
             glpk.glp_set_row_bnds(self.lp, 1 + i, var_type, lb, ub)
         constr_coefs_sparse = scipy.sparse.coo_matrix(constr_coefs)
@@ -251,11 +250,11 @@ class LpProblem:
             glpk.glp_scale_prob(self.lp, opt)
         # sjj = [glpk.glp_get_sjj(self.lp, 1+i) for i in range(glpk.glp_get_num_cols(self.lp))]
         # rii = [glpk.glp_get_rii(self.lp, 1+i) for i in range(glpk.glp_get_num_rows(self.lp))]
-        self.simplex_result = glpk.glp_simplex(self.lp, None)
-        return self.results(short_result)
+        simplex_result = glpk.glp_simplex(self.lp, None)
+        return self.results(simplex_result, short_result)
 
     def get_row_value(self, row):
-        """Retrieve right hand side of a constraint (auxiliary variable)
+        """Retrieve right-hand side of a constraint (auxiliary variable)
 
         :param row: index into constraints
         :type row: int
@@ -265,11 +264,13 @@ class LpProblem:
         assert row <= glpk.glp_get_num_rows(self.lp)
         return glpk.glp_get_row_prim(self.lp, row)
 
-    def results(self, short_result=False):
+    def results(self, simplex_result, short_result=False):
         """Collect results for the optimization
 
         Alternatively, reduced results, e.g. when doing pFBA
 
+        :param simplex_result: result from the LP solver
+        :type simplex_result: int
         :param short_result: short result (only objective value and status) or long result
         :type short_result: bool (default:False)
         :return: result from the FBA optimization
@@ -277,13 +278,14 @@ class LpProblem:
         """
         self.res = {
             'fun': glpk.glp_get_obj_val(self.lp),
-            'simplex_status': self.simplex_result,
-            'simplex_message': simplex_status[self.simplex_result],
+            'success': glpk.glp_get_status(self.lp) == glpk.GLP_OPT,
         }
         if short_result is False:
             self.res['x'] = np.array([glpk.glp_get_col_prim(self.lp, 1 + i) for i in range(self.ncols)])
             self.res['reduced_costs'] = np.array([glpk.glp_get_col_dual(self.lp, 1 + i) for i in range(self.ncols)])
             self.res['shadow_prices'] = np.array([glpk.glp_get_row_dual(self.lp, 1 + i) for i in range(self.nrows)])
+            self.res['simplex_result'] = simplex_result
+            self.res['simplex_message'] = simplex_status[simplex_result],
             self.res['generic_status'] = generic_status[glpk.glp_get_status(self.lp)]
             self.res['primal_status'] = prim_status[glpk.glp_get_prim_stat(self.lp)]
             self.res['dual_status'] = dual_status[glpk.glp_get_dual_stat(self.lp)]

xbanalysis/problems/rba_problem.py

+"""Implementation of GBA Problem class.
+
+   - several process machines supported (yet to be validated)
+   - non-enzymatic metabolic reactions supported (yet to be validated)
+
+   Yet to implement:
+   - enzyme degradation
+   - enzyme transitions
+   - several compartments constraints
+   - target concentrations and target fluxes
+
+Peter Schubert, HHU Duesseldorf, June 2022
+"""
+
+import numpy as np
+
+from xbanalysis.problems.lp_problem import LpProblem
+
+
+class RbaProblem:
+
+    def __init__(self, xba_model):
+        """
+
+        :param xba_model:
+        """
+        # metabolic reactions and enzyme transitions (no FBA reactions, no protein synthesis, no degradation)
+        self.mr_rids = [r.id for r in xba_model.reactions.values() if r.sboterm.is_in_branch('SBO:0000167')
+                        and r.sboterm.is_in_branch('SBO:0000179') is False]
+        self.mrid2enz = {rid: xba_model.reactions[rid].enzyme for rid in self.mr_rids}
+        self.enz_sids = [enz_sid for enz_sid in self.mrid2enz.values() if enz_sid is not None]
+        self.rev_enz_sids = [enz_sid for mr_rid, enz_sid in self.mrid2enz.items()
+                             if enz_sid is not None and xba_model.reactions[mr_rid].reversible]
+        self.processes = {s.rba_process_machine: sid for sid, s in xba_model.species.items()
+                          if hasattr(s, 'rba_process_machine')}
+        self.psenz_rids = [xba_model.species[enz_sid].ps_rid for enz_sid in self.enz_sids]
+        self.pspm_rids = [xba_model.species[enz_sid].ps_rid for enz_sid in self.processes.values()]
+        self.densities = {cid: c.rba_frac_gcdw for cid, c in xba_model.compartments.items()
+                          if hasattr(c, 'rba_frac_gcdw')}
+
+        # metabolites and enzymes used in metabolic reactions and enzyme transitions
+        sids = set()
+        for rid in self.mr_rids:
+            for sid in xba_model.reactions[rid].reactants:
+                sids.add(sid)
+            for sid in xba_model.reactions[rid].products:
+                sids.add(sid)
+        used_sids = list(sids)
+        self.sids = [sid for sid in used_sids if xba_model.species[sid].constant is False]
+        self.ext_conc = {sid: xba_model.species[sid].initial_conc
+                         for sid in used_sids if xba_model.species[sid].constant is True}
+
+        self.rba_cols = self.mr_rids + self.enz_sids + list(self.processes)
+        self.rba_rows = (self.sids + list(self.processes)
+                         + [sid + '_feff' for sid in self.enz_sids]
+                         + [sid + '_reff' for sid in self.rev_enz_sids]
+                         + list(self.densities))
+
+        self.is_rba_model = self.check_rba_model(xba_model)
+        if self.is_rba_model is False:
+            print('Error: not a GBA model - missing parameters')
+            return
+
+        fix_m, var_m, row_bnds_m = self.construct_mass_balance(xba_model)
+        fix_pc, var_pc, row_bnds_pc = self.construct_process_capacities(xba_model)
+        fix_eff, var_eff, row_bnds_eff = self.construct_efficencies(xba_model)
+        fix_cd, var_cd, row_bnds_cd = self.construct_densities(xba_model)
+        self.fix_cd = fix_cd
+        self.fix_mat = np.vstack((fix_m, fix_pc, fix_eff, fix_cd))
+        self.var_mat = np.vstack((var_m, var_pc, var_eff, var_cd))
+        self.fix_row_bnds = np.vstack((np.zeros_like(row_bnds_m), np.zeros_like(row_bnds_pc),
+                                       row_bnds_eff, row_bnds_cd))
+        self.var_row_bnds = np.vstack((row_bnds_m, row_bnds_pc,
+                                       np.zeros_like(row_bnds_eff), np.zeros_like(row_bnds_cd)))
+
+        # in RBA we only check feasibility.
+        # Here we maximize enzyme amount and check that objective value e.g. > 1e-10 mmol/gcdw
+        self.obj_dir = 'maximize'
+        self.obj_coefs = np.hstack((np.zeros(len(self.mr_rids)),
+                                    np.ones(len(self.enz_sids) + len(self.processes))))
+        col_lbs = np.hstack(([xba_model.reactions[rid].fbc_lb for rid in self.mr_rids],
+                             np.zeros(len(self.enz_sids) + len(self.processes))))
+        col_ubs = np.hstack(([xba_model.reactions[rid].fbc_ub for rid in self.mr_rids],
+                             np.nan * np.ones(len(self.enz_sids) + len(self.processes))))
+        self.col_bnds = np.vstack((col_lbs, col_ubs))
+        self.scaling = '2N'
+        self.max_iter = 20
+
+    def check_rba_model(self, xba_model):
+        """Check if rba related data available in the XBA model
+
+        :param xba_model: Xba model with data extracted from SBML file
+        :type xba_model: XbaModel
+        :return: indicator if model contains data required for RBA problem formulation
+        :rtype: bool
+        """
+        is_rba_model = True
+
+        r_check = {'flux lower bound': 'fbc_lb', 'flux upper bound': 'fbc_ub', 'kapp forward': 'rba_kapp_f'}
+        for check, param in r_check.items():
+            missing = [rid for rid in self.mr_rids if hasattr(xba_model.reactions[rid], param) is False]
+            if len(missing) > 0:
+                print(f'missing {check} in:', missing)
+                is_rba_model = False
+
+        missing = [rid for rid in self.mr_rids if xba_model.reactions[rid].reversible
+                   and hasattr(xba_model.reactions[rid], 'rba_kapp_r') is False]
+        if len(missing) > 0:
+            print('missing kapp reverse in:', missing)
+            is_rba_model = False
+
+        missing = [sid for sid in self.sids + self.enz_sids + list(self.processes.values())
+                   if hasattr(xba_model.species[sid], 'mw') is False]
+        if len(missing) > 0:
+            print('missing molecular weight in:', missing)
+            is_rba_model = False
+
+        missing = [sid for sid in list(self.processes.values())
+                   if hasattr(xba_model.species[sid], 'rba_process_capacity') is False]
+        if len(missing) > 0:
+            print('missing process capacity in:', missing)
+            is_rba_model = False
+
+        if len(self.densities) == 0:
+            print('missing a density constraint')
+            is_rba_model = False
+        return is_rba_model
+
+    def rba_optimize(self):
+        """
+
+        :return:
+        """
+        if self.is_rba_model is False:
+            print('Error: not a RBA model - missing parameters')
+            res = {'success': False, 'fun': 0.0, 'message': 'not an RBA model'}
+            return res
+
+        res = {}
+        # Preanalysis: find range of growth rates
+        gr = {'lb': 0.0, 'ub': np.nan, 'try': 1.0}
+        for i in range(self.max_iter):
+            constr_coefs = self.fix_mat + gr['try'] * self.var_mat
+            # row_bnds = self.fix_row_bnds + gr['try'] * self.var_row_bnds
+            lp = LpProblem()
+            lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
+            res = lp.solve(scaling=self.scaling)
+            del lp
+            if res['success'] and res['fun'] > 1e-10:
+                gr['lb'] = gr['try']
+                gr['try'] *= 10
+            else:
+                if gr['try'] <= 1e-4:
+                    break
+                gr['ub'] = gr['try']
+                gr['try'] /= 10
+            if gr['lb'] > 0.0 and np.isfinite(gr['ub']):
+                break
+
+        if gr['lb'] == 0.0 or np.isnan(gr['ub']):
+            res = {'success': False, 'fun': 0.0, 'message': 'no solution found'}
+            return res
+
+        # find exact growth rate
+        for i in range(self.max_iter):
+            gr['try'] = 0.5 * (gr['lb'] + gr['ub'])
+            constr_coefs = self.fix_mat + gr['try'] * self.var_mat
+            # row_bnds = self.fix_row_bnds + gr['try'] * self.var_row_bnds
+            lp = LpProblem()
+            lp.configure(self.obj_coefs, constr_coefs, self.fix_row_bnds, self.col_bnds, self.obj_dir)
+            res = lp.solve(scaling=self.scaling)
+            del lp
+            if res['success'] and res['fun'] > 1e-10:
+                gr['lb'] = gr['try']
+            else:
+                gr['ub'] = gr['try']
+            if (gr['ub'] - gr['lb']) < gr['try'] * 1e-5:
+                break
+        res['fun'] = gr['try']
+        return res
+
+    def construct_mass_balance(self, xba_model):
+
+        fix_m_x_mr = xba_model.get_stoic_matrix(self.sids, self.mr_rids)
+        var_m_x_enz = xba_model.get_stoic_matrix(self.sids, self.psenz_rids)
+        var_m_x_pm = xba_model.get_stoic_matrix(self.sids, self.pspm_rids)
+        fix_m = np.hstack((fix_m_x_mr,
+                           np.zeros((len(self.sids), len(self.enz_sids)+len(self.processes)))))
+        var_m = np.hstack((np.zeros((len(self.sids), len(self.mr_rids))),
+                           var_m_x_enz,
+                           var_m_x_pm))
+        targets = np.array([getattr(xba_model.species[sid], 'rba_target', 0.0) for sid in self.sids])
+        row_bnds_m = np.vstack((targets, targets)).T
+        return fix_m, var_m, row_bnds_m
+
+    def construct_process_capacities(self, xba_model):
+
+        var_pc_x_enz = np.zeros((len(self.processes), len(self.psenz_rids)))
+        for idx, pm in enumerate(self.processes):
+            var_pc_x_enz[idx] = [xba_model.species[enz_sid].rba_process_costs.get(pm, 0.0)
+                                 for enz_sid in self.enz_sids]
+        var_pc_x_pm = np.zeros((len(self.processes), len(self.processes)))
+        for idx, pm in enumerate(self.processes):
+            var_pc_x_pm[idx] = [xba_model.species[enz_sid].rba_process_costs.get(pm, 0.0)
+                                for enz_sid in self.processes.values()]
+        fix_pc_x_pm = -np.diag([xba_model.species[enz_sid].rba_process_capacity
+                                for enz_sid in self.processes.values()]) * 3600
+        fix_pc = np.hstack((np.zeros((len(self.processes), len(self.mr_rids) + len(self.enz_sids))),
+                            fix_pc_x_pm))
+        var_pc = np.hstack((np.zeros((len(self.processes), len(self.mr_rids))),
+                            var_pc_x_enz, var_pc_x_pm))
+        row_bnds_pc = np.zeros((len(self.processes), 2))
+        return fix_pc, var_pc, row_bnds_pc
+
+    def get_efficiencies(self, xba_model):
+
+        efficiencies = np.zeros((2, len(self.enz_sids)))
+        idx = 0
+        for mr_rid, enz_sid in self.mrid2enz.items():
+            if enz_sid is not None:
+                r = xba_model.reactions[mr_rid]
+                saturation = 1.0
+                if hasattr(r, 'rba_km_ext'):
+                    for sid in r.reactants:
+                        if sid in self.ext_conc:
+                            saturation *= self.ext_conc[sid] / (self.ext_conc[sid] + r.rba_km_ext)
+                efficiencies[0, idx] = r.rba_kapp_f * 3600 * saturation
+                saturation = 1.0
+                if r.reversible:
+                    if hasattr(r, 'rba_km_ext'):
+                        for sid in r.products:
+                            if sid in self.ext_conc:
+                                saturation *= self.ext_conc[sid] / (self.ext_conc[sid] + r.rba_km_ext)
+                    efficiencies[1, idx] = r.rba_kapp_r * 3600 * saturation
+                idx += 1
+        return efficiencies
+
+    def construct_efficencies(self, xba_model):
+
+        efficiencies = self.get_efficiencies(xba_model)
+        mask = np.array([enz_sid is not None for enz_sid in self.mrid2enz.values()])
+        fix_feff_x_mr = np.diag(np.ones(len(self.mrid2enz)))[mask]
+        fix_feff_x_enz = np.diag(-efficiencies[0])
+        fix_feff_x_pm = np.zeros((len(self.enz_sids), len(self.processes)))
+        fix_feff = np.hstack((fix_feff_x_mr, fix_feff_x_enz, fix_feff_x_pm))
+        row_bnds_feff = np.vstack((np.ones(len(self.enz_sids)) * np.nan, np.zeros(len(self.enz_sids)))).T
+
+        mask = [enz_sid is not None and xba_model.reactions[mr_rid].reversible
+                for mr_rid, enz_sid in self.mrid2enz.items()]
+        fix_reff_x_mr = np.diag(-np.ones(len(self.mrid2enz)))[mask]
+        mask = [enz_sid in self.rev_enz_sids for enz_sid in self.enz_sids]
+        fix_reff_x_enz = np.diag(-efficiencies[1])[mask]
+        fix_reff_x_pm = np.zeros((len(self.rev_enz_sids), len(self.processes)))
+        fix_reff = np.hstack((fix_reff_x_mr, fix_reff_x_enz, fix_reff_x_pm))
+        row_bnds_reff = np.vstack((np.ones(len(self.rev_enz_sids)) * np.nan, np.zeros(len(self.rev_enz_sids)))).T
+        fix_eff = np.vstack((fix_feff, fix_reff))
+        var_eff = np.zeros_like(fix_eff)
+        row_bnds_eff = np.vstack((row_bnds_feff, row_bnds_reff))
+        return fix_eff, var_eff, row_bnds_eff
+
+    def construct_densities(self, xba_model):
+        fix_cd_x_mr = np.zeros((len(self.densities), len(self.mr_rids)))
+        fix_cd_x_enz = np.zeros((len(self.densities), len(self.enz_sids)))
+        for row_idx, cid in enumerate(self.densities):
+            for col_idx, sid in enumerate(self.enz_sids):
+                if xba_model.species[sid].compartment == cid:
+                    fix_cd_x_enz[row_idx, col_idx] = xba_model.species[sid].mw / 1000
+        fix_cd_x_pm = np.zeros((len(self.densities), len(self.processes)))
+        for row_idx, cid in enumerate(self.densities):
+            for col_idx, sid in enumerate(self.processes.values()):
+                if xba_model.species[sid].compartment == cid:
+                    fix_cd_x_pm[row_idx, col_idx] = xba_model.species[sid].mw / 1000
+        fix_cd = np.hstack((fix_cd_x_mr, fix_cd_x_enz, fix_cd_x_pm))
+
+        var_cd = np.zeros_like(fix_cd)
+        row_bnds_cd = np.vstack((np.zeros(len(self.densities)),
+                                 [val for val in self.densities.values()])).T
+        return fix_cd, var_cd, row_bnds_cd

xbanalysis/utils/utils.py

@@ -6,6 +6,9 @@ Peter Schubert, HHU Duesseldorf, May 2022
 import sys
 import re
 
+xml_gba_ns = 'http://www.hhu.de/ccb/gba/ns'
+xml_rba_ns = 'http://www.hhu.de/ccb/rba/ns'
+
 
 def get_sbml_ids(df, sbo_id):
     """Get SBML IDs with specified sboterm-id.