SboTerm class implemented together with branch support

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

@@ -6,6 +6,7 @@ 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:
@@ -23,3 +24,5 @@ class XbaCompartment:
             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_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']
\ No newline at end of file
+            self.sboterm = SboTerm(s_parameter['sboterm'])

xbanalysis/model/xba_reaction.py

@@ -7,6 +7,7 @@ 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:
@@ -15,14 +16,14 @@ class XbaReaction:
         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.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 != 'SBO:0000247':
+            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'])
@@ -59,36 +60,36 @@ class XbaReaction:
 
     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

@@ -5,6 +5,7 @@ Peter Schubert, HHU Duesseldorf, May 2022
 
 import sbmlxdf
 from xbanalysis.utils.utils import xml_gba_ns, xml_rba_ns
+from xbanalysis.model.sbo_term import SboTerm
 
 
 class XbaSpecies:
@@ -13,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']
@@ -26,10 +27,13 @@ class XbaSpecies:
             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 len(attrs) > 0:
-                if ('process_machine' in attrs) and ('process_capacity_per_s' in attrs):
+            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

xbanalysis/problems/fba_problem.py

@@ -37,12 +37,22 @@ class FbaProblem:
         """
         self.xba_model = xba_model
 
-        self.metab_sids = [s.id for s in xba_model.species.values() if s.sboterm == 'SBO:0000247']
-        # identify reactions that do not use enzymes in reactants and products
-        self.mr_rids = [r.id for r in xba_model.reactions.values() if r.metabs_only]
+        # 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])
-        self.s_mat = xba_model.get_stoic_matrix(self.metab_sids, self.mr_rids)
+
+        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:

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

@@ -91,16 +91,16 @@ 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, constr_coefs, row_bnds, col_bnds, direction='maximize'):

xbanalysis/problems/rba_problem.py

@@ -24,13 +24,9 @@ class RbaProblem:
 
         :param xba_model:
         """
-        # metabolites and enzymes in the model
-        metab_sids = [s.id for s in xba_model.species.values() if s.sboterm == 'SBO:0000247']
-        self.m_sids = [sid for sid in metab_sids if xba_model.species[sid].constant is False]
-        self.ext_conc = {sid: xba_model.species[sid].initial_conc
-                         for sid in metab_sids if xba_model.species[sid].constant is True}
-        self.mr_rids = [r.id for r in xba_model.reactions.values()
-                        if r.sboterm in ('SBO:0000167', 'SBO:0000179', 'SBO:0000182') and r.metabs_only]
+        # 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()
@@ -42,8 +38,20 @@ class RbaProblem:
         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.m_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))
@@ -60,7 +68,10 @@ class RbaProblem:
         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.row_bnds = np.vstack((row_bnds_m, row_bnds_pc, row_bnds_eff, row_bnds_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
@@ -98,7 +109,7 @@ class RbaProblem:
             print('missing kapp reverse in:', missing)
             is_rba_model = False
 
-        missing = [sid for sid in self.m_sids + self.enz_sids + list(self.processes.values())
+        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)
@@ -130,8 +141,9 @@ class RbaProblem:
         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.row_bnds, self.col_bnds, self.obj_dir)
+            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:
@@ -145,16 +157,17 @@ class RbaProblem:
             if gr['lb'] > 0.0 and np.isfinite(gr['ub']):
                 break
 
-        if gr['lb'] == 0.0 or np.isfinite(gr['ub']) is False:
-            res['success'] = False
+        if gr['lb'] == 0.0 or np.isnan(gr['ub']):
+            res = {'success': False, 'fun': 0.0, 'message': 'no solution found'}
             return res
 
-        # find eaxt growth rate
+        # 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.row_bnds, self.col_bnds, self.obj_dir)
+            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:
@@ -168,15 +181,16 @@ class RbaProblem:
 
     def construct_mass_balance(self, xba_model):
 
-        fix_m_x_mr = xba_model.get_stoic_matrix(self.m_sids, self.mr_rids)
-        var_m_x_enz = xba_model.get_stoic_matrix(self.m_sids, self.psenz_rids)
-        var_m_x_pm = xba_model.get_stoic_matrix(self.m_sids, self.pspm_rids)
+        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.m_sids), len(self.enz_sids)+len(self.processes)))))
-        var_m = np.hstack((np.zeros((len(self.m_sids), len(self.mr_rids))),
+                           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))
-        row_bnds_m = np.zeros((len(self.m_sids), 2))
+        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):