new gpr parsing strategy. still gives infeasible models

Scripts/protein_costs_analysis.py

@@ -1,9 +1,8 @@
 import os
 import pandas as pd
 
-DATA_DIR = os.path.join(os.path.split(os.getcwd())[0], 'Data')
-PAM_DATA_FILE_PATH = os.path.join(DATA_DIR, 'proteinAllocationModel_iML1515_EnzymaticData_py.xls')
-RESULT_DATA_FILE = os.path.join(os.path.split(os.getcwd())[0], 'Results', 'protein_costs_glycolysis_tca.xlsx')
+PAM_DATA_FILE_PATH = os.path.join('Data', 'proteinAllocationModel_iML1515_EnzymaticData_py.xls')
+RESULT_DATA_FILE = os.path.join('Results', 'protein_costs_glycolysis_tca.xlsx')
 
 active_enzyme_info = pd.read_excel(PAM_DATA_FILE_PATH, sheet_name='ActiveEnzymes')
 
@@ -15,7 +14,9 @@
 protein_efficiency_df = active_enzyme_info.filter(['rxnID','molMass','kcat'])
 protein_efficiency_df = protein_efficiency_df.assign(rxnName = lambda x: x.rxnID.str.rstrip('_fb'))
 protein_efficiency_glyc_tca = protein_efficiency_df[protein_efficiency_df.rxnName.isin(tca_glycolysis_reactions)]
-protein_efficiency_glyc_tca  = protein_efficiency_glyc_tca.assign(relEfficiency = lambda x: x['molMass']/x['kcat']).sort_values(by = 'relEfficiency', ascending=False)
+protein_efficiency_glyc_tca  = protein_efficiency_glyc_tca.assign(kcat = lambda x:x['kcat']*3600,
+                                                                  relEfficiency = lambda x: x['molMass']/x['kcat']
+                                                                  ).sort_values(by = 'relEfficiency', ascending=False)
 print(protein_efficiency_glyc_tca)
 
 protein_efficiency_glyc_tca.to_excel(RESULT_DATA_FILE, sheet_name='protein_costs')

Scripts/toy_ec_pam.py

@@ -4,9 +4,9 @@
 import numpy as np
 
 #importing the tools from the PAModelpy package
-from PAModelpy.EnzymeSectors import ActiveEnzymeSector, UnusedEnzymeSector, TransEnzymeSector
-from PAModelpy.PAModel import PAModel
-from PAModelpy.configuration import Config
+from src.PAModelpy.EnzymeSectors import ActiveEnzymeSector, UnusedEnzymeSector, TransEnzymeSector
+from src.PAModelpy.PAModel import PAModel
+from src.PAModelpy.configuration import Config
 
 Config.BIOMASS_REACTION = 'R7'
 Config.GLUCOSE_EXCHANGE_RXNID = 'R1'
@@ -104,14 +104,21 @@ def build_toy_gem():
 def build_active_enzyme_sector(Config):
     kcat_fwd = [1, 0.5, 1, 1, 0.5 ,0.45, 1.5]  # High turnover for exchange reactions
     kcat_rev = [kcat for kcat in kcat_fwd]
+    protein2gene = {}
     rxn2kcat = {}
     for i in range(n-3): # all reactions have an enzyme, except excretion reactions
         rxn_id = f'R{i+1}'
         # 1e-6 to correct for the unit transformation in the model (meant to make the calculations preciser for different unit dimensions)
         #dummy molmass like in MATLAB script
-        rxn2kcat = {**rxn2kcat, **{rxn_id: {f'E{i+1}':{'f': kcat_fwd[i]/(3600*1e-6), 'b': kcat_rev[i]/(3600*1e-6), 'molmass': 1e6}}}}
-
-    return ActiveEnzymeSector(rxn2protein = rxn2kcat, configuration=Config)
+        rxn2kcat = {**rxn2kcat, **{rxn_id:
+                                       {f'E{i+1}':
+                                            {'f': kcat_fwd[i]/(3600*1e-6),
+                                             'b': kcat_rev[i]/(3600*1e-6),
+                                             'molmass': 1e6,
+                                             'protein_reaction_association':[[f'E{i+1}']]}}}}
+        protein2gene = {**protein2gene,**{f'E{i+1}':[[f'g{i+1}']]}}
+
+    return ActiveEnzymeSector(rxn2protein = rxn2kcat, protein2gene = protein2gene, configuration=Config)
 
 def build_unused_protein_sector(Config):
     return UnusedEnzymeSector(id_list = ['R1'], ups_mu=[-0.01*1e-3], ups_0=[0.1*1e-3], mol_mass= [1], configuration=Config)

src/PAModelpy/CatalyticEvent.py

@@ -4,14 +4,15 @@
 (e.g. one conversion of substrate to product, can be catalyzed by multiple enzymes)
 """
 import cobra
-from cobra import DictList, Object
+from cobra import DictList, Object, Reaction
 from cobra.exceptions import OptimizationError
 from cobra.util.solver import check_solver_status
 from optlang.symbolics import Zero
-from typing import Optional, Dict
+from typing import Optional, Dict, Union
 from warnings import warn
 from copy import copy, deepcopy
 
+
 class CatalyticEvent(Object):
     """
         CatalyticEvent is a class for holding information regarding the
@@ -49,6 +50,7 @@ def __init__(
         #relation to reaction
         self.rxn_id = rxn_id
         self.rxn = None
+        self.catalytic_reactions = DictList()
 
         #relation to enzymes
         self.kcats = kcats2enzymes
@@ -155,11 +157,22 @@ def model(self, model):
             rxn = cobra.Reaction(id = self.rxn_id)
             self._model.add_reactions([rxn])
             self.rxn = rxn
+        #add reaction constraint
+        ce_constraint = self._model.problem.Constraint(Zero, name= self.id, lb = 0, ub=0)
+        self._model.add_cons_vars([ce_constraint])
+        self._model.constraints[self.id].set_linear_coefficients({
+            self.rxn.forward_variable: 1,
+            self.rxn.reverse_variable: -1,
+        })
+        self.constraints[self.id] = ce_constraint
+
         #add enzymes to the model if they are not there already
         for enzyme in self.enzymes:
             if enzyme in self._model.enzymes:
                 self.constraints ={**self.constraints, **enzyme._constraints}
+                enzyme._constraints[self.id] = ce_constraint
                 enzyme_model = self._model.enzymes.get_by_id(enzyme.id)
+
                 if self.rxn_id in enzyme_model.rxn2kcat.keys():
                     if self not in enzyme_model.catalytic_events:
                         enzyme_model.add_catalytic_event(self, kcats= {enzyme: enzyme.rxn2kcat[self.rxn_id]})
@@ -168,6 +181,11 @@ def model(self, model):
             else:
                 self._model.add_enzymes([enzyme])
 
+            # add enzyme-reaction association
+            self.add_enzyme_reaction_association(enzyme)
+
+
+
     def add_enzymes(self,enzyme_kcat_dict: dict):
         """
         Add enzymes to the catalytic event and create bindings to the related model.
@@ -184,10 +202,8 @@ def add_enzymes(self,enzyme_kcat_dict: dict):
         for enzyme, kcat in enzyme_kcat_dict.items():
             # check if the enzyme is already associated to the catalytic event
             if enzyme in self.enzymes:
-                # print(enzyme)
-                # warn(f'Enzyme {enzyme.id} is already associated with catalytic event {self.id}. This enzyme will be updated')
-                self.change_kcat_values({enzyme.id: enzyme.get_kcat_values(rxn_ids = [self.rxn_id])})
-                # print(self.enzymes)
+                self.add_enzyme_reaction_association(enzyme)
+                self.change_kcat_values({enzyme.id: enzyme.get_kcat_values(rxn_ids = [self.rxn_id+"_"+enzyme.id])})
                 continue
 
             self.enzymes.append(enzyme)
@@ -213,23 +229,48 @@ def add_enzymes(self,enzyme_kcat_dict: dict):
             #add constraints to the catalytic event
             self.constraints = {**self.constraints, **enzyme._constraints}
 
+            #add reaction-protein association
+            self.add_enzyme_reaction_association(enzyme)
+
             #connect the enzyme variable to the enzyme in the model and the reaction
             for direction, kcatvalue in kcat.items():
-                coeff = kcatvalue * 3600 * 1e-6
                 #add enzyme to the associated reaction with kinetic constants
                 #and relate enzyme to the catalytic event
                 if direction == 'f':
                     self.constraints[f'EC_{enzyme.id}_{direction}'].set_linear_coefficients({
-                        self.rxn.forward_variable: 1/coeff,
                         enzyme_var.forward_variable: -1
                         })
 
                 elif direction == 'b':
                     self.constraints[f'EC_{enzyme.id}_{direction}'].set_linear_coefficients({
-                        self.rxn.reverse_variable: 1/coeff,
                         enzyme_var.reverse_variable: -1
                         })
 
+    def add_enzyme_reaction_association(self, enzyme):
+        catalytic_reaction = Reaction(self.id + "_" + enzyme.id)
+        self._model.add_reactions([catalytic_reaction])
+        self.catalytic_reactions.append(catalytic_reaction)
+        self._model.constraints[self.id].set_linear_coefficients({
+            catalytic_reaction.forward_variable: -1,
+            catalytic_reaction.reverse_variable: 1,
+        })
+        self.add_catalytic_reaction_to_enzyme_constraint(catalytic_reaction, enzyme)
+
+    def add_catalytic_reaction_to_enzyme_constraint(self, catalytic_reaction:Reaction,
+                                                    enzyme):
+        #remove relation to metabolic reaction
+        for dir in ['f', 'b']:
+            self._model._change_kcat_in_enzyme_constraint(self.rxn, enzyme.id,
+                                                      dir, 0)
+        kcat_dict = self.kcats[enzyme]
+        for direction, kcatvalue in kcat_dict.items():
+            self._model._change_kcat_in_enzyme_constraint(catalytic_reaction, enzyme.id,
+                                                          direction, kcatvalue)
+
+        #change rxn2kcat dict for correct referencing
+        del enzyme.rxn2kcat[self.rxn_id]
+        enzyme.rxn2kcat[catalytic_reaction.id] = kcat_dict
+
 
     def remove_enzymes(self, enzyme_list: list):
         """
@@ -268,18 +309,6 @@ def remove_enzymes(self, enzyme_list: list):
                 self.enzyme_variables.remove(enzyme_var)
             #set coefficient in constraint to 0
             for constraint in set([cons for name, cons in self.constraints.items() if enz.id in name]):
-                # self.constraints[constraint.name] = constraint
-                # coeff = 0
-                # #set coefficients to 0
-                # if constraint.name[-1] == 'f':
-                #     constraint.set_linear_coefficients({
-                #         self.rxn.forward_variable: coeff
-                #         })
-                #
-                # elif constraint.name[-1] == 'b':
-                #     constraint.set_linear_coefficients({
-                #         self.rxn.reverse_variable: coeff
-                #         })
                 if constraint in self._model.constraints.values():
                     self._model.remove_cons_vars([constraint])
                 #remove constraint from list of r=constraints
@@ -330,25 +359,12 @@ def change_kcat_values(self, enzyme_kcat_dict : dict):
 
             enzyme_var = self._model.enzyme_variables.get_by_id(enzyme)
             if enzyme_var not in self.enzyme_variables: self.enzyme_variables.add(enzyme_var)
-            enzyme_obj = self._model.enzymes.get_by_id(enzyme)
-            enzyme_var.change_kcat_values({self.rxn: kcat_dict})
+            ce_reaction = self._model.reactions.get_by_id(self.rxn_id+"_"+enzyme.id)
+            enzyme_var.change_kcat_values({ce_reaction: kcat_dict})
             for direction, kcat in kcats_change.items():
                 #change enzyme variable
-                enzyme_var.kcats[self.rxn_id][direction] = kcat
-                # get constraint
-                constraint_id = f'EC_{enzyme}_{direction}'
-                constraint = enzyme_obj._constraints[constraint_id]
-                # change kcat value in the constraint
-                coeff =  kcat * 3600 * 1e-6
-                if direction == 'f':
-                    self._model.constraints[constraint_id].set_linear_coefficients({
-                        self.rxn.forward_variable: 1/coeff
-                        })
-                elif direction == 'b':
-                    self._model.constraints[constraint_id].set_linear_coefficients({
-                        self.rxn.reverse_variable: 1/coeff
-                        })
-            self._model.solver.update()
+                enzyme_var.kcats[ce_reaction.id][direction] = kcat
+                self._model._change_kcat_in_enzyme_constraint(ce_reaction, enzyme,direction, kcat)
 
     def __copy__(self) -> 'CatalyticEvent':
         """ Copy the CatalyticEvent