From 0553b9f322fe5cb9fda5be76851c339920a1a987 Mon Sep 17 00:00:00 2001 From: Phil Date: Wed, 5 Aug 2026 13:58:34 -0400 Subject: [PATCH 1/2] fix: collect essential knock-in identifiers, not their cost values `essential_kis` was built from `cmp_ki_cost[er]`, which yields the intervention cost of each essential reaction rather than the reaction identifier. On e_coli_core with every reaction made a knock-in candidate the set came out as {1.0, 2.0, 3.0} -- the distinct lump costs. `SDProblem.__init__` then tested `r in self.essential_kis` with `r` a cobra Reaction, so the membership check could never succeed: 0 of 56 binaries were pinned, and reactions essential to a module stayed optional. Knock-in problems could therefore return solutions that omit a reaction the module cannot do without. Collect identifiers instead and compare against `r.id`. Measured on e_coli_core with a PROTECT module and all reactions as knock-in candidates: forced-on binaries go from 0 to 11. The default knockout path is unaffected -- `essential_kis` is only populated from reactions carrying a knock-in cost, so it stays empty when `ki_cost` is unused. MCS on e_coli_core (growth >= 0.001, max_cost 3) returns the identical 353 designs before and after, and the test suite passes (368 passed, 2 skipped, performance tests deselected). Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01RFtof9nZXFvCNoXz19po8C --- straindesign/compute_strain_designs.py | 2 +- straindesign/strainDesignProblem.py | 9 +++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/straindesign/compute_strain_designs.py b/straindesign/compute_strain_designs.py index 4f5c167..217e6bf 100644 --- a/straindesign/compute_strain_designs.py +++ b/straindesign/compute_strain_designs.py @@ -767,7 +767,7 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions: # remove ko-costs (and thus knockability) of essential reactions [cmp_ko_cost.pop(er) for er in essential_reacs if er in cmp_ko_cost] - essential_kis = set(cmp_ki_cost[er] for er in essential_reacs if er in cmp_ki_cost) + essential_kis = set(er for er in essential_reacs if er in cmp_ki_cost) # Build MILP kwargs1 = kwargs kwargs1[KOCOST] = cmp_ko_cost diff --git a/straindesign/strainDesignProblem.py b/straindesign/strainDesignProblem.py index 6bfe031..626a581 100644 --- a/straindesign/strainDesignProblem.py +++ b/straindesign/strainDesignProblem.py @@ -86,9 +86,10 @@ class SDProblem: the big-M method by default (with COBRA standard M=1000). M should be chosen 'sufficiently large' to avoid computational artifacts and 'sufficiently small' to avoid numerical issues. - essential_kis (optional (set)): - A set of reactions that are marked as addable and that are essential for at least one of the - strain design modules. Providing such "essential knock-ins" may speed up the strain design computation. + essential_kis (optional (set of str)): + Reaction identifiers that are marked as addable and that are essential for at least one of the + strain design modules. Their intervention binaries are fixed to 1, since a solution that omits + them cannot satisfy the module they are essential for. Returns: (SDProblem): @@ -167,7 +168,7 @@ def __init__(self, model: Model, sd_modules: List[SDModule], *args, **kwargs): else: self.b_ineq = [0.0, float(self.max_cost), np.inf] self.z_map_constr_ineq = sparse.csc_matrix((numr, 3)) - self.lb = [1.0 if r in self.essential_kis else 0.0 for r in model.reactions] + self.lb = [1.0 if r.id in self.essential_kis else 0.0 for r in model.reactions] self.ub = [1.0 - float(i) for i in self.z_non_targetable] self.idx_z = [i for i in range(0, numr)] self.c = [0.0] * numr From 91ac88e1f2f18058a55bf792b5c973677cb43695 Mon Sep 17 00:00:00 2001 From: Phil Date: Wed, 5 Aug 2026 14:51:00 -0400 Subject: [PATCH 2/2] fix: reject designs that zero a variable whose bounds exclude zero, in verify_sd verify_sd rebuilds the continuous part of the problem as a plain LP and decides feasibility from it. Variables an intervention touches are removed by dropping their columns, which equals v=0 only while every constraint that should still bite survives. That does not hold in general. prevent_boundary_knockouts turns a bound that excludes zero, such as ATPM >= 3.15, into a row with no z-mapping precisely so a knockout contradicts it rather than deleting it. For PROTECT modules, reassign_lb_ub_from_ineq then folds single-variable rows back into variable bounds, and those bounds disappear together with the column. The check saw an empty system and reported it feasible, which also made the "strain already meets the requirements" shortcut fire on knock-in problems and skip the MILP entirely. Test that case before building the LP: zeroing a variable whose bounds exclude zero contradicts the bound, so the region is empty whatever the rest of the system does. Column removal is kept for every other variable, since absence owes nothing to feasibility tolerances and is a stronger statement than an interval of [0, 0]. e_coli_core MCS (growth >= 0.001, max_cost 3) returns the identical 353 designs and still reports {ATPM} as a size-1 MCS. EFV enumeration via a PROTECT module with all reactions as knock-in candidates (max_cost 45) returns 607 designs, all members of the minimal-support antichain, with the two size-26 supports set-identical to those obtained independently as the Berge dual of the model's minimal cut sets. Test suite: 368 passed, 2 skipped. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01RFtof9nZXFvCNoXz19po8C --- straindesign/strainDesignMILP.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/straindesign/strainDesignMILP.py b/straindesign/strainDesignMILP.py index f7c603a..fdc6064 100644 --- a/straindesign/strainDesignMILP.py +++ b/straindesign/strainDesignMILP.py @@ -275,6 +275,17 @@ def verify_sd(self, sols) -> List: if np.logical_xor(sol[0,z_i],sense==-1) ] active_eqs = [i for i in range(self.cont_MILP.z_map_constr_eq.shape[1]) if i not in inactive_eqs] + # Zeroing a variable whose bounds exclude zero contradicts that bound, so the + # region is empty whatever the remaining system does. This has to be tested + # here rather than left to the LP: prevent_boundary_knockouts keeps such a + # bound (e.g. ATPM >= 3.15) as a row with no z-mapping so that a knockout + # contradicts it, but reassign_lb_ub_from_ineq later folds single-variable + # rows back into variable bounds, and those vanish together with the column. + if any(self.cont_MILP.lb[j] > 0.0 or self.cont_MILP.ub[j] < 0.0 for j in inactive_vars): + valid[i] = False + continue + # Otherwise drop the columns outright. Absence is a stronger statement than an + # interval of [0, 0], since it owes nothing to feasibility tolerances. lp = MILP_LP(A_ineq=self.cont_MILP.A_ineq[active_ineqs, :][:, active_vars], b_ineq=[self.cont_MILP.b_ineq[i] for i in active_ineqs], A_eq=self.cont_MILP.A_eq[active_eqs, :][:, active_vars],