improving linear pathways

modelpruner/graph/linear_pathway.py

-"""
+"""Implementation of LinearPathway Class.
 
-Peter Schubert, October 2022
+Peter Schubert, HHU Duesseldorf, October 2022
 """
 
 
 import numpy as np
 
 from modelpruner.graph.pathway_segment import PathwaySegment
+from modelpruner.graph.pathway_segment import ReactionLeg
 
 
 class LinearPathway:
@@ -21,11 +22,48 @@ class LinearPathway:
         self.init_sid = init_sid
         self.balance_c2_sids = set(balance_c2_sids)
         self.ng = network_graph
+        self.net_reactants = {}
+        self.net_products = {}
 
+        # construct the linear pathway
         self.segments = [self.get_initial_segment()]
         self.add_downstream_segments()
         self.add_upstream_segments()
 
+        # collect pathway data
+        self.sids = []
+        self.rids = [self.segments[0].producing.rid]
+        self.reversible = self.segments[0].reversible
+        scale = 1.0
+        direction = 1.0 if self.segments[0].producing.forward is True else -1.0
+        self.rids_scale = {self.segments[0].producing.rid: scale * direction}
+        for pwy_seg in self.segments:
+            self.sids.append(pwy_seg.sid)
+            self.rids.append(pwy_seg.consuming.rid)
+            scale *= pwy_seg.producing.stoic / pwy_seg.consuming.stoic
+            direction = 1.0 if pwy_seg.consuming.forward is True else -1.0
+            self.rids_scale[pwy_seg.consuming.rid] = scale * direction
+
+        net_stoic = self.get_net_stoichiometry()
+        self.net_reactants = {sid: -stoic for sid, stoic in net_stoic.items() if stoic < 0.0}
+        self.net_products = {sid: stoic for sid, stoic in net_stoic.items() if stoic > 0.0}
+
+    def get_net_stoichiometry(self):
+
+        substrates = set()
+        for rid in self.rids_scale:
+            substrates |= set(self.ng.edges[rid].reactants)
+            substrates |= set(self.ng.edges[rid].products)
+        net_stoic = {substrate: 0.0 for substrate in substrates}
+
+        for rid, scale in self.rids_scale.items():
+            for sid, stoic in self.ng.edges[rid].reactants.items():
+                net_stoic[sid] -= stoic * scale
+            for sid, stoic in self.ng.edges[rid].products.items():
+                net_stoic[sid] += stoic * scale
+
+        return net_stoic
+
     def get_initial_segment(self):
         """get initial pathway segment for inital node.
 
@@ -37,64 +75,56 @@ class LinearPathway:
         :rtype: PathwaySegment
         """
         self.balance_c2_sids.remove(self.init_sid)
-
         connectivity = self.get_us_ds_neighbors(self.init_sid)
-        pwy_seg = PathwaySegment(connectivity['sid'])
-        pwy_seg.set_upstream(connectivity['p_rid'], connectivity['nb_us'])
-        pwy_seg.set_downstream(connectivity['c_rid'], connectivity['nb_ds'])
-        pwy_seg.set_reversible(connectivity['reversible'])
-        return pwy_seg
+        return PathwaySegment(connectivity)
 
     def add_downstream_segments(self):
         """add downstream segments """
 
         lcount = 1
         upstream_seg = self.segments[-1]
-        while upstream_seg.nb_ds is not None and lcount < 50:
+        while ((upstream_seg.consuming is not None) and
+               (upstream_seg.consuming.nb is not None) and
+               (lcount < 50)):
             lcount += 1
             pwy_rev = upstream_seg.reversible
-            self.balance_c2_sids.remove(upstream_seg.nb_ds)
+            self.balance_c2_sids.remove(upstream_seg.consuming.nb)
             connectivity = self.get_ds_neighbor(upstream_seg)
-
-            if 'c_rid' in connectivity:
-                pwy_seg = PathwaySegment(connectivity['sid'])
-                pwy_seg.set_upstream(upstream_seg.c_rid, upstream_seg.sid)
-                pwy_seg.set_downstream(connectivity['c_rid'], connectivity['nb_ds'])
+            pwy_seg = PathwaySegment(connectivity)
             pwy_seg.set_reversible(pwy_rev and connectivity['reversible'])
             self.segments.append(pwy_seg)
 
-                if connectivity['nb_ds'] is not None:
-                    # treat a reversible pathway that becomes directional (i.e. continue ds)
-                    if (pwy_rev is True) and (connectivity['reversible'] is False):
+            if connectivity['consuming'] is None:
+                # termination of the pathway at a c1 metabolite (blocked node)
+                break
+
+            if connectivity['consuming'].nb is None:
+                # termination of the pathway with connectivity towards the network
+                if pwy_rev is True:
+                    if connectivity['reversible'] is False:
+                        # irreversible termination in a reversible pathway
                         for pwy_seg in self.segments:
                             pwy_seg.set_reversible(False)
-                    # in other cases, just continue ds
-                # cases of path termination
                     else:
-                    # pathway termination of a reversible pathway (reverse pwy_segs and continue ds)
-                    if (pwy_rev is True) and (connectivity['reversible'] is True):
+                        # reversible termination in a reversible pathway
                         for pwy_seg in self.segments:
-                            pwy_seg.set_reversible(False)
                             pwy_seg.reverse_direction()
                         self.segments.reverse()
-                    # in other cases, the path termination terminates the while loop
 
-            # i.e. no consuming node, but only one producing node
-            #  we reverse tha pwy_segments and their direction and continue ds
+            if connectivity['consuming'].nb is not upstream_seg.sid:
+                # downstream connectivity to a c2 node
+                if (pwy_rev is True) and (connectivity['reversible'] is False):
+                    for pwy_seg in self.segments:
+                        pwy_seg.set_reversible(False)
             elif pwy_rev is True:
-                assert 'p_rid' in connectivity
+                # pathway reversal situation
                 assert connectivity['reversible'] is False
-                pwy_seg = PathwaySegment(connectivity['sid'])
-                pwy_seg.set_upstream(upstream_seg.c_rid, upstream_seg.sid)
-                pwy_seg.set_downstream(connectivity['p_rid'], connectivity['nb_us'])
-                self.segments.append(pwy_seg)
-                for pwy_seg in self.segments:
+                for pwy_seg in self.segments[:-1]:
                     pwy_seg.set_reversible(False)
                     pwy_seg.reverse_direction()
                 self.segments.reverse()
             else:
-                # conditions when directional pathway has producing reaction
-                print(f'inconsistent configuration, starting at {self.init_sid}')
+                print(f'inconsistent configuration at {pwy_seg.sid}, starting from {self.init_sid}')
                 break
             upstream_seg = self.segments[-1]
 
@@ -102,37 +132,46 @@ class LinearPathway:
 
         lcount = 1
         downstream_seg = self.segments[0]
-        while downstream_seg.nb_us in self.balance_c2_sids and lcount < 50:
+        while ((downstream_seg.producing is not None) and
+               (downstream_seg.producing.nb in self.balance_c2_sids) and
+               (lcount < 50)):
             lcount += 1
             pwy_rev = downstream_seg.reversible
-            self.balance_c2_sids.remove(downstream_seg.nb_us)
+            self.balance_c2_sids.remove(downstream_seg.producing.nb)
             connectivity = self.get_us_neighbor(downstream_seg)
-
-            if 'p_rid' in connectivity:
-                pwy_seg = PathwaySegment(connectivity['sid'])
-                pwy_seg.set_downstream(downstream_seg.p_rid, downstream_seg.sid)
-                pwy_seg.set_upstream(connectivity['p_rid'], connectivity['nb_us'])
+            pwy_seg = PathwaySegment(connectivity)
             pwy_seg.set_reversible(pwy_rev and connectivity['reversible'])
             self.segments.insert(0, pwy_seg)
-                if (connectivity['reversible'] is False) and (pwy_rev is True):
+
+            if connectivity['producing'] is None:
+                # termination of the pathway at a c1 metabolite (blocked node)
+                break
+
+            if connectivity['producing'].nb is None:
+                # termination of the pathway with connectivity towards the network
+                if (pwy_rev is True) and (connectivity['reversible'] is False):
+                    for pwy_seg in self.segments:
+                        pwy_seg.set_reversible(False)
+                break
+
+            if connectivity['producing'].nb is not downstream_seg.sid:
+                # upstream connectivity to a c2 node
+                if (pwy_rev is True) and (connectivity['reversible'] is False):
                     for pwy_seg in self.segments:
                         pwy_seg.set_reversible(False)
             elif pwy_rev is True:
-                assert 'c_rid' in connectivity
+                # pathway reversal situation
                 assert connectivity['reversible'] is False
-                pwy_seg = PathwaySegment(connectivity['sid'])
-                pwy_seg.set_downstream(downstream_seg.c_rid, downstream_seg.sid)
-                pwy_seg.set_upstream(connectivity['c_rid'], connectivity['nb_ds'])
-                self.segments.insert(0, pwy_seg)
-                for pwy_seg in self.segments:
+                for pwy_seg in self.segments[:-1]:
+                    pwy_seg.set_reversible(False)
                     pwy_seg.reverse_direction()
                 self.segments.reverse()
             else:
-                print(f'inconsistent configuration, starting at {self.init_sid}')
+                print(f'inconsistent configuration at {pwy_seg.sid}, starting from {self.init_sid}')
                 break
             downstream_seg = self.segments[0]
 
-    def get_neighbor_node(self, substrate_sids):
+    def get_neighbor_node(self, substrates):
         """Get a c2 node from reactants/product sids
 
         If none of the reactant/products is a c2 node,
@@ -141,12 +180,12 @@ class LinearPathway:
         In case of several c2 nodes, return the first node
         after sorting.
 
-        :param substrate_sids: metabolite ids (reactants or products)
-        :type substrate_sids: list-like of str
+        :param substrates: substrates with stoichiometry
+        :type substrates: dict (key: sid - str, val: stoichiometry - float)
         :return: metabolite id (sid) for neighbor node
         :rtype: str or None
         """
-        nbs = substrate_sids.intersection(self.balance_c2_sids)
+        nbs = set(substrates.keys()).intersection(self.balance_c2_sids)
         if len(nbs) == 0:
             nb = None
         else:
@@ -164,49 +203,57 @@ class LinearPathway:
         cur_node = self.ng.nodes[init_sid]
 
         rids = [rid for rid in (cur_node.consuming_rids | cur_node.producing_rids)]
-        revs = [rev for rev in (cur_node.consuming_rids | cur_node.producing_rids).values()]
+        revs = [val['rev'] for val in (cur_node.consuming_rids | cur_node.producing_rids).values()]
         rev = bool(np.all(np.array(revs)))
-        p_rid = None
-        c_rid = None
-        nb_us = None
-        nb_ds = None
+        producing_leg = None
+        consuming_leg = None
 
         if cur_node.connectivity == 1:
             # c1 nodes are actually blocked nodes where a pathway ends
             rid = rids[0]
 
             if rid in cur_node.consuming_rids:
-                c_rid = rid
-                nb_ds = self.get_neighbor_node(self.ng.edges[c_rid].products)
+                stoic = cur_node.consuming_rids[rid]['stoic']
+                nb = self.get_neighbor_node(self.ng.edges[rid].products)
+                consuming_leg = ReactionLeg(rid, stoic, nb, True)
             else:
-                p_rid = rid
-                nb_us = self.get_neighbor_node(self.ng.edges[p_rid].reactants)
-
-            # in case of a reversible reaction, configure consuming path
-            if (rev is True) and (c_rid is None):
-                c_rid = p_rid
-                nb_ds = nb_us
-                p_rid = None
-                nb_us = None
+                stoic = cur_node.producing_rids[rid]['stoic']
+                nb = self.get_neighbor_node(self.ng.edges[rid].reactants)
+                if rev is False:
+                    producing_leg = ReactionLeg(rid, stoic, nb, True)
+                else:
+                    # in case of a reversible producing reaction, configure consuming path
+                    consuming_leg = ReactionLeg(rid, stoic, nb, False)
 
         if cur_node.connectivity == 2:
             if rev is True:
                 # process a node with reversible fluxes on both sides
                 #  one being considered as producing, one as consuming reaction
-                nb_sids = [None, None]
-                for idx, rid in enumerate(rids):
-                    if rid in cur_node.consuming_rids:
-                        nb_sids[idx] = self.get_neighbor_node(self.ng.edges[rid].products)
+                nbs = [None, None]
+                stoics = [0.0, 0.0]
+                fwds = [True, True]
+
+                # arbitrarily select first reaction as producing and second as consuming
+                if rids[0] in cur_node.consuming_rids:
+                    fwds[0] = False
+                if rids[1] in cur_node.producing_rids:
+                    fwds[1] = False
+
+                for idx in range(2):
+                    if rids[idx] in cur_node.consuming_rids:
+                        stoics[idx] = cur_node.consuming_rids[rids[idx]]['stoic']
+                        nbs[idx] = self.get_neighbor_node(self.ng.edges[rids[idx]].products)
                     else:
-                        nb_sids[idx] = self.get_neighbor_node(self.ng.edges[rid].reactants)
+                        stoics[idx] = cur_node.producing_rids[rids[idx]]['stoic']
+                        nbs[idx] = self.get_neighbor_node(self.ng.edges[rids[idx]].reactants)
 
-                # if one of the reversible reactions connects to network, let it be the producing side
-                if nb_sids[0] is None:
-                    p_rid, c_rid = rids
-                    nb_us, nb_ds = nb_sids
+                # preferably the producing leg connects to the network (nb = None)
+                if nbs[0] is None:
+                    producing_leg = ReactionLeg(rids[0], stoics[0], nbs[0], fwds[0])
+                    consuming_leg = ReactionLeg(rids[1], stoics[1], nbs[1], fwds[1])
                 else:
-                    c_rid, p_rid = rids
-                    nb_ds, nb_us = nb_sids
+                    producing_leg = ReactionLeg(rids[1], stoics[1], nbs[1], not fwds[1])
+                    consuming_leg = ReactionLeg(rids[0], stoics[0], nbs[0], not fwds[0])
 
             # get neighbors of a node in a directional pathway segment
             else:
@@ -221,29 +268,44 @@ class LinearPathway:
                 if drid in cur_node.consuming_rids:
                     c_rid = drid
                     p_rid = orid
+                    c_stoic = cur_node.consuming_rids[c_rid]['stoic']
                     nb_ds = self.get_neighbor_node(self.ng.edges[c_rid].products)
+                    c_fwd = True
 
+                    # the 'producing' reaction can be a producing or a reversible consuming reaction
                     if p_rid in cur_node.producing_rids:
+                        p_stoic = cur_node.producing_rids[p_rid]['stoic']
                         nb_us = self.get_neighbor_node(self.ng.edges[p_rid].reactants)
-                    # in case the 2nd reaction is also a consuming reaction, it must be reversible
+                        p_fwd = True
                     else:
-                        assert cur_node.consuming_rids[p_rid] is True
+                        # case of a consuming reaction run in reverse
+                        assert cur_node.consuming_rids[p_rid]['rev'] is True
+                        p_stoic = cur_node.consuming_rids[p_rid]['stoic']
                         nb_us = self.get_neighbor_node(self.ng.edges[p_rid].products)
+                        p_fwd = False
 
                 # case, where the selected directional reaction is a producing reaction
                 else:
                     p_rid = drid
                     c_rid = orid
+                    p_stoic = cur_node.producing_rids[p_rid]['stoic']
                     nb_us = self.get_neighbor_node(self.ng.edges[p_rid].reactants)
+                    p_fwd = True
 
                     if c_rid in cur_node.consuming_rids:
+                        c_stoic = cur_node.consuming_rids[c_rid]['stoic']
                         nb_ds = self.get_neighbor_node(self.ng.edges[c_rid].products)
+                        c_fwd = True
                     else:
-                        assert cur_node.producing_rids[c_rid] is True
+                        assert cur_node.producing_rids[c_rid]['rev'] is True
+                        c_stoic = cur_node.producing_rids[c_rid]['stoic']
                         nb_ds = self.get_neighbor_node(self.ng.edges[c_rid].reactants)
+                        c_fwd = False
 
-        return {'sid': init_sid, 'p_rid': p_rid, 'nb_us': nb_us,
-                'c_rid': c_rid, 'nb_ds': nb_ds, 'reversible': rev}
+                producing_leg = ReactionLeg(p_rid, p_stoic, nb_us, p_fwd)
+                consuming_leg = ReactionLeg(c_rid, c_stoic, nb_ds, c_fwd)
+
+        return {'sid': init_sid, 'producing': producing_leg, 'consuming': consuming_leg, 'reversible': rev}
 
     def get_ds_neighbor(self, us_pwy_seg):
         """get downstream neighbor at us_sid coming from p_rid.
@@ -263,51 +325,92 @@ class LinearPathway:
           the upstream pathway is reversible (leading to a directional pathway
           in opposite direction)
         """
-        sid = us_pwy_seg.nb_ds
+        sid = us_pwy_seg.consuming.nb
         cur_node = self.ng.nodes[sid]
-        revs = [rev for rev in (cur_node.consuming_rids | cur_node.producing_rids).values()]
+        revs = [val['rev'] for val in (cur_node.consuming_rids | cur_node.producing_rids).values()]
         rev = bool(np.all(np.array(revs)))
-        c_rid = None
-        nb_ds = None
+        producing_leg = None
+        consuming_leg = None
 
-        if cur_node.connectivity == 2:
+        # collect information for upstream leg
+        p_rid = us_pwy_seg.consuming.rid
+        us_nb = us_pwy_seg.sid
+        if p_rid in cur_node.producing_rids:
+            p_stoic = cur_node.producing_rids[p_rid]['stoic']
+            p_fwd = True
+        else:
+            assert cur_node.consuming_rids[p_rid]['rev'] is True
+            p_stoic = cur_node.consuming_rids[p_rid]['stoic']
+            p_fwd = False
 
-            p_rid = us_pwy_seg.c_rid
-            c_rid = [rid for rid in (cur_node.consuming_rids | cur_node.producing_rids) if rid != p_rid][0]
+        if cur_node.connectivity == 1:
+            producing_leg = ReactionLeg(p_rid, p_stoic, us_nb, p_fwd)
 
+        if cur_node.connectivity == 2:
+            c_rid = [rid for rid in (cur_node.consuming_rids | cur_node.producing_rids) if rid != p_rid][0]
             if c_rid in cur_node.consuming_rids:
-                nb_ds = self.get_neighbor_node(self.ng.edges[c_rid].products)
+                c_stoic = cur_node.consuming_rids[c_rid]['stoic']
+                ds_nb = self.get_neighbor_node(self.ng.edges[c_rid].products)
+                c_fwd = True
             else:
-                nb_ds = self.get_neighbor_node(self.ng.edges[c_rid].reactants)
-                if cur_node.producing_rids[c_rid] is False:
-                    # flow reversal in case the producing reaction is directional
-                    return {'sid': sid, 'p_rid': c_rid, 'nb_us': nb_ds, 'reversible': rev}
+                # a 'producing' reaction consuming the metabolite
+                c_stoic = cur_node.producing_rids[c_rid]['stoic']
+                ds_nb = self.get_neighbor_node(self.ng.edges[c_rid].reactants)
+                c_fwd = False
+
+            # check for reversal in case the producing reaction is directional
+            if (c_rid in cur_node.consuming_rids) or (cur_node.producing_rids[c_rid]['rev'] is True):
+                producing_leg = ReactionLeg(p_rid, p_stoic, us_nb, p_fwd)
+                consuming_leg = ReactionLeg(c_rid, c_stoic, ds_nb, c_fwd)
+            else:
+                producing_leg = ReactionLeg(c_rid, c_stoic, ds_nb, not c_fwd)
+                consuming_leg = ReactionLeg(p_rid, p_stoic, us_nb, not p_fwd)
 
-        return {'sid': sid, 'c_rid': c_rid, 'nb_ds': nb_ds, 'reversible': rev}
+        return {'sid': sid, 'producing': producing_leg, 'consuming': consuming_leg, 'reversible': rev}
 
     def get_us_neighbor(self, ds_pwy_seg):
 
-        sid = ds_pwy_seg.nb_us
+        sid = ds_pwy_seg.producing.nb
         cur_node = self.ng.nodes[sid]
-        revs = [rev for rev in (cur_node.consuming_rids | cur_node.producing_rids).values()]
+        revs = [val['rev'] for val in (cur_node.consuming_rids | cur_node.producing_rids).values()]
         rev = bool(np.all(np.array(revs)))
-        p_rid = None
-        nb_us = None
+        producing_leg = None
+        consuming_leg = None
 
-        if cur_node.connectivity == 2:
+        # collect information for downstream leg
+        c_rid = ds_pwy_seg.producing.rid
+        ds_nb = ds_pwy_seg.sid
+        if c_rid in cur_node.consuming_rids:
+            c_stoic = cur_node.consuming_rids[c_rid]['stoic']
+            c_fwd = True
+        else:
+            assert cur_node.producing_rids[c_rid]['rev'] is True, f'{sid},{c_rid}'
+            c_stoic = cur_node.producing_rids[c_rid]['stoic']
+            c_fwd = False
 
-            c_rid = ds_pwy_seg.p_rid
-            p_rid = [rid for rid in (cur_node.consuming_rids | cur_node.producing_rids) if rid != c_rid][0]
+        if cur_node.connectivity == 1:
+            consuming_leg = ReactionLeg(c_rid, c_stoic, ds_nb, c_fwd)
 
+        if cur_node.connectivity == 2:
+            p_rid = [rid for rid in (cur_node.consuming_rids | cur_node.producing_rids) if rid != c_rid][0]
             if p_rid in cur_node.producing_rids:
-                nb_us = self.get_neighbor_node(self.ng.edges[p_rid].reactants)
+                p_stoic = cur_node.producing_rids[p_rid]['stoic']
+                us_nb = self.get_neighbor_node(self.ng.edges[p_rid].reactants)
+                p_fwd = True
             else:
-                nb_us = self.get_neighbor_node(self.ng.edges[p_rid].products)
-                if cur_node.consuming_rids[p_rid] is False:
-                    # flow reversal in case the consuming reaction is directional
-                    return {'sid': sid, 'c_rid': p_rid, 'nb_ds': nb_us, 'reversible': rev}
+                p_stoic = cur_node.consuming_rids[p_rid]['stoic']
+                us_nb = self.get_neighbor_node(self.ng.edges[p_rid].products)
+                p_fwd = False
+
+            # check for reversal in case the producing reaction is directional
+            if (p_rid in cur_node.producing_rids) or (cur_node.consuming_rids[p_rid]['rev'] is True):
+                producing_leg = ReactionLeg(p_rid, p_stoic, us_nb, p_fwd)
+                consuming_leg = ReactionLeg(c_rid, c_stoic, ds_nb, c_fwd)
+            else:
+                producing_leg = ReactionLeg(c_rid, c_stoic, ds_nb, not c_fwd)
+                consuming_leg = ReactionLeg(p_rid, p_stoic, us_nb, not p_fwd)
 
-        return {'sid': sid, 'p_rid': p_rid, 'nb_us': nb_us, 'reversible': rev}
+        return {'sid': sid, 'producing': producing_leg, 'consuming': consuming_leg, 'reversible': rev}
 
     def get_reactions(self):
 

modelpruner/graph/linear_pathways.py

-"""
+"""Implementation of LinearPathways Class.
 
-Peter Schubert, October 2022
+Peter Schubert, HHU Duesseldorf, October 2022
 """
 
 from modelpruner.graph.linear_pathway import LinearPathway
@@ -8,7 +8,14 @@ from modelpruner.graph.linear_pathway import LinearPathway
 
 class LinearPathways:
 
-    def __init__(self, network_graph, biomass_rid=None):
+    # TODO: rework again to first get lin_pathways from each of the nodes (excluding biomass precursors)
+    #  idea to get as long pathways as possible
+    #  sort for longest pathways
+    #  remove duplicated pathways
+    #  accept longest pathways first and drop pathways that contain already consumed reactions
+    #  check if c2_nodes remain, and for these request again pathways
+
+    def __init__(self, network_graph, protected_rids=None):
         """
 
         Construct linear pathways from metabolites with a connectivity of 2
@@ -19,12 +26,17 @@ class LinearPathways:
 
         self.ng = network_graph
 
-        if biomass_rid in self.ng.edges:
-            biomass_sids = self.ng.edges[biomass_rid].reactants | self.ng.edges[biomass_rid].products
-        else:
-            biomass_sids = set()
-        self.c2_sids = {node.id: None for node in self.ng.nodes.values()
-                        if node.id not in biomass_sids and node.is_c2_node()}
+        # TODO: biomass protected, external metabolites protected
+        # metabolites connected to protected reactions get protected
+        protected_sids = set()
+        if protected_rids is not None:
+            for rid in protected_rids:
+                if rid in self.ng.edges:
+                    protected_sids |= set(self.ng.edges[rid].reactants.keys())
+                    protected_sids |= set(self.ng.edges[rid].products.keys())
+
+        self.c2_sids = {node.sid: None for node in self.ng.nodes.values()
+                        if node.sid not in protected_sids and node.is_c2_node()}
         self.pwy_rids = {}
         self.pathways = []
 
@@ -45,8 +57,8 @@ class LinearPathways:
     def pwy_substrates(self, lin_pwy):
         sids = set()
         for pwy_seg in lin_pwy.segments:
-            sids |= set(self.ng.edges[pwy_seg.p_rid].reactants)
-            sids |= set(self.ng.edges[pwy_seg.p_rid].products)
-            sids |= set(self.ng.edges[pwy_seg.c_rid].reactants)
-            sids |= set(self.ng.edges[pwy_seg.c_rid].products)
+            sids |= set(self.ng.edges[pwy_seg.producing.rid].reactants.keys())
+            sids |= set(self.ng.edges[pwy_seg.producing.rid].products.keys())
+            sids |= set(self.ng.edges[pwy_seg.consuming.rid].reactants.keys())
+            sids |= set(self.ng.edges[pwy_seg.consuming.rid].products.keys())
         return sids

modelpruner/graph/metabolite_node.py

+"""Implementation of MetaboliteNode Class.
+
+Peter Schubert, HHU Duesseldorf, October 2022
+"""
 
 import numpy as np
 
@@ -5,17 +9,17 @@ import numpy as np
 class MetaboliteNode:
 
     def __init__(self, sid):
-        self.id = sid
+        self.sid = sid
         self.producing_rids = {}
         self.consuming_rids = {}
         self.connectivity = 0
 
-    def add_producing_rid(self, rid, reversible):
-        self.producing_rids[rid] = reversible
+    def add_producing_rid(self, rid, stoic, reversible):
+        self.producing_rids[rid] = {'stoic': stoic, 'rev': reversible}
         self.connectivity = len(self.producing_rids) + len(self.consuming_rids)
 
-    def add_consuming_rid(self, rid, reversible):
-        self.consuming_rids[rid] = reversible
+    def add_consuming_rid(self, rid, stoic, reversible):
+        self.consuming_rids[rid] = {'stoic': stoic, 'rev': reversible}
         self.connectivity = len(self.producing_rids) + len(self.consuming_rids)
 
     def is_c2_node(self):
@@ -23,10 +27,10 @@ class MetaboliteNode:
             if len(self.producing_rids) == 1 and len(self.consuming_rids) == 1:
                 return True
             if len(self.producing_rids) == 2:
-                if bool(np.any(np.array([val for val in self.producing_rids.values()]))) is True:
+                if bool(np.any(np.array([val['rev'] for val in self.producing_rids.values()]))) is True:
                     return True
             if len(self.consuming_rids) == 2:
-                if bool(np.any(np.array([val for val in self.consuming_rids.values()]))) is True:
+                if bool(np.any(np.array([val['rev'] for val in self.consuming_rids.values()]))) is True:
                     return True
         elif self.connectivity == 1:  # blocked metabolite
             return True

modelpruner/graph/network_graph.py

-"""
+"""Implementation of NetworkGraph Class.
 
-Peter Schubert, October 2022
+Peter Schubert, HHU Duesseldorf, October 2022
 """
 
 from modelpruner.graph.metabolite_node import MetaboliteNode
@@ -25,8 +25,8 @@ class NetworkGraph:
             rid = fba_model.rids[col]
             reversible = bool(fba_model.reversible[col])
             if stoic > 0:
-                self.edges[rid].add_product(sid)
-                self.nodes[sid].add_producing_rid(rid, reversible)
+                self.edges[rid].add_product(sid, stoic)
+                self.nodes[sid].add_producing_rid(rid, stoic, reversible)
             else:
-                self.edges[rid].add_reactant(sid)
-                self.nodes[sid].add_consuming_rid(rid, reversible)
+                self.edges[rid].add_reactant(sid, -stoic)
+                self.nodes[sid].add_consuming_rid(rid, -stoic, reversible)

modelpruner/graph/pathway_segment.py

+"""Implementation of PathwaySegment Class.
+
+Peter Schubert, HHU Duesseldorf, October 2022
 """
 
 
-Peter Schubert, October 2022
-"""
+class ReactionLeg:
+
+    def __init__(self, rid, stoic, nb, forward):
+        self.rid = rid
+        self.stoic = stoic
+        self.nb = nb
+        self.forward = forward
 
 
 class PathwaySegment:
 
-    def __init__(self, sid):
-        self.sid = sid
-        self.reversible = False
-        self.p_rid = None
-        self.nb_us = None
-        self.c_rid = None
-        self.nb_ds = None
+    def __init__(self, connectivity):
+        self.sid = connectivity['sid']
+        self.reversible = connectivity.get('reversible', False)
+        self.producing = connectivity.get('producing')
+        self.consuming = connectivity.get('consuming')
 
-    def set_upstream(self, rid, sid):
-        self.p_rid = rid
-        self.nb_us = sid
+    def set_producing(self, reaction_leg):
+        """Set the producing reaction leg.
 
-    def set_downstream(self, rid, sid):
-        self.c_rid = rid
-        self.nb_ds = sid
+        A reversible 'consuming' reaction can be a producing reaction
+
+        :param reaction_leg: producing reaction leg with stoic,
+        :type reaction_leg: ReactionLeg
+        """
+        self.producing = reaction_leg
+
+    def set_consuming(self, reaction_leg):
+        self.consuming = reaction_leg
 
     def set_reversible(self, reversible):
         """Set reversible flag of pathway segment.
@@ -37,12 +48,12 @@ class PathwaySegment:
         self.reversible = reversible
 
     def reverse_direction(self):
-        temp_rid = self.p_rid
-        temp_sid = self.nb_us
-        self.p_rid = self.c_rid
-        self.nb_us = self.nb_ds
-        self.c_rid = temp_rid
-        self.nb_ds = temp_sid
+
+        tmp_leg = self.producing
+        self.producing = self.consuming
+        self.consuming = tmp_leg
+        self.producing.forward = not self.producing.forward
+        self.consuming.forward = not self.consuming.forward
 
     def get_reactions(self):
         """Get the reactions connecting the pathway segments
@@ -51,8 +62,8 @@ class PathwaySegment:
         :rtype: set of str
         """
         rids = set()
-        if self.p_rid is not None:
-            rids.add(self.p_rid)
-        if self.c_rid is not None:
-            rids.add(self.c_rid)
+        if type(self.producing) is ReactionLeg:
+            rids.add(self.producing.rid)
+        if type(self.consuming) is ReactionLeg:
+            rids.add(self.consuming.rid)
         return rids

modelpruner/graph/reaction_edge.py

+"""Implementation of ReactionEdge Class.
+
+Peter Schubert, HHU Duesseldorf, October 2022
+"""
+
 
 class ReactionEdge:
 
     def __init__(self, rid, reversible):
 
-        self.id = rid
+        self.rid = rid
         self.reversible = reversible
-        self.reactants = set()
-        self.products = set()
+        self.reactants = {}
+        self.products = {}
 
-    def add_product(self, sid):
-        self.products.add(sid)
+    def add_product(self, sid, stoic):
+        self.products[sid] = stoic
 
-    def add_reactant(self, sid):
-        self.reactants.add(sid)
+    def add_reactant(self, sid, stoic):
+        self.reactants[sid] = stoic