kChip CommScores simulation and CommKinetics updates

modelseedpy/core/msminimalmedia.py

@@ -4,7 +4,6 @@
 from modelseedpy.core.msmodelutl import MSModelUtil
 from modelseedpy.fbapkg.basefbapkg import BaseFBAPkg
 from itertools import combinations, permutations, chain
-from optlang import Variable, Constraint
 from cobra.medium import minimal_medium
 from optlang.symbolics import Zero
 from math import isclose, inf, factorial
@@ -123,10 +122,10 @@ def _min_consumption_objective(model_util, interacting):
                 print(f"The {var_name} variable already exists in "
                       f"{model_util.model.id} and thus is skipped.\n")
                 continue
-            vars[rxn.id] = Variable(var_name, lb=0, ub=1, type="binary")
+            vars[rxn.id] = model_util.model.problem.Variable(var_name, lb=0, ub=1, type="binary")
             model_util.add_cons_vars([vars[rxn.id]])
             # bin_flux: {rxn_bin}*1000 >= {rxn_rev_flux}
-            model_util.create_constraint(Constraint(Zero, lb=0, ub=None, name=cons_name),
+            model_util.create_constraint(model_util.model.problem.Constraint(Zero, lb=0, ub=None, name=cons_name),
                                          coef={vars[rxn.id]: 1000, rxn.reverse_variable: -1})
         return vars
 
@@ -181,7 +180,7 @@ def minimize_components(org_model, min_growth=0.1, environment=None,
             sol_media = _exchange_solution(sol_rxns_dict)
             min_media = sol_media if len(sol_media) < len(min_media) else min_media
             ## omit the solution from future searches
-            model_util.create_constraint(Constraint(
+            model_util.create_constraint(model_util.model.problem.Constraint(
                 Zero, lb=None, ub=len(sol_dict)-1, name=f"exclude_sol{len(solution_dicts)}"), sol_dict)
 
             # search the permutation space by omitting previously investigated solution_dicts
@@ -219,8 +218,7 @@ def _knockout(org_model, rxnVar, variables, sol_dict, sol_index, interacting):
                          for rxnVar2 in sol_dict if (
                                  rxnVar != rxnVar2 and any(["_e0" in met.id for met in rxnVar2.metabolites]))
                          })
-        knocked_model_utl.create_constraint(Constraint(Zero, lb=0.1, ub=None,
-                                                       name=f"{rxnVar.name}-sol{sol_index}"), coef)
+        knocked_model_utl.create_constraint(org_model.problem.Constraint(Zero, lb=0.1, ub=None, name=f"{rxnVar.name}-sol{sol_index}"), coef)
         return knocked_model_utl.optimize()
 
     @staticmethod

modelseedpy/core/msmodelutl.py

@@ -13,7 +13,6 @@
 from modelseedpy.core.fbahelper import FBAHelper
 from itertools import chain
 from optlang.symbolics import Zero
-from optlang import Constraint, Objective
 from math import isclose
 from multiprocessing import Value
 
@@ -650,7 +649,7 @@ def remove_cons_vars(self, vars_cons):
         self.model.solver.update()
 
     def add_objective(self, objective, direction="max", coef=None):
-        self.model.objective = Objective(objective, direction=direction)
+        self.model.objective = self.model.problem.Objective(objective, direction=direction)
         self.model.solver.update()
         if coef:
             self.model.objective.set_linear_coefficients(coef)
@@ -668,16 +667,16 @@ def biomass_expression(self):
                 # returns the biomass expression of the lowest cytoplasmic compartment
                 return met.constraint.expression
 
-    def add_minimal_objective_cons(self, min_value=0.1, objective_expr=None):
-        if "min_value" not in self.model.constraints:
+    def add_minimal_objective_cons(self, min_value=0.1, objective_expr=None, name="min_value"):
+        if name not in self.model.constraints:
             objective_expr = objective_expr or self.model.objective.expression
-            self.create_constraint(Constraint(objective_expr, lb=min_value, ub=None, name="min_value"))
+            self.create_constraint(self.model.problem.Constraint(objective_expr, lb=min_value, ub=None, name=name))
             # print(self.model.constraints["min_value"])
         else:
-            print(f"The min_value constraint already exists in {self.model.id}, "
+            print(f"The {name} constraint already exists in {self.model.id}, "
                   f"hence the lb is simply updated from"
-                  f" {self.model.constraints['min_value'].lb} to {min_value}.\n")
-            self.model.constraints["min_value"].lb = min_value
+                  f" {self.model.constraints[name].lb} to {min_value}.\n")
+            self.model.constraints[name].lb = min_value
 
     def add_exchange_to_model(self, cpd, rxnID):
         self.model.add_boundary(metabolite=Metabolite(id=cpd.id, name=cpd.name, compartment="e0"),
@@ -1107,7 +1106,7 @@ def resource_balance_constraint(self, flux_limit=140):
         for rxn in self.model.reactions:
             if "EX_" not in rxn.id:
                 vars_coef[rxn.forward_variable] = vars_coef[rxn.reverse_variable] = 1
-        self.create_constraint(Constraint(Zero, lb=0, ub=flux_limit, name="resource_balance_limit"), coef=vars_coef)
+        self.create_constraint(self.model.problem.Constraint(Zero, lb=0, ub=flux_limit, name="resource_balance_limit"), coef=vars_coef)
 
     def apply_test_condition(self, condition, model=None):
         """Applies constraints and objective of specified condition to model
@@ -1646,7 +1645,7 @@ def add_atp_hydrolysis(self, compartment):
     def costless_excreta(self, pfba=False):
         # the double optimization is intended to truly find the maximal biomass growth
         original_objective = self.model.objective
-        minObj_cons = Constraint(self.model.objective.expression, lb=self.model.slim_optimize(), name="minObj")
+        minObj_cons = self.model.problem.Constraint(self.model.objective.expression, lb=self.model.slim_optimize(), name="minObj")
         self.add_cons_vars([minObj_cons])
         if pfba:
             self.model.problem.constraints.minObj_cons.lb = self.model.slim_optimize()

modelseedpy/fbapkg/basefbapkg.py

@@ -3,15 +3,15 @@
 from __future__ import absolute_import
 
 import logging
-import re  # !!! import is never used
-from optlang.symbolics import Zero, add  # !!! add is never used
-import json as _json  # !!! import is never used
+import re
+from optlang.symbolics import Zero, add
+import json as _json
 from cobra.core import (
     Gene,
     Metabolite,
     Model,
     Reaction,
-)  # !!! none of these imports are used
+)
 from modelseedpy.fbapkg.mspackagemanager import MSPackageManager
 from modelseedpy.core.msmodelutl import MSModelUtil
 from modelseedpy.core.exceptions import FeasibilityError
@@ -30,9 +30,7 @@ class BaseFBAPkg:
     Base class for FBA packages
     """
 
-    def __init__(
-        self, model, name, variable_types={}, constraint_types={}, reaction_types={}
-    ):
+    def __init__(self, model, name, variable_types={}, constraint_types={}, reaction_types={}):
         if isinstance(model, MSModelUtil):
             self.model = model.model
             self.modelutl = model
@@ -75,9 +73,7 @@ def clear(self):
             self.constraints[obj_type] = {}
         self.model.remove_cons_vars(cobra_objs)
 
-    def build_variable(
-        self, obj_type, lower_bound, upper_bound, vartype, cobra_obj=None
-    ):
+    def build_variable(self, obj_type, lower_bound, upper_bound, vartype, cobra_obj=None):
         name = None
         if self.variable_types[obj_type] == "none":
             count = len(self.variables[obj_type])
@@ -107,8 +103,7 @@ def build_constraint(
         if name in self.constraints[obj_type]:
             self.model.remove_cons_vars(self.constraints[obj_type][name])
         self.constraints[obj_type][name] = self.model.problem.Constraint(
-            Zero, lb=lower_bound, ub=upper_bound, name=name + "_" + obj_type
-        )
+            Zero, lb=lower_bound, ub=upper_bound, name=name + "_" + obj_type)
         self.model.add_cons_vars(self.constraints[obj_type][name])
         self.model.solver.update()
         if len(coef) > 0:

modelseedpy/fbapkg/commkineticpkg.py

@@ -7,19 +7,24 @@
 from modelseedpy.core.fbahelper import FBAHelper
 
 # Base class for FBA packages
-# TODO This code may need to refactored to not require an MSCommunity object, which is circular logic
 class CommKineticPkg(BaseFBAPkg):
     def __init__(self, model):
-        BaseFBAPkg.__init__(self, model, "community kinetics", {}, {"commkin": "string"})
+        BaseFBAPkg.__init__(self, model, "community kinetics", {}, {"commKin": "string"})
 
-    def build_package(self, kinetic_coef, community_model):
+    def build_package(self, kinetic_coef, community_model, probs=None):
         self.validate_parameters({}, [], {"kinetic_coef": kinetic_coef, "community": community_model})
-        for member in self.parameters["community"].members:
-            self.build_constraint(member)
+        for species in self.parameters["community"].members:
+            if species.id+"_commKin" in self.model.constraints:
+                print(f"Removing {species.id+'_commKin'} from {self.model.id}")
+                self.model.remove_cons_vars(self.model.constraints[species.id+"_commKin"])
+            self.build_constraint(species, probs)
 
-    def build_constraint(self, member):
-        bioRXN = member.primary_biomass
-        coef = {bioRXN.forward_variable: -self.parameters["kinetic_coef"], bioRXN.reverse_variable: self.parameters["kinetic_coef"]}
-        for reaction in member.reactions:
-            if reaction != bioRXN:   coef[reaction.forward_variable] = coef[reaction.reverse_variable] = 1
-        return BaseFBAPkg.build_constraint(self, "commkin", None, 0, coef, "Species" + str(member.index))
+    def build_constraint(self, species, probs):
+        coef = {species.primary_biomass.forward_variable: -1 * self.parameters["kinetic_coef"],
+                species.primary_biomass.reverse_variable: self.parameters["kinetic_coef"]}
+        for rxn in self.model.reactions:
+            rxnIndex = int(FBAHelper.rxn_compartment(rxn)[1:])
+            if (rxnIndex == species.index and rxn != species.primary_biomass):
+                val = 1 if not isinstance(probs, dict) else probs.get(rxn.id, 1)
+                coef[rxn.forward_variable] = coef[rxn.reverse_variable] = val
+        return BaseFBAPkg.build_constraint(self, "commKin", None, 0, coef, species.id)