StrainDesign Developer’s Guide

Reference paper: Schneider P., Bekiaris P. S., von Kamp A., Klamt S. — StrainDesign: a comprehensive Python package for computational design of metabolic networks. Bioinformatics, btac632 (2022). DOI: 10.1093/bioinformatics/btac632

Table of Contents

Architecture Overview
Repository Structure
Public API & Entry Points
Preprocessing Pipeline
- 4.1 Model Cleaning
- 4.2 Compression Pass 1
- 4.3 FVA
- 4.4 GPR Integration
- 4.5 Compression Pass 2
- 4.6 Network Compression — Algorithm Details
- 4.7 Regulatory Interventions
MILP Construction
- 5.1 Overview & z-map Matrices
- 5.2 Primal LP from COBRApy Model
- 5.3 SUPPRESS Module — Farkas’ Lemma Dual
- 5.4 PROTECT Module
- 5.5 OPTKNOCK — LP Duality for Bilevel Optimization
- 5.6 ROBUSTKNOCK — Three-Level Duality
- 5.7 OPTCOUPLE — Growth Coupling Potential
- 5.8 Linking Binary Variables: link_z
- 5.9 Big-M Bounding
- 5.10 Indicator Constraints vs. Big-M
Solver Backends
- 6.1 MILP_LP Unified Interface
- 6.2 CPLEX
- 6.3 Gurobi
- 6.4 SCIP
- 6.5 GLPK
Computation Approaches — SDMILP
- 7.1 ANY
- 7.2 BEST
- 7.3 POPULATE
- 7.4 Exclusion Constraints
- 7.5 Solution Verification
Post-processing & Result Container
- 8.1 Solution Decompression
- 8.2 GPR Translation
- 8.3 SDSolutions
Constraint & Expression Parsing
Supporting Utilities
- 10.1 LP Suppression (networktools.py)
Known Issues, Urgent Actions & Future Work
MILP Performance Notes
Testing

1. Architecture Overview

StrainDesign solves metabolic strain design problems: given a genome-scale metabolic model (COBRApy Model) and a set of biological objectives, find the minimal sets of gene/reaction knockouts (and additions) that force the network into a desired flux state.

All methods reduce to Mixed-Integer Linear Programs (MILP). The key intellectual contribution is the translation of nested/bilevel biological objectives into a flat MILP via LP duality and Farkas’ lemma, so a single MILP call suffices to find strain designs. The pipeline is:

User (model + SDModules)
        │
        ▼
   Preprocessing        networktools.py, compression.py, lptools.py
   ─ FVA cleanup
   ─ GPR extension
   ─ Network compression
        │
        ▼
   MILP Construction    strainDesignProblem.py  (SDProblem)
   ─ build_primal_from_cbm
   ─ farkas_dualize / LP_dualize
   ─ link_z (bounding + indicator constraints)
        │
        ▼
   Solver Execution     strainDesignMILP.py  (SDMILP)
   ─ ANY / BEST / POPULATE loops
   ─ verify_sd
   ─ exclusion constraints
        │
        ▼
   Post-processing      compute_strain_designs.py, strainDesignSolutions.py
   ─ expand_sd (decompress)
   ─ GPR translation
   ─ SDSolutions

All methods (MCS/suppress/protect, OptKnock, RobustKnock, OptCouple) use the same pipeline — only how addModule builds the LP blocks differs.

2. Repository Structure

straindesign/
├── __init__.py                     # Package exports & avail_solvers detection
├── names.py                        # All string constants (solver names, module types, etc.)
│
├── compute_strain_designs.py       # Main user-facing orchestration function
├── strainDesignModule.py           # SDModule: problem specification (a dict subclass)
├── strainDesignProblem.py          # SDProblem: translates model+modules → MILP matrices
├── strainDesignMILP.py             # SDMILP: solves the MILP, manages solution loop
├── strainDesignSolutions.py        # SDSolutions: result container, GPR translation
│
├── compression.py                  # Network compression (RREF, RationalMatrix, nullspace)
├── networktools.py                 # GPR extension, regulatory extension, compress wrappers, LP suppression
├── lptools.py                      # FVA, flux space plotting, solver selection
├── parse_constr.py                 # Constraint/expression string → matrix conversion
├── indicatorConstraints.py         # IndicatorConstraints data class
├── solver_interface.py             # MILP_LP factory: instantiates the right solver
├── cplex_interface.py              # CPLEX backend (Cplex_MILP_LP)
├── gurobi_interface.py             # Gurobi backend (Gurobi_MILP_LP)
├── scip_interface.py               # SCIP backend (SCIP_MILP_LP)
├── glpk_interface.py               # GLPK backend (GLPK_MILP_LP)
├── efmtool_cmp_interface.py        # EFMtool JAR interface (legacy compression backend)
├── pool.py                         # SDPool: Windows-compatible multiprocessing pool
└── efmtool.jar                     # Bundled EFMtool binary

Tests are in:

tests/
└── *.py   (pytest test files, one per feature area)

Configuration / packaging:

setup.py           # python_requires=">=3.7"; classifiers include 3.9–3.13
conda-recipe/meta.yaml
.github/workflows/CI-test.yml

3. Public API & Entry Points

Primary function

from straindesign import compute_strain_designs, SDModule

sols = compute_strain_designs(
    model,
    sd_modules   = [sd_module1, sd_module2],
    solver       = 'gurobi',      # 'cplex' | 'gurobi' | 'scip' | 'glpk'
    max_cost     = 3,
    max_solutions= 10,
    solution_approach = 'best',   # 'any' | 'best' | 'populate'
    ko_cost      = {'r1': 1, 'r2': 2},   # reaction knockouts
    ki_cost      = {'r3': 1},             # reaction additions
    gko_cost     = {'geneA': 1},          # gene knockouts
    gki_cost     = {'geneB': 1},          # gene additions
    time_limit   = 3600,
)

SDModule specification

# MCS-style: suppress undesired subspace, protect desired one
SDModule(model, 'suppress', constraints='EX_product_e >= 0.1')
SDModule(model, 'protect',  constraints='BIOMASS >= 0.05')

# With inner optimization (force to inner optimum first)
SDModule(model, 'suppress',
         constraints='EX_product_e >= 0.1',
         inner_objective='BIOMASS', inner_opt_sense='maximize')

# OptKnock
SDModule(model, 'optknock',
         inner_objective='BIOMASS',
         outer_objective='EX_product_e')

# RobustKnock
SDModule(model, 'robustknock',
         inner_objective='BIOMASS',
         outer_objective='EX_product_e')

# OptCouple
SDModule(model, 'optcouple',
         inner_objective='BIOMASS',
         prod_id='EX_product_e',
         min_gcp=0.1)

Result access

sols.status           # 'optimal', 'infeasible', 'time_limit', ...
sols.reaction_sd      # list of dicts: {reac_id: -1 (KO) | +1 (KI) | 0 (not added KI)}
sols.gene_sd          # list of dicts: gene-level interventions (if gene_kos=True)
sols.get_strain_designs()          # same as reaction_sd
sols.get_gene_reac_sd_assoc()      # (reac_sd, assoc, gene_sd) with GPR linkage
sols.get_strain_design_costs()     # list of floats
sols.get_reaction_sd_bnds()        # intervention expressed as flux bounds

Lower-level access

from straindesign import SDProblem, SDMILP, MILP_LP
from straindesign import fva, select_solver
from straindesign.compression import compress_cobra_model

4. Preprocessing Pipeline

compute_strain_designs (compute_strain_designs.py) orchestrates all preprocessing before constructing the MILP. The pipeline uses two-pass compression: the model is compressed once before GPR extension (so that FVA and GPR processing operate on a smaller model), then again after GPR extension (to compress the newly added gene pseudoreactions). This roughly halves model size at each downstream step and yields a ~4x end-to-end speedup on genome-scale models like iML1515.

4.1 Model Cleaning

Functions: remove_ext_mets, remove_dummy_bounds (in networktools.py, re-exported from compression.py).

Steps:

Remove dummy bounds — COBRApy uses ±1000 as a proxy for ±∞. Anything above a bound_thres (from cobra.Configuration) is replaced with ±inf so big-M bounding works correctly.
Deep-copy the model so the user’s model is never mutated.
Remove external metabolites — metabolites with no stoichiometric coupling to the rest of the network (boundary metabolites already handled by COBRApy exchange reactions, but orphaned ones are pruned).

4.2 Compression Pass 1 (before GPR extension)

File: compression.py Orchestrated by: networktools.compress_model

Before FVA or GPR processing, the model is compressed with propagate_gpr=True. This merges coupled and parallel reactions while propagating GPR rules through the compression map: when reactions are lumped, their GPR rules are AND-combined (using sympy And/Or constructors for flattening and deduplication). The result is a much smaller model that retains correct GPR annotations.

After compress #1:

compress_modules(sd_modules, cmp_mapReac_1) — rewrites constraint/objective reaction IDs.
compress_ki_ko_cost(ko_cost, ki_cost, cmp_mapReac_1) — maps reaction and regulatory costs (gene costs are not yet present and are handled separately after GPR extension).

4.3 FVA

Function: bound_blocked_or_irrevers_fva, fva (from lptools.py), speedy_fva (from speedy_fva.py)

Now running on the compressed model (roughly half the reactions of the original), FVA identifies:

Blocked reactions (FVA min = FVA max = 0): forced to zero, removed from intervention candidates.
Irreversible reactions (FVA min >= 0 or FVA max <= 0): bounds tightened.
Essential reactions (cannot be zero under module constraints): excluded from knockout candidates.

fva() in lptools.py delegates to speedy_fva(), which replaces the legacy implementation with a two-phase approach that reduces LP count by 70-80% on genome-scale models.

speedy_fva algorithm

Phase 1 — Scan LPs resolve the majority of bounds cheaply:

v=0 feasibility (Phase 1a): If lb[j] <= 0 <= ub[j], then zero flux is feasible, so min x_j = 0 (or max x_j = 0 when the bound is one-sided). No LP needed.
min(sum|x|) scan LP (Phase 1b): Builds an extended LP with splitting variables for reversible reactions to minimize total absolute flux. The optimal vertex pushes many reactions to zero or their bounds. A vectorized bound scan marks any reaction whose flux equals its variable bound (within tolerance).
Push-to-bounds iteration (Phase 1c): Using the same LP with dual simplex warm-start, iteratively sets directed objectives that push unresolved reactions toward their upper/lower bounds. Each round: set c[j] = -1 for unresolved-max reactions, solve, scan; then c[j] = +1 for unresolved-min reactions, solve, scan. Stops when a round resolves fewer than 5 new bounds. Warm-started dual simplex makes each re-solve near-instant (the basis stays primal feasible after objective changes).

Phase 2 — Individual LPs for remaining unresolved objectives:

Sequential mode (default for small/medium models): warm-started dual simplex with periodic LP rebuild (every 200 solves) to limit degeneration. Each LP vertex is checked for co-optimization: a bound scan on the full solution vector that may resolve other unresolved directions for free.
Parallel mode (for large models, configurable via threads): process-based SDPool with per-worker LP instances.

Optional compression: speedy_fva can compress the model internally (auto-enabled for n >= 200). This uses single-pass nullspace compression + conservation removal, with a fast model copy that avoids expensive deepcopy(solver).

Benchmarks (Gurobi, iML1515, compressed, single-threaded):

Configuration	Standard FVA	speedy_fva	LP count	Speedup
Unconstrained	90s	15s	71% fewer	6.0x
Biomass >= 0.1	17s	13s	41% fewer	1.3x

Phase 1 resolves 70-80% of objectives in the unconstrained case, but only ~40% when constraints limit the feasible space (fewer reactions can achieve their variable bounds).

Design decisions and rejected alternatives

KKT dual certification (investigated, not pursued): We explored using KKT optimality conditions to certify co-optimal reactions without solving additional LPs. After solving min x_j, we obtain a vertex solution and its basis. For another reaction k, optimality at the same vertex requires dual variables satisfying S^T lambda + mu_ub - mu_lb = e_k with correct complementarity. This can be checked by:

Extracting the basis B from the solver
Computing reduced costs for objective e_k: pi = B^{-T} e_{pos(k)}, then r_i = -A_i^T pi for nonbasic variables
Verifying sign conditions: r_i >= 0 for nonbasic-at-lb, r_i <= 0 for nonbasic-at-ub

The problem is that bound scanning already handles nonbasic variables (they are at their bounds by definition in a simplex vertex). KKT certification only adds value for basic variables (at interior points), which requires a sparse LU factorization of B per LP solve (1-5ms for m=900) — comparable to or exceeding the cost of the warm-started dual simplex LP itself (~0.5ms). The overhead exceeds the savings.

An earlier implementation (DualChecker) used null-space projection to test dual feasibility via sparse LU. On compressed models, this produced zero certifications because compression creates linearly dependent column subsets (coupled reactions with diagonal values of 2/3 instead of 1), causing the factorization check to reject all candidates. On uncompressed models certifications were possible but the null-space dimension was too large (d=39-304 on iMLcore) for efficient LP-based certification. The approach was removed in favor of the simpler bound scan, which captures the same information for nonbasic variables at zero cost.

Intelligent objective ordering (considered, marginal benefit): After each LP solve, the basis determines which variables are basic (interior) vs nonbasic (at bounds). Solving next for a basic variable should require fewer dual simplex pivots than jumping to a distant nonbasic variable. However, after bound scanning, the remaining unresolved reactions are precisely those that were basic (interior) in every solution seen — these “hard” reactions resist being pushed to bounds regardless of ordering. The objective change is always a rank-1 update, and pivot count depends more on problem structure than ordering. Warm-started dual simplex is already so fast per-LP (~0.5ms) that ordering heuristics cannot recover their own overhead.

4.4 GPR Integration

File: networktools.py Key functions: reduce_gpr, extend_model_gpr

When gene-level interventions are requested (gko_cost or gki_cost provided), the model is extended with artificial reactions that represent gene-level knockouts. This allows the MILP’s binary variables z to operate on genes rather than reactions. Because the model has already been compressed, there are fewer reactions to extend and fewer GPR rules to process.

Step 1 — Reduce GPR rules (`reduce_gpr`)

Uses AST-based GPR evaluation (evaluate_gpr_ast):

Each reaction’s GPR rule is parsed to an AST (ast.parse).
Genes are removed if their knockout/addition cannot affect any non-blocked, non-essential reaction.
This is determined by evaluating the DNF/CNF GPR expression with gene states {True, False, None (undetermined)}.
Reduces the gene intervention space before extension.

Step 2 — Extend model with artificial gene reactions (`extend_model_gpr`)

For each gene/isoenzyme group that participates in the GPR rules, artificial reactions and metabolites are added to encode the boolean logic:

AND gate: If a reaction requires g1 AND g2, an artificial metabolite and two artificial reactions are added. Knocking out either gene removes the pathway.
OR gate: If g1 OR g2, parallel paths are created. Both must be knocked out to abolish the reaction.

The function returns a reac_map dict mapping original reactions to the new artificial ones in the extended model. This map is used later to:

Translate gko_cost/gki_cost into reaction costs on the extended model.
Expand computed solutions back to gene-level interventions.

Important: Gene names starting with digits are prefixed with 'g' (cobra.manipulation.rename_genes) because Python AST cannot parse identifiers starting with digits.

Step 3 — GPR solution translation (`SDSolutions.init`)

After computation, the artificial reactions in solutions are translated back to gene-level interventions using gpr_eval. This evaluates each solution’s intervention set against the DNF GPR rules to determine which reactions are truly knocked out (or added).

gpr_eval(cj_terms, interv):

cj_terms: the GPR rule as a list of “conjunctive clauses” (AND-groups within OR expression).
interv: current gene intervention state.
Returns True/False/nan (active/inactive/undetermined).

4.5 Compression Pass 2 (after GPR extension)

After GPR extension adds gene pseudoreactions, the model is compressed a second time. This pass merges newly coupled gene pseudoreactions and parallel reactions introduced by the GPR encoding. The no_par_compress_reacs set is rebuilt from the (now rewritten) sd_modules to protect constraint/objective reactions.

The final compression map is the concatenation of both passes: cmp_mapReac = cmp_mapReac_1 + cmp_mapReac_2. When decompressing solutions, expand_sd processes this list in reverse order (pass 2 first, then pass 1).

A final FVA pass identifies essential reactions in the fully compressed model and removes them from the knockout candidates.

4.6 Network Compression — Algorithm Details

File: compression.py Orchestrated by: networktools.compress_model

Network compression reduces model size while preserving the strain design problem. It is lossless — all solutions in the compressed space correspond to solutions in the original, and vice versa (via the cmp_mapReac mapping).

RationalMatrix

The core data structure for numerically exact RREF computation. It uses sparse dual-integer storage (numerator matrix + denominator matrix, both scipy.sparse.csr_matrix). Operations are performed in rational arithmetic using fractions.Fraction and sympy.Rational. The model stays in exact rational arithmetic throughout both compression passes; no float conversion is performed between passes.

Key operations: row operations, scaling, finding pivots, extracting nullspace.

Floating-point coefficients are converted via float_to_rational with configurable precision (detect_max_precision scans the model’s stoichiometry).

Compression Algorithm (StoichMatrixCompressor)

The main class StoichMatrixCompressor performs the following in each cycle, repeating until no further reduction is possible:

Parallel reactions (compress_model_parallel): Reactions with identical stoichiometry (same column in S) are merged via hash-based O(nnz) comparison. This is cheap and runs first to shrink the model before the expensive RREF step. The resulting “super-reaction” carries the combined cost, and the z-variable controls all of them simultaneously. Hash matches are verified with exact index/data comparison to prevent collisions.
Conservation relation removal (remove_conservation_relations): Compute the left nullspace of the stoichiometric matrix S over the rationals (homogeneous dependencies among metabolite rows). Linearly dependent rows are removed. This reduces the row count of S for the subsequent RREF computation.
Exit check: If the loop has completed at least one full cycle (including a coupled step) and neither parallel nor coupled compression reduced the reaction count, exit. This ensures the model is always fully cleaned (conservation relations removed) on exit.
Coupled reactions (compress_model_coupled): Two reactions r_i and r_j are coupled if there exists a scalar alpha such that in every feasible steady state, v_i = alpha * v_j. These can be lumped into one reaction. Coupling is detected by computing the right nullspace (RREF) and identifying proportional columns. Lumpable reactions are merged and their costs combined. When propagate_gpr=True, GPR rules are AND-combined during merging. Contradicting groups (where bounds intersection yields lb > ub or lb = ub = 0) are detected and removed.

Blocked reactions are removed once before the loop via remove_blocked_reactions.

Compression Result

cmp_mapReac: a list of dicts, one per compression pass. Each dict maps compressed reaction IDs to their original reactions with scaling factors.

This is used by:

compress_modules(sd_modules, cmp_mapReac) — rewrites constraints/objectives in terms of compressed reaction IDs.
compress_ki_ko_cost(cmp_ko_cost, cmp_ki_cost, cmp_mapReac) — maps intervention costs.
expand_sd(sd, cmp_mapReac) — after solving, maps binary solutions back to original reactions (processes passes in reverse order).

Backend Selection

The compression backend is chosen by the compression_backend parameter:

'sparse_rref' (default): Pure Python, uses RationalMatrix for exact arithmetic. No Java dependency.
'efmtool_rref' (legacy): Uses the bundled efmtool.jar via JPype for RREF. Requires pip install straindesign[java]. Available via compress_model_efmtool.

4.7 Regulatory Interventions

Function: extend_model_regulatory (networktools.py)

Regulatory interventions represent changes that activate or silence a pathway without directly knocking out a gene. They are encoded similarly to gene interventions: artificial reactions are added to the model, and their “presence” or “absence” is controlled by binary z-variables.

Regulatory constraints are split by parsability: constraints that reference only existing reaction IDs (e.g., r6 >= 4.5) are applied before compression pass 1; constraints that reference gene IDs (e.g., g4 <= 0.4) are deferred until after extend_model_gpr creates the corresponding pseudoreactions.

After solving, postprocess_reg_sd in compute_strain_designs.py converts regulatory intervention values from numeric (1/0) to boolean (True/False) in the solution dict.

5. MILP Construction

File: strainDesignProblem.py Classes: SDProblem, ContMILP Key functions: build_primal_from_cbm, LP_dualize, farkas_dualize, reassign_lb_ub_from_ineq, prevent_boundary_knockouts, link_z

5.1 Overview & z-map Matrices

The MILP has the form:

minimize   c' * x
subject to:
  A_ineq * x  <=  b_ineq
  A_eq   * x   =  b_eq
  lb  <=  x  <=  ub
  x[0..num_z-1]  ∈ {0,1}         (binary z-variables: interventions)
  x[num_z..]     ∈ ℝ             (continuous: flux variables + dual variables)
  indicator constraints or big-M

The first num_z variables are binary: z_i = 1 means reaction i is knocked out (or added, if z_inverted[i]).

Three z-map matrices track which z-variables affect which parts of the LP:

z_map_vars[i, j]: z-variable i controls continuous variable j (knockout/addition of that reaction flux).
z_map_constr_ineq[i, j]: z-variable i controls inequality constraint j.
z_map_constr_eq[i, j]: z-variable i controls equality constraint j.

A value of +1 means “z=1 removes this constraint/variable” (knockout); -1 means “z=1 adds this” (knock-in).

These maps are propagated through all dualization steps (see below) and ultimately used in link_z to wire up big-M coefficients or indicator constraints.

MILP header rows (fixed positions):

Row 0 (idx_row_maxcost): −cost' z ≤ 0 (ensures non-negative cost)
Row 1 (idx_row_mincost): cost' z ≤ max_cost
Row 2 (idx_row_obj): objective constraint (used by fixObjective)

5.2 Primal LP from COBRApy Model

Function: build_primal_from_cbm(model, V_ineq, v_ineq, V_eq, v_eq, c=None)

Constructs the standard LP from a COBRApy model:

A_eq  * v  =  b_eq    (stoichiometric matrix S stacked with V_eq)
A_ineq* v  <= b_ineq  (V_ineq constraints, e.g. flux ratio constraints)
lb <= v <= ub         (reaction bounds from model)

A_eq = [S; V_eq], b_eq = [0; v_eq]
z_map_vars = I (identity — each z directly maps to its reaction variable)
z_map_constr_ineq = 0, z_map_constr_eq = 0 (in the primal, knockouts affect variables, not constraints)

prevent_boundary_knockouts: If a variable has lb < 0 or ub > 0, knocking it out to zero requires converting these non-zero bounds into inequality constraints. This prevents situations where a “zero flux” constraint on a reaction with lb = −1000 would be applied incorrectly.

5.3 SUPPRESS Module — Farkas’ Lemma Dual

Function: farkas_dualize(A_ineq_p, b_ineq_p, A_eq_p, b_eq_p, lb_p, ub_p, ...)

The SUPPRESS module wants to make a certain flux subspace infeasible. Farkas’ lemma provides a certificate of infeasibility:

Farkas’ Lemma: The system {A_ineq x ≤ b_ineq, A_eq x = b_eq, lb ≤ x ≤ ub} is infeasible if and only if there exist multipliers y such that:

A_ineq' y_ineq + A_eq' y_eq = 0   (dual feasibility)
b_ineq' y_ineq + b_eq' y_eq < 0   (certificate: b'y < 0)
y_ineq ≥ 0

The Farkas dual is constructed by calling LP_dualize with zero objective (c = 0), yielding dual variables y. Then the constraint c_dual' y ≤ −1 is appended (equivalently: b_ineq' y_ineq + b_eq' y_eq ≤ −1), which enforces the Farkas infeasibility certificate.

When an intervention (z) removes a primal constraint or variable, the corresponding dual variable is also removed — this is the mechanism by which knockouts “activate” the SUPPRESS condition. When the dual system becomes feasible with the ≤ −1 constraint, the corresponding primal flux subspace is infeasible, i.e., the targeted flux state is suppressed.

Known limitation: For the special case A x = b, x ∈ ℝ, b ≈ 0 (homogeneous equality with unrestricted variable), the Farkas certificate requires b'y ≠ 0 rather than < 0. This edge case is not currently implemented (see section 11).

5.4 PROTECT Module

The PROTECT module keeps a flux subspace feasible. It uses the primal LP directly (build_primal_from_cbm) without dualization, then calls reassign_lb_ub_from_ineq to tighten bounds. As long as the primal system is feasible, the protected state is reachable. The objective is set to zero (feasibility only).

5.5 OPTKNOCK — LP Duality for Bilevel Optimization

OptKnock solves:

maximize  outer_objective(v)
s.t.      v ∈ argmax { inner_objective(v) : v feasible }

This bilevel program is converted to a single-level MILP using LP duality:

LP Duality Theorem: For a minimization LP {min c'x : Ax ≤ b, x ≥ 0}, the dual is {max b'y : A'y ≤ c, y ≥ 0}. At optimality, primal and dual objectives are equal: c'x* = b'y*.

LP_dualize(A_ineq_p, b_ineq_p, A_eq_p, b_eq_p, lb_p, ub_p, c_p, ...)

Constructs the dual LP. The primal constraints become dual variables and vice versa:

Primal	→	Dual
`x ≥ 0`		`y ≤ c_i` (inequality constraint)
`x ∈ ℝ`		`y = c_i` (equality constraint)
`x ≤ 0`		`y ≥ c_i` (flipped inequality)
`A_ineq row`		`y_ineq ≥ 0`
`A_eq row`		`y_eq ∈ ℝ`
Objective `b'y`		`max → min`, `b` becomes new objective

The OptKnock MILP block is assembled as:

# Build: [primal (with constraints) | dual (without constraints)]
# Connect: c_inner_primal == c_inner_dual (strong duality)
A_eq_combined = block_diag(A_eq_v, A_eq_dual) stacked with
                hstack(c_primal, c_dual_objective) == 0

The outer objective is set as the MILP objective. The strong duality equality guarantees that v is at the inner optimum.

5.6 ROBUSTKNOCK — Three-Level Duality

RobustKnock solves a min-max problem:

maximize  min_v { outer(v) : v ∈ argmax inner(v) }

This requires two applications of LP duality:

First, the inner bilevel problem {inner(v) at max} is dualized into primal+dual pair (like OptKnock).
The combined primal+dual system is treated as a new primal, and dualized again to convert the outer minimization to a maximization.

The result is a three-layer LP construction assembled into a single MILP block. The outer objective (minimized production, maximized over worst case) becomes the MILP objective.

This is the most complex module and the most sensitive to numerical issues.

5.7 OPTCOUPLE — Growth Coupling Potential

OptCouple maximizes the Growth Coupling Potential (GCP):

GCP = max_v { inner(v) } − max_v { inner(v) | prod = 0 }

The two terms are the max growth with and without production. The objective is their difference.

Construction:

Build primal for the “with production” LP → (A_ineq_p, ..., c_p).
Build primal for the “baseline” LP (no production constraint) → (A_ineq_b, ..., c_b).
Dualize the baseline LP → c_dual_b (dual objective coefficients).
Assemble combined block with strong duality for both LPs.
Final objective: c_p - c_b (growth with production minus baseline).

min_gcp parameter sets a lower bound on GCP in the MILP.

5.8 Linking Binary Variables: `link_z`

Function: link_z in SDProblem (called at end of __init__).

This is the step that physically wires the binary variables z into the MILP. It proceeds in 5 steps:

Step 1 — Split knockable equalities into pairs of inequalities: An equality A_eq[j] x = b that is controlled by z is split into A_ineq x ≤ b and −A_ineq x ≤ −b. Both inequalities are individually linked to z, so when z is activated, both are activated together (enforcing zero).

Step 2 — Convert variable knockouts to inequality knockouts: A variable knockout (controlled by z_map_vars) is converted to: x ≤ 0 and −x ≤ 0. These bound constraints are added to A_ineq and linked to z via z_map_constr_ineq.

Step 3 — Compute big-M values via LP bounding (parallel): For each knockable constraint A[i] x ≤ b[i], compute:

M_i = max { A[i] x : x in most-relaxed feasible LP }

The “most-relaxed LP” removes all knockable constraints (allowing all knockouts) to get an upper bound on A[i] x. This gives the tightest valid big-M for each constraint.

Parallelized via SDPool when there are > 1000 constraints. Each worker (worker_init/worker_compute) holds a private LP object.

Step 4 — Link z via big-M for bounded constraints (finite M): For constraints where a finite M_i was found:

z=1 knocks out A[i] x ≤ b: add coefficient −M_i in column z_i of row i: A[i] x − M_i * z ≤ b → when z=1: A[i] x ≤ b + M_i (always satisfied if M_i tight enough).
z=0 knocks out: A[i] x + (M_i − b) * z ≤ M_i.

Step 5 — Create indicator constraints for remaining knockables: Constraints where M_i = ∞ (unbounded LP subspace) cannot use big-M. These become indicator constraints passed to the solver backend:

z_i = indicval  →  A[i] x ≤ b[i]   (sense 'L')
z_i = indicval  →  A[i] x  = b[i]  (sense 'E')

5.9 Big-M Bounding

The quality of big-M values is critical:

Too small → valid solutions may be cut off (incorrect results).
Too large → numerical instability, especially for GLPK.
M = ∞ → the constraint must use an indicator constraint (solver-dependent, slower).

StrainDesign computes per-constraint tight bounds via LP. The LP maximizes A[i] x subject to the static (non-knockable) constraints, giving the tightest valid M_i for each knockable constraint.

For GLPK (which has no indicator constraint support), M defaults to the COBRApy bound threshold (default: 1000). This is conservative and may cause numerical issues on models with large flux ranges. Setting M manually via the M parameter allows tuning.

Round-up: computed Ms are ceiled to 5 significant decimal digits.

5.10 Indicator Constraints vs. Big-M

Feature	Indicator Constraints	Big-M
Solver support	CPLEX, Gurobi, SCIP	All (required for GLPK)
Numerical stability	Excellent (no large coefficients)	Can be poor with large M
Solver speed	Typically faster (native branching)	Depends on M tightness
When used	Unbounded constraints (M = ∞)	Bounded constraints (finite M)
Forced by user	`M=some_value` forces big-M everywhere	Default for GLPK

The IndicatorConstraints class (indicatorConstraints.py) stores:

binv: indices of binary variables that trigger each constraint
A: coefficient matrix for the triggered constraints
b: RHS
sense: 'L' (≤), 'E' (=), 'G' (≥)
indicval: which value of z (0 or 1) activates the constraint

Each solver backend translates IndicatorConstraints into its native format. For GLPK, the solver interface translates them to big-M constraints using the provided M.

Sign convention in z-maps (subtle!):

z_map_constr_ineq[i, j] = +1: z_i = 1 → constraint j removed (knockout)
z_map_constr_ineq[i, j] = −1: z_i = 1 → constraint j active (knock-in)

This sign is inverted when KIs vs. KOs are assigned (in z_kos_kis diagonal matrix at end of SDProblem.__init__), so the final indicator constraint indicval correctly maps to 0 (KO: constraint active when z=0) or 1 (KI: constraint active when z=1).

6. Solver Backends

6.1 MILP_LP Unified Interface

File: solver_interface.py Class: MILP_LP

MILP_LP is a factory/dispatcher: its __init__ detects the requested solver and instantiates the appropriate backend class (Cplex_MILP_LP, Gurobi_MILP_LP, etc.). All backends expose the same interface:

class MILP_LP:
    solve()                       # → (x, opt, status)
    slim_solve()                  # → float (objective value only, fast)
    populate(n)                   # → (x_list, opt, status)  [CPLEX/Gurobi only]
    set_objective(c)
    set_ub(ub_list)               # [(idx, val), ...]
    set_time_limit(t)
    add_ineq_constraints(A, b)
    add_eq_constraints(A, b)
    set_ineq_constraint(row, c, b) # Overwrite a specific row
    clear_objective()
    set_objective_idx(idx_val_list)

slim_solve() is used heavily for verification — it runs the LP and returns just the objective value (or nan for infeasible). This avoids extracting the full solution vector.

6.2 CPLEX

File: cplex_interface.py Class: Cplex_MILP_LP

Native indicator constraint support (cplex.indicator_constraints.add)
Solution pool (populate_solution_pool) for enumerating multiple solutions
Thread control via cplex.parameters.threads
Seed via cplex.parameters.randomseed
Status mapped from CPLEX codes to StrainDesign names (OPTIMAL, INFEASIBLE, etc.)
Advantages: Best solution pool support, fastest for large MILPs in practice, population-based enumeration.
Limitation: Licensed (academic/commercial). No indicator constraints → automatic big-M fallback.

6.3 Gurobi

File: gurobi_interface.py Class: Gurobi_MILP_LP

Native indicator constraint support (gp.Model.addGenConstrIndicator)
Solution pool via Gurobi SolCount / PoolSolutions
Seed via Seed parameter
Status mapped from Gurobi Status codes
Advantages: Often fastest solver, excellent solution pool. MIPFocus parameter available for tuning.
Note: Requires a Gurobi license (free for academic use).

6.4 SCIP

File: scip_interface.py Class: SCIP_MILP_LP

Indicator constraints via pyscipopt addConsIndicator
No native solution pool → falls back to iterative compute_optimal
Seed via setSeed
Advantages: Open-source, no license needed, handles large models well.
Limitation: No populate support; POPULATE mode falls back to BEST.

6.5 GLPK

File: glpk_interface.py Class: GLPK_MILP_LP

No indicator constraint support — all indicator constraints are translated to big-M at construction time.
Uses swiglpk (bundled with COBRApy).
No solution pool → iterative enumeration only.
Advantages: Always available without additional installation.
Limitation: Poorest performance on large MILPs. Sensitive to M value choice. Default M = 1000 (COBRApy bound threshold).

7. Computation Approaches — SDMILP

File: strainDesignMILP.py Class: SDMILP (inherits from SDProblem and MILP_LP)

SDMILP.__init__ calls SDProblem.__init__ to build the MILP matrices, then MILP_LP.__init__ to set up the solver backend.

7.1 ANY

Method: SDMILP.compute(**kwargs)

Finds feasible solutions without enforcing global optimality:

Solve MILP as-is (no cost-minimization objective, just feasibility).
Verify each solution.
Add exclusion constraint.
Repeat until max_solutions or infeasibility.

For MCS problems, z_sum ≤ max_cost and z_sum ≥ 0 are the only cost constraints. Results may have suboptimal (higher) cardinality compared to BEST.

Warning: For OptKnock/RobustKnock, ANY may return the wild-type (z=0) immediately if that satisfies the outer objective. Use BEST for these methods.

7.2 BEST

Method: SDMILP.compute_optimal(**kwargs)

Finds globally optimal solutions:

Solve MILP with full objective (minimize cost for MCS; maximize outer for OptKnock, etc.).
Verify solution.
Add exclusion constraint (excludes solution and supersets for MCS; only the exact solution for OptKnock via add_exclusion_constraints_ineq).
Fix objective: fixObjective(c_bu, opt) — constrain next solution to have at least the same objective value.
Among solutions with equal objective, minimize intervention cost: setMinIntvCostObjective().
Repeat until max_solutions, infeasibility, or time limit.

For MCS computations, the MILP first checks if the wild-type already satisfies all modules (via verify_sd(sparse.csr_matrix((1, num_z)))).

7.3 POPULATE

Method: SDMILP.enumerate(**kwargs)

Uses the solver’s native solution pool (CPLEX/Gurobi):

Call populateZ(max_solutions) to generate a pool.
Verify all solutions in the pool.
Add exclusion constraints for invalid solutions.
Repeat.

Falls back to compute_optimal for SCIP and GLPK (which lack native pool support).

Gurobi pool notes: Pool quality/diversity can be tuned via Gurobi’s PoolSearchMode and PoolGap parameters, but this is not currently exposed in the interface.

7.4 Exclusion Constraints

Methods: add_exclusion_constraints(z), add_exclusion_constraints_ineq(z)

For a solution with binary vector z:

Case 1 — z is all zeros (wild-type): Add sum(z_i) ≥ 1 → forces at least one intervention next time. Implemented as 1 ≥ −1 (infeasible if all z = 0).

Case 2 — z has exactly one non-zero entry (single intervention): Set ub[idx] = 0 → permanently disable that intervention.

Case 3 — z has multiple non-zero entries: Integer cut: sum(z_i : i active) ≤ sum(z_i_current) − 1 This excludes the exact solution and all its supersets (for MCS: if removing r1 and r2 works, we don’t want r1+r2+r3 as a “new” solution since it’s a non-minimal superset).

add_exclusion_constraints_ineq uses a different cut that excludes only the exact solution, not its supersets (used for OptKnock where we want all optima).

7.5 Solution Verification

Method: SDMILP.verify_sd(sols)

After each MILP solution, the binary vector z is verified by building and solving a continuous LP with the corresponding constraint/variable removal applied:

From cont_MILP (saved before link_z): the original LP matrices with z-map information.
For each z_i in the solution: determine which variables (z_map_vars) and constraints (z_map_constr_ineq, z_map_constr_eq) are deactivated.
Build the restricted LP: remove deactivated constraints/variables.
Solve via slim_solve().
If slim_solve() returns nan → infeasible → solution is valid (for SUPPRESS) or invalid (for PROTECT).

This catches numerical artifacts where the MILP returns a slightly non-integer solution that passes rounding but does not actually satisfy the biological constraint.

cont_MILP is the continuous pre-link_z version of the problem. It preserves the structural relationship between z-variables and constraints/variables before big-M coefficients obscure it. This is why the continuous validation is numerically robust even if M values were imprecise.

8. Post-processing & Result Container

8.1 Solution Decompression

Function: expand_sd(sd, cmp_mapReac) in networktools.py

After solving on the compressed model, solutions refer to compressed reaction IDs. expand_sd maps each compressed intervention back to the set of original reactions it represents (using cmp_mapReac).

A single compressed KO may expand to multiple original reaction KOs (coupled/parallel reactions).
The expanded solution preserves the cost structure.

Function: filter_sd_maxcost(sd, max_cost, uncmp_ko_cost, uncmp_ki_cost) After expansion, the cost of a solution in original space may differ from compressed space (due to parallel reactions with non-unit costs). Solutions exceeding max_cost in original space are filtered out.

8.2 GPR Translation

SDSolutions translates gene interventions to reaction phenotypes using gpr_eval.

For each gene-level solution gene_sd:

Collect all affected reactions from the model’s GPR rules.
Evaluate each reaction’s GPR with the gene states using gpr_eval.
Reactions whose GPR evaluates to False (given knockouts) are collected as reaction knockouts.
Reactions requiring knocked-in genes that are True (given additions) are collected as reaction additions.

get_gene_reac_sd_assoc() returns three parallel lists:

reac_sd: reaction-level intervention dicts
assoc: association dict showing which genes caused which reaction change
gene_sd: the original gene interventions

8.3 SDSolutions

File: strainDesignSolutions.py

SDSolutions stores:

reaction_sd: list of dicts {reac_id: -1 | 0 | +1} for each solution
gene_sd: list of gene-level dicts (if gene interventions were used)
status: solver status string
sd_setup: full setup dict (model_id, modules, costs, solver, etc.) — can be serialized to JSON and reloaded

Serialization: sd_setup is JSON-compatible, allowing CNApy to save .sd files that can be reloaded into compute_strain_designs via the sd_setup dict parameter.

9. Constraint & Expression Parsing

File: parse_constr.py

Handles user-supplied constraints and objectives in multiple formats:

Input formats (all equivalent):

# String
"r1 + 3*r2 = 0.3, -5*r3 - r4 <= -0.5"

# List of strings
["r1 + 3*r2 = 0.3", "-5*r3 - r4 <= -0.5"]

# List of structured constraints
[[{'r1': 1.0, 'r2': 3.0}, '=', 0.3], [{'r3': -5.0, 'r4': -1.0}, '<=', -0.5]]

Key functions:

Function	Input → Output
`parse_constraints(constr, reaction_ids)`	Any format → `[[dict, sense, rhs], ...]`
`parse_linexpr(expr, reaction_ids)`	Any format → `[dict, ...]`
`lineq2mat(equations, reaction_ids)`	String list → `(A_ineq, b_ineq, A_eq, b_eq)` sparse matrices
`lineqlist2mat(eqlist, reaction_ids)`	Structured list → sparse matrices
`linexpr2dict(expr, reaction_ids)`	String → `{reac_id: coeff}`
`linexprdict2mat(expr_dict, reaction_ids)`	Dict → sparse row vector

The parse_constraints function uses regex splitting on ',' and '\n' to separate multiple constraints in a single string, then dispatches to lineq2list which uses further regex for coefficient/identifier/operator extraction.

10. Supporting Utilities

`lptools.py`

fva(model, variables, constraints, solver): Flux variability analysis. Parallelized via SDPool using fva_worker_init/fva_worker_compute. Returns a pandas.DataFrame with min/max flux values.
select_solver(solver, model): Selects solver based on availability, preference, and model settings. Priority: CPLEX > Gurobi > SCIP > GLPK.
plot_flux_space(model, ...): 2D/3D convex hull visualization using scipy.spatial.ConvexHull and Matplotlib.
DisableLogger: Context manager to suppress logging during FVA.

`networktools.py` — LP Suppression

Problem

COBRApy’s Model is tightly coupled to its solver backend (optlang). Every mutation — setting a reaction’s bounds, renaming a reaction, adding/removing metabolites — immediately updates the solver’s LP representation via set_bounds, set_linear_coefficients, _populate_solver, etc. During compression and preprocessing, hundreds of such mutations occur on a model copy. Synchronizing each one to the solver is pure overhead: the LP is rebuilt from scratch at the end anyway.

Early approaches patched optlang classes (Constraint.set_linear_coefficients, Variable.set_bounds, Objective.set_linear_coefficients) at the class level during suppression. This worked for optlang-level calls but missed cobra-level code paths: Reaction.add_metabolites accesses self.forward_variable, self.reverse_variable, and model.constraints[met.id] to build solver constraint updates. These paths reached the solver through cobra internals, not optlang’s public API.

Architecture: class-level cobra patching

suppress_lp_context(model) is a context manager that patches both optlang and cobra classes at the class level for the duration of the block. All solver-bound operations become no-ops, and cobra methods that would access the solver are redirected to safe stubs.

Patches applied by _suppress_lp_updates:

Target	Original	Suppressed replacement
`Reaction._set_id_with_model`	Renames solver variable	`_suppressed_set_id`: writes `_id` + updates `DictList._dict` index only
`Reaction.lower_bound` (setter)	Calls `update_variable_bounds` → solver	`_suppressed_set_lb`: writes `_lower_bound` only
`Reaction.upper_bound` (setter)	Calls `update_variable_bounds` → solver	`_suppressed_set_ub`: writes `_upper_bound` only
`Reaction.update_variable_bounds`	Accesses `forward_variable.set_bounds(...)`	No-op
`Model._populate_solver`	Full solver rebuild	No-op
`Model.remove_reactions`	Removes solver variables/constraints	`_suppressed_remove_reactions`: removes from model lists only
`Model.remove_metabolites`	Removes solver constraints	`_suppressed_remove_metabolites`: removes from model lists only
`Container.__getitem__`	Raises `KeyError` for missing keys	`_permissive_container_getitem`: returns `_SolverStub` for missing keys
optlang `Constraint.set_linear_coefficients`	Updates solver constraint	No-op
optlang `Variable.set_bounds`	Updates solver variable bounds	No-op
optlang `Objective.set_linear_coefficients`	Updates solver objective	No-op

_SolverStub is a lightweight null object returned by the permissive Container.__getitem__ when cobra code tries to access a solver variable/constraint that doesn’t exist (because _populate_solver is suppressed). It is hashable (used as a dict key in set_linear_coefficients calls) and has no-op set_bounds/set_linear_coefficients methods. This single class eliminates the need to patch every individual cobra method that might touch the solver.

_suppressed_set_id and DictList invariant: cobra’s DictList._dict maps id → int index (not id → object). The suppressed setter must maintain this invariant: it pops the old key and inserts the new key with the same integer index. Violating this causes ValueError in DictList.__setitem__ due to its duplicate-check logic.

Objective preservation through compression

When suppress_lp_context enters, it captures the model’s objective coefficients into model._suppressed_obj (a dict mapping reaction_id → float coefficient). This dict is updated by compression code as reactions are merged:

Coupled compression (compression.py, _apply_compression_to_model): when reactions are merged, their objective contributions are weighted by compression factors and summed into the surviving reaction’s entry.
Parallel compression (compression.py, compress_model_parallel): when identical-stoichiometry reactions are merged, their objective contributions are summed directly.
Legacy Java compression (efmtool_cmp_interface.py, compress_model_java): same logic using jBigFraction2sympyRat conversion factors.

On exit, if the model was modified (reaction IDs changed), the context manager rebuilds the solver from the model’s current state and restores the objective from _suppressed_obj. If the model was not modified (e.g., compute_strain_designs only modifies a copy), no rebuild occurs — this is critical because an unconditional rebuild would produce a solver with the wrong objective direction for methods like RobustKnock that have inner/outer objectives.

# Simplified exit logic
current_ids = {r.id for r in model.reactions}
if current_ids != pre_ids:
    model._solver = solver_interface.Model()
    model._populate_solver(model.reactions, model.metabolites)
    for rxn in model.reactions:
        c = final_obj.get(rxn.id, 0.0)
        if c != 0:
            rxn.objective_coefficient = c

Design principles

All plumbing in the context manager. Compression code, GPR extension, and other preprocessing functions are backend-agnostic — they manipulate the cobra model without knowing whether LP updates are suppressed. The context manager handles all solver synchronization.
Class-level patching, not instance-level. Patching at the class level (e.g., Reaction.lower_bound = property(fget, suppressed_fset)) ensures that ALL reactions are covered, including ones created during compression. Instance-level patches would miss newly created objects.
Conditional rebuild. The solver is only rebuilt when the model was actually modified. This prevents regressions in methods like RobustKnock/OptKnock/OptCouple that have complex objective structures.

Performance impact

On iML1515, LP suppression reduces preprocessing time from ~131s to ~70s (46% reduction). The bulk of the savings comes from avoiding ~2000 set_linear_coefficients calls during compression, each of which would otherwise rebuild solver constraint data structures.

`pool.py` — SDPool

Windows-compatible multiprocessing.Pool subclass. On Windows, the initializer function is pickled to a temporary file and loaded by worker processes, avoiding a performance regression from passing the initializer directly (COBRApy issue #997). Uses 'spawn' context universally to avoid fork-related instability.

`names.py`

All string constants. Importing with from straindesign.names import * is used throughout the codebase. Constants include module type strings, solver names, status codes, and parameter keys. This single-file approach ensures consistent naming.

Notable: PROTECT = 'mcs_lin' and SUPPRESS = 'mcs_bilvl' are defined first (legacy names), then immediately overwritten with 'protect' and 'suppress'. This is vestigial code — the legacy names are not used and could be removed.

`init.py`

Detects available solvers at import time:

avail_solvers = set()
try: import cplex; avail_solvers.add('cplex')
try: import gurobipy; avail_solvers.add('gurobi')
try: import pyscipopt; avail_solvers.add('scip')
try: import swiglpk; avail_solvers.add('glpk')

Exports: all public classes, functions, and the avail_solvers set.

Also calls _start_jvm() from efmtool_cmp_interface.py — see next section.

`efmtool_cmp_interface.py` — JPype/JVM Initialization

JPype and Java are optional dependencies. StrainDesign works without them; the default compression backend is sparse_rref (pure Python). The efmtool Java backend (compression_backend='efmtool_rref') requires jpype1 and a JVM.

Conditional eager JVM startup

__init__.py calls _start_jvm() at package import time. This function:

Checks if jpype1 is installed (via find_spec). If not → immediate return, no-op.
Adds efmtool.jar to the classpath.
Starts the JVM with jpype.startJVM("--enable-native-access=ALL-UNNAMED").
Loads all Java classes via import jpype.imports + Python import statements.
Registers atexit handler for clean JVM shutdown.

If any step fails (no Java installed, jar missing, etc.), the function returns silently. Users who never use efmtool_rref are unaffected. When _init_java() is later called (on first efmtool use), it checks _JAVA_INITIALIZED and falls back to JClass() loading.

Why eager, not lazy?

This design is critical for stability. In v1.14, from .efmtool import * in __init__.py started the JVM at import time. When we switched to lazy initialization (JVM started on first efmtool use), JPype began crashing with SIGBUS/SIGSEGV on Linux and macOS CI runners when processing larger matrices (iMLcore, 586 reactions).

Root causes identified:

JVM + OpenBLAS pthread conflict (jpype#808): The JVM modifies pthread stack allocation. If started after NumPy/OpenBLAS has already spawned worker threads, subsequent JNI calls can trigger SIGSEGV.
GC finalization race (jpype#934): Python’s garbage collector can attempt to finalize JPype proxy objects during a JNI call, causing Bus error. Mitigated with gc.disable()/gc.enable() around heavy Java calls in basic_columns_rat_java and compress_model_java.
import jpype.imports vs JClass(): The Python import-style class loading (import ch.javasoft...) sets up JNI references differently from JClass(). The import style (used in v1.14) is more stable on CI runners.

CI considerations

Java 21 (Temurin) is used on all CI platforms. The --enable-native-access=ALL-UNNAMED flag suppresses Java 17+ warnings about JPype’s System.load() calls. Despite all mitigations, JPype JNI crashes remain non-deterministic at ~1-in-20 frequency on some Linux runners (jpype#934). The iMLcore efmtool parity tests are marked --large (skipped on CI) since the e_coli_core tests exercise the same code path on a smaller, more reliable matrix size.

11. Known Issues, Urgent Actions & Future Work

11.1 Urgent Action Items

⚠ Farkas special case (homogeneous equality): In farkas_dualize, there is a documented comment noting the unimplemented edge case:

“In the case of (1) A x = b, (2) x ∈ ℝ, (3) b ≈ 0, Farkas’ lemma is special, because b'y ≠ 0 is required…”

This could cause silent incorrect results for SUPPRESS modules with equality constraints where b ≈ 0 (e.g., free reactions with zero RHS). Impact: depends on whether the compressed, preprocessed models trigger this case. A defensive check + explicit splitting into two inequalities would fix it.

⚠ Big-M for GLPK with large flux ranges: GLPK uses M = cobra_bound_threshold = 1000 by default. For models with non-standard bounds (e.g., ATPM = 8.39, but oxygen intake up to 10,000), this can produce incorrect strain designs. Models should always have remove_dummy_bounds applied, but verifying this is enforced is important.

⚠ add_exclusion_constraints_ineq has a bug:

A_ineq = [1.0 if z[j, i] else -1.0 for i in self.idx_z]
A_ineq.resize((1, self.A_ineq.shape[1]))  # ← list.resize doesn't exist

list.resize is a NumPy array method, not a list method. This would raise an AttributeError at runtime. This method is not used in the current main computation paths (OptKnock uses add_exclusion_constraints), but should be fixed. The fix: convert to np.array before calling .resize.

⚠ Legacy name double-assignment in names.py:

PROTECT = 'mcs_lin'
SUPPRESS = 'mcs_bilvl'
PROTECT = 'protect'      # overwrites silently
SUPPRESS = 'suppress'    # overwrites silently

The first assignments are dead code and should be removed to avoid confusion.

11.2 Low-Hanging Fruits

Solution space plotting issues (reported): lptools.plot_flux_space uses scipy.spatial.ConvexHull and Matplotlib. Known issues include:

3D plots with degenerate hulls (lower-dimensional polytopes) cause QhullError.
Color scaling in plot_flux_space_colored may fail for single-point solutions.
Fix: add degeneracy checks and fallback to scatter plots for ≤ 2 unique points.

Populate mode for SCIP/GLPK: Currently silently falls back to BEST. Should log a clear warning and document that SCIP/GLPK do not support POPULATE.

Gurobi pool tuning: The populate method in gurobi_interface.py could expose PoolSearchMode=2 (focused enumeration) for better solution diversity. Currently uses default settings.

RobustKnock + POPULATE: RobustKnock with POPULATE is not well-tested. The three-level dual structure is particularly sensitive to solver pool heuristics. Add explicit test cases.

Parallel FVA compression for large models: For genome-scale models (> 2000 reactions), compression can be slow. The RREF step in RationalMatrix has O(n²) row operations. Potential improvement: use block-sparse RREF or exploit model structure (e.g., subsystems).

FVA Phase 2 parallelization for non-Gurobi solvers: Gurobi releases the GIL during model.optimize(), enabling ThreadPoolExecutor-based parallelism with batched warm-start (4.3x speedup measured on iML1515). CPLEX, GLPK, and SCIP hold the GIL, so they must use process-based SDPool with its ~4-6s spawn overhead. Investigating ctypes/Cython wrappers that release the GIL for these solvers could enable thread-based parallelism across all backends.

Essential KI detection: essential_kis is computed but may not always be passed correctly through the compression/GPR pipeline. Verify that the set persists through compress_ki_ko_cost.

Opposite-direction parallel merging: Two irreversible reactions with opposite stoichiometry and complementary bounds (e.g., A→B (0,inf) and B→A (0,inf)) could be merged into one reversible reaction (-inf,inf). This would help eliminate flux cycles. However, a computed strain design might need to knock out one direction specifically (the other being essential), which limits the merger to non-knockable reactions. Few reactions are essential in both directions, so the practical benefit is small.

11.3 Future Development Directions

Nullspace-based primal/dual formalization: The referenced paper (DOI: 10.1093/bioinformatics/btz393, von Kamp & Klamt 2019) presents a nullspace-based approach to computing minimal cut sets. Integrating this formalization could:

Avoid explicit dualization for certain model types.
Give a more compact MILP representation.
Improve numerical conditioning.

A key idea: the primal flux space has a right nullspace (internal cycles). The Farkas dual operates in the left nullspace. Exploiting this duality directly, rather than constructing full dual matrices, could reduce MILP size significantly.

Variable splitting: Reversible reactions split into forward/reverse can sometimes improve MILP performance (tighter LP relaxation). This is a pre-processing step that should be optional and tested for impact on different model types.

On the mathematical basis of big-M values (investigated, not pursued):

This section documents a thorough investigation into computing theoretically optimal big-M values for the Farkas dual MILP. The investigation confirmed the mathematical foundations but revealed that the problem is inherently as hard as MCS enumeration itself. No implementation was made; this section exists to prevent re-derivation.

Where M comes from — the primal view:

In the Farkas dual formulation (section 5.3), each knockable constraint i is linked to a binary z_i via either big-M (A[i]x - M_i*z_i ≤ b[i]) or an indicator constraint (z_i = 0 → A[i]x ≤ b[i]). The big-M value M_i must be large enough that when z_i = 1, the constraint A[i]x ≤ b[i] + M_i is never binding — i.e., M_i ≥ max{A[i]x - b[i]} over all feasible x in every MCS that knocks out constraint i.

For SUPPRESS modules, the Farkas dual constraints encode the primal flux space through duality. Each knockable dual constraint corresponds to a primal reaction. The M value for reaction j’s dual constraint is determined by the primal flux behavior:

M_j = 1 / min|v_j| taken over all Elementary Flux Vectors (EFVs) in which reaction j participates.

Here, an EFV is a support-minimal feasible flux distribution satisfying the target constraint. MCS are the minimal hitting sets of EFV supports. For reaction j to appear in any MCS of cardinality > 1, there must exist an EFV where j carries flux — and the smallest such flux across all EFVs determines the tightest valid M.

Why single-KO M is trivial:

For reactions that are individually essential (size-1 MCS candidates), M_j = 1/FVA_min(|v_j|). This is a single LP: minimize |v_j| subject to Sv = 0, bounds, and the target constraint. This is already computed during preprocessing FVA and is the basis for the current M-bounding in step 3 of link_z.

Why multi-KO M is hard:

For reactions that only appear in MCS of size ≥ 2, M_j depends on the combinatorial structure of which other reactions are simultaneously knocked out. The minimum |v_j| is not taken over the full flux space (where j can always carry large flux via alternative pathways) but over the restricted space where a specific set of other reactions are forced to zero — and j must be essential in that restricted space.

Formally, computing the optimal M_j requires solving:

M_j = max over all MCS K containing j of:
        1 / min{ |v_j| : Sv = 0, lb ≤ v ≤ ub, v_target ≥ t,
                          v_i = 0 for i in K\{j} }

This is a max-min over an exponential number of MCS. Even computing the inner minimization for a single known MCS K is just an LP, but enumerating which K to consider requires MCS enumeration — the very problem M values are meant to help solve.

Combined primal/Farkas MILP — formulation and results:

We attempted to bypass MCS enumeration by embedding the combinatorial search in a single MILP. The formulation combines primal feasibility with a Farkas certificate:

Primal: S*v = 0, lb_i*(1-k_i) ≤ v_i ≤ ub_i*(1-k_i) (big-M knockout), v_target ≥ target_lb. Binary k_i = 1 knocks out reaction i.
Essential participation: Sign binary forces v_j ≥ δ or v_j ≤ -δ.
Farkas certificate: Same k_i binaries control dual constraints via indicators. When k_i = 0 (active): S[:,i]'u - r_i + s_i = δ_{i=target}. When k_i = 1 (knocked): r_i = s_i = 0. Farkas sum: Σ(ub_i*s_i - lb_i*r_i) ≤ target_lb - ε.
Objective: minimize |v_j| = t, giving M_j = 1/t.

The Farkas certificate proves that removing j from the active set (while keeping the knocked reactions knocked) makes the system infeasible — i.e., j is essential in the knockout context defined by the k_i values. The optimizer simultaneously searches for the knockout set and the flux distribution where j’s participation is minimized.

Practical obstacle — primal leakage:

The big-M knockout constraints v_i ≤ ub_i*(1-k_i) allow leakage when k_i is within solver tolerance of 1 but not exactly 1. With Gurobi’s default IntFeasTol = 10⁻⁵, a variable with k_i = 0.9999963 (accepted as integer 1) gives leakage v_i ≤ ub_i × 3.7×10⁻⁶. With default bounds ±1000, this is 0.0037 per reaction. Across 15+ knocked reactions, the cumulative leakage invalidates the primal solution entirely — the model appears feasible when it should be infeasible.

This was diagnosed by inspecting k/v values at full precision (tests/check_k_precision.py): PIt2r showed k = 0.9999963, v = 0.003679, meaning the “knocked” reaction carried measurable flux that kept the system alive.

Mitigation — FVA-tight bounds + solver tolerances:

Two measures reduce leakage to negligible levels:

FVA-derived bounds (replace ±1000 with actual flux ranges, typically O(1)–O(10))
Tight solver tolerances (IntFeasTol = 10⁻⁹, FeasibilityTol = 10⁻⁹, NumericFocus = 3)

With both, per-reaction leakage drops to O(bound × 10⁻⁹) ≈ 10⁻⁸.

Validation on E. coli core (95 reactions, 353 MCS up to size 3):

Ground truth M values were computed independently: enumerate all MCS via indicator-based MILP (no big-M), then for each MCS, solve one LP per knocked reaction to find the minimum |v_j| in the Farkas certificate. M_gt[j] = max over all MCS containing j. (See tests/ground_truth_M.py.)

With tight tolerances and FVA bounds, the combined MILP matched ground truth exactly for 6 of 15 tested reactions, spanning 4 orders of magnitude:

Reaction	M_gt	M_milp	Ratio	#KO	Time
TALA	413,793	413,793	1.0000	15	~120s
TKT1	413,793	413,818	1.0001	16	~120s
NADH16	68,966	68,966	1.0000	6	~35s
EX_etoh_e	68,966	68,966	1.0000	23	~120s
ETOHt2r	68,966	68,966	1.0000	9	~120s
ATPM	0.12	0.12	1.0000	36	0.8s

Reactions with intermediate M values (PFL: 329, NADTRHD: 0.91) did not converge within 120s, returning inflated values. The combinatorial search space is too large for the solver to explore efficiently.

Why we did not pursue this further:

The investigation confirmed that the formulation is mathematically correct — it produces exact M values when the solver converges. However, it also confirmed the fundamental insight: computing the optimal M_j for multi-reaction MCS is a combinatorial problem of the same complexity as MCS enumeration itself. There is no shortcut:

For size-1 MCS: M = 1/FVA_min — already computed cheaply via LP.
For size-k MCS (k ≥ 2): the optimal M depends on which other reactions are simultaneously knocked out, requiring exploration of the space of knockout combinations. This is exactly what the MCS MILP does.
The “FVA over the full Farkas polyhedron” approach (computing bounds with all constraints active) gives a lower bound on M — it tells you the range of dual variables when no constraints are relaxed. But in an actual MCS, some constraints ARE relaxed (z = 1), and the remaining dual variables may need larger values to construct a valid Farkas certificate. The FVA approach is a natural progression from the single-KO case but does not capture the combinatorial interaction that determines the true M for multi-reaction MCS.

In summary: the relationship M_j = 1/min|v_j| over EFVs containing j is exact but not efficiently computable. Any finite M chosen without enumerating MCS risks being too small (cutting off valid MCS) or too large (degrading solver performance). Indicator constraints avoid the M problem entirely at the cost of solver performance (lazy enforcement, weaker LP relaxation). This is the fundamental trade-off; there is no free lunch.

Investigation scripts were removed after documenting findings. The formulations above can be reconstructed from the mathematical descriptions if needed.

Reaction scaling for better-conditioned M values: Trace metabolites (e.g., MoYB at 7×10⁻⁷ mmol/gDW/h) create small inhomogeneous bounds in the primal, which after dualization produce small coefficients in the Farkas objective, forcing dual variables to be very large to satisfy b^T u ≤ -1. This leads to M values spanning many orders of magnitude (irreducible condition number ∝ κ(S)), harming LP relaxation quality.

Column scaling (per reaction) addresses this: multiply flux variable v_j by a power-of-10 factor α_j so that inhomogeneous bounds move into a reasonable range (e.g., [0.1, 10]):

Stoichiometry: S[:, j] /= α_j (coefficients shrink)
Bounds: lb_j *= α_j, ub_j *= α_j (bounds grow)
Mass balance: (S[:,j]/α_j) × (α_j × v_j) = S[:,j] × v_j (unchanged)

Row scaling (per metabolite) rebalances after column scaling: multiply metabolite row i by β_i (power of 10) to bring the max coefficient back to O(1). This does not change the solution space (just scales the mass balance equation).

Implementation would go between build_primal_from_cbm and dualization in addModule. The scaling factors must be tracked for unscaling continuous solution values in non-MCS problem types (OptKnock, etc.). For MCS, only binary z-variables are returned, so unscaling is not needed.

The FVA bounds from bound_blocked_or_irrevers_fva (already computed in preprocessing) provide the tight flux ranges needed to determine α_j. A per-module FVA considering module constraints could provide even tighter ranges.

Size-1 MCS via FVA on SUPPRESS modules: The current preprocessing (line 346–356 in compute_strain_designs.py) runs FVA on PROTECT (and OPTKNOCK etc.) modules to find essential reactions — reactions that must carry flux in the protected region — and marks them non-knockable. SUPPRESS modules are explicitly skipped (if m[MODULE_TYPE] != SUPPRESS).

The idea: also run FVA on SUPPRESS modules (without inner objective) to identify reactions essential for the target (undesired) region. Combining the two FVAs yields size-1 MCS cheaply and further reduces the knockable reaction set:

Step 1 — FVA on PROTECT (already implemented): Reactions where min(|flux|) > 0 under PROTECT constraints are essential for the protected region. These are already marked non-knockable.

Step 2 — FVA on SUPPRESS (new): Reactions where min(|flux|) > 0 under SUPPRESS constraints are essential for the target region. Knocking out any single one of these makes the target infeasible — each is a size-1 MCS candidate.

Step 3 — Cross-check: Size-1 candidates from step 2 that are NOT in the essential set from step 1 are true size-1 MCS: they disrupt the target without being required for the protected region. Add them to the output.

Step 4 — Mark non-knockable: All reactions identified in steps 1 and 2 can be removed from the knockable set for the main MILP:

Step 1 essentials: knocking them out violates protection (already handled).
Step 2 essentials not in step 1: they are size-1 MCS, already found. Since any larger MCS containing a size-1 MCS is not minimal, these reactions cannot participate in any new MCS. Safe to exclude.
Step 2 essentials also in step 1: essential for both regions, obviously non-knockable.

Benefits:

Size-1 MCS are found for free (just FVA, no MILP). This is the cheapest possible way to discover them.
Fewer knockable reactions in the main MILP → fewer z-variables → smaller B&B tree.
Reactions with tiny flux bounds (e.g., trace metabolite carriers at 1e-7) that are essential in the SUPPRESS region get excluded early, improving the conditioning of the main MILP’s M-value computation.
The model may simplify further (fewer knockables may enable additional compression passes or remove degenerate constraints), though the extent of this benefit is unclear.

Caveats:

Only finds size-1 MCS. Size ≥ 2 MCS still require the full MILP.
The inner objective (if present in the SUPPRESS module) is ignored for this FVA. This is correct for standard MCS but may miss reactions that are only essential under the bilevel structure.
Adds O(2n) LP solves per SUPPRESS module (one FVA). For compressed models this is fast.

MILP warm-starting: When BEST generates multiple solutions iteratively, each call starts cold. Warm-starting from the previous solution’s LP relaxation (available in CPLEX/Gurobi) could speed up subsequent solves significantly.

Improving RobustKnock numerics: The three-level dual construction in ROBUSTKNOCK creates very large matrices with potential for numerical blow-up. Applying scaling (row/column normalization) to the LP before and after dualization could help. Alternatively, use iterative refinement at the inner LP level.

SOS1/SOS2 constraints for bilevel formulations: Gurobi, CPLEX, and SCIP support SOS1 (Special Ordered Sets type 1) and SOS2 constraints natively. These could be useful for connecting variables of primal and KKT parts of multi-level optimization problems, encoding complementarity conditions (primal-dual coupling) more efficiently, or as an alternative to big-M for disjunctive constraints in bilevel formulations. SOS1 constraints enforce that at most one variable in a set is nonzero, which is exactly the complementarity condition in KKT-based reformulations.

Sparse matrix format conversions: The computation pipeline contains many sparse matrix format conversions (tocsr, tocsc, todok, tocoo) that may not all be necessary. A profiling pass could identify redundant conversions. Key areas: SDProblem.__init__ performs multiple tocsc/tocsr conversions during link_z; addModule builds matrices in one format then converts for block assembly; compression.py converts between formats during RREF and coupled-reaction detection. Eliminating unnecessary conversions could reduce preprocessing time, especially for large models.

Integration of regulatory networks beyond simple T/F interventions: The current regulatory intervention model (active/inactive) is binary. Future work could integrate thermodynamic constraints or kinetic feasibility checks.

SCIP native solution enumeration: SCIP supports solution enumeration beyond the incidental pool (limits/maxsol, default 100). Since version 2.0, the cons_countsols constraint handler can count and enumerate all feasible solutions of a constraint integer program via the count command. Key details:

constraints/countsols/collect = TRUE stores detected solutions (default FALSE = count only).
To enumerate all optimal solutions: solve to optimality to get c*, add the objective as a constraint with both bounds equal to c*, then run count.
Restarts must be turned off during counting (use SCIPsetParamsCountsols()).
SCIP uses unrestricted subtree detection, which can detect several solutions at once, so a soft solution limit may be exceeded before SCIP stops.
Collected solutions are stored with respect to active (non-presolved) variables only; they are lifted back into the original variable space when written to file.
This is not equivalent to Gurobi’s ranked near-optimal pool (PoolSearchMode). SCIP enumerates feasibility/optimality, not a ranked set of near-optimal solutions.

Currently, StrainDesign’s SCIP populate method uses a custom workaround (iterative solve with solution exclusion constraints). Replacing this with native cons_countsols enumeration could improve performance and correctness. However, pyscipopt does not yet expose the count command or cons_countsols parameters — this would require upstream pyscipopt changes or direct C API calls.

pyscipopt gaps for LP method and basis control: Several SoPlex and SCIP features that would benefit StrainDesign are not exposed through the pyscipopt Python bindings:

LP basis get/set: SCIP internally maintains LP bases via SoPlex, but pyscipopt does not expose SCIPlpGetBasisInd, SCIPgetLPBasisInd, or the SoPlex getBasis/setBasis methods. This prevents explicit basis warm-starting for LPs solved through the SCIP/SoPlex interface.
LP algorithm selection for SCIP_LP (pure LP via pyscipopt.LP): The pyscipopt.LP class wraps SoPlex directly but does not expose SoPlex’s setIntParam(ALGORITHM, ...) for selecting primal vs dual simplex. The SCIP_MILP path supports this via lp/initalgorithm and lp/resolvealgorithm parameters, but the direct SoPlex LP wrapper does not.
Solution enumeration: The count command and cons_countsols constraint handler (see above) are not accessible from pyscipopt.

These gaps are candidates for upstream pyscipopt issues or contributions.

Improved CNApy integration: The sd_setup dict format is designed for CNApy interoperability. Ensuring full round-trip serialization (JSON → compute_strain_designs → JSON) with all parameter types (modules, costs, constraints) is an ongoing compatibility concern.

12. MILP Performance Notes

The following design choices have a measurable impact on MILP enumeration performance. They are listed roughly by impact, based on benchmarking with e_coli_core (95 rxns, 353 MCS) and iML1515 (1920 compressed rxns, 393 MCS).

Unbounded reactions save constraints and variables. Reactions with infinite bounds (-inf / inf) do not generate big-M or indicator constraints because no finite bound needs enforcement. Calling remove_dummy_bounds early (replacing ±1000 with ±inf) is therefore not just cosmetic — it directly reduces MILP size. Tightening bounds to FVA ranges has the opposite effect: more finite bounds means more constraints.

Omit non-knockable z-variables. The SD formulation creates one binary z per reaction, but only knockable reactions (non-zero cost, non-essential) actually need one. Non-knockable z-variables have ub=0 and are fixed to zero; with solver presolve disabled they become dead weight. _trim_z_variables removes them before the solver sees the problem. On iML1515: 462 → 205 binaries after trimming.

Aggressive compression. Two-pass network compression (before and after GPR extension) reduces the stoichiometric matrix substantially. On iML1515: 2712 → 1920 reactions, 1728 → 974 metabolites. Every removed reaction is one fewer z-variable, one fewer column in all constraint matrices, and fewer indicator constraints.

GPR propagation through compression. Merging GPR rules during compression (via propagate_gpr=True) and only running extend_model_gpr on the compressed model avoids creating gene pseudoreactions for reactions that were already merged away. This reduces both preprocessing time (FVA on ~190 vs ~2700 reactions) and final model size.

Condensed formulation (no reaction splitting). The SD “condensed” formulation handles reversible reactions without splitting them into forward/reverse pairs. This produces fewer z-variables and fewer indicator constraints than the “split” variant. Benchmarks show condensed is consistently faster (4.0s vs 4.1s on iMLcore, but the gap widens on larger models).

Auxiliary-variable indicators do not help. Replacing SD’s complex indicators (z=0 → A·x = b) with FLB-style auxiliary variables (A·x + s = b always, z=0 → s = 0) was tested. The extra continuous variables and equality constraints outweigh the benefit of simpler indicator propagation (~15% slower on iMLcore).

Solver settings matter. Gurobi’s presolve is critical for indicator-constraint MILPs. With Presolve=0 (which was temporarily needed for a Gurobi 13 bug), solve times were ~1.6× slower across all formulations. On iML1515 the full pipeline went from ~12 minutes (presolve off) to ~7 minutes (presolve on).

13. Testing

Tests are in tests/. They are run with:

pytest tests -v --log-cli-level=INFO --junit-xml=test-results.xml

CI matrix (.github/workflows/CI-test.yml):

OS: ubuntu-latest, windows-latest
Python: 3.10, 3.11, 3.12, 3.13
Package managers: pip, conda
CPLEX excluded for Python 3.13 (max supported: 3.12)
Segfault on Ubuntu (JPype JVM shutdown) handled via JUnit XML exit-code check

Key test areas to verify after changes:

Constraint parsing — any change to parse_constr.py must be verified against all input formats.
MILP construction — after changes to link_z, build_primal_from_cbm, or dualization functions: run with a small toy model and verify solutions against known correct answers.
Compression — changes to compression.py or networktools.py must verify that expand_sd correctly reconstructs original-space solutions.
Solver backends — changes to any *_interface.py file require testing with the specific solver installed.
GPR translation — changes to networktools.extend_model_gpr or SDSolutions.gpr_eval must be tested with models that have AND/OR GPR logic (e.g., iJO1366 for E. coli).

Adding new tests:

Use a small toy model (3-5 reactions) for unit tests of MILP construction.
Use e_coli_core for integration tests of the full pipeline.
Use iJO1366 for performance regression tests (track solve time).
Always test with at least GLPK (always available) and one commercial solver if possible.

Performance profiling:

import cProfile
cProfile.run("compute_strain_designs(model, sd_modules=[...], solver='glpk')", 'profile_out')
import pstats
pstats.Stats('profile_out').sort_stats('cumulative').print_stats(30)

The hot spots are typically: link_z (LP bounding), fva (preprocessing), and the solver’s solve() loop.

Quick Reference: Data Flow Through Key Functions

compute_strain_designs(model, **kwargs)
    │
    │── remove_dummy_bounds(model)
    │── cmp_model = model.copy()
    │── remove_ext_mets(cmp_model)
    │── [if reg] split regulatory constraints:
    │       ├── reaction-based → extend_model_regulatory(cmp_model) now
    │       └── gene-based → deferred until after GPR extension
    │
    │── ─── COMPRESS #1 (before GPR extension) ───
    │── compress_model(cmp_model, propagate_gpr=True) → cmp_mapReac_1
    │── compress_modules(sd_modules, cmp_mapReac_1)
    │── compress_ki_ko_cost(rxn_ko, rxn_ki, cmp_mapReac_1)
    │
    │── ─── FVA (on compressed model) ───
    │── bound_blocked_or_irrevers_fva(cmp_model, ...)
    │── fva(cmp_model, ...) → essential_reacs
    │
    │── ─── GPR extension (on compressed model) ───
    │── [if gene_kos] rename_genes(cmp_model, ...)
    │── [if gene_kos] reduce_gpr(cmp_model, essential_reacs, gkis, gkos)
    │── [if gene_kos] extend_model_gpr(cmp_model, ...) → reac_map
    │── [if gene_kos] apply deferred gene-based regulatory constraints
    │── [if gene_kos] merge gene costs into cmp_ko_cost / cmp_ki_cost
    │
    │── ─── COMPRESS #2 (after GPR extension) ───
    │── compress_model(cmp_model) → cmp_mapReac_2
    │── compress_modules(sd_modules, cmp_mapReac_2)
    │── compress_ki_ko_cost(all_ko, all_ki, cmp_mapReac_2)
    │── cmp_mapReac = cmp_mapReac_1 + cmp_mapReac_2
    │
    │── fva(cmp_model, ...) → final essential_reacs
    │
    │── ─── MILP construction & solving ───
    │── SDMILP(cmp_model, sd_modules, ko_cost, ki_cost, solver, M, ...)
    │       ├── SDProblem.__init__(...)
    │       │       ├── [for each module] addModule(sd_module)
    │       │       │       ├── build_primal_from_cbm(model, V_ineq, v_ineq, ...)
    │       │       │       ├── [SUPPRESS] farkas_dualize(...)
    │       │       │       ├── [PROTECT]  reassign_lb_ub_from_ineq(...)
    │       │       │       ├── [OPTKNOCK] LP_dualize(...) + block_diag assembly
    │       │       │       ├── [ROBUSTKNOCK] LP_dualize x2
    │       │       │       └── [OPTCOUPLE] LP_dualize + GCP objective
    │       │       └── link_z()
    │       │               ├── split equalities → inequalities
    │       │               ├── variable KOs → inequality KOs
    │       │               ├── parallel LP bounding → big-M values
    │       │               ├── insert big-M into A_ineq
    │       │               ├── lump duplicate pairs → equalities
    │       │               └── remaining → IndicatorConstraints
    │       └── MILP_LP.__init__(c, A_ineq, ..., solver=solver)
    │               └── instantiate Cplex/Gurobi/SCIP/GLPK backend
    ├── [BEST]     sdmilp.compute_optimal(max_solutions, time_limit)
    ├── [ANY]      sdmilp.compute(max_solutions, time_limit)
    ├── [POPULATE] sdmilp.enumerate(max_solutions, time_limit)
    │       each iteration:
    │       ├── solveZ() / populateZ(n)
    │       ├── verify_sd(z) → valid[]
    │       └── add_exclusion_constraints(z)
    ├── expand_sd(sd, cmp_mapReac)          ← processes pass 2, then pass 1
    ├── filter_sd_maxcost(sd, max_cost, ...)
    ├── postprocess_reg_sd(sd, ...)
    └── SDSolutions(model, sd, status, sd_setup)
            └── [if gene_kos] translate gene_sd → reaction_sd via gpr_eval