diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml
index 2cf9fdc..6de21e4 100644
--- a/.github/workflows/python-publish.yml
+++ b/.github/workflows/python-publish.yml
@@ -38,6 +38,20 @@ jobs:
         pip install yapf
     - name: format repository with yapf in google-python format
       run: yapf -i -r --style ./.github/style.yapf .
+    # Build and validate the distributions before anything is published or pushed:
+    # a failure here leaves the repository untouched, with no tag, release or
+    # version-bump commit to unwind.
+    - name: Install build and twine
+      run: |
+        pip install --upgrade build twine
+    - name: Clean PyPi build directories
+      run: rm -rf dist
+    - name: Build for PyPi
+      run: |
+        python -m build
+    - name: Check the built distributions
+      run: |
+        twine check dist/*
     - name: retrieve name and E-Mail configuration
       run: |
         git config --global user.name ${{ secrets.NAME_GITHUB }}
@@ -52,6 +66,13 @@ jobs:
       run: git fetch
     - name: push code to main
       run: git push
+    - name: Publish to PyPI
+      uses: pypa/gh-action-pypi-publish@release/v1
+      with:
+        password: ${{ secrets.PYPI_API_TOKEN }}
+    # Tag and release only once PyPI has accepted the upload. A version cannot be
+    # re-uploaded to PyPI even after deletion, so the release must not advertise a
+    # version that never made it there.
     - name: Create Release
       env:
         GITHUB_TOKEN: ${{ github.token }}
@@ -59,21 +80,6 @@ jobs:
         gh release create "v${{ github.event.inputs.version }}" \
           --title "v${{ github.event.inputs.version }}" \
           --notes "${{ github.event.inputs.description }}"
-    - name: Install build and twine
-      run: |
-        pip install --upgrade build twine
-    - name: Clean PyPi build directories
-      run: rm -rf dist
-    - name: Build for PyPi
-      run: |
-        python -m build
-    - name: Check the built distributions
-      run: |
-        twine check dist/*
-    - name: Publish to PyPI
-      uses: pypa/gh-action-pypi-publish@release/v1
-      with:
-        password: ${{ secrets.PYPI_API_TOKEN }}
     - name: Install conda dependencies
       run: conda install anaconda-client conda-build
     - name: Clean Anaconda build directories
diff --git a/straindesign/compression.py b/straindesign/compression.py
index 5c08692..a1ba305 100644
--- a/straindesign/compression.py
+++ b/straindesign/compression.py
@@ -101,7 +101,6 @@ def _lcm_list(numbers: List[int]) -> int:
 
 # Rational Matrix with Sparse Storage
 
-
 _INT64_MAX = (1 << 63) - 1
 
 
@@ -372,8 +371,9 @@ def scale_column(self, col: int, scalar_num: int, scalar_den: int) -> None:
         num_lil, den_lil = self._num_sparse, self._den_sparse
         num_csc = num_lil.tocsc() if num_lil.format != 'csc' else num_lil
         den_csc = den_lil.tocsc() if den_lil.format != 'csc' else den_lil
-        entries = [(num_csc.indices[i], int(num_csc.data[i]), int(den_csc.data[i]))
-                   for i in range(num_csc.indptr[col], num_csc.indptr[col + 1])]
+        entries = [
+            (num_csc.indices[i], int(num_csc.data[i]), int(den_csc.data[i])) for i in range(num_csc.indptr[col], num_csc.indptr[col + 1])
+        ]
         for row, cur_num, cur_den in entries:
             if cur_num == 0:
                 continue
@@ -407,10 +407,9 @@ def to_sparse_csr(self) -> Tuple[csr_matrix, int]:
         Returns numerator matrix scaled by LCM of denominators, plus the LCM.
         """
         if self.is_bigint():
-            raise OverflowError(
-                "coefficients exceed int64 and cannot be stored in a scipy sparse matrix; "
-                "use to_coo_exact() / to_sparse_pattern(), or the public sparse_nullspace() helper "
-                "which returns an ExactCOO in that case.")
+            raise OverflowError("coefficients exceed int64 and cannot be stored in a scipy sparse matrix; "
+                                "use to_coo_exact() / to_sparse_pattern(), or the public sparse_nullspace() helper "
+                                "which returns an ExactCOO in that case.")
         # Empty matrix (e.g. a 0-dimensional nullspace)
         if self._den_sparse is None:
             return csr_matrix((self._rows, self._cols), dtype=np.int64), 1
@@ -421,8 +420,7 @@ def to_sparse_csr(self) -> Tuple[csr_matrix, int]:
         common_denom = _lcm_list([int(d) for d in dens if d != 0])
         coo_num = self._num_sparse.tocoo()
         coo_den = self._den_sparse.tocoo()
-        scaled_data = [int(num) * (common_denom // int(den)) if num != 0 else 0
-                       for num, den in zip(coo_num.data, coo_den.data)]
+        scaled_data = [int(num) * (common_denom // int(den)) if num != 0 else 0 for num, den in zip(coo_num.data, coo_den.data)]
         scaled = csr_matrix((scaled_data, (coo_num.row, coo_num.col)), shape=(self._rows, self._cols), dtype=np.int64)
         return scaled, common_denom
 
@@ -443,13 +441,17 @@ def to_coo_exact(self) -> 'ExactCOO':
         if self._dict_frac is not None:
             for r, cd in self._dict_frac.items():
                 for c, f in cd.items():
-                    rows.append(int(r)); cols.append(int(c)); fracs.append(f)
+                    rows.append(int(r))
+                    cols.append(int(c))
+                    fracs.append(f)
         elif self._num_sparse is not None:
             coo_num = self._num_sparse.tocoo()
             coo_den = self._den_sparse.tocoo()
             for r, c, n, d in zip(coo_num.row, coo_num.col, coo_num.data, coo_den.data):
                 if n != 0:
-                    rows.append(int(r)); cols.append(int(c)); fracs.append(Fraction(int(n), int(d)))
+                    rows.append(int(r))
+                    cols.append(int(c))
+                    fracs.append(Fraction(int(n), int(d)))
         denom = _lcm_list([f.denominator for f in fracs]) if fracs else 1
         data = [int(f.numerator) * (denom // f.denominator) for f in fracs]
         return ExactCOO(rows, cols, data, (self._rows, self._cols), denom)
@@ -470,7 +472,9 @@ def to_sparse_pattern(self) -> Tuple[csr_matrix, Dict[int, Dict[int, Fraction]]]
                 rd = {}
                 for c, f in cd.items():
                     if f != 0:
-                        ai.append(int(r)); aj.append(int(c)); rd[int(c)] = f
+                        ai.append(int(r))
+                        aj.append(int(c))
+                        rd[int(c)] = f
                 if rd:
                     row_data[int(r)] = rd
             pattern = csr_matrix(([1] * len(ai), (ai, aj)), shape=(self._rows, self._cols), dtype=np.int8)
@@ -695,7 +699,7 @@ def _eliminate(prd, pivot_val, pivot_col, targets, index):
         _eliminate(prd, pval, pcol, targets, False)
         holders = pivcol_holders.get(pcol)
         if holders:
-            holders.difference_update(above)   # pcol is now cleared from those rows
+            holders.difference_update(above)  # pcol is now cleared from those rows
 
     # Final GCD reduction of the pivot rows (insurance; rows are already reduced per step).
     for key in pivot_keys:
@@ -706,14 +710,12 @@ def _eliminate(prd, pivot_val, pivot_col, targets, index):
                 row_data[c] //= row_gcd
 
     # Translate results back to original column space, keyed by pivot index (rref_data[i] = pivot i).
-    original_data = {i: {col_order[sc]: v for sc, v in data[key].items()}
-                     for i, key in enumerate(pivot_keys)}
+    original_data = {i: {col_order[sc]: v for sc, v in data[key].items()} for i, key in enumerate(pivot_keys)}
     pivot_cols_original = [col_order[p] for p in pivot_cols_sorted]
 
     return original_data, rank, pivot_cols_original
 
 
-
 def _nullspace_sparse(matrix: RationalMatrix) -> RationalMatrix:
     """Compute nullspace using integer RREF with row scaling.
 
@@ -940,7 +942,10 @@ def __init__(self, metas: int, reacs: int):
 class _WorkRecord:
     """Mutable state during compression algorithm."""
 
-    def __init__(self, stoich: RationalMatrix, meta_names: List[str], reac_names: List[str],
+    def __init__(self,
+                 stoich: RationalMatrix,
+                 meta_names: List[str],
+                 reac_names: List[str],
                  bounds: Optional[List[Tuple[float, float]]] = None):
         rows, cols = stoich.get_row_count(), stoich.get_column_count()
         self.pre = RationalMatrix.identity(rows)
@@ -1346,10 +1351,9 @@ def _small_bound(r):
                     return min(fin) if fin else None
 
                 def _lam(r):
-                    return 1.0 if ratios[r] is None else float(abs(ratios[r]))   # master's own ratio is 1
+                    return 1.0 if ratios[r] is None else float(abs(ratios[r]))  # master's own ratio is 1
 
-                bounded = [(b, r) for r in group for b in [_small_bound(r)]
-                           if b is not None and 0.1 <= _lam(r) <= 10]
+                bounded = [(b, r) for r in group for b in [_small_bound(r)] if b is not None and 0.1 <= _lam(r) <= 10]
                 keep = min(bounded)[1] if bounded else max(group, key=lambda r: nnz[r])
                 self._restore_group_scale(work, group, ratios, keep)
 
@@ -1363,8 +1367,7 @@ def _lam(r):
 
         return contradicting_removed
 
-    def _restore_group_scale(self, work: _WorkRecord, group: List[int],
-                             ratios: List[Optional[Fraction]], keep: int) -> None:
+    def _restore_group_scale(self, work: _WorkRecord, group: List[int], ratios: List[Optional[Fraction]], keep: int) -> None:
         """Express a merged group in the units of one of its members.
 
         A lump's ratios are fixed but its overall scale is free, and merging into ``group[0]`` can
@@ -1376,7 +1379,7 @@ def _restore_group_scale(self, work: _WorkRecord, group: List[int],
         master = group[0]
         if keep == master:
             return
-        lam = abs(ratios[keep])                     # |.| so the reaction keeps its orientation
+        lam = abs(ratios[keep])  # |.| so the reaction keeps its orientation
         if lam == 0 or lam == 1:
             return
         # v_master = ratios[keep] * v_keep, so re-expressing the lump in v_keep multiplies the
@@ -1509,11 +1512,12 @@ def remove_conservation_relations(model) -> None:
         model.remove_metabolites(dependent_mets)
 
 
-def compress_cobra_model(model,
-                         methods: Optional[List[Union[str, CompressionMethod]]] = None,
-                         in_place: bool = True,
-                         suppressed_reactions: Set[str] = set(),
-                         protected_reactions: Set[str] = set()) -> CompressionResult:
+def compress_cobra_model(
+    model,
+    methods: Optional[List[Union[str, CompressionMethod]]] = None,
+    in_place: bool = True,
+    suppressed_reactions: Set[str] = set(),
+    protected_reactions: Set[str] = set()) -> CompressionResult:
     """
     Compress a COBRA model using nullspace-based coupling detection.
 
@@ -1561,7 +1565,8 @@ def compress_cobra_model(model,
     # Run compression
     compressor = StoichMatrixCompressor(*compression_methods)
     bounds = [(float(r.lower_bound), float(r.upper_bound)) for r in model.reactions]
-    compression_record = compressor.compress(stoich_matrix, metabolite_names, reaction_names, suppressed_reactions, bounds, protected_reactions)
+    compression_record = compressor.compress(stoich_matrix, metabolite_names, reaction_names, suppressed_reactions, bounds,
+                                             protected_reactions)
 
     # Apply to model (uses direct manipulation, bypasses solver)
     reaction_map = _apply_compression_to_model(model, compression_record, reaction_names)
@@ -1691,10 +1696,7 @@ def _apply_compression_to_model(model, compression_record, original_reaction_nam
         # Update stored objective through compression factors
         obj_dict = getattr(model, '_suppressed_obj', None)
         if obj_dict is not None:
-            merged_obj = sum(
-                obj_dict.pop(original_reaction_names[idx], 0.0) * float(coeff)
-                for idx, coeff in contributing
-            )
+            merged_obj = sum(obj_dict.pop(original_reaction_names[idx], 0.0) * float(coeff) for idx, coeff in contributing)
             if merged_obj != 0:
                 obj_dict[main_rxn.id] = merged_obj
 
@@ -1784,12 +1786,12 @@ def stoichmat_coeff_to_fraction(model) -> None:
     for rxn in model.reactions:
         for met, coeff in rxn._metabolites.items():
             if isinstance(coeff, Fraction):
-                continue                                        # already exact
+                continue  # already exact
             elif isinstance(coeff, (float, int)):
-                rxn._metabolites[met] = float_to_fraction(coeff)          # -> Fraction
-            elif hasattr(coeff, 'p'):                           # sympy.Rational -> Fraction
+                rxn._metabolites[met] = float_to_fraction(coeff)  # -> Fraction
+            elif hasattr(coeff, 'p'):  # sympy.Rational -> Fraction
                 rxn._metabolites[met] = Fraction(int(coeff.p), int(coeff.q))
-            elif hasattr(coeff, 'numerator'):                   # other Rational -> Fraction
+            elif hasattr(coeff, 'numerator'):  # other Rational -> Fraction
                 rxn._metabolites[met] = Fraction(coeff.numerator, coeff.denominator)
             else:
                 raise TypeError(f"Unsupported coefficient type: {type(coeff)}")
@@ -1848,8 +1850,10 @@ def _expr_to_gpr_string(expr):
         return expr
     op, args = expr
     other = 'or' if op == 'and' else 'and'
-    parts = [f'({s})' if isinstance(a, tuple) and a[0] == other else s
-             for a, s in sorted(((a, _expr_to_gpr_string(a)) for a in args), key=lambda p: p[1])]
+    parts = [
+        f'({s})' if isinstance(a, tuple) and a[0] == other else s
+        for a, s in sorted(((a, _expr_to_gpr_string(a)) for a in args), key=lambda p: p[1])
+    ]
     return f' {op} '.join(parts)
 
 
@@ -1868,7 +1872,6 @@ def _combine_gprs(gpr_bodies, op):
 
 # High-Level Compression API
 
-
 # Monotone (positive-unate) GPR-rule simplification
 #
 # Pipeline:  parse -> minimal SOP (DNF + absorption) -> algebraic factoring.
@@ -1896,35 +1899,63 @@ def _gpr_parse(s):
     Accepts both ``and``/``or`` and ``*``/``+`` operators, and is robust to any gene id,
     including digit-leading or dotted names.
     """
-    toks = list(_gpr_tokenize(s)); pos = 0
-    def peek(): return toks[pos] if pos < len(toks) else None
+    toks = list(_gpr_tokenize(s))
+    pos = 0
+
+    def peek():
+        return toks[pos] if pos < len(toks) else None
+
     def eat():
-        nonlocal pos; t = toks[pos]; pos += 1; return t
+        nonlocal pos
+        t = toks[pos]
+        pos += 1
+        return t
+
     def p_or():
         n = [p_and()]
-        while peek() in ('or', '+'): eat(); n.append(p_and())
+        while peek() in ('or', '+'):
+            eat()
+            n.append(p_and())
         return ('OR', n) if len(n) > 1 else n[0]
+
     def p_and():
         n = [p_atom()]
-        while peek() in ('and', '*'): eat(); n.append(p_atom())
+        while peek() in ('and', '*'):
+            eat()
+            n.append(p_atom())
         return ('AND', n) if len(n) > 1 else n[0]
+
     def p_atom():
-        if peek() == '(': eat(); e = p_or(); eat(); return e
+        if peek() == '(':
+            eat()
+            e = p_or()
+            eat()
+            return e
         return ('VAR', eat())
+
     return p_or()
 
 
 # variable <-> bit mapping (reset per rule via simplify_gpr_string)
-_GPR_VMAP = {}; _GPR_VINV = []
+_GPR_VMAP = {}
+_GPR_VINV = []
+
+
 def _gpr_bit(v):
     i = _GPR_VMAP.get(v)
     if i is None:
-        i = len(_GPR_VINV); _GPR_VMAP[v] = i; _GPR_VINV.append(v)
+        i = len(_GPR_VINV)
+        _GPR_VMAP[v] = i
+        _GPR_VINV.append(v)
     return 1 << i
+
+
 def _gpr_lits_of(mask):
     out = []
     while mask:
-        l = mask & -mask; out.append(('VAR', _GPR_VINV[l.bit_length() - 1])); mask ^= l
+        l = mask & -mask
+        out.append(('VAR', _GPR_VINV[l.bit_length() - 1]))
+        mask ^= l
     return out
 
 
@@ -1945,11 +1976,12 @@ def _gpr_absorb(cubes):
 
 def _gpr_to_dnf(node):
     t = node[0]
-    if t == 'VAR':   return [_gpr_bit(node[1])]
+    if t == 'VAR': return [_gpr_bit(node[1])]
     if t == 'CONST': return [] if not node[1] else [0]
     if t == 'OR':
         cov = []
-        for ch in node[1]: cov += _gpr_to_dnf(ch)
+        for ch in node[1]:
+            cov += _gpr_to_dnf(ch)
         return _gpr_absorb(cov)
     if t == 'AND':
         cov = [0]
@@ -1961,8 +1993,10 @@ def _gpr_to_dnf(node):
 
 
 def _gpr_common(cubes):
-    it = iter(cubes); c = next(it)
-    for x in it: c &= x
+    it = iter(cubes)
+    c = next(it)
+    for x in it:
+        c &= x
     return c
 
 
@@ -1971,7 +2005,9 @@ def _gpr_lit_counts(F):
     for c in F:
         m = c
         while m:
-            l = m & -m; cnt[l] = cnt.get(l, 0) + 1; m ^= l
+            l = m & -m
+            cnt[l] = cnt.get(l, 0) + 1
+            m ^= l
     return cnt
 
 
@@ -1991,18 +2027,21 @@ def _gpr_candidate_divisors(F):
     if len(F) < 2: return []
     cnt = _gpr_lit_counts(F)
     reps = sorted((x for x, n in cnt.items() if n >= 2), key=lambda x: -cnt[x])
-    seen = set(); out = []
+    seen = set()
+    out = []
     for l in reps:
         K = tuple(sorted(_gpr_one_kernel(F, l)))
         if len(K) >= 2 and K not in seen:
-            seen.add(K); out.append(list(K))
+            seen.add(K)
+            out.append(list(K))
     return out
 
 
 def _gpr_divide(F, D):
     """Exact algebraic division: (Q, R) with D*Q disjoint-union R == F (correctness guaranteed
     regardless of divisor quality -- a quotient cube is accepted only if D*Q stays inside F)."""
-    Fs = set(F); quo = None
+    Fs = set(F)
+    quo = None
     for d in D:
         vd = {c & ~d for c in F if (c & d) == d}
         quo = vd if quo is None else (quo & vd)
@@ -2015,8 +2054,8 @@ def _gpr_divide(F, D):
 
 def _gpr_factor(F):
     F = _gpr_absorb(F)
-    if not F:     return ('CONST', False)
-    if F == [0]:  return ('CONST', True)
+    if not F: return ('CONST', False)
+    if F == [0]: return ('CONST', True)
     if len(F) == 1:
         lits = _gpr_lits_of(F[0])
         return lits[0] if len(lits) == 1 else ('AND', lits)
@@ -2044,9 +2083,9 @@ def _gpr_factor(F):
 def _gpr_est_cubes(node):
     """Upper bound on DNF cube count (product across ANDs, sum across ORs); cheap, no expansion."""
     t = node[0]
-    if t == 'VAR':   return 1
+    if t == 'VAR': return 1
     if t == 'CONST': return 1
-    if t == 'OR':    return sum(_gpr_est_cubes(c) for c in node[1])
+    if t == 'OR': return sum(_gpr_est_cubes(c) for c in node[1])
     if t == 'AND':
         p = 1
         for c in node[1]:
@@ -2056,6 +2095,8 @@ def _gpr_est_cubes(node):
 
 
 _GPR_WARN = []
+
+
 def _gpr_factor_auto(node, budget=50000):
     """Global factoring within budget; AND-split above it. Never splits an OR unless one single
     OR-block alone exceeds budget (logged as a last resort -- raise the budget to avoid)."""
@@ -2083,7 +2124,9 @@ def simplify_gpr_string(rule, budget=50000):
     """Return a leaf-minimized, boolean-equivalent monotone GPR string ('' passes through)."""
     if not rule or not rule.strip():
         return rule
-    _GPR_VMAP.clear(); _GPR_VINV.clear(); _GPR_WARN.clear()
+    _GPR_VMAP.clear()
+    _GPR_VINV.clear()
+    _GPR_WARN.clear()
     return _gpr_to_string(_gpr_factor_auto(_gpr_parse(rule), budget))
 
 
@@ -2102,14 +2145,14 @@ def simplify_model_gprs(model, budget=50000):
         try:
             new = simplify_gpr_string(s, budget)
             if new and new != s:
-                r.gene_reaction_rule = new; nchg += 1
+                r.gene_reaction_rule = new
+                nchg += 1
         except Exception as e:
             logging.warning('gpr_simplify: kept original GPR for %s (%s)' % (r.id, type(e).__name__))
     logging.info('  GPR rule simplification: %d rules, %d rewritten.' % (n, nchg))
 
 
-def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False,
-                   no_coupled_compress_reacs=set()):
+def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False, no_coupled_compress_reacs=set()):
     """Compress a metabolic model using multiple techniques.
 
     Performs blocked reaction removal, conservation relation removal, and
@@ -2148,8 +2191,7 @@ def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False,
 
             # 1. Parallel (cheap — hash-based, no RREF)
             LOG.info(f'  Compression {run}: Lumping parallel reactions.')
-            reac_map_exp = compress_model_parallel(model, no_par_compress_reacs,
-                                                    propagate_gpr=propagate_gpr)
+            reac_map_exp = compress_model_parallel(model, no_par_compress_reacs, propagate_gpr=propagate_gpr)
             parallel_changed = numr > len(reac_map_exp)
             if parallel_changed:
                 LOG.info(f'  Reduced to {len(reac_map_exp)} reactions.')
@@ -2169,9 +2211,7 @@ def compress_model(model, no_par_compress_reacs=set(), propagate_gpr=False,
             # 4. Coupled (expensive — nullspace/RREF)
             numr_pre = len(model.reactions)
             LOG.info(f'  Compression {run}: Lumping coupled reactions.')
-            reac_map_exp = compress_model_coupled(model,
-                                                  propagate_gpr=propagate_gpr,
-                                                  protected_reactions=no_coupled_compress_reacs)
+            reac_map_exp = compress_model_coupled(model, propagate_gpr=propagate_gpr, protected_reactions=no_coupled_compress_reacs)
             for new_reac, old_reac_val in reac_map_exp.items():
                 old_reacs = [r for r in no_par_compress_reacs if r in old_reac_val]
                 if old_reacs:
@@ -2228,8 +2268,7 @@ def compress_model_coupled(model, propagate_gpr=False, protected_reactions=set()
         for r in model.reactions:
             r.gene_reaction_rule = ''
 
-        result = compress_cobra_model(model, methods=CompressionMethod.standard(), in_place=True,
-                                      protected_reactions=protected_reactions)
+        result = compress_cobra_model(model, methods=CompressionMethod.standard(), in_place=True, protected_reactions=protected_reactions)
         reaction_map = result.reaction_map
 
         # Propagate GPR rules: AND-combine contributing reactions' GPR ASTs
@@ -2270,8 +2309,7 @@ def compress_model_parallel(model, protected_rxns=set(), propagate_gpr=False):
     ub = [float(r.upper_bound) for r in model.reactions]
     fwd = [1 if (isinf(u) and f > 0 or isinf(l) and f < 0) else 0 for f, l, u in zip(factor, lb, ub)]
     rev = [1 if (isinf(l) and f > 0 or isinf(u) and f < 0) else 0 for f, l, u in zip(factor, lb, ub)]
-    inh = [i + 1 if not ((isinf(ub[i]) or ub[i] == 0) and (isinf(lb[i]) or lb[i] == 0)) else 0
-           for i in range(len(model.reactions))]
+    inh = [i + 1 if not ((isinf(ub[i]) or ub[i] == 0) and (isinf(lb[i]) or lb[i] == 0)) else 0 for i in range(len(model.reactions))]
 
     # Canonical scale-invariant key per reaction: normalize the stoichiometry row by its first
     # nonzero coefficient in exact rational arithmetic, so reactions parallel up to any rational
@@ -2350,10 +2388,7 @@ def _parallel_key(i):
         else:
             scales = [abs(factor[j]) for j in group]
             total = sum(Fraction(s).limit_denominator(1000) for s in scales)
-            rational_map[model.reactions[i].id] = {
-                old_reac_ids[j]: Fraction(abs(factor[j])).limit_denominator(1000) / total
-                for j in group
-            }
+            rational_map[model.reactions[i].id] = {old_reac_ids[j]: Fraction(abs(factor[j])).limit_denominator(1000) / total for j in group}
 
     return rational_map
 
diff --git a/straindesign/compute_strain_designs.py b/straindesign/compute_strain_designs.py
index 6db62c6..4f5c167 100644
--- a/straindesign/compute_strain_designs.py
+++ b/straindesign/compute_strain_designs.py
@@ -67,8 +67,9 @@ def _essentials_from_limits(flux_limits):
     ``flux_limits`` is an FVA result over the module's constrained polytope. Both bounds must share
     a sign and stay clear of zero, so the reaction carries flux in every point of the module.
     """
-    return {reac_id for reac_id, limits in flux_limits.iterrows()
-            if np.min(abs(limits)) > _ESSENTIALITY_TOL and np.prod(np.sign(limits)) > 0}
+    return {
+        reac_id for reac_id, limits in flux_limits.iterrows() if np.min(abs(limits)) > _ESSENTIALITY_TOL and np.prod(np.sign(limits)) > 0
+    }
 
 
 def reduce_model_gprs(model, essential_reacs, gkis, gkos):
@@ -265,7 +266,6 @@ def is_gene_essential_to_reaction_ast(reaction, gene_id):
 
 
 @with_suppressed_lp
-
 def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
     """Computes strain designs for a user-defined strain design problem
 
@@ -577,20 +577,21 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
         for r in cmp_model.reactions:
             can_fwd, can_rev = _rev[r.id]
             if not can_fwd and float(r._upper_bound) > 0.0:
-                r._upper_bound = min(0.0, float(r._upper_bound)); _n_tight += 1
+                r._upper_bound = min(0.0, float(r._upper_bound))
+                _n_tight += 1
             if not can_rev and float(r._lower_bound) < 0.0:
-                r._lower_bound = max(0.0, float(r._lower_bound)); _n_tight += 1
-        logging.info('  Reversibility pre-tightening fixed %d reaction directions (%.1fs).'
-                     % (_n_tight, time.time() - t0))
+                r._lower_bound = max(0.0, float(r._lower_bound))
+                _n_tight += 1
+        logging.info('  Reversibility pre-tightening fixed %d reaction directions (%.1fs).' % (_n_tight, time.time() - t0))
         logging.info('Compressing Network (' + str(len(cmp_model.reactions)) + ' reactions).')
         t0 = time.time()
-        cmp_mapReac_1 = compress_model(cmp_model, no_par_compress_reacs,
-                                        propagate_gpr=True,
-                                        no_coupled_compress_reacs=no_coupled_compress_reacs)
+        cmp_mapReac_1 = compress_model(cmp_model,
+                                       no_par_compress_reacs,
+                                       propagate_gpr=True,
+                                       no_coupled_compress_reacs=no_coupled_compress_reacs)
         sd_modules = compress_modules(sd_modules, cmp_mapReac_1)
         # Compress reaction + regulatory costs only (gene costs not yet added)
-        cmp_ko_cost, cmp_ki_cost, cmp_mapReac_1 = compress_ki_ko_cost(
-            uncmp_ko_cost, uncmp_ki_cost, cmp_mapReac_1)
+        cmp_ko_cost, cmp_ki_cost, cmp_mapReac_1 = compress_ki_ko_cost(uncmp_ko_cost, uncmp_ki_cost, cmp_mapReac_1)
         logging.info('  Compressed to ' + str(len(cmp_model.reactions)) + ' reactions (%.1fs).' % (time.time() - t0))
     else:
         cmp_mapReac_1 = []
@@ -627,8 +628,7 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
                              len([True for r in cmp_model.reactions if r.gene_reaction_rule]) < num_gpr):
             num_genes = len(cmp_model.genes)
             num_gpr = len([True for r in cmp_model.reactions if r.gene_reaction_rule])
-            logging.info('  Simplified to ' + str(num_genes) + ' genes and ' +
-                str(num_gpr) + ' gpr rules (%.1fs).' % (time.time() - t_gpr))
+            logging.info('  Simplified to ' + str(num_genes) + ' genes and ' + str(num_gpr) + ' gpr rules (%.1fs).' % (time.time() - t_gpr))
         logging.info('  Extending metabolic network with gpr associations.')
         reac_map = extend_model_gpr(cmp_model, has_gene_names)
         for i, m in enumerate(sd_modules):
@@ -668,11 +668,12 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
         logging.info('Compressing after GPR extension (' + str(len(cmp_model.reactions)) + ' reactions).')
         t0 = time.time()
         no_par_compress_reacs = _collect_no_par_compress_reacs(sd_modules)
-        cmp_mapReac_2 = compress_model(cmp_model, no_par_compress_reacs,
-)
+        cmp_mapReac_2 = compress_model(
+            cmp_model,
+            no_par_compress_reacs,
+        )
         sd_modules = compress_modules(sd_modules, cmp_mapReac_2)
-        cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2 = compress_ki_ko_cost(
-            cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2)
+        cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2 = compress_ki_ko_cost(cmp_ko_cost, cmp_ki_cost, cmp_mapReac_2)
         cmp_mapReac = cmp_mapReac_1 + cmp_mapReac_2
         logging.info('  Compressed to ' + str(len(cmp_model.reactions)) + ' reactions (%.1fs).' % (time.time() - t0))
     else:
@@ -684,14 +685,15 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
     t0 = time.time()
     # Save pre-FVA bounds for dump_preprocessed (bound config experiments)
     pre_fva_bounds = {r.id: (r.lower_bound, r.upper_bound) for r in cmp_model.reactions}
-    # Subset-scope the bound-stripping FVA off reactions already at (0,+inf). 
-    # FVA on such a reaction can only ever tighten a genuinely blocked one 
-    # to (0,0); leaving it at (0,+inf) changes no feasible flux (the stoichiometry 
+    # Subset-scope the bound-stripping FVA off reactions already at (0,+inf).
+    # FVA on such a reaction can only ever tighten a genuinely blocked one
+    # to (0,0); leaving it at (0,+inf) changes no feasible flux (the stoichiometry
     # already forces it to 0) and a blocked knockable can never belong to a minimal cut set.
-    _fva_scope = [r.id for r in cmp_model.reactions
-                  if not (float(r.lower_bound) == 0.0
-                          and np.isinf(float(r.upper_bound))
-                          and float(r.upper_bound) > 0)]
+    _fva_scope = [
+        r.id
+        for r in cmp_model.reactions
+        if not (float(r.lower_bound) == 0.0 and np.isinf(float(r.upper_bound)) and float(r.upper_bound) > 0)
+    ]
     essential_reacs = set()
     suppress_essential = set()
     cmp_size1_mcs = []
@@ -700,19 +702,17 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
     # With exactly one classical module, one FVA over the constrained module polytope can serve both
     # model-bound tightening and module essentiality. This is only sound for a single module: applying
     # one module's tighter ranges to the shared model could otherwise alter another module's polytope.
-    fold_module_fva = (
-        len(sd_modules) == 1
-        and sd_modules[0][MODULE_TYPE] in [SUPPRESS, PROTECT]
-        and sd_modules[0][INNER_OBJECTIVE] is None
-    )
+    fold_module_fva = (len(sd_modules) == 1 and sd_modules[0][MODULE_TYPE] in [SUPPRESS, PROTECT] and
+                       sd_modules[0][INNER_OBJECTIVE] is None)
     if fold_module_fva:
         module = sd_modules[0]
         fold_scope = sorted(set(_fva_scope) | set(knockable_ids))
-        flux_limits = bound_blocked_or_irrevers_fva(
-            cmp_model, solver=kwargs[SOLVER], constraints=module[CONSTRAINTS],
-            compress=False, reaction_list=fold_scope)
-        module_limits = flux_limits.loc[
-            [reac_id for reac_id in knockable_ids if reac_id in flux_limits.index]]
+        flux_limits = bound_blocked_or_irrevers_fva(cmp_model,
+                                                    solver=kwargs[SOLVER],
+                                                    constraints=module[CONSTRAINTS],
+                                                    compress=False,
+                                                    reaction_list=fold_scope)
+        module_limits = flux_limits.loc[[reac_id for reac_id in knockable_ids if reac_id in flux_limits.index]]
         module['fva_bounds'] = module_limits
         essentials_in_module = _essentials_from_limits(module_limits)
         if module[MODULE_TYPE] == SUPPRESS:
@@ -721,8 +721,7 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
             essential_reacs.update(essentials_in_module)
         logging.info('  Folded model/module FVA done (%.1fs).' % (time.time() - t0))
     else:
-        bound_blocked_or_irrevers_fva(
-            cmp_model, solver=kwargs[SOLVER], compress=False, reaction_list=_fva_scope)
+        bound_blocked_or_irrevers_fva(cmp_model, solver=kwargs[SOLVER], compress=False, reaction_list=_fva_scope)
         logging.info('  FVA done (%.1fs).' % (time.time() - t0))
 
         # FVA to identify essential reactions and size-1 MCS before building MILP
@@ -733,8 +732,11 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
         # module and flows to SDMILP via sd_modules. Scoping to knockable reactions keeps
         # the LP count down (and only knockable reactions carry z-links to tighten anyway).
         for module in sd_modules:
-            flux_limits = fva(cmp_model, solver=kwargs[SOLVER], constraints=module[CONSTRAINTS],
-                              compress=False, reaction_list=knockable_ids)
+            flux_limits = fva(cmp_model,
+                              solver=kwargs[SOLVER],
+                              constraints=module[CONSTRAINTS],
+                              compress=False,
+                              reaction_list=knockable_ids)
             module['fva_bounds'] = flux_limits
             essentials_in_module = _essentials_from_limits(flux_limits)
             if module[MODULE_TYPE] == SUPPRESS:
@@ -814,34 +816,35 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
         import os, pickle as _pickle
         dump_path = dump_preprocessed
         with open(dump_path, 'wb') as f:
-            _pickle.dump({
-                'cmp_model': cmp_model,
-                'sd_modules': sd_modules,
-                'kwargs_milp': kwargs_milp,
-                'kwargs_computation': kwargs_computation,
-                'solution_approach': solution_approach,
-                'cmp_mapReac': cmp_mapReac,
-                # Expansion/filtering data
-                'uncmp_ko_cost': uncmp_ko_cost,
-                'uncmp_ki_cost': uncmp_ki_cost,
-                'uncmp_reg_cost': uncmp_reg_cost,
-                'orig_model': orig_model,
-                'orig_sd_modules': orig_sd_modules,
-                'orig_ko_cost': orig_ko_cost,
-                'orig_ki_cost': orig_ki_cost,
-                'orig_reg_cost': orig_reg_cost,
-                'gene_kos': kwargs['gene_kos'],
-                'orig_gko_cost': locals().get('orig_gko_cost'),
-                'orig_gki_cost': locals().get('orig_gki_cost'),
-                'max_cost': kwargs[MAX_COST],
-                'cmp_size1_mcs': cmp_size1_mcs,
-                'pre_fva_bounds': pre_fva_bounds,
-            }, f)
+            _pickle.dump(
+                {
+                    'cmp_model': cmp_model,
+                    'sd_modules': sd_modules,
+                    'kwargs_milp': kwargs_milp,
+                    'kwargs_computation': kwargs_computation,
+                    'solution_approach': solution_approach,
+                    'cmp_mapReac': cmp_mapReac,
+                    # Expansion/filtering data
+                    'uncmp_ko_cost': uncmp_ko_cost,
+                    'uncmp_ki_cost': uncmp_ki_cost,
+                    'uncmp_reg_cost': uncmp_reg_cost,
+                    'orig_model': orig_model,
+                    'orig_sd_modules': orig_sd_modules,
+                    'orig_ko_cost': orig_ko_cost,
+                    'orig_ki_cost': orig_ki_cost,
+                    'orig_reg_cost': orig_reg_cost,
+                    'gene_kos': kwargs['gene_kos'],
+                    'orig_gko_cost': locals().get('orig_gko_cost'),
+                    'orig_gki_cost': locals().get('orig_gki_cost'),
+                    'max_cost': kwargs[MAX_COST],
+                    'cmp_size1_mcs': cmp_size1_mcs,
+                    'pre_fva_bounds': pre_fva_bounds,
+                },
+                f)
         logging.info('Preprocessed data saved to ' + dump_path)
         logging.info('  Resume with:')
         logging.info('    from straindesign import compute_strain_designs_from_preprocessed')
-        logging.info("    sol = compute_strain_designs_from_preprocessed('%s', seed=42)" %
-                     dump_path.replace('\\', '\\\\'))
+        logging.info("    sol = compute_strain_designs_from_preprocessed('%s', seed=42)" % dump_path.replace('\\', '\\\\'))
 
         # Return early with size-1 MCS only (or empty)
         setup = deepcopy(cmp_sd_solution.sd_setup) if 'cmp_sd_solution' in dir() else {MODEL_ID: orig_model.id}
@@ -885,15 +888,14 @@ def compute_strain_designs(model: Model, **kwargs: dict) -> SDSolutions:
     if kwargs['gene_kos']:
         setup.update({GKOCOST: orig_gko_cost, GKICOST: orig_gki_cost})
 
-    sd_solutions = _decompress_solutions(
-        cmp_sd_solution, cmp_mapReac, cmp_size1_mcs,
-        kwargs[MAX_COST], uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost,
-        orig_model, setup, kwargs['gene_kos'],
-        locals().get('orig_gko_cost'), locals().get('orig_gki_cost'))
+    sd_solutions = _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, kwargs[MAX_COST], uncmp_ko_cost, uncmp_ki_cost,
+                                         uncmp_reg_cost, orig_model, setup, kwargs['gene_kos'],
+                                         locals().get('orig_gko_cost'),
+                                         locals().get('orig_gki_cost'))
     sd_solutions._cmp_model = cmp_model
-    logging.info(str(sd_solutions.get_num_materialized()) + ' solutions found'
-                 + (' (lazy, estimated %d total).' % sd_solutions.get_num_sols()
-                    if sd_solutions.is_lazy else '.'))
+    logging.info(
+        str(sd_solutions.get_num_materialized()) + ' solutions found' +
+        (' (lazy, estimated %d total).' % sd_solutions.get_num_sols() if sd_solutions.is_lazy else '.'))
 
     return sd_solutions
 
@@ -915,9 +917,8 @@ def postprocess_reg_sd(reg_cost, sd):
 LAZY_EXPANSION_THRESHOLD = 100_000
 
 
-def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs,
-                          max_cost, uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost,
-                          orig_model, setup, gene_kos, orig_gko_cost, orig_gki_cost):
+def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, max_cost, uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost, orig_model,
+                          setup, gene_kos, orig_gko_cost, orig_gki_cost):
     """Decompress MILP solutions, using lazy expansion if estimated count exceeds threshold."""
     logging.info('  Decompressing.')
 
@@ -933,9 +934,8 @@ def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs,
 
     if estimated > LAZY_EXPANSION_THRESHOLD:
         logging.info('  Estimated %d expanded solutions - using lazy expansion.' % estimated)
-        sd, group_map, compressed_sd = _build_lazy_representatives(
-            cmp_sds, cmp_size1_mcs, cmp_mapReac, max_cost,
-            uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost)
+        sd, group_map, compressed_sd = _build_lazy_representatives(cmp_sds, cmp_size1_mcs, cmp_mapReac, max_cost, uncmp_ko_cost,
+                                                                   uncmp_ki_cost, uncmp_reg_cost)
 
         status = cmp_sd_solution.status
         if status not in [OPTIMAL, TIME_LIMIT_W_SOL] and sd:
@@ -995,8 +995,7 @@ def _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs,
     return sd_solutions
 
 
-def _build_lazy_representatives(cmp_sds, cmp_size1_mcs, cmp_mapReac, max_cost,
-                                uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost):
+def _build_lazy_representatives(cmp_sds, cmp_size1_mcs, cmp_mapReac, max_cost, uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost):
     """Build one representative expanded solution per compressed group.
 
     Returns (sd, group_map, compressed_sd).
@@ -1033,9 +1032,7 @@ def _build_lazy_representatives(cmp_sds, cmp_size1_mcs, cmp_mapReac, max_cost,
     return sd, group_map, compressed_sd
 
 
-def compute_strain_designs_from_preprocessed(dump, seed=None, solver=None,
-                                             solution_approach=None, max_solutions=None,
-                                             time_limit=None):
+def compute_strain_designs_from_preprocessed(dump, seed=None, solver=None, solution_approach=None, max_solutions=None, time_limit=None):
     """Load preprocessed model and run MILP solve with optional overrides.
 
     Args:
@@ -1094,8 +1091,7 @@ def compute_strain_designs_from_preprocessed(dump, seed=None, solver=None,
     from straindesign.networktools import suppress_lp_context
     with suppress_lp_context(cmp_model):
         logging.info('Loading preprocessed data from ' + (dump if isinstance(dump, str) else 'dict input.'))
-        logging.info('  Seed: %s, Solver: %s, Approach: %s' % (
-            kwargs_milp.get(SEED), kwargs_milp.get(SOLVER), sol_approach))
+        logging.info('  Seed: %s, Solver: %s, Approach: %s' % (kwargs_milp.get(SEED), kwargs_milp.get(SOLVER), sol_approach))
 
         t0 = time.time()
         sd_milp = SDMILP(cmp_model, sd_modules, **kwargs_milp)
@@ -1116,13 +1112,11 @@ def compute_strain_designs_from_preprocessed(dump, seed=None, solver=None,
     if gene_kos:
         setup.update({GKOCOST: orig_gko_cost, GKICOST: orig_gki_cost})
 
-    sd_solutions = _decompress_solutions(
-        cmp_sd_solution, cmp_mapReac, cmp_size1_mcs,
-        max_cost, uncmp_ko_cost, uncmp_ki_cost, uncmp_reg_cost,
-        orig_model, setup, gene_kos, orig_gko_cost, orig_gki_cost)
+    sd_solutions = _decompress_solutions(cmp_sd_solution, cmp_mapReac, cmp_size1_mcs, max_cost, uncmp_ko_cost, uncmp_ki_cost,
+                                         uncmp_reg_cost, orig_model, setup, gene_kos, orig_gko_cost, orig_gki_cost)
     sd_solutions._cmp_model = cmp_model
-    logging.info(str(sd_solutions.get_num_materialized()) + ' solutions found'
-                 + (' (lazy, estimated %d total).' % sd_solutions.get_num_sols()
-                    if sd_solutions.is_lazy else '.'))
+    logging.info(
+        str(sd_solutions.get_num_materialized()) + ' solutions found' +
+        (' (lazy, estimated %d total).' % sd_solutions.get_num_sols() if sd_solutions.is_lazy else '.'))
 
     return sd_solutions
diff --git a/straindesign/cplex_interface.py b/straindesign/cplex_interface.py
index 11453d4..e064b1e 100644
--- a/straindesign/cplex_interface.py
+++ b/straindesign/cplex_interface.py
@@ -95,7 +95,18 @@ class Cplex_MILP_LP(Cplex):
             A CPLEX MILP/LP interface class.
     """
 
-    def __init__(self, c=None, A_ineq=None, b_ineq=None, A_eq=None, b_eq=None, lb=None, ub=None, vtype=None, indic_constr=None, seed=None, milp_threads=None):
+    def __init__(self,
+                 c=None,
+                 A_ineq=None,
+                 b_ineq=None,
+                 A_eq=None,
+                 b_eq=None,
+                 lb=None,
+                 ub=None,
+                 vtype=None,
+                 indic_constr=None,
+                 seed=None,
+                 milp_threads=None):
         super().__init__()
         self.objective.set_sense(self.objective.sense.minimize)
         try:
@@ -331,14 +342,12 @@ def set_lp_method(self, method):
             method: LP_METHOD_AUTO, LP_METHOD_PRIMAL, LP_METHOD_DUAL, or LP_METHOD_BARRIER
         """
         # CPLEX lpmethod: 0=auto, 1=primal, 2=dual, 4=barrier
-        _map = {LP_METHOD_AUTO: 0, LP_METHOD_PRIMAL: 1,
-                LP_METHOD_DUAL: 2, LP_METHOD_BARRIER: 4}
+        _map = {LP_METHOD_AUTO: 0, LP_METHOD_PRIMAL: 1, LP_METHOD_DUAL: 2, LP_METHOD_BARRIER: 4}
         self.parameters.lpmethod.set(_map.get(method, 0))
 
     def get_lp_method(self):
         """Return the current LP method as a solver-neutral string."""
-        _rmap = {0: LP_METHOD_AUTO, 1: LP_METHOD_PRIMAL,
-                 2: LP_METHOD_DUAL, 4: LP_METHOD_BARRIER}
+        _rmap = {0: LP_METHOD_AUTO, 1: LP_METHOD_PRIMAL, 2: LP_METHOD_DUAL, 4: LP_METHOD_BARRIER}
         return _rmap.get(self.parameters.lpmethod.get(), LP_METHOD_AUTO)
 
     def get_basis(self):
@@ -362,9 +371,7 @@ def set_basis(self, basis):
         """
         if basis is None:
             return
-        self.start.set_start(
-            col_status=basis['vbasis'], row_status=basis['cbasis'],
-            col_primal=[], row_primal=[], col_dual=[], row_dual=[])
+        self.start.set_start(col_status=basis['vbasis'], row_status=basis['cbasis'], col_primal=[], row_primal=[], col_dual=[], row_dual=[])
 
     def set_time_limit(self, t):
         """Set the computation time limit (in seconds)"""
diff --git a/straindesign/glpk_interface.py b/straindesign/glpk_interface.py
index 90ebaab..80a3e1c 100644
--- a/straindesign/glpk_interface.py
+++ b/straindesign/glpk_interface.py
@@ -387,8 +387,7 @@ def set_lp_method(self, method):
         # GLPK meth: 1=primal, 2=dual, 3=dual+pricing
         if method == LP_METHOD_BARRIER:
             logging.warning('GLPK does not support barrier method, falling back to dual simplex.')
-        _map = {LP_METHOD_AUTO: 1, LP_METHOD_PRIMAL: 1,
-                LP_METHOD_DUAL: 2, LP_METHOD_BARRIER: 2}
+        _map = {LP_METHOD_AUTO: 1, LP_METHOD_PRIMAL: 1, LP_METHOD_DUAL: 2, LP_METHOD_BARRIER: 2}
         self.lp_params.meth = _map.get(method, 1)
 
     def get_lp_method(self):
diff --git a/straindesign/gurobi_interface.py b/straindesign/gurobi_interface.py
index a9aa730..3087d09 100644
--- a/straindesign/gurobi_interface.py
+++ b/straindesign/gurobi_interface.py
@@ -31,6 +31,7 @@
 # Shared quiet environment — avoids creating a new license session per model
 _quiet_env = None
 
+
 def _get_quiet_env():
     """Return a shared Gurobi environment with OutputFlag=0."""
     global _quiet_env
@@ -107,7 +108,18 @@ class Gurobi_MILP_LP(gp.Model):
             A Gurobi MILP/LP interface class.
     """
 
-    def __init__(self, c=None, A_ineq=None, b_ineq=None, A_eq=None, b_eq=None, lb=None, ub=None, vtype=None, indic_constr=None, seed=None, milp_threads=None):
+    def __init__(self,
+                 c=None,
+                 A_ineq=None,
+                 b_ineq=None,
+                 A_eq=None,
+                 b_eq=None,
+                 lb=None,
+                 ub=None,
+                 vtype=None,
+                 indic_constr=None,
+                 seed=None,
+                 milp_threads=None):
         super().__init__(env=_get_quiet_env())
         try:
             numvars = A_ineq.shape[1]
@@ -142,8 +154,8 @@ def __init__(self, c=None, A_ineq=None, b_ineq=None, A_eq=None, b_eq=None, lb=No
                 cols = A_csr.indices[start:end]
                 vals = A_csr.data[start:end]
                 lhs = gp.quicksum(float(val) * x[int(col)] for col, val in zip(cols, vals))
-                self.addGenConstrIndicator(x[indic_constr.binv[i]], bool(indic_constr.indicval[i]),
-                                           lhs, '=' if indic_constr.sense[i] == 'E' else '<', indic_constr.b[i])
+                self.addGenConstrIndicator(x[indic_constr.binv[i]], bool(indic_constr.indicval[i]), lhs,
+                                           '=' if indic_constr.sense[i] == 'E' else '<', indic_constr.b[i])
 
         # set parameters
         self.params.OutputFlag = 0
@@ -379,14 +391,12 @@ def set_lp_method(self, method):
         Args:
             method: LP_METHOD_AUTO, LP_METHOD_PRIMAL, LP_METHOD_DUAL, or LP_METHOD_BARRIER
         """
-        _map = {LP_METHOD_AUTO: -1, LP_METHOD_PRIMAL: 0,
-                LP_METHOD_DUAL: 1, LP_METHOD_BARRIER: 2}
+        _map = {LP_METHOD_AUTO: -1, LP_METHOD_PRIMAL: 0, LP_METHOD_DUAL: 1, LP_METHOD_BARRIER: 2}
         self.params.Method = _map.get(method, -1)
 
     def get_lp_method(self):
         """Return the current LP method as a solver-neutral string."""
-        _rmap = {-1: LP_METHOD_AUTO, 0: LP_METHOD_PRIMAL,
-                 1: LP_METHOD_DUAL, 2: LP_METHOD_BARRIER}
+        _rmap = {-1: LP_METHOD_AUTO, 0: LP_METHOD_PRIMAL, 1: LP_METHOD_DUAL, 2: LP_METHOD_BARRIER}
         return _rmap.get(self.params.Method, LP_METHOD_AUTO)
 
     def get_basis(self):
diff --git a/straindesign/lptools.py b/straindesign/lptools.py
index 1ba04fe..9da5494 100644
--- a/straindesign/lptools.py
+++ b/straindesign/lptools.py
@@ -319,7 +319,10 @@ def fva_legacy(model, **kwargs) -> DataFrame:
     if status not in [OPTIMAL, UNBOUNDED]:
         logging.info('FVA problem not feasible.')
         return DataFrame(
-            {"minimum": [nan] * numr, "maximum": [nan] * numr},
+            {
+                "minimum": [nan] * numr,
+                "maximum": [nan] * numr
+            },
             index=reaction_ids,
         )
 
@@ -328,18 +331,14 @@ def fva_legacy(model, **kwargs) -> DataFrame:
     x = [nan] * 2 * numr
 
     if processes > 1 and numr > 300:
-        with SDPool(processes, initializer=fva_worker_init,
-                    initargs=(A_ineq, b_ineq, A_eq, b_eq, lb, ub, solver)) as pool:
+        with SDPool(processes, initializer=fva_worker_init, initargs=(A_ineq, b_ineq, A_eq, b_eq, lb, ub, solver)) as pool:
             chunk_size = len(reaction_ids) // processes
-            for i, value in pool.imap_unordered(fva_worker_compute, range(2 * numr),
-                                                chunksize=chunk_size):
+            for i, value in pool.imap_unordered(fva_worker_compute, range(2 * numr), chunksize=chunk_size):
                 x[i] = value
     elif processes > 1 and numr > 500 and solver == GLPK:
-        with SDPool(processes, initializer=fva_worker_init_glpk,
-                    initargs=(A_ineq, b_ineq, A_eq, b_eq, lb, ub)) as pool:
+        with SDPool(processes, initializer=fva_worker_init_glpk, initargs=(A_ineq, b_ineq, A_eq, b_eq, lb, ub)) as pool:
             chunk_size = len(reaction_ids) // processes
-            for i, value in pool.imap_unordered(fva_worker_compute_glpk, range(2 * numr),
-                                                chunksize=chunk_size):
+            for i, value in pool.imap_unordered(fva_worker_compute_glpk, range(2 * numr), chunksize=chunk_size):
                 x[i] = value
     else:
         fva_worker_init(A_ineq, b_ineq, A_eq, b_eq, lb, ub, solver)
@@ -352,8 +351,7 @@ def fva_legacy(model, **kwargs) -> DataFrame:
         logging.warning(f'FVA: {len(nan_remaining)}/{2*numr} LP solves returned NaN, re-solving.')
         _BATCH = 50
         while nan_remaining:
-            lp_retry = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq,
-                               lb=lb, ub=ub, solver=solver)
+            lp_retry = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb, ub=ub, solver=solver)
             prev_retry = 0
             for i in nan_remaining[:_BATCH]:
                 C = idx2c(i, prev_retry)
@@ -374,8 +372,7 @@ def fva_legacy(model, **kwargs) -> DataFrame:
                 sig = sign(mod(i, 2) - 0.5)
                 c_vec = [0.0] * numr
                 c_vec[col] = sig
-                lp_last = MILP_LP(c=c_vec, A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq,
-                                  lb=lb, ub=ub, solver=solver)
+                lp_last = MILP_LP(c=c_vec, A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb, ub=ub, solver=solver)
                 x[i] = lp_last.slim_solve()
             nan_remaining = [i for i in nan_remaining if isnan(x[i])]
         if nan_remaining:
@@ -383,8 +380,10 @@ def fva_legacy(model, **kwargs) -> DataFrame:
 
     x = [v if abs(v) >= 1e-11 else 0.0 for v in x]
     return DataFrame(
-        {"minimum": [x[i] for i in range(1, 2 * numr, 2)],
-         "maximum": [-x[i] for i in range(0, 2 * numr, 2)]},
+        {
+            "minimum": [x[i] for i in range(1, 2 * numr, 2)],
+            "maximum": [-x[i] for i in range(0, 2 * numr, 2)]
+        },
         index=reaction_ids,
     )
 
@@ -623,8 +622,7 @@ def slim_fba_via_cmp(model, cmp_model, cmp_map, **kwargs) -> float:
         if isinstance(obj, str):
             obj = linexpr2dict(obj, orig_reaction_ids)
     else:
-        obj = {r.id: r.objective_coefficient for r in model.reactions
-               if r.objective_coefficient != 0}
+        obj = {r.id: r.objective_coefficient for r in model.reactions if r.objective_coefficient != 0}
 
     # Trace each objective reaction through compression
     obj_cmp = {}
@@ -673,8 +671,7 @@ def slim_fba_via_cmp(model, cmp_model, cmp_map, **kwargs) -> float:
     lb = [float(v.lower_bound) for v in cmp_model.reactions]
     ub = [float(v.upper_bound) for v in cmp_model.reactions]
 
-    fba_prob = MILP_LP(c=c, A_ineq=A_ineq, b_ineq=b_ineq,
-                       A_eq=A_eq, b_eq=b_eq, lb=lb, ub=ub, solver=solver)
+    fba_prob = MILP_LP(c=c, A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb, ub=ub, solver=solver)
     opt_cx = fba_prob.slim_solve()
 
     if isnan(opt_cx):
@@ -906,6 +903,7 @@ def yopt(model, **kwargs) -> Solution:
     else:
         status = INFEASIBLE
 
+
 def expand_fluxes(fluxes_cmp, cmp_map, orig_reaction_ids):
     """Expand a compressed flux vector to the full (uncompressed) model.
 
@@ -944,6 +942,7 @@ def expand_fluxes(fluxes_cmp, cmp_map, orig_reaction_ids):
             fluxes[rid] = 0.0
     return fluxes
 
+
 def _make_fix_constraint(axes, ax_idx, ax_type, value):
     """Create an equality constraint fixing axis ax_idx to value."""
     if ax_type == 'rate':
@@ -960,8 +959,7 @@ def _optimize_axis(model, ax_idx, axes, ax_type, constraints, solver, sense):
     if ax_type == 'rate':
         sol = fba(model, obj=axes[ax_idx][0], constraints=constraints, solver=solver, obj_sense=sense)
     else:  # yield
-        sol = yopt(model, obj_num=axes[ax_idx][0], obj_den=axes[ax_idx][1],
-                   constraints=constraints, solver=solver, obj_sense=sense)
+        sol = yopt(model, obj_num=axes[ax_idx][0], obj_den=axes[ax_idx][1], constraints=constraints, solver=solver, obj_sense=sense)
     return sol
 
 
@@ -1000,8 +998,7 @@ def _fba_project(obj_dict):
 
     # Step 1: Find 4 extremes
     extremes = []
-    for coeff, sense in [(ax0_coeff, 'maximize'), (ax0_coeff, 'minimize'),
-                         (ax1_coeff, 'maximize'), (ax1_coeff, 'minimize')]:
+    for coeff, sense in [(ax0_coeff, 'maximize'), (ax0_coeff, 'minimize'), (ax1_coeff, 'maximize'), (ax1_coeff, 'minimize')]:
         sol = fba(model, obj=coeff, constraints=constraints, solver=solver, obj_sense=sense)
         if sol.status == OPTIMAL:
             x0 = sum(c * sol.fluxes.get(r, 0) for r, c in ax0_coeff.items())
@@ -1027,16 +1024,13 @@ def _fba_project(obj_dict):
     unique_pts.sort(key=lambda p: arctan2(p[1] - cy, p[0] - cx))
 
     # Step 4: Recursive edge refinement
-    diameter = max(
-        ((a[0]-b[0])**2 + (a[1]-b[1])**2)**0.5
-        for a in unique_pts for b in unique_pts
-    )
+    diameter = max(((a[0] - b[0])**2 + (a[1] - b[1])**2)**0.5 for a in unique_pts for b in unique_pts)
     min_edge = 1e-10 * diameter if diameter > 0 else 1e-10
 
     def _refine(vi, vj, depth):
         if depth > 50:
             return [vi]
-        edge_len = ((vi[0]-vj[0])**2 + (vi[1]-vj[1])**2)**0.5
+        edge_len = ((vi[0] - vj[0])**2 + (vi[1] - vj[1])**2)**0.5
         if edge_len < min_edge:
             return [vi]
         # Outward normal (perpendicular to edge, pointing outward from centroid)
@@ -1063,8 +1057,8 @@ def _refine(vi, vj, depth):
             dist /= norm_len
         if dist > tol and edge_len > min_edge:
             # Check p_new is not a duplicate of vi or vj
-            if ((abs(p_new[0]-vi[0]) < tol and abs(p_new[1]-vi[1]) < tol) or
-                (abs(p_new[0]-vj[0]) < tol and abs(p_new[1]-vj[1]) < tol)):
+            if ((abs(p_new[0] - vi[0]) < tol and abs(p_new[1] - vi[1]) < tol) or
+                (abs(p_new[0] - vj[0]) < tol and abs(p_new[1] - vj[1]) < tol)):
                 return [vi]
             left = _refine(vi, p_new, depth + 1)
             right = _refine(p_new, vj, depth + 1)
@@ -1091,6 +1085,7 @@ def _trace_boundary_adaptive(model, axes, ax_types, constraints, solver, max_dep
     Returns (upper_boundary, lower_boundary) as sorted lists of (x, y) tuples.
     Uses recursive midpoint refinement where linear interpolation error exceeds tolerance.
     """
+
     def _fix_and_opt(x_val, sense):
         constr = constraints.copy()
         constr.append(_make_fix_constraint(axes, 0, ax_types[0], x_val))
@@ -1100,11 +1095,8 @@ def _fix_and_opt(x_val, sense):
         return nan
 
     # Step 1: axis-0 range endpoints (already known from val_limits, but we need y values)
-    x_min, x_max = ceil_dec(
-        _optimize_axis(model, 0, axes, ax_types[0], constraints, solver, 'minimize').objective_value, 8
-    ), floor_dec(
-        _optimize_axis(model, 0, axes, ax_types[0], constraints, solver, 'maximize').objective_value, 8
-    )
+    x_min, x_max = ceil_dec(_optimize_axis(model, 0, axes, ax_types[0], constraints, solver, 'minimize').objective_value,
+                            8), floor_dec(_optimize_axis(model, 0, axes, ax_types[0], constraints, solver, 'maximize').objective_value, 8)
 
     # Step 2: y values at endpoints
     y_min_at_xmin = _fix_and_opt(x_min, 'minimize')
@@ -1168,10 +1160,7 @@ def _trace_polytope_3d_rate(model, axes, constraints, solver):
     tol = 1e-8
 
     def _project(sol):
-        return tuple(
-            ceil_dec(sum(c * sol.fluxes.get(r, 0) for r, c in coeffs[j].items()), 9)
-            for j in range(3)
-        )
+        return tuple(ceil_dec(sum(c * sol.fluxes.get(r, 0) for r, c in coeffs[j].items()), 9) for j in range(3))
 
     def _is_dup(pt, pts):
         return any(all(abs(pt[k] - q[k]) < tol for k in range(3)) for q in pts)
@@ -1248,10 +1237,8 @@ def _hull_face_polygons(hull):
         # Build orthonormal basis in the face plane
         ref = array([1, 0, 0]) if abs(normal[0]) < 0.9 else array([0, 1, 0])
         u = ref - normal * normal.dot(ref)
-        u = u / (u.dot(u) ** 0.5)
-        v = array([normal[1]*u[2] - normal[2]*u[1],
-                    normal[2]*u[0] - normal[0]*u[2],
-                    normal[0]*u[1] - normal[1]*u[0]])
+        u = u / (u.dot(u)**0.5)
+        v = array([normal[1] * u[2] - normal[2] * u[1], normal[2] * u[0] - normal[0] * u[2], normal[0] * u[1] - normal[1] * u[0]])
         # Sort by angle in face plane
         angles = [arctan2((pt - centroid).dot(v), (pt - centroid).dot(u)) for pt in pts]
         order = sorted(range(len(verts)), key=lambda i: angles[i])
@@ -1420,14 +1407,13 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]:
             variable contains information about which datapoints need to be connected in triangles to
             render a closed surface. The last variable contains the matplotlib object.
     """
-    
+
     cmp_model = kwargs.pop('cmp_model', None)
     cmp_map = kwargs.pop('cmp_map', None)
 
     _orig_axes = None  # store original axis names for labelling
     if cmp_model is not None and cmp_map is not None:
-        from straindesign.networktools import (resolve_gene_constraints,
-            compress_constraints, _build_cmp_reverse_map)
+        from straindesign.networktools import (resolve_gene_constraints, compress_constraints, _build_cmp_reverse_map)
         # Resolve gene constraints on the original model, then compress
         if CONSTRAINTS in kwargs and kwargs[CONSTRAINTS]:
             kwargs[CONSTRAINTS] = resolve_gene_constraints(model, kwargs[CONSTRAINTS])
@@ -1455,7 +1441,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]:
                 axes[i] = [reverse.get(a, a) if isinstance(a, str) else a for a in ax]
         # Switch to compressed model
         model = cmp_model
-      
+
     reaction_ids = model.reactions.list_attr("id")
 
     if CONSTRAINTS in kwargs:
@@ -1592,8 +1578,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]:
             vertices = _trace_polygon_rate_rate(model, axes, kwargs[CONSTRAINTS], solver)
         else:
             adapt_depth = max(5, int(log2(max(points, 2))))
-            upper, lower = _trace_boundary_adaptive(
-                model, axes, ax_type, kwargs[CONSTRAINTS], solver, max_depth=adapt_depth)
+            upper, lower = _trace_boundary_adaptive(model, axes, ax_type, kwargs[CONSTRAINTS], solver, max_depth=adapt_depth)
             # Build polygon from upper (left-to-right) + reversed lower (right-to-left)
             if upper and lower:
                 vertices = upper + list(reversed(lower))
@@ -1634,9 +1619,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]:
             x0, y0 = unique_verts[0]
             dx, dy = unique_verts[1][0] - x0, unique_verts[1][1] - y0
             span = max(abs(dx), abs(dy), 1e-12)
-            is_collinear = all(
-                abs((v[0] - x0) * dy - (v[1] - y0) * dx) / span < 1e-6
-                for v in unique_verts[2:])
+            is_collinear = all(abs((v[0] - x0) * dy - (v[1] - y0) * dx) / span < 1e-6 for v in unique_verts[2:])
 
         if is_collinear and len(unique_verts) <= 1:
             # Collapsed to a point
@@ -1786,8 +1769,8 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]:
 
             elif n_yields == 1:
                 # 1 yield + 2 rate: slice along yield, trace rate-rate polygon per slice
-                datapoints, triang, slice_outlines = _trace_3d_slice_polygon(
-                    model, axes, ax_type, val_limits, kwargs[CONSTRAINTS], solver, points)
+                datapoints, triang, slice_outlines = _trace_3d_slice_polygon(model, axes, ax_type, val_limits, kwargs[CONSTRAINTS], solver,
+                                                                             points)
 
             else:
                 # 2+ yields: grid-based scanning (fallback)
@@ -1807,8 +1790,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]:
                     if abs(y_hi - y_lo) < 1e-10:
                         y_space = [y_lo]
                     else:
-                        n_pts = max(3, int(points * abs(y_hi - y_lo) / max(1e-10,
-                            max(abs(val_limits[1][1] - val_limits[1][0]), 1e-10))))
+                        n_pts = max(3, int(points * abs(y_hi - y_lo) / max(1e-10, max(abs(val_limits[1][1] - val_limits[1][0]), 1e-10))))
                         n_pts = min(n_pts, points)
                         y_space = linspace(y_lo, y_hi, n_pts).tolist()
                     upper_slice = []
@@ -1843,8 +1825,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]:
                 triang = []
                 for s in range(len(slices) - 1):
                     _triangulate_strips(datapoints_top[s], datapoints_top[s + 1], datapoints, triang)
-                    _triangulate_strips(datapoints_bottom[s], datapoints_bottom[s + 1], datapoints, triang,
-                                        flip_winding=True)
+                    _triangulate_strips(datapoints_bottom[s], datapoints_bottom[s + 1], datapoints, triang, flip_winding=True)
                 front_top = [t[0] for t in datapoints_top]
                 front_bot = [b[0] for b in datapoints_bottom]
                 _triangulate_strips(front_top, front_bot, datapoints, triang, flip_winding=True)
@@ -1852,8 +1833,7 @@ def plot_flux_space(model, axes, **kwargs) -> Tuple[list, list, list]:
                 back_bot = [b[-1] for b in datapoints_bottom]
                 _triangulate_strips(back_top, back_bot, datapoints, triang)
                 _triangulate_strips(datapoints_top[0], datapoints_bottom[0], datapoints, triang)
-                _triangulate_strips(datapoints_top[-1], datapoints_bottom[-1], datapoints, triang,
-                                    flip_winding=True)
+                _triangulate_strips(datapoints_top[-1], datapoints_bottom[-1], datapoints, triang, flip_winding=True)
 
         if not datapoints:
             raise Exception('No feasible points found. Problem may be infeasible.')
@@ -1882,10 +1862,8 @@ def _normal_color(face_pts):
                 return 0.5
             e1 = pts[1] - pts[0]
             e2 = pts[2] - pts[0]
-            normal = array([e1[1]*e2[2] - e1[2]*e2[1],
-                            e1[2]*e2[0] - e1[0]*e2[2],
-                            e1[0]*e2[1] - e1[1]*e2[0]])
-            length = (normal.dot(normal)) ** 0.5
+            normal = array([e1[1] * e2[2] - e1[2] * e2[1], e1[2] * e2[0] - e1[0] * e2[2], e1[0] * e2[1] - e1[1] * e2[0]])
+            length = (normal.dot(normal))**0.5
             if length > 0:
                 normal = normal / length
             # Map normal direction to scalar: use spherical angles
@@ -1899,8 +1877,7 @@ def _normal_color(face_pts):
             rng = mx - mn if mx > mn else 1
             face_colors = plt.get_cmap(cmap)([(c - mn) / rng for c in color_vals])
             face_colors[:, 3] = 1.0
-            collection = Poly3DCollection(poly_verts, facecolors=face_colors,
-                                          edgecolors='black', linewidths=1.0)
+            collection = Poly3DCollection(poly_verts, facecolors=face_colors, edgecolors='black', linewidths=1.0)
             ax3.add_collection3d(collection)
             plot1 = collection
         elif triang:
@@ -1911,16 +1888,14 @@ def _normal_color(face_pts):
             rng = mx - mn if mx > mn else 1
             face_colors = plt.get_cmap(cmap)([(c - mn) / rng for c in color_vals])
             face_colors[:, 3] = 1.0
-            collection = Poly3DCollection(tri_verts, facecolors=face_colors,
-                                          edgecolors='none', linewidths=0)
+            collection = Poly3DCollection(tri_verts, facecolors=face_colors, edgecolors='none', linewidths=0)
             ax3.add_collection3d(collection)
             # Draw slice polygon outlines (exact contour at each yield level)
             if slice_outlines:
                 for idx_list in slice_outlines:
                     pts = array([datapoints[i] for i in idx_list])
                     pts = array(list(pts) + [pts[0]])  # close the loop
-                    ax3.plot(pts[:, 0], pts[:, 1], pts[:, 2],
-                             color='gray', linewidth=0.4)
+                    ax3.plot(pts[:, 0], pts[:, 1], pts[:, 2], color='gray', linewidth=0.4)
             plot1 = collection
         else:
             plot1 = ax3.scatter(x, y, z, s=20)
diff --git a/straindesign/networktools.py b/straindesign/networktools.py
index 060dc88..0c5027e 100644
--- a/straindesign/networktools.py
+++ b/straindesign/networktools.py
@@ -53,6 +53,7 @@ def _sb_noop(self, lb, ub):
 
 # -- Cobra-level suppression replacements ------------------------------------
 
+
 def _suppressed_set_id(self, value):
     """Bypass solver variable rename: write _id and update DictList index."""
     old_id = self._id
@@ -158,7 +159,7 @@ def _suppressed_copy(model):
     extension's ``add_metabolites`` calls run without building (or pushing constraints into) a live
     solver.
     """
-    iface = model.solver.interface          # captured before stubbing
+    iface = model.solver.interface  # captured before stubbing
     saved = model._solver
     orig_copy = next(o for cls, attr, o in _ORIG_COBRA if cls is Model and attr == 'copy')
     try:
@@ -169,6 +170,7 @@ def _suppressed_copy(model):
     new._solver = _CarrierSolver(iface)
     return new
 
+
 _ORIG_CONTAINER_GETITEM = None  # saved Container.__getitem__
 
 
@@ -272,10 +274,10 @@ def _suppressed_add_metabolites(self, metabolite_list):
 
 # -- Saved originals (None = not suppressed) ----------------------------------
 
-_ORIG_SLC = None   # (cls, method) for Constraint.set_linear_coefficients
-_ORIG_SB = None    # (cls, method) for Variable.set_bounds
+_ORIG_SLC = None  # (cls, method) for Constraint.set_linear_coefficients
+_ORIG_SB = None  # (cls, method) for Variable.set_bounds
 _ORIG_OSLC = None  # (cls, method) for Objective.set_linear_coefficients
-_ORIG_COBRA = []   # list of (cls, attr_name, original) for cobra-level patches
+_ORIG_COBRA = []  # list of (cls, attr_name, original) for cobra-level patches
 
 
 def _suppress_lp_updates(model):
@@ -438,8 +440,7 @@ def suppress_lp_context(model):
                 del model._suppressed_obj
             if current_ids != _pre_ids:
                 if model.groups:
-                    kept = {c.id for c in
-                            list(model.reactions) + list(model.metabolites) + list(model.genes)}
+                    kept = {c.id for c in list(model.reactions) + list(model.metabolites) + list(model.genes)}
                     for grp in model.groups:
                         stale = [m for m in grp.members if m.id not in kept]
                         if stale:
@@ -461,11 +462,13 @@ def with_suppressed_lp(func):
 
     The first positional argument must be a cobra Model.
     """
+
     @wraps(func)
     def wrapper(*args, **kwargs):
         model = args[0]
         with suppress_lp_context(model):
             return func(*args, **kwargs)
+
     return wrapper
 
 
@@ -473,6 +476,7 @@ def wrapper(*args, **kwargs):
 # I/O Suppression
 # =============================================================================
 
+
 @contextmanager
 def _silent_io():
     """Suppress stdout, stderr and logging."""
@@ -608,6 +612,7 @@ def gene_kos_to_constraints(model, gene_kos):
 
     return [[{r: 1}, '=', 0] for r in sorted(knocked_out_reactions)]
 
+
 def _build_cmp_reverse_map(cmp_map):
     """Build reverse lookup: original reaction ID -> final compressed ID.
 
@@ -668,11 +673,10 @@ def compress_constraints(constraints, cmp_mapReac):
                 coeff_dict = c[0]
                 lumped = [k for k in coeff_dict if k in old_reac_val]
                 if lumped:
-                    coeff_dict[new_reac] = sum(
-                        coeff_dict.pop(k) * old_reac_val[k] for k in lumped
-                    )
+                    coeff_dict[new_reac] = sum(coeff_dict.pop(k) * old_reac_val[k] for k in lumped)
     return constraints
 
+
 def resolve_gene_constraints(model, constraints):
     """Scan constraints for gene IDs/names and replace with reaction constraints.
 
@@ -1358,8 +1362,7 @@ def filter_sd_maxcost(sd, max_cost, kocost, kicost):
     # non-made KIs are marked by 0.0 and non-made KOs don't appear.
     # We count costs of interventions made, which are marked by v != 0.
     if max_cost:
-        costs = [np.sum([(kocost[k] if k in kocost else kicost.get(k, 0)) if v != 0 else 0
-                         for k, v in m.items()]) for m in sd]
+        costs = [np.sum([(kocost[k] if k in kocost else kicost.get(k, 0)) if v != 0 else 0 for k, v in m.items()]) for m in sd]
         sd = [sd[i] for i in range(len(sd)) if costs[i] <= max_cost + 1e-8]
         # sort strain designs by intervention costs
         [s.update({'**cost**': c}) for s, c in zip(sd, costs)]
@@ -1425,14 +1428,12 @@ def bound_blocked_or_irrevers_fva(model, **kwargs):
 from operator import attrgetter as _attrgetter
 from fractions import Fraction as _Fraction
 import math as _math
-from cobra.io.dict import (_metabolite_to_dict, _metabolite_from_dict,
-                           _gene_to_dict, gene_from_dict, _fix_type, _update_optional,
-                           _OPTIONAL_REACTION_ATTRIBUTES, _ORDERED_OPTIONAL_REACTION_KEYS,
-                           _OPTIONAL_MODEL_ATTRIBUTES, _ORDERED_OPTIONAL_MODEL_KEYS)
+from cobra.io.dict import (_metabolite_to_dict, _metabolite_from_dict, _gene_to_dict, gene_from_dict, _fix_type, _update_optional,
+                           _OPTIONAL_REACTION_ATTRIBUTES, _ORDERED_OPTIONAL_REACTION_KEYS, _OPTIONAL_MODEL_ATTRIBUTES,
+                           _ORDERED_OPTIONAL_MODEL_KEYS)
 from cobra.util.solver import set_objective as _set_objective
 
-_INF_TOKENS = {'inf': _math.inf, 'infinity': _math.inf, '+inf': _math.inf,
-               '-inf': -_math.inf, '-infinity': -_math.inf}
+_INF_TOKENS = {'inf': _math.inf, 'infinity': _math.inf, '+inf': _math.inf, '-inf': -_math.inf, '-infinity': -_math.inf}
 
 
 def _num_to_json(v):
@@ -1442,13 +1443,13 @@ def _num_to_json(v):
     if isinstance(v, bool):
         return v
     if isinstance(v, _Fraction):
-        return str(v)                                    # exact: '1/3', '10', '-1/2'
+        return str(v)  # exact: '1/3', '10', '-1/2'
     num, den = getattr(v, 'numerator', None), getattr(v, 'denominator', None)
     if num is not None and den is not None and not isinstance(v, int):
-        return '%d/%d' % (int(num), int(den))            # sympy Rational etc.
-    v = _fix_type(v)                                     # numpy -> native python
+        return '%d/%d' % (int(num), int(den))  # sympy Rational etc.
+    v = _fix_type(v)  # numpy -> native python
     if isinstance(v, float) and (_math.isinf(v) or _math.isnan(v)):
-        return str(v)                                    # 'inf', '-inf', 'nan'
+        return str(v)  # 'inf', '-inf', 'nan'
     return v
 
 
@@ -1461,7 +1462,7 @@ def _num_from_json(v):
             return _INF_TOKENS[low]
         if low == 'nan':
             return _math.nan
-        return _Fraction(v)                              # 'num/den' or integer string
+        return _Fraction(v)  # 'num/den' or integer string
     return v
 
 
@@ -1473,8 +1474,7 @@ def _reaction_to_json(reaction, obj_coeff=0):
     new['id'] = reaction.id
     new['name'] = reaction.name
     new['metabolites'] = _OrderedDict(
-        (str(met), _num_to_json(reaction.metabolites[met]))
-        for met in sorted(reaction.metabolites, key=_attrgetter('id')))
+        (str(met), _num_to_json(reaction.metabolites[met])) for met in sorted(reaction.metabolites, key=_attrgetter('id')))
     new['lower_bound'] = _num_to_json(reaction.lower_bound)
     new['upper_bound'] = _num_to_json(reaction.upper_bound)
     new['gene_reaction_rule'] = reaction.gene_reaction_rule
@@ -1495,9 +1495,8 @@ def _reaction_from_json(reaction, model):
         if k in {'objective_coefficient', 'reversibility', 'reaction'}:
             continue
         elif k == 'metabolites':
-            new_reaction.add_metabolites(_OrderedDict(
-                (model.metabolites.get_by_id(str(met)), _num_from_json(coeff))
-                for met, coeff in v.items()))
+            new_reaction.add_metabolites(
+                _OrderedDict((model.metabolites.get_by_id(str(met)), _num_from_json(coeff)) for met, coeff in v.items()))
         elif k in {'lower_bound', 'upper_bound'}:
             setattr(new_reaction, k, _num_from_json(v))
         else:
@@ -1552,10 +1551,8 @@ def model_from_dict(obj):
     model.add_metabolites([_metabolite_from_dict(m) for m in obj['metabolites']])
     model.genes.extend([gene_from_dict(g) for g in obj['genes']])
     model.add_reactions([_reaction_from_json(r, model) for r in obj['reactions']])
-    objective_reactions = [r for r in obj['reactions']
-                           if _num_from_json(r.get('objective_coefficient', 0)) != 0]
-    coefficients = {model.reactions.get_by_id(r['id']): _num_from_json(r['objective_coefficient'])
-                    for r in objective_reactions}
+    objective_reactions = [r for r in obj['reactions'] if _num_from_json(r.get('objective_coefficient', 0)) != 0]
+    coefficients = {model.reactions.get_by_id(r['id']): _num_from_json(r['objective_coefficient']) for r in objective_reactions}
     _set_objective(model, coefficients)
     for k, v in obj.items():
         if k in {'id', 'name', 'notes', 'compartments', 'annotation'}:
diff --git a/straindesign/parse_constr.py b/straindesign/parse_constr.py
index ce628e2..81843f7 100644
--- a/straindesign/parse_constr.py
+++ b/straindesign/parse_constr.py
@@ -327,4 +327,3 @@ def linexprdict2str(D):
         return expr
     else:
         return ""
-
diff --git a/straindesign/scip_interface.py b/straindesign/scip_interface.py
index 8798918..7a7fc8c 100644
--- a/straindesign/scip_interface.py
+++ b/straindesign/scip_interface.py
@@ -94,7 +94,18 @@ class SCIP_MILP(pso.Model):
             A SCIP MILP interface class.
     """
 
-    def __init__(self, c=None, A_ineq=None, b_ineq=None, A_eq=None, b_eq=None, lb=None, ub=None, vtype=None, indic_constr=None, seed=None, milp_threads=None):
+    def __init__(self,
+                 c=None,
+                 A_ineq=None,
+                 b_ineq=None,
+                 A_eq=None,
+                 b_eq=None,
+                 lb=None,
+                 ub=None,
+                 vtype=None,
+                 indic_constr=None,
+                 seed=None,
+                 milp_threads=None):
         super().__init__()
         # uncomment to forward SCIP output to python terminal
         # self.redirectOutput()
@@ -371,16 +382,14 @@ def set_lp_method(self, method):
             method: LP_METHOD_AUTO, LP_METHOD_PRIMAL, LP_METHOD_DUAL, or LP_METHOD_BARRIER
         """
         # SCIP lp/initalgorithm + lp/resolvealgorithm: 's'=auto, 'p'=primal, 'd'=dual, 'b'=barrier
-        _map = {LP_METHOD_AUTO: 's', LP_METHOD_PRIMAL: 'p',
-                LP_METHOD_DUAL: 'd', LP_METHOD_BARRIER: 'b'}
+        _map = {LP_METHOD_AUTO: 's', LP_METHOD_PRIMAL: 'p', LP_METHOD_DUAL: 'd', LP_METHOD_BARRIER: 'b'}
         algo = _map.get(method, 's')
         self.setParam('lp/initalgorithm', algo)
         self.setParam('lp/resolvealgorithm', algo)
 
     def get_lp_method(self):
         """Return the current LP method as a solver-neutral string."""
-        _rmap = {'s': LP_METHOD_AUTO, 'p': LP_METHOD_PRIMAL,
-                 'd': LP_METHOD_DUAL, 'b': LP_METHOD_BARRIER}
+        _rmap = {'s': LP_METHOD_AUTO, 'p': LP_METHOD_PRIMAL, 'd': LP_METHOD_DUAL, 'b': LP_METHOD_BARRIER}
         return _rmap.get(self.getParam('lp/initalgorithm'), LP_METHOD_AUTO)
 
     def get_basis(self):
diff --git a/straindesign/solver_interface.py b/straindesign/solver_interface.py
index b89a074..0462cec 100644
--- a/straindesign/solver_interface.py
+++ b/straindesign/solver_interface.py
@@ -102,7 +102,8 @@ class MILP_LP(object):
 
     def __init__(self, **kwargs):
         allowed_keys = {
-            'c', 'A_ineq', 'b_ineq', 'A_eq', 'b_eq', 'lb', 'ub', 'vtype', 'indic_constr', 'M', SOLVER, 'skip_checks', 'tlim', SEED, MILP_THREADS
+            'c', 'A_ineq', 'b_ineq', 'A_eq', 'b_eq', 'lb', 'ub', 'vtype', 'indic_constr', 'M', SOLVER, 'skip_checks', 'tlim', SEED,
+            MILP_THREADS
         }
         # set all keys passed in kwargs
         for key, value in kwargs.items():
diff --git a/straindesign/speedy_fva.py b/straindesign/speedy_fva.py
index 50be8f0..b4190f2 100644
--- a/straindesign/speedy_fva.py
+++ b/straindesign/speedy_fva.py
@@ -47,8 +47,12 @@
 from cobra import Configuration
 
 from straindesign.lptools import (
-    select_solver, idx2c, fva_worker_init, fva_worker_compute,
-    fva_worker_init_glpk, fva_worker_compute_glpk,
+    select_solver,
+    idx2c,
+    fva_worker_init,
+    fva_worker_compute,
+    fva_worker_init_glpk,
+    fva_worker_compute_glpk,
 )
 from straindesign.solver_interface import MILP_LP
 from straindesign.pool import SDPool
@@ -56,21 +60,24 @@
 from straindesign.names import CONSTRAINTS, SOLVER, OPTIMAL, UNBOUNDED, INFEASIBLE, GLPK, LP_METHOD_DUAL
 from straindesign.networktools import suppress_lp_context
 from straindesign.compression import (
-    compress_cobra_model, CompressionMethod, remove_conservation_relations,
-    stoichmat_coeff_to_fraction, stoichmat_coeff2float, remove_blocked_reactions,
+    compress_cobra_model,
+    CompressionMethod,
+    remove_conservation_relations,
+    stoichmat_coeff_to_fraction,
+    stoichmat_coeff2float,
+    remove_blocked_reactions,
 )
 
-
 # Tunable threshold at which parallel FVA kicks in. Is compared against the number of
 # LPs that are left to be solved by Phase 2. Empirically, parallel is only worthwhile
 # at high LP counts.
 _PARALLEL_PHASE2_MIN = 4000
 
-
 # ---------------------------------------------------------------------------
 # Compression helpers
 # ---------------------------------------------------------------------------
 
+
 def _compress_for_fva(model):
     """Copy and compress model for FVA (single-pass coupled compression + conservation removal).
 
@@ -95,9 +102,7 @@ def _compress_for_fva(model):
         # Single-pass coupled compression (NULLSPACE only, no RECURSIVE iteration)
         for r in cmp_model.reactions:
             r.gene_reaction_rule = ''
-        result = compress_cobra_model(cmp_model,
-                                      methods=[CompressionMethod.NULLSPACE],
-                                      in_place=True)
+        result = compress_cobra_model(cmp_model, methods=[CompressionMethod.NULLSPACE], in_place=True)
         rmap = result.reaction_map
         if len(rmap) < n_before:
             cmp_maps.append(rmap)
@@ -120,9 +125,7 @@ def _map_constraints(parsed_constraints, cmp_maps, cmp_reaction_ids):
                 coeff_dict = constraint[0]
                 lumped = [k for k in coeff_dict if k in orig_map]
                 if lumped:
-                    coeff_dict[cmp_id] = sum(
-                        coeff_dict.pop(k) * float(orig_map[k]) for k in lumped
-                    )
+                    coeff_dict[cmp_id] = sum(coeff_dict.pop(k) * float(orig_map[k]) for k in lumped)
     # Remove references to reactions not in compressed model (e.g., blocked)
     for constraint in parsed_constraints:
         for k in list(constraint[0].keys()):
@@ -166,8 +169,10 @@ def _expand_fva(fva_cmp, cmp_maps, orig_reaction_ids):
             result_max[rxn_id] = 0.0
 
     return DataFrame(
-        {"minimum": [result_min[r] for r in orig_reaction_ids],
-         "maximum": [result_max[r] for r in orig_reaction_ids]},
+        {
+            "minimum": [result_min[r] for r in orig_reaction_ids],
+            "maximum": [result_max[r] for r in orig_reaction_ids]
+        },
         index=orig_reaction_ids,
     )
 
@@ -176,6 +181,7 @@ def _expand_fva(fva_cmp, cmp_maps, orig_reaction_ids):
 # Global scan LP helper
 # ---------------------------------------------------------------------------
 
+
 def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0):
     """Build LP for min sum(|x|) via variable splitting.
 
@@ -193,9 +199,9 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0):
     tol = 1e-9
 
     # Classify reactions
-    fwd = lb >= -tol          # forward-only or fixed
-    bwd = ub <= tol           # backward-only or fixed
-    rev = (~fwd) & (~bwd)     # truly reversible
+    fwd = lb >= -tol  # forward-only or fixed
+    bwd = ub <= tol  # backward-only or fixed
+    rev = (~fwd) & (~bwd)  # truly reversible
     n_rev = int(rev.sum())
     rev_idx = np.where(rev)[0]
 
@@ -210,9 +216,9 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0):
         elif bwd[j]:
             c[j] = -1.0  # |x_j| = -x_j
         else:
-            c[j] = 1.0   # |x_j| = x_j (fwd or fixed)
+            c[j] = 1.0  # |x_j| = x_j (fwd or fixed)
     for k, j in enumerate(rev_idx):
-        c[n + k] = 1.0       # p_k
+        c[n + k] = 1.0  # p_k
         c[n + n_rev + k] = 1.0  # n_k
 
     # Equality constraints: original S*x = 0 (+ extras) + splitting equalities
@@ -233,8 +239,7 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0):
         rows_all = np.concatenate([rows_split, rows_split, rows_split])
         cols_all = np.concatenate([cols_x, cols_p, cols_n])
         data_all = np.concatenate([data_x, data_p, data_n])
-        A_split = sparse.csr_matrix((data_all, (rows_all, cols_all)),
-                                    shape=(n_rev, n_ext))
+        A_split = sparse.csr_matrix((data_all, (rows_all, cols_all)), shape=(n_rev, n_ext))
 
         # Extend original equalities to n_ext columns
         S_ext = sparse.hstack([S_eq, sparse.csr_matrix((S_eq.shape[0], 2 * n_rev))])
@@ -270,9 +275,14 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0):
         lb_ext[n + n_rev + k] = 0.0
         ub_ext[n + n_rev + k] = min(-lb[j], BIG)
 
-    lp = MILP_LP(c=c.tolist(), A_ineq=A_ineq_ext, b_ineq=list(b_ineq),
-                 A_eq=A_eq_full, b_eq=b_eq_full,
-                 lb=lb_ext.tolist(), ub=ub_ext.tolist(), solver=solver)
+    lp = MILP_LP(c=c.tolist(),
+                 A_ineq=A_ineq_ext,
+                 b_ineq=list(b_ineq),
+                 A_eq=A_eq_full,
+                 b_eq=b_eq_full,
+                 lb=lb_ext.tolist(),
+                 ub=ub_ext.tolist(),
+                 solver=solver)
     return lp, n
 
 
@@ -280,6 +290,7 @@ def _build_abssum_lp(S_eq, b_eq, A_ineq, b_ineq, lb, ub, solver, BIG=1000.0):
 # Main entry point
 # ---------------------------------------------------------------------------
 
+
 def speedy_fva(model, **kwargs):
     """Accelerated FVA using global scan LPs and KKT-based optimality propagation.
 
@@ -345,10 +356,8 @@ def speedy_fva(model, **kwargs):
         kwargs[CONSTRAINTS] = resolve_gene_constraints(orig_model, kwargs[CONSTRAINTS])
         kwargs[CONSTRAINTS] = parse_constraints(kwargs[CONSTRAINTS], orig_reaction_ids)
         if cmp_maps:
-            kwargs[CONSTRAINTS] = _map_constraints(
-                kwargs[CONSTRAINTS], cmp_maps, reaction_ids)
-        A_ineq_extra, b_ineq_extra, A_eq_extra, b_eq_extra = lineqlist2mat(
-            kwargs[CONSTRAINTS], reaction_ids)
+            kwargs[CONSTRAINTS] = _map_constraints(kwargs[CONSTRAINTS], cmp_maps, reaction_ids)
+        A_ineq_extra, b_ineq_extra, A_eq_extra, b_eq_extra = lineqlist2mat(kwargs[CONSTRAINTS], reaction_ids)
 
     if SOLVER not in kwargs:
         kwargs[SOLVER] = None
@@ -376,14 +385,16 @@ def speedy_fva(model, **kwargs):
     lb = np.array([v.lower_bound for v in model.reactions], dtype=np.float64)
     ub = np.array([v.upper_bound for v in model.reactions], dtype=np.float64)
 
-    lp = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq,
-                 lb=lb.tolist(), ub=ub.tolist(), solver=solver)
+    lp = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb.tolist(), ub=ub.tolist(), solver=solver)
 
     _, _, status = lp.solve()
     if status not in [OPTIMAL, UNBOUNDED]:
         logging.info('FVA problem not feasible.')
         return DataFrame(
-            {"minimum": [nan] * n_orig, "maximum": [nan] * n_orig},
+            {
+                "minimum": [nan] * n_orig,
+                "maximum": [nan] * n_orig
+            },
             index=reaction_ids,
         )
 
@@ -438,8 +449,7 @@ def _bound_scan(x_local):
     # ------------------------------------------------------------------
     _tp1 = _time.perf_counter()
     # 1a: v=0 feasibility check — free resolutions, no LP needed
-    v0_feasible = (not np.any(lb > tol_bound) and not np.any(ub < -tol_bound)
-                   and not has_constraints)
+    v0_feasible = (not np.any(lb > tol_bound) and not np.any(ub < -tol_bound) and not has_constraints)
     if v0_feasible:
         zero_lb = np.abs(lb) < tol_bound
         newly_min = zero_lb & (~res_min)
@@ -456,8 +466,7 @@ def _bound_scan(x_local):
     if precheck:
         # 1b: min(sum(|x|)) scan — pushes reactions toward zero
         # Effective for resolving lb=0 / ub=0 bounds in one shot.
-        scan_lp, n_scan = _build_abssum_lp(
-            A_eq, b_eq, A_ineq, b_ineq, lb, ub, solver)
+        scan_lp, n_scan = _build_abssum_lp(A_eq, b_eq, A_ineq, b_ineq, lb, ub, solver)
         scan_lp.set_lp_method(LP_METHOD_DUAL)
         n_ext = len(scan_lp.c)
 
@@ -478,9 +487,8 @@ def _bound_scan(x_local):
 
         if verbose:
             n_done_iter = int(res_max.sum() + res_min.sum())
-            logging.debug(
-                f"  Phase 1 min|x|: +{resolved_absmin} "
-                f"({n_done_iter}/{2*n_orig} resolved)")
+            logging.debug(f"  Phase 1 min|x|: +{resolved_absmin} "
+                          f"({n_done_iter}/{2*n_orig} resolved)")
 
         # 1c: Iterative push-to-bounds — directed per-reaction objectives
         # Push unresolved-max reactions toward ub, unresolved-min toward lb.
@@ -510,7 +518,7 @@ def _bound_scan(x_local):
                     x_scan = np.array(x_list_scan[:n_scan], dtype=np.float64)
                     before = int(res_max.sum() + res_min.sum())
                     _bound_scan(x_scan)
-        
+
                     np.maximum(incumbent_max, x_scan, out=incumbent_max)
                     np.minimum(incumbent_min, x_scan, out=incumbent_min)
                     resolved_push_ub = int(res_max.sum() + res_min.sum()) - before
@@ -535,7 +543,7 @@ def _bound_scan(x_local):
                     x_scan = np.array(x_list_scan[:n_scan], dtype=np.float64)
                     before = int(res_max.sum() + res_min.sum())
                     _bound_scan(x_scan)
-        
+
                     np.maximum(incumbent_max, x_scan, out=incumbent_max)
                     np.minimum(incumbent_min, x_scan, out=incumbent_min)
                     resolved_push_lb = int(res_max.sum() + res_min.sum()) - before
@@ -543,10 +551,9 @@ def _bound_scan(x_local):
 
             if verbose:
                 n_done_iter = int(res_max.sum() + res_min.sum())
-                logging.debug(
-                    f"  Phase 1 push {push_iter}: "
-                    f"ub +{resolved_push_ub}, lb +{resolved_push_lb} "
-                    f"({n_done_iter}/{2*n_orig} resolved)")
+                logging.debug(f"  Phase 1 push {push_iter}: "
+                              f"ub +{resolved_push_ub}, lb +{resolved_push_lb} "
+                              f"({n_done_iter}/{2*n_orig} resolved)")
 
             push_iter += 1
             if resolved_this_round < 5:
@@ -570,29 +577,22 @@ def _bound_scan(x_local):
         unresolved = []
         for j in range(n_orig):
             if not res_max[j]:
-                unresolved.append(2 * j)      # even = max
+                unresolved.append(2 * j)  # even = max
             if not res_min[j]:
                 unresolved.append(2 * j + 1)  # odd = min
 
         x_par = [nan] * (2 * n_orig)
         t0 = _time.perf_counter()
         if solver == GLPK:
-            with SDPool(threads, initializer=fva_worker_init_glpk,
-                        initargs=(A_ineq, b_ineq, A_eq, b_eq,
-                                  lb.tolist(), ub.tolist())) as pool:
+            with SDPool(threads, initializer=fva_worker_init_glpk, initargs=(A_ineq, b_ineq, A_eq, b_eq, lb.tolist(), ub.tolist())) as pool:
                 chunk_size = max(1, len(unresolved) // threads)
-                for i, value in pool.imap_unordered(
-                        fva_worker_compute_glpk, unresolved,
-                        chunksize=chunk_size):
+                for i, value in pool.imap_unordered(fva_worker_compute_glpk, unresolved, chunksize=chunk_size):
                     x_par[i] = value
         else:
             with SDPool(threads, initializer=fva_worker_init,
-                        initargs=(A_ineq, b_ineq, A_eq, b_eq,
-                                  lb.tolist(), ub.tolist(), solver)) as pool:
+                        initargs=(A_ineq, b_ineq, A_eq, b_eq, lb.tolist(), ub.tolist(), solver)) as pool:
                 chunk_size = max(1, len(unresolved) // threads)
-                for i, value in pool.imap_unordered(
-                        fva_worker_compute, unresolved,
-                        chunksize=chunk_size):
+                for i, value in pool.imap_unordered(fva_worker_compute, unresolved, chunksize=chunk_size):
                     x_par[i] = value
         t_solve += _time.perf_counter() - t0
         lps_solved += len(unresolved)
@@ -602,9 +602,7 @@ def _bound_scan(x_local):
         if nan_idx:
             _BATCH = 50
             while nan_idx:
-                lp_retry = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq,
-                                   A_eq=A_eq, b_eq=b_eq,
-                                   lb=lb.tolist(), ub=ub.tolist(), solver=solver)
+                lp_retry = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb.tolist(), ub=ub.tolist(), solver=solver)
                 prev_retry = 0
                 for i in nan_idx[:_BATCH]:
                     C = idx2c(i, prev_retry)
@@ -638,8 +636,7 @@ def _bound_scan(x_local):
 
         def _rebuild_lp():
             nonlocal lp, prev_col
-            lp = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq,
-                         lb=lb.tolist(), ub=ub.tolist(), solver=solver)
+            lp = MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb.tolist(), ub=ub.tolist(), solver=solver)
             prev_col = -1
 
         _rebuild_lp()
@@ -738,7 +735,10 @@ def _rebuild_lp():
     incumbent_min[np.abs(incumbent_min) < 1e-11] = 0.0
 
     fva_result = DataFrame(
-        {"minimum": incumbent_min, "maximum": incumbent_max},
+        {
+            "minimum": incumbent_min,
+            "maximum": incumbent_max
+        },
         index=reaction_ids,
     )
 
@@ -746,12 +746,10 @@ def _rebuild_lp():
         cmp_msg = ""
         if cmp_maps:
             cmp_msg = f" (compressed {n_original}→{n_orig} rxns)"
-        logging.debug(
-            f"  speedy_fva done{cmp_msg}: {lps_solved} LPs, "
-            f"{total_bound_resolved} bound-resolved, "
-            f"{2*n_orig} total objectives")
-        logging.debug(
-            f"  timing: solve={t_solve:.2f}s, threads={threads}")
+        logging.debug(f"  speedy_fva done{cmp_msg}: {lps_solved} LPs, "
+                      f"{total_bound_resolved} bound-resolved, "
+                      f"{2*n_orig} total objectives")
+        logging.debug(f"  timing: solve={t_solve:.2f}s, threads={threads}")
 
     t_phase['phase2'] = _time.perf_counter() - _tp2
     t_phase['total'] = _time.perf_counter() - _tp0
@@ -777,10 +775,10 @@ def _rebuild_lp():
 # Fast exact reversibility (sign-only FVA) for pre-compression tightening
 # ---------------------------------------------------------------------------
 
-_REV_SCAN_TOL = 1e-3     # co-option certifies only on flux comfortably above solver noise
+_REV_SCAN_TOL = 1e-3  # co-option certifies only on flux comfortably above solver noise
 _REV_REBUILD_EVERY = 200
-_ZERO_SNAP = 1e-11        # |flux| below this is solver noise, not a direction; 0 disables snapping
-_DEGEN_TOL = 1e-6         # warm-start guard: fresh optimum must not fall below a known-achievable incumbent
+_ZERO_SNAP = 1e-11  # |flux| below this is solver noise, not a direction; 0 disables snapping
+_DEGEN_TOL = 1e-6  # warm-start guard: fresh optimum must not fall below a known-achievable incumbent
 
 
 def _rev_structural_sweep(model):
@@ -797,22 +795,28 @@ def _rev_structural_sweep(model):
         for e in met_rx.values():
             prod = [(i, 'f') for i, c in e if c > 0 and af[i]] + [(i, 'r') for i, c in e if c < 0 and ar[i]]
             cons = [(i, 'f') for i, c in e if c < 0 and af[i]] + [(i, 'r') for i, c in e if c > 0 and ar[i]]
-            prx = {i for i, _ in prod}; crx = {i for i, _ in cons}
+            prx = {i for i, _ in prod}
+            crx = {i for i, _ in cons}
 
             def kill(i, d):
                 nonlocal changed
                 if d == 'f' and af[i]:
-                    af[i] = False; changed = True
+                    af[i] = False
+                    changed = True
                 if d == 'r' and ar[i]:
-                    ar[i] = False; changed = True
+                    ar[i] = False
+                    changed = True
+
             if not prod:
-                for i, d in cons: kill(i, d)
+                for i, d in cons:
+                    kill(i, d)
             elif len(prx) == 1:
                 s = next(iter(prx))
                 for i, d in cons:
                     if i == s: kill(i, d)
             if not cons:
-                for i, d in prod: kill(i, d)
+                for i, d in prod:
+                    kill(i, d)
             elif len(crx) == 1:
                 s = next(iter(crx))
                 for i, d in prod:
@@ -860,33 +864,42 @@ def fast_reversibility(model, solver=None, compress=True):
     lb = np.array([float(r.lower_bound) for r in m.reactions])
     ub = np.array([float(r.upper_bound) for r in m.reactions])
     S = sparse.csr_matrix(create_stoichiometric_matrix(m))
-    A_ineq = sparse.csr_matrix((0, n)); b_ineq = []
-    A_eq = S; b_eq = [0.0] * S.shape[0]
+    A_ineq = sparse.csr_matrix((0, n))
+    b_ineq = []
+    A_eq = S
+    b_eq = [0.0] * S.shape[0]
 
     def build():
-        return MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq,
-                       lb=lb.tolist(), ub=ub.tolist(), solver=solver)
+        return MILP_LP(A_ineq=A_ineq, b_ineq=b_ineq, A_eq=A_eq, b_eq=b_eq, lb=lb.tolist(), ub=ub.tolist(), solver=solver)
+
     lp = build()
 
-    incumbent_max = np.full(n, -np.inf); incumbent_min = np.full(n, np.inf)
+    incumbent_max = np.full(n, -np.inf)
+    incumbent_min = np.full(n, np.inf)
     # Bound signs are exact model data, not solver output, so they decide on the exact sign: a
     # direction counts as blocked only if the bound itself forbids it. A tolerance here would
     # discard a direction whose achievable flux is merely small (a max of 1e-8 is still forward).
-    res_max = ub <= 0.0                   # fwd already blocked by bounds (sweep/original)
+    res_max = ub <= 0.0  # fwd already blocked by bounds (sweep/original)
     res_min = lb >= 0.0
     incumbent_max[res_max] = ub[res_max]
     incumbent_min[res_min] = lb[res_min]
     fixed = ub == lb
-    res_max[fixed] = True; res_min[fixed] = True
-    incumbent_max[fixed] = ub[fixed]; incumbent_min[fixed] = lb[fixed]
+    res_max[fixed] = True
+    res_min[fixed] = True
+    incumbent_max[fixed] = ub[fixed]
+    incumbent_min[fixed] = lb[fixed]
 
     def scan(x):
-        nm = (~res_max) & (x > _REV_SCAN_TOL); res_max[nm] = True
+        nm = (~res_max) & (x > _REV_SCAN_TOL)
+        res_max[nm] = True
         np.maximum(incumbent_max, x, out=incumbent_max)
-        nn = (~res_min) & (x < -_REV_SCAN_TOL); res_min[nn] = True
+        nn = (~res_min) & (x < -_REV_SCAN_TOL)
+        res_min[nn] = True
         np.minimum(incumbent_min, x, out=incumbent_min)
 
-    n_lp = 0; prev_col = -1; seq = 0
+    n_lp = 0
+    prev_col = -1
+    seq = 0
 
     def solve_dir(j, direction):
         nonlocal prev_col, seq, n_lp
@@ -897,14 +910,17 @@ def solve_dir(j, direction):
         else:
             lp.set_objective_idx(C)
         prev_col = j
-        r = lp.solve(); n_lp += 1; seq += 1
+        r = lp.solve()
+        n_lp += 1
+        seq += 1
         return r
 
     # Feasibility preflight: one zero-objective solve proves the polytope is non-empty (so a later
     # infeasible status can only come from the objective change, not the model) and its flux vector
     # seeds every incumbent, so no subsequent warm-start optimum can silently contradict a flux
     # already witnessed as achievable.
-    x_feas, _, status_feas = lp.solve(); n_lp += 1
+    x_feas, _, status_feas = lp.solve()
+    n_lp += 1
     if status_feas == INFEASIBLE:
         raise ValueError('fast_reversibility: the model has no steady-state flux distribution.')
     if status_feas == OPTIMAL and x_feas:
@@ -914,15 +930,20 @@ def unknown(j, direction):
         """Record 'direction not determined': the incumbent goes to +/-inf so the direction is
         reported as achievable. Tightening on an unproven bound could delete a feasible pathway;
         reporting a spurious direction only forgoes tightening."""
-        if direction == 1: res_max[j] = True; incumbent_max[j] = np.inf
-        else: res_min[j] = True; incumbent_min[j] = -np.inf
+        if direction == 1:
+            res_max[j] = True
+            incumbent_max[j] = np.inf
+        else:
+            res_min[j] = True
+            incumbent_min[j] = -np.inf
 
     for j in range(n):
         for direction in (1, -1):
             if (direction == 1 and res_max[j]) or (direction == -1 and res_min[j]):
                 continue
             if seq > 0 and seq % _REV_REBUILD_EVERY == 0:
-                lp = build(); prev_col = -1
+                lp = build()
+                prev_col = -1
             x_list, obj_val, status = solve_dir(j, direction)
             if status != OPTIMAL:
                 # UNBOUNDED is a proven infinite direction, every other nonoptimal status (time
@@ -934,14 +955,19 @@ def unknown(j, direction):
             degen = (direction == 1 and np.isfinite(inc) and val < inc - _DEGEN_TOL * (1 + abs(inc))) or \
                     (direction == -1 and np.isfinite(inc) and val > inc + _DEGEN_TOL * (1 + abs(inc)))
             if degen:
-                lp = build(); prev_col = -1
+                lp = build()
+                prev_col = -1
                 x_list, obj_val, status = solve_dir(j, direction)
                 if status != OPTIMAL:
                     unknown(j, direction)
                     continue
                 val = -obj_val if direction == 1 else obj_val
-            if direction == 1: res_max[j] = True; incumbent_max[j] = max(incumbent_max[j], val)
-            else: res_min[j] = True; incumbent_min[j] = min(incumbent_min[j], val)
+            if direction == 1:
+                res_max[j] = True
+                incumbent_max[j] = max(incumbent_max[j], val)
+            else:
+                res_min[j] = True
+                incumbent_min[j] = min(incumbent_min[j], val)
             scan(np.array(x_list[:n], dtype=np.float64))
 
     # (4) expand compressed min/max back to original reactions
@@ -952,5 +978,4 @@ def unknown(j, direction):
     df = DataFrame({"minimum": incumbent_min, "maximum": incumbent_max}, index=cmp_rid)
     if cmp_maps:
         df = _expand_fva(df, cmp_maps, orig_rid)
-    return {r: (float(df.at[r, 'maximum']) > 0.0, float(df.at[r, 'minimum']) < 0.0)
-            for r in orig_rid}
+    return {r: (float(df.at[r, 'maximum']) > 0.0, float(df.at[r, 'minimum']) < 0.0) for r in orig_rid}
diff --git a/straindesign/strainDesignMILP.py b/straindesign/strainDesignMILP.py
index 325e77e..f7c603a 100644
--- a/straindesign/strainDesignMILP.py
+++ b/straindesign/strainDesignMILP.py
@@ -116,7 +116,7 @@ def _trim_z_variables(self):
             self._z_orig_indices = None  # no trimming needed
             return
 
-        self._z_orig_indices = keep_z       # trimmed_idx -> orig_idx
+        self._z_orig_indices = keep_z  # trimmed_idx -> orig_idx
         self._orig_num_z = self.num_z
 
         n_cont = len(self.c) - self.num_z
diff --git a/straindesign/strainDesignModule.py b/straindesign/strainDesignModule.py
index df24d3f..e7cd16a 100644
--- a/straindesign/strainDesignModule.py
+++ b/straindesign/strainDesignModule.py
@@ -231,8 +231,8 @@ def __init__(self, model, module_type, *args, **kwargs):
         self[MODEL_ID] = model.id
         self[MODULE_TYPE] = module_type
         allowed_keys = {
-            CONSTRAINTS, INNER_OBJECTIVE, INNER_OPT_SENSE, OUTER_OBJECTIVE, OUTER_OPT_SENSE,
-            INNER_OPT_TOL, OUTER_OPT_TOL, PROD_ID, 'skip_checks', MIN_GCP, 'reac_ids'
+            CONSTRAINTS, INNER_OBJECTIVE, INNER_OPT_SENSE, OUTER_OBJECTIVE, OUTER_OPT_SENSE, INNER_OPT_TOL, OUTER_OPT_TOL, PROD_ID,
+            'skip_checks', MIN_GCP, 'reac_ids'
         }
         # set all keys passed in kwargs as properties of the SD_Module object
         for key, value in kwargs.items():
@@ -262,8 +262,8 @@ def __init__(self, model, module_type, *args, **kwargs):
             elif self[INNER_OPT_SENSE] not in [MINIMIZE, MAXIMIZE] or self[OUTER_OPT_SENSE] not in [MINIMIZE, MAXIMIZE]:
                 raise Exception('Inner and outer optimization sense must be "' + MINIMIZE + '" or "' + MAXIMIZE + '" (default).')
             if ((self[INNER_OBJECTIVE] == None) or (self[OUTER_OBJECTIVE] == None)):
-                raise Exception('When module type is "' + OPTKNOCK + '", "' + ROBUSTKNOCK +
-                                '" or "' + DOUBLEOPT + '", an inner and outer objective function must be provided.')
+                raise Exception('When module type is "' + OPTKNOCK + '", "' + ROBUSTKNOCK + '" or "' + DOUBLEOPT +
+                                '", an inner and outer objective function must be provided.')
         elif (self[MODULE_TYPE] == OPTCOUPLE):
             if self[INNER_OPT_SENSE] is None:
                 self[INNER_OPT_SENSE] = MAXIMIZE
@@ -277,8 +277,7 @@ def __init__(self, model, module_type, *args, **kwargs):
                 raise Exception('When module type is "' + OPTCOUPLE + '", the production reaction id must be provided.')
         if self[MODULE_TYPE] in [PROTECT, SUPPRESS] and self[OUTER_OBJECTIVE] is not None:
             if self[INNER_OBJECTIVE] is None:
-                raise Exception('When outer_objective is set for "' + self[MODULE_TYPE] +
-                                '", inner_objective must also be provided.')
+                raise Exception('When outer_objective is set for "' + self[MODULE_TYPE] + '", inner_objective must also be provided.')
             if self[OUTER_OPT_SENSE] is None:
                 self[OUTER_OPT_SENSE] = MAXIMIZE
             if self[OUTER_OPT_SENSE] not in [MINIMIZE, MAXIMIZE]:
diff --git a/straindesign/strainDesignProblem.py b/straindesign/strainDesignProblem.py
index d994c5e..6bfe031 100644
--- a/straindesign/strainDesignProblem.py
+++ b/straindesign/strainDesignProblem.py
@@ -278,9 +278,9 @@ def _module_bound_override(self, sd_module):
                 logging.warning('  Module FVA range for %s is inconsistent (min %g > max %g), no bound '
                                 'override applied.' % (rid, lim.minimum, lim.maximum))
                 continue
-            lo = 0.0 if lim.minimum >= tol else None       # convincingly non-negative in-module
-            hi = 0.0 if lim.maximum <= -tol else None      # convincingly non-positive in-module
-            if lo is not None or hi is not None:           # both sides fixed == blocked in-module
+            lo = 0.0 if lim.minimum >= tol else None  # convincingly non-negative in-module
+            hi = 0.0 if lim.maximum <= -tol else None  # convincingly non-positive in-module
+            if lo is not None or hi is not None:  # both sides fixed == blocked in-module
                 override[rid] = (lo, hi)
         return override
 
@@ -335,23 +335,19 @@ def addModule(self, sd_module):
                 # Inequality: c_v · x_v + tol · c_dual · d <= 0  (actual >= tol * optimal)
                 n_v, n_ref, n_d = len(c_v), len(c_inner), len(c_inner_dual)
                 c_dual_scaled = [inner_opt_tol * c for c in c_inner_dual]
-                ineq_link = sparse.hstack((sparse.csr_matrix(c_v), sparse.csr_matrix((1, n_ref)),
-                                           sparse.csr_matrix(c_dual_scaled)))
-                eq_link = sparse.hstack((sparse.csr_matrix((1, n_v)), sparse.csr_matrix(c_inner),
-                                          sparse.csr_matrix(c_inner_dual)))
-                A_ineq_p = sparse.vstack((sparse.block_diag((A_ineq_v, A_ineq_inner, A_ineq_dual)),
-                                          ineq_link)).tocsr()
+                ineq_link = sparse.hstack((sparse.csr_matrix(c_v), sparse.csr_matrix((1, n_ref)), sparse.csr_matrix(c_dual_scaled)))
+                eq_link = sparse.hstack((sparse.csr_matrix((1, n_v)), sparse.csr_matrix(c_inner), sparse.csr_matrix(c_inner_dual)))
+                A_ineq_p = sparse.vstack((sparse.block_diag((A_ineq_v, A_ineq_inner, A_ineq_dual)), ineq_link)).tocsr()
                 b_ineq_p = b_ineq_v + b_ineq_inner + b_ineq_dual + [0.0]
-                A_eq_p = sparse.vstack((sparse.block_diag((A_eq_v, A_eq_inner, A_eq_dual)),
-                                        eq_link)).tocsr()
+                A_eq_p = sparse.vstack((sparse.block_diag((A_eq_v, A_eq_inner, A_eq_dual)), eq_link)).tocsr()
                 b_eq_p = b_eq_v + b_eq_inner + b_eq_dual + [0.0]
                 lb_p = lb_v + lb_inner + lb_dual
                 ub_p = ub_v + ub_inner + ub_dual
                 z_map_vars_p = sparse.hstack((z_map_vars_v, z_map_vars_inner, z_map_vars_dual))
-                z_map_constr_ineq_p = sparse.hstack((z_map_constr_ineq_v, z_map_constr_ineq_inner,
-                    z_map_constr_ineq_dual, sparse.csc_matrix((self.num_z, 1))))
-                z_map_constr_eq_p = sparse.hstack((z_map_constr_eq_v, z_map_constr_eq_inner,
-                    z_map_constr_eq_dual, sparse.csc_matrix((self.num_z, 1))))
+                z_map_constr_ineq_p = sparse.hstack(
+                    (z_map_constr_ineq_v, z_map_constr_ineq_inner, z_map_constr_ineq_dual, sparse.csc_matrix((self.num_z, 1))))
+                z_map_constr_eq_p = sparse.hstack(
+                    (z_map_constr_eq_v, z_map_constr_eq_inner, z_map_constr_eq_dual, sparse.csc_matrix((self.num_z, 1))))
             else:
                 # Exact inner optimality (original code)
                 A_ineq_p = sparse.block_diag((A_ineq_v, A_ineq_dual)).tocsr()
@@ -679,12 +675,10 @@ def addModule(self, sd_module):
                 A_eq_base = sparse.block_diag((A_eq_p, A_eq_p, A_eq_dl))
                 b_eq_ext = b_eq_p + b_eq_p + b_eq_dl
                 # Equality anchor: c_out · x_ref + c_dl · d_out = 0
-                eq_link = sparse.hstack((sparse.csr_matrix((1, n_p)),
-                    sparse.csr_matrix(c_out_in_p), sparse.csr_matrix(c_dl)))
+                eq_link = sparse.hstack((sparse.csr_matrix((1, n_p)), sparse.csr_matrix(c_out_in_p), sparse.csr_matrix(c_dl)))
                 # Relaxed inequality: c_out · x_actual + tol * c_dl · d_out <= 0
                 c_dl_scaled = [outer_opt_tol * v for v in c_dl]
-                iq_link = sparse.hstack((sparse.csr_matrix(c_out_in_p),
-                    sparse.csr_matrix((1, n_p)), sparse.csr_matrix(c_dl_scaled)))
+                iq_link = sparse.hstack((sparse.csr_matrix(c_out_in_p), sparse.csr_matrix((1, n_p)), sparse.csr_matrix(c_dl_scaled)))
                 A_eq_p = sparse.vstack((A_eq_base, eq_link)).tocsr()
                 b_eq_p = b_eq_ext + [0.0]
                 A_ineq_p = sparse.vstack((A_ineq_ext, iq_link)).tocsr()
@@ -692,16 +686,16 @@ def addModule(self, sd_module):
                 lb_p = lb_p + lb_p + lb_dl
                 ub_p = ub_p + ub_p + ub_dl
                 z_map_vars_p = sparse.hstack((z_map_vars_p, z_map_vars_p, z_map_vars_dl))
-                z_map_constr_ineq_p = sparse.hstack((z_map_constr_ineq_p, z_map_constr_ineq_p,
-                    z_map_constr_ineq_dl, sparse.csc_matrix((self.num_z, 1))))
-                z_map_constr_eq_p = sparse.hstack((z_map_constr_eq_p, z_map_constr_eq_p,
-                    z_map_constr_eq_dl, sparse.csc_matrix((self.num_z, 1))))
+                z_map_constr_ineq_p = sparse.hstack(
+                    (z_map_constr_ineq_p, z_map_constr_ineq_p, z_map_constr_ineq_dl, sparse.csc_matrix((self.num_z, 1))))
+                z_map_constr_eq_p = sparse.hstack(
+                    (z_map_constr_eq_p, z_map_constr_eq_p, z_map_constr_eq_dl, sparse.csc_matrix((self.num_z, 1))))
             else:
                 # Exact outer optimality (original code)
                 A_ineq_p = sparse.block_diag((A_ineq_p, A_ineq_dl)).tocsr()
                 b_ineq_p = b_ineq_p + b_ineq_dl
-                A_eq_p = sparse.vstack((sparse.block_diag((A_eq_p, A_eq_dl)),
-                    sparse.hstack((sparse.csr_matrix(c_out_in_p), sparse.csr_matrix(c_dl))))).tocsr()
+                A_eq_p = sparse.vstack((sparse.block_diag(
+                    (A_eq_p, A_eq_dl)), sparse.hstack((sparse.csr_matrix(c_out_in_p), sparse.csr_matrix(c_dl))))).tocsr()
                 b_eq_p = b_eq_p + b_eq_dl + [0.0]
                 lb_p = lb_p + lb_dl
                 ub_p = ub_p + ub_dl
@@ -817,7 +811,7 @@ def link_z(self):
         #    constraint (gurobi/cplex) or the constant self.M (glpk/user-M).
         knockable_constr_ineq = np.unique(self.z_map_constr_ineq.nonzero()[1])
 
-        _idxz = set(self.idx_z)                                          # O(1) membership in the scan below
+        _idxz = set(self.idx_z)  # O(1) membership in the scan below
         cont_vars = [i not in _idxz for i in range(0, numvars)]
         M_lb = [self.lb[i] for i in np.nonzero(cont_vars)[0]]
         M_ub = [self.ub[i] for i in np.nonzero(cont_vars)[0]]
@@ -849,7 +843,7 @@ def link_z(self):
                 max_Ax[i] = np.inf
                 n_multi += 1
         logging.info('  Bounding MILP: %d constraints (%d zero, %d single-var, %d multi-var->indicator/M).' %
-                      (num_Ms, n_zero, n_single, n_multi))
+                     (num_Ms, n_zero, n_single, n_multi))
 
         # round Ms up to 5 digits
         Ms = [np.ceil(M * 1e5) / 1e5 if not isinf(M) else self.M for M in max_Ax]
@@ -1038,9 +1032,9 @@ def link_z(self):
                 continue
             col = Aic.getcol(z).tocoo()
             if any((r not in budget_rows) and v != 0 for r, v in zip(col.row.tolist(), col.data.tolist())):
-                continue                       # finite-M z-link -> controls a row
+                continue  # finite-M z-link -> controls a row
             if Aec is not None and Aec.getcol(z).nnz:
-                continue                       # equality z-link
+                continue  # equality z-link
             self.ub[z] = 0.0
             n_free += 1
         if n_free:
@@ -1126,7 +1120,7 @@ def build_primal_from_cbm(model, V_ineq=None, v_ineq=None, V_eq=None, v_eq=None,
                     lb[i] = max(lb[i], float(lo))
                 if hi is not None:
                     ub[i] = min(ub[i], float(hi))
-                if lb[i] > ub[i]:            # numeric guard: never emit an inconsistent block
+                if lb[i] > ub[i]:  # numeric guard: never emit an inconsistent block
                     lb[i], ub[i] = float(lo), float(hi)
     z_map_vars = sparse.identity(numr, 'd', format="csc")
     z_map_constr_eq = sparse.csc_matrix((numr, A_eq.shape[0]))
diff --git a/straindesign/strainDesignSolutions.py b/straindesign/strainDesignSolutions.py
index 2d86e51..2889f27 100644
--- a/straindesign/strainDesignSolutions.py
+++ b/straindesign/strainDesignSolutions.py
@@ -185,13 +185,13 @@ def _translate_genes_to_reactions(sd_list, model):
                     candidate_reacs.update(r.id for r in model.genes.get_by_id(g).reactions)
             for r in candidate_reacs:
                 gpr_r = rxn_gpr[r]
-                if gpr_r.eval(ko_off):             # reaction still possible under the interventions
-                    if not gpr_r.eval(all_off):    # ... only because of a knock-in
+                if gpr_r.eval(ko_off):  # reaction still possible under the interventions
+                    if not gpr_r.eval(all_off):  # ... only because of a knock-in
                         reac_ki.add(r)
-                else:                              # reaction dead under the interventions
-                    if gpr_r.eval(noki_off):       # ... the knock-out is what killed it
+                else:  # reaction dead under the interventions
+                    if gpr_r.eval(noki_off):  # ... the knock-out is what killed it
                         reac_ko.add(r)
-                    else:                          # ... dead regardless (e.g. an un-made knock-in)
+                    else:  # ... dead regardless (e.g. an un-made knock-in)
                         reac_no_ki.add(r)
             reaction_sd[i].update({k: -1.0 for k in reac_ko})
             reaction_sd[i].update({k: 1.0 for k in reac_ki})
@@ -461,8 +461,7 @@ def expand_group(self, grp_idx):
 
         # Expand
         expanded = expand_sd([compressed_sd[grp_idx].copy()], cmp_mapReac)
-        expanded = filter_sd_maxcost(expanded, meta['max_cost'],
-                                     meta['uncmp_ko_cost'], meta['uncmp_ki_cost'])
+        expanded = filter_sd_maxcost(expanded, meta['max_cost'], meta['uncmp_ko_cost'], meta['uncmp_ki_cost'])
         # Postprocess regulatory interventions (inline to avoid circular import)
         reg_cost = meta.get('uncmp_reg_cost', {})
         for s in expanded:
@@ -475,12 +474,11 @@ def expand_group(self, grp_idx):
 
         # GPR translation + costs/bounds
         if self._model is None:
-            raise RuntimeError(
-                'This SDSolutions was loaded without a model, so compressed '
-                'groups cannot be expanded. Reload with '
-                'SDSolutions.load(file, model=True) to rebuild the embedded '
-                'model, or call attach_model(model) with the original model '
-                'that was passed to compute_strain_designs.')
+            raise RuntimeError('This SDSolutions was loaded without a model, so compressed '
+                               'groups cannot be expanded. Reload with '
+                               'SDSolutions.load(file, model=True) to rebuild the embedded '
+                               'model, or call attach_model(model) with the original model '
+                               'that was passed to compute_strain_designs.')
         model = self._model
         if self.is_gene_sd:
             reaction_sd_exp, gene_sd_exp = self._translate_genes_to_reactions(expanded, model)
@@ -490,8 +488,7 @@ def expand_group(self, grp_idx):
             gene_sd_exp = None
             cost_sd = expanded
 
-        sd_cost_exp, itv_bounds_exp, has_complex = self._compute_costs_and_bounds(
-            cost_sd, reaction_sd_exp, model, self.sd_setup)
+        sd_cost_exp, itv_bounds_exp, has_complex = self._compute_costs_and_bounds(cost_sd, reaction_sd_exp, model, self.sd_setup)
         if has_complex:
             self.has_complex_regul_itv = True
 
@@ -676,11 +673,11 @@ def load(cls, filename, model=None, cmp_model=None):
             obj = pickle.load(f)
 
         def _resolve(arg, embedded):
-            if arg is True:                       # rebuild the embedded snapshot
+            if arg is True:  # rebuild the embedded snapshot
                 return model_from_dict(embedded) if embedded is not None else None
-            if arg is None or arg is False:       # attach nothing
+            if arg is None or arg is False:  # attach nothing
                 return None
-            return arg                            # an explicit cobra model
+            return arg  # an explicit cobra model
 
         obj._model = _resolve(model, getattr(obj, '_embedded_model_dict', None))
         obj._cmp_model = _resolve(cmp_model, getattr(obj, '_embedded_cmp_model_dict', None))
diff --git a/tests/conftest.py b/tests/conftest.py
index f4556fc..23d3ff5 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -3,15 +3,15 @@
 from cobra import Configuration
 from straindesign.names import *
 
-
 # ---------------------------------------------------------------------------
 # Custom CLI flags for test_performance.py tiered benchmarks
 # ---------------------------------------------------------------------------
 
+
 def pytest_addoption(parser):
     for name, help_text in [
         ("--medium", "Run iMLcore genome-scale benchmarks (~4 min total)."),
-        ("--large",  "Run iML1515 large-model benchmarks (several min/solver)."),
+        ("--large", "Run iML1515 large-model benchmarks (several min/solver)."),
     ]:
         try:
             parser.addoption(name, action="store_true", default=False, help=help_text)
@@ -22,14 +22,12 @@ def pytest_addoption(parser):
 def pytest_configure(config):
     for marker, desc in [
         ("medium", "genome-scale benchmark; enable with --medium"),
-        ("large",  "large-model benchmark; enable with --large"),
+        ("large", "large-model benchmark; enable with --large"),
     ]:
         config.addinivalue_line("markers", f"{marker}: {desc}")
     # Suppress known third-party warnings
-    config.addinivalue_line("filterwarnings",
-        "ignore:FigureCanvasTemplate is non-interactive:UserWarning")
-    config.addinivalue_line("filterwarnings",
-        "ignore:builtin type Swig.*has no __module__ attribute:DeprecationWarning")
+    config.addinivalue_line("filterwarnings", "ignore:FigureCanvasTemplate is non-interactive:UserWarning")
+    config.addinivalue_line("filterwarnings", "ignore:builtin type Swig.*has no __module__ attribute:DeprecationWarning")
 
 
 def pytest_collection_modifyitems(config, items):
@@ -40,6 +38,7 @@ def pytest_collection_modifyitems(config, items):
                 if marker in item.keywords:
                     item.add_marker(skip)
 
+
 cobra_conf = Configuration()
 bound_thres = max((abs(cobra_conf.lower_bound), abs(cobra_conf.upper_bound)))
 
diff --git a/tests/test_03_plots.py b/tests/test_03_plots.py
index 9c1e1b5..0805974 100644
--- a/tests/test_03_plots.py
+++ b/tests/test_03_plots.py
@@ -28,8 +28,7 @@ def test_plot_3d_space(curr_solver, model_weak_coupling):
 def test_plot_2d_point(curr_solver, model_weak_coupling):
     """0D: fix both axes to single values, should show a dot."""
     constr = ['r4 = 0', 'r7 = 0', 'r9 = 0', 'r_P = 2', 'r_BM = 4']
-    dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM'),
-                                    constraints=constr, plt_backend='template')
+    dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM'), constraints=constr, plt_backend='template')
     assert len(dp) == 1
 
 
@@ -37,8 +36,7 @@ def test_plot_2d_point(curr_solver, model_weak_coupling):
 def test_plot_2d_line(curr_solver, model_weak_coupling):
     """1D: fix one axis, line along the other."""
     constr = ['r4 = 0', 'r7 = 0', 'r9 = 0', 'r_P = 2']
-    dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM'),
-                                    constraints=constr, plt_backend='template')
+    dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM'), constraints=constr, plt_backend='template')
     assert len(dp) >= 2
 
 
@@ -46,6 +44,5 @@ def test_plot_2d_line(curr_solver, model_weak_coupling):
 def test_plot_3d_point(curr_solver, model_weak_coupling):
     """0D in 3D: all axes fixed, should show a dot."""
     constr = ['r4 = 0', 'r7 = 0', 'r9 = 0', 'r_BM = 2', 'r_P = 2', 'r_Q = 2']
-    dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM', 'r_Q'),
-                                    constraints=constr, plt_backend='template')
+    dp, tri, _ = sd.plot_flux_space(model_weak_coupling, ('r_P', 'r_BM', 'r_Q'), constraints=constr, plt_backend='template')
     assert len(dp) == 1
diff --git a/tests/test_04_preprocessing.py b/tests/test_04_preprocessing.py
index a0f60d9..005e7a5 100644
--- a/tests/test_04_preprocessing.py
+++ b/tests/test_04_preprocessing.py
@@ -20,9 +20,9 @@
 from cobra.core.gene import GPR
 from sympy import simplify_logic
 
-
 # ── GPR extension + compression ──────────────────────────────────────
 
+
 @pytest.mark.timeout(15)
 def test_gpr_extension_compression1(model_gpr):
     gkocost = {
@@ -61,33 +61,33 @@ def test_gpr_extension_compression2(model_gpr):
 
 # ── GPR propagation helper unit tests ────────────────────────────────
 
+
 class TestGprAstToExpr:
     # None is special-cased (identity check); the rest share one assertion shape.
     def test_none_returns_none(self):
         assert _gpr_ast_to_expr(None) is None
 
-    @pytest.mark.parametrize("node,expected", [
-        (ast.Name(id='g1'), 'g1'),
-        (ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]),
-         ('and', ('g1', 'g2'))),
-        (ast.BoolOp(op=ast.Or(), values=[ast.Name(id='g1'), ast.Name(id='g2')]),
-         ('or', ('g1', 'g2'))),
-        # (g1 and g2) or g3
-        (ast.BoolOp(op=ast.Or(), values=[
-            ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]),
-            ast.Name(id='g3')]),
-         ('or', (('and', ('g1', 'g2')), 'g3'))),
-        # g1 and (g2 and g3) -> flattened
-        (ast.BoolOp(op=ast.And(), values=[
-            ast.Name(id='g1'),
-            ast.BoolOp(op=ast.And(), values=[ast.Name(id='g2'), ast.Name(id='g3')])]),
-         ('and', ('g1', 'g2', 'g3'))),
-    ], ids=['single_gene', 'and', 'or', 'nested', 'nested_same_op_flattened'])
+    @pytest.mark.parametrize(
+        "node,expected",
+        [
+            (ast.Name(id='g1'), 'g1'),
+            (ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), ('and', ('g1', 'g2'))),
+            (ast.BoolOp(op=ast.Or(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), ('or', ('g1', 'g2'))),
+            # (g1 and g2) or g3
+            (ast.BoolOp(op=ast.Or(), values=[ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]),
+                                             ast.Name(id='g3')]), ('or', (('and', ('g1', 'g2')), 'g3'))),
+            # g1 and (g2 and g3) -> flattened
+            (ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'),
+                                              ast.BoolOp(op=ast.And(), values=[ast.Name(id='g2'), ast.Name(id='g3')])]),
+             ('and', ('g1', 'g2', 'g3'))),
+        ],
+        ids=['single_gene', 'and', 'or', 'nested', 'nested_same_op_flattened'])
     def test_gpr_ast_to_expr(self, node, expected):
         assert _gpr_ast_to_expr(node) == expected
 
 
 class TestExprToGprString:
+
     @pytest.mark.parametrize("expr,expected", [
         (None, ''),
         ('g1', 'g1'),
@@ -95,7 +95,8 @@ class TestExprToGprString:
         (('or', ['g1', 'g2']), 'g1 or g2'),
         (('or', [('and', ['g1', 'g2']), 'g3']), '(g1 and g2) or g3'),
         (('and', [('or', ['g1', 'g2']), 'g3']), '(g1 or g2) and g3'),
-    ], ids=['none', 'single_gene', 'and', 'or', 'nested_and_in_or', 'nested_or_in_and'])
+    ],
+                             ids=['none', 'single_gene', 'and', 'or', 'nested_and_in_or', 'nested_or_in_and'])
     def test_expr_to_gpr_string(self, expr, expected):
         assert _expr_to_gpr_string(expr) == expected
 
@@ -112,38 +113,42 @@ def test_roundtrips_through_cobra(self):
 
 
 class TestCombineGprAnd:
-    @pytest.mark.parametrize("nodes,expected", [
-        ([None, None], ''),
-        ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')])],
-         'g1 and g2'),
-        # Empty GPR (None) should be skipped in AND combination.
-        ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), None, None],
-         'g1 and g2'),
-        ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]),
-          ast.Name(id='g3')],
-         'g1 and g2 and g3'),
-        # AND of overlapping expressions should simplify: (g1 and g2) AND (g1 and g3).
-        ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]),
-          ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g3')])],
-         'g1 and g2 and g3'),
-        ([], ''),
-    ], ids=['all_empty', 'single_non_empty', 'skip_empty', 'two_non_empty',
-            'simplification', 'empty_list'])
+
+    @pytest.mark.parametrize(
+        "nodes,expected",
+        [
+            ([None, None], ''),
+            ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')])], 'g1 and g2'),
+            # Empty GPR (None) should be skipped in AND combination.
+            ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]), None, None], 'g1 and g2'),
+            ([ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]),
+              ast.Name(id='g3')], 'g1 and g2 and g3'),
+            # AND of overlapping expressions should simplify: (g1 and g2) AND (g1 and g3).
+            ([
+                ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g2')]),
+                ast.BoolOp(op=ast.And(), values=[ast.Name(id='g1'), ast.Name(id='g3')])
+            ], 'g1 and g2 and g3'),
+            ([], ''),
+        ],
+        ids=['all_empty', 'single_non_empty', 'skip_empty', 'two_non_empty', 'simplification', 'empty_list'])
     def test_combine_gprs_and(self, nodes, expected):
         assert _combine_gprs(nodes, 'and') == expected
 
 
 class TestCombineGprOr:
-    @pytest.mark.parametrize("nodes,expected", [
-        # If any reaction has empty GPR (always active), result is empty.
-        ([ast.Name(id='g1'), None], ''),
-        ([None, None], ''),
-        ([ast.Name(id='g1'), ast.Name(id='g2')], 'g1 or g2'),
-        # OR with duplicate terms should deduplicate: g1 OR g1 -> g1.
-        ([ast.Name(id='g1'), ast.Name(id='g1')], 'g1'),
-        ([], ''),
-    ], ids=['any_empty_returns_empty', 'all_empty', 'two_non_empty',
-            'deduplication', 'empty_list'])
+
+    @pytest.mark.parametrize(
+        "nodes,expected",
+        [
+            # If any reaction has empty GPR (always active), result is empty.
+            ([ast.Name(id='g1'), None], ''),
+            ([None, None], ''),
+            ([ast.Name(id='g1'), ast.Name(id='g2')], 'g1 or g2'),
+            # OR with duplicate terms should deduplicate: g1 OR g1 -> g1.
+            ([ast.Name(id='g1'), ast.Name(id='g1')], 'g1'),
+            ([], ''),
+        ],
+        ids=['any_empty_returns_empty', 'all_empty', 'two_non_empty', 'deduplication', 'empty_list'])
     def test_combine_gprs_or(self, nodes, expected):
         assert _combine_gprs(nodes, 'or') == expected
 
@@ -158,12 +163,14 @@ def test_no_absorption(self):
 
 # ── GPR propagation integration tests (model_gpr.xml) ────────────────
 
+
 @pytest.fixture
 def gpr_model():
     return read_sbml_model('tests/model_gpr.xml')
 
 
 class TestModelGprCompression:
+
     def test_coupled_compression_propagates_gpr(self, gpr_model):
         """Coupled compression should AND-combine GPR rules, skipping empty ones."""
         remove_blocked_reactions(gpr_model)
@@ -234,8 +241,7 @@ def test_coupled_group_r4_r5_r6_rdex(self, gpr_model):
 
         from cobra.core.gene import GPR
         result_sympy = GPR.from_string(target_rxn.gene_reaction_rule).as_symbolic()
-        expected = GPR.from_string(
-            '(g1 and g4 and g7 and g8) or (g1 and g4 and g5 and g8 and g9)').as_symbolic()
+        expected = GPR.from_string('(g1 and g4 and g7 and g8) or (g1 and g4 and g5 and g8 and g9)').as_symbolic()
 
         assert simplify_logic(result_sympy ^ expected) == False, \
             f"GPR mismatch. Got: {result_sympy}, expected: {expected}"
@@ -267,4 +273,3 @@ def test_coupled_group_r3_rpex(self, gpr_model):
 
         assert simplify_logic(result_sympy ^ expected) == False, \
             f"GPR mismatch. Got: {result_sympy}, expected: {expected}"
-
diff --git a/tests/test_05_straindesign.py b/tests/test_05_straindesign.py
index e3381a5..f24ea1f 100644
--- a/tests/test_05_straindesign.py
+++ b/tests/test_05_straindesign.py
@@ -259,32 +259,47 @@ def test_doubleopt(curr_solver, model_doubleopt):
     need to achieve 60% of their optima.
     """
     # Knockable: organism-internal reactions + shared sinks (enable cross-feeding KOs)
-    ko_reacs = [r.id for r in model_doubleopt.reactions if r.id.startswith('A_R') or r.id.startswith('B_R')
-                or r.id in ['A_10', 'B_10', 'shared_R_D', 'shared_R_C']]
+    ko_reacs = [
+        r.id
+        for r in model_doubleopt.reactions
+        if r.id.startswith('A_R') or r.id.startswith('B_R') or r.id in ['A_10', 'B_10', 'shared_R_D', 'shared_R_C']
+    ]
     kocost = {r: 1 for r in ko_reacs}
 
     # --- Exact DOUBLEOPT: tight coupling via shared sink KOs ---
-    modules_exact = [sd.SDModule(model_doubleopt, DOUBLEOPT,
-                                 inner_objective='A_BM',
-                                 outer_objective='B_BM',
-                                 constraints=['A_BM >= 0.1', 'B_BM >= 0.1'])]
+    modules_exact = [
+        sd.SDModule(model_doubleopt, DOUBLEOPT, inner_objective='A_BM', outer_objective='B_BM', constraints=['A_BM >= 0.1', 'B_BM >= 0.1'])
+    ]
     sol_exact = sd.compute_strain_designs(model_doubleopt.copy(),
-        sd_modules=modules_exact, max_cost=6, max_solutions=inf,
-        solution_approach=POPULATE, ko_cost=kocost, solver=curr_solver, compress=False)
+                                          sd_modules=modules_exact,
+                                          max_cost=6,
+                                          max_solutions=inf,
+                                          solution_approach=POPULATE,
+                                          ko_cost=kocost,
+                                          solver=curr_solver,
+                                          compress=False)
     assert len(sol_exact.reaction_sd) > 0, \
         "Exact DOUBLEOPT should find solutions when shared sinks are knockable"
     min_cost_exact = min(sum(abs(v) for v in s.values()) for s in sol_exact.reaction_sd)
 
     # --- Relaxed DOUBLEOPT: 60% optimality tolerance finds cheaper solutions ---
-    modules_relaxed = [sd.SDModule(model_doubleopt, DOUBLEOPT,
-                                   inner_objective='A_BM',
-                                   outer_objective='B_BM',
-                                   inner_opt_tol=0.6,
-                                   outer_opt_tol=0.6,
-                                   constraints=['A_BM >= 0.1', 'B_BM >= 0.1'])]
+    modules_relaxed = [
+        sd.SDModule(model_doubleopt,
+                    DOUBLEOPT,
+                    inner_objective='A_BM',
+                    outer_objective='B_BM',
+                    inner_opt_tol=0.6,
+                    outer_opt_tol=0.6,
+                    constraints=['A_BM >= 0.1', 'B_BM >= 0.1'])
+    ]
     sol_relaxed = sd.compute_strain_designs(model_doubleopt.copy(),
-        sd_modules=modules_relaxed, max_cost=3, max_solutions=inf,
-        solution_approach=POPULATE, ko_cost=kocost, solver=curr_solver, compress=False)
+                                            sd_modules=modules_relaxed,
+                                            max_cost=3,
+                                            max_solutions=inf,
+                                            solution_approach=POPULATE,
+                                            ko_cost=kocost,
+                                            solver=curr_solver,
+                                            compress=False)
     assert len(sol_relaxed.reaction_sd) > 0, \
         "Relaxed DOUBLEOPT should find solutions with opt_tol=0.6"
     min_cost_relaxed = min(sum(abs(v) for v in s.values()) for s in sol_relaxed.reaction_sd)
@@ -348,9 +363,14 @@ def test_dump_preprocessed(model_small_example, tmp_path):
 
     # Step 1: Dump preprocessed data (should return early without solving)
     sol_dump = sd.compute_strain_designs(model_small_example,
-        sd_modules=modules, max_cost=inf, max_solutions=inf,
-        solution_approach='any', ki_cost=kicost, solver=solver,
-        compress=True, dump_preprocessed=dump_path)
+                                         sd_modules=modules,
+                                         max_cost=inf,
+                                         max_solutions=inf,
+                                         solution_approach='any',
+                                         ki_cost=kicost,
+                                         solver=solver,
+                                         compress=True,
+                                         dump_preprocessed=dump_path)
     import os
     assert os.path.exists(dump_path), "Dump file should exist"
 
@@ -387,9 +407,13 @@ def test_lazy_expansion(model_small_example):
     csd.LAZY_EXPANSION_THRESHOLD = 1  # Force lazy mode
     try:
         sol = sd.compute_strain_designs(model_small_example,
-            sd_modules=modules, max_cost=inf, max_solutions=inf,
-            solution_approach='any', ki_cost=kicost, solver=solver,
-            compress=True)
+                                        sd_modules=modules,
+                                        max_cost=inf,
+                                        max_solutions=inf,
+                                        solution_approach='any',
+                                        ki_cost=kicost,
+                                        solver=solver,
+                                        compress=True)
         if sol.is_lazy:
             assert sol.get_num_materialized() < sol.get_num_sols(), \
                 "Lazy mode: materialized < estimated total"
diff --git a/tests/test_07_compression.py b/tests/test_07_compression.py
index a60961b..d5d8aae 100644
--- a/tests/test_07_compression.py
+++ b/tests/test_07_compression.py
@@ -150,11 +150,9 @@ def test_fba_optimum_recovered_through_map():
     model = load_model("e_coli_core")
     cmp_maps = nt.compress_model(model)
     cmp_id, factor = _trace_lump(cmp_maps, biomass)
-    assert cmp_id in [r.id for r in model.reactions], (
-        f"compression map names {cmp_id}, which is not in the compressed model")
+    assert cmp_id in [r.id for r in model.reactions], (f"compression map names {cmp_id}, which is not in the compressed model")
     val = sd.fba(model, obj={cmp_id: 1}, obj_sense='maximize').objective_value
-    assert abs(factor * val - ref) < 1e-6, (
-        f"recovered optimum {factor * val} != uncompressed {ref}")
+    assert abs(factor * val - ref) < 1e-6, (f"recovered optimum {factor * val} != uncompressed {ref}")
 
 
 def test_cobra_optimize_after_compression():
@@ -188,8 +186,7 @@ def test_cobra_optimize_after_compression():
     sol = model_cmp.optimize()
     assert sol.status == 'optimal', f"Expected optimal solution, got {sol.status}"
     val_expanded = sol.objective_value * biomass_coeff
-    assert abs(val_orig - val_expanded) < 1e-6, (
-        f"Expanded objective mismatch: original={val_orig}, expanded={val_expanded}")
+    assert abs(val_orig - val_expanded) < 1e-6, (f"Expanded objective mismatch: original={val_orig}, expanded={val_expanded}")
 
 
 # =============================================================================
@@ -252,13 +249,12 @@ def test_mcs_e_coli_core():
     solver = SCIP if SCIP in strong_solvers else next(iter(strong_solvers))
     model = load_model('e_coli_core')
     modules = [sd.SDModule(model, SUPPRESS, constraints='BIOMASS_Ecoli_core_w_GAM >= 0.001')]
-    sols = sd.compute_strain_designs(model,
-                                     sd_modules=modules,
-                                     solution_approach=POPULATE,
-                                     max_cost=3,
-                                     gene_kos=True,
-                                     solver=solver,
-)
-    assert len(sols.reaction_sd) == 455, (
-        f"Expected 455 MCS for e_coli_core, got {len(sols.reaction_sd)}")
-
+    sols = sd.compute_strain_designs(
+        model,
+        sd_modules=modules,
+        solution_approach=POPULATE,
+        max_cost=3,
+        gene_kos=True,
+        solver=solver,
+    )
+    assert len(sols.reaction_sd) == 455, (f"Expected 455 MCS for e_coli_core, got {len(sols.reaction_sd)}")
diff --git a/tests/test_08_gene_ko.py b/tests/test_08_gene_ko.py
index b28d1f6..b1925b1 100644
--- a/tests/test_08_gene_ko.py
+++ b/tests/test_08_gene_ko.py
@@ -17,7 +17,9 @@ def wt_growth(ecoli_core):
 
 # ── gene_kos_to_constraints unit tests ───────────────────────────────
 
+
 class TestGeneKosToConstraints:
+
     def test_single_gene_by_id(self, ecoli_core):
         """b0727 (sucB) is in AND rule for AKGDH → single KO kills AKGDH."""
         c = gene_kos_to_constraints(ecoli_core, ['b0727'])
@@ -62,7 +64,9 @@ def test_empty_input(self, ecoli_core):
 
 # ── resolve_gene_constraints unit tests ──────────────────────────────
 
+
 class TestResolveGeneConstraints:
+
     def test_string_format(self, ecoli_core):
         """Gene KO as string 'b0727 = 0' should produce reaction constraints."""
         c = resolve_gene_constraints(ecoli_core, 'b0727 = 0')
@@ -116,7 +120,9 @@ def test_sd_grammar_mixed(self, ecoli_core):
 
 # ── FBA integration tests ────────────────────────────────────────────
 
+
 class TestFbaWithGeneKO:
+
     def test_nonessential_single_ko(self, ecoli_core, wt_growth):
         """b0727 (sucB) KO reduces growth slightly but is not lethal."""
         sol = sd.fba(ecoli_core, constraints='b0727 = 0')
@@ -160,7 +166,9 @@ def test_sd_grammar_minus_one_in_fba(self, ecoli_core, wt_growth):
 
 # ── FVA integration tests ────────────────────────────────────────────
 
+
 class TestFvaWithGeneKO:
+
     def test_knocked_reaction_fixed_at_zero(self, ecoli_core):
         """FVA after b0727 KO: AKGDH must have min=max=0."""
         fva_r = sd.fva(ecoli_core, constraints='b0727 = 0')
diff --git a/tests/test_09_performance.py b/tests/test_09_performance.py
index 715052d..9e6252c 100644
--- a/tests/test_09_performance.py
+++ b/tests/test_09_performance.py
@@ -88,7 +88,8 @@ def _git_sha() -> str:
     try:
         return subprocess.check_output(
             ["git", "rev-parse", "--short", "HEAD"],
-            cwd=TESTS_DIR, stderr=subprocess.DEVNULL,
+            cwd=TESTS_DIR,
+            stderr=subprocess.DEVNULL,
         ).decode().strip()
     except Exception:
         return "unknown"
@@ -107,8 +108,7 @@ def _solver_ver(solver: str) -> str:
     return "n/a"
 
 
-def record(name: str, solver: str, model_id: str,
-           elapsed: float, n_sol: int, status: str) -> None:
+def record(name: str, solver: str, model_id: str, elapsed: float, n_sol: int, status: str) -> None:
     entry = {
         "test": name,
         "solver": solver,
@@ -128,11 +128,11 @@ def record(name: str, solver: str, model_id: str,
 
 STRONG_SOLVERS = [s for s in [CPLEX, GUROBI] if s in sd.avail_solvers]
 
-
 # ---------------------------------------------------------------------------
 # Model fixtures  (session-scoped → loaded once per run)
 # ---------------------------------------------------------------------------
 
+
 @pytest.fixture(scope="session")
 def model_core():
     return load_model("e_coli_core")
@@ -152,6 +152,7 @@ def model_imlcore():
 # Quick suite — e_coli_core  (MCS correctness + speed, ~8 s / solver)
 # ===========================================================================
 
+
 @pytest.mark.parametrize("solver", STRONG_SOLVERS)
 @pytest.mark.timeout(180)
 def test_ecoli_core_mcs_455(solver, model_core):
@@ -164,18 +165,15 @@ def test_ecoli_core_mcs_455(solver, model_core):
     t0 = time.perf_counter()
     sol = sd.compute_strain_designs(
         m,
-        sd_modules=[sd.SDModule(m, SUPPRESS,
-                                constraints="BIOMASS_Ecoli_core_w_GAM >= 0.001")],
+        sd_modules=[sd.SDModule(m, SUPPRESS, constraints="BIOMASS_Ecoli_core_w_GAM >= 0.001")],
         solution_approach=POPULATE,
         max_cost=3,
         gene_kos=True,
         solver=solver,
     )
     elapsed = time.perf_counter() - t0
-    record("mcs_455", solver, "e_coli_core", elapsed,
-           len(sol.reaction_sd), sol.status)
-    assert len(sol.reaction_sd) == 455, (
-        f"[{solver}] Expected 455 MCS, got {len(sol.reaction_sd)}")
+    record("mcs_455", solver, "e_coli_core", elapsed, len(sol.reaction_sd), sol.status)
+    assert len(sol.reaction_sd) == 455, (f"[{solver}] Expected 455 MCS, got {len(sol.reaction_sd)}")
 
 
 # ===========================================================================
@@ -189,8 +187,17 @@ def test_ecoli_core_mcs_455(solver, model_core):
 
 # Shared KO / KI costs for model_weak_coupling  (from test_05_straindesign.py)
 _WEAK_KO = {
-    'r1': 1, 'r2': 1, 'r4': 1.1, 'r5': 0.75, 'r7': 0.8, 'r8': 1,
-    'r9': 1, 'r_S': 1.0, 'r_P': 1, 'r_BM': 1, 'r_Q': 1.5,
+    'r1': 1,
+    'r2': 1,
+    'r4': 1.1,
+    'r5': 0.75,
+    'r7': 0.8,
+    'r8': 1,
+    'r9': 1,
+    'r_S': 1.0,
+    'r_P': 1,
+    'r_BM': 1,
+    'r_Q': 1.5,
 }
 _WEAK_KI = {'r3': 0.6, 'r6': 1.0}
 _WEAK_REG = {'r6 >= 4.5': 1.2}
@@ -208,9 +215,7 @@ def test_weak_mcs_wgcp(solver, model_weak):
     """
     m = model_weak.copy()
     modules = [
-        sd.SDModule(m, SUPPRESS,
-                    inner_objective="r_BM",
-                    constraints=["r_P - 0.4 r_S <= 0", "r_S >= 0.1"]),
+        sd.SDModule(m, SUPPRESS, inner_objective="r_BM", constraints=["r_P - 0.4 r_S <= 0", "r_S >= 0.1"]),
         sd.SDModule(m, PROTECT, constraints=["r_BM >= 0.2"]),
     ]
     t0 = time.perf_counter()
@@ -222,14 +227,12 @@ def test_weak_mcs_wgcp(solver, model_weak):
         max_solutions=inf,
         ko_cost=_WEAK_KO,
         ki_cost=_WEAK_KI,
-        reg_cost=dict(_WEAK_REG),   # fresh copy: extend_model_regulatory mutates its arg
+        reg_cost=dict(_WEAK_REG),  # fresh copy: extend_model_regulatory mutates its arg
         solver=solver,
     )
     elapsed = time.perf_counter() - t0
-    record("mcs_wgcp", solver, "weak_coupling", elapsed,
-           len(sol.reaction_sd), sol.status)
-    assert len(sol.reaction_sd) == 3, (
-        f"[{solver}] Expected 3 wGCP MCS solutions, got {len(sol.reaction_sd)}")
+    record("mcs_wgcp", solver, "weak_coupling", elapsed, len(sol.reaction_sd), sol.status)
+    assert len(sol.reaction_sd) == 3, (f"[{solver}] Expected 3 wGCP MCS solutions, got {len(sol.reaction_sd)}")
 
 
 @pytest.mark.parametrize("solver", STRONG_SOLVERS)
@@ -242,10 +245,7 @@ def test_weak_optknock(solver, model_weak):
     """
     m = model_weak.copy()
     modules = [
-        sd.SDModule(m, OPTKNOCK,
-                    outer_objective="r_P",
-                    inner_objective="r_BM",
-                    constraints="r_BM >= 1"),
+        sd.SDModule(m, OPTKNOCK, outer_objective="r_P", inner_objective="r_BM", constraints="r_BM >= 1"),
     ]
     t0 = time.perf_counter()
     sol = sd.compute_strain_designs(
@@ -256,14 +256,12 @@ def test_weak_optknock(solver, model_weak):
         max_solutions=3,
         ko_cost=_WEAK_KO,
         ki_cost=_WEAK_KI,
-        reg_cost=dict(_WEAK_REG),   # fresh copy: extend_model_regulatory mutates its arg
+        reg_cost=dict(_WEAK_REG),  # fresh copy: extend_model_regulatory mutates its arg
         solver=solver,
     )
     elapsed = time.perf_counter() - t0
-    record("optknock", solver, "weak_coupling", elapsed,
-           len(sol.reaction_sd), sol.status)
-    assert len(sol.reaction_sd) == 3, (
-        f"[{solver}] Expected 3 OptKnock solutions, got {len(sol.reaction_sd)}")
+    record("optknock", solver, "weak_coupling", elapsed, len(sol.reaction_sd), sol.status)
+    assert len(sol.reaction_sd) == 3, (f"[{solver}] Expected 3 OptKnock solutions, got {len(sol.reaction_sd)}")
 
 
 @pytest.mark.parametrize("solver", STRONG_SOLVERS)
@@ -276,10 +274,9 @@ def test_weak_robustknock(solver, model_weak):
     """
     m = model_weak.copy()
     modules = [
-        sd.SDModule(m, ROBUSTKNOCK,
-                    outer_objective="r_P",
-                    inner_objective="r_BM",
-                    constraints=[[{"r_BM": 1.0}, ">=", 1.0]]),
+        sd.SDModule(m, ROBUSTKNOCK, outer_objective="r_P", inner_objective="r_BM", constraints=[[{
+            "r_BM": 1.0
+        }, ">=", 1.0]]),
     ]
     t0 = time.perf_counter()
     sol = sd.compute_strain_designs(
@@ -290,14 +287,12 @@ def test_weak_robustknock(solver, model_weak):
         max_solutions=3,
         ko_cost=_WEAK_KO,
         ki_cost=_WEAK_KI,
-        reg_cost=dict(_WEAK_REG),   # fresh copy: extend_model_regulatory mutates its arg
+        reg_cost=dict(_WEAK_REG),  # fresh copy: extend_model_regulatory mutates its arg
         solver=solver,
     )
     elapsed = time.perf_counter() - t0
-    record("robustknock", solver, "weak_coupling", elapsed,
-           len(sol.reaction_sd), sol.status)
-    assert len(sol.reaction_sd) >= 2, (
-        f"[{solver}] Expected ≥2 RobustKnock solutions, got {len(sol.reaction_sd)}")
+    record("robustknock", solver, "weak_coupling", elapsed, len(sol.reaction_sd), sol.status)
+    assert len(sol.reaction_sd) >= 2, (f"[{solver}] Expected ≥2 RobustKnock solutions, got {len(sol.reaction_sd)}")
 
 
 @pytest.mark.parametrize("solver", STRONG_SOLVERS)
@@ -310,10 +305,7 @@ def test_weak_optcouple(solver, model_weak):
     """
     m = model_weak.copy()
     modules = [
-        sd.SDModule(m, OPTCOUPLE,
-                    prod_id="r_P",
-                    inner_objective="r_BM",
-                    min_gcp=1.0),
+        sd.SDModule(m, OPTCOUPLE, prod_id="r_P", inner_objective="r_BM", min_gcp=1.0),
     ]
     t0 = time.perf_counter()
     sol = sd.compute_strain_designs(
@@ -324,20 +316,19 @@ def test_weak_optcouple(solver, model_weak):
         max_solutions=3,
         ko_cost=_WEAK_KO,
         ki_cost=_WEAK_KI,
-        reg_cost=dict(_WEAK_REG),   # fresh copy: extend_model_regulatory mutates its arg
+        reg_cost=dict(_WEAK_REG),  # fresh copy: extend_model_regulatory mutates its arg
         solver=solver,
     )
     elapsed = time.perf_counter() - t0
-    record("optcouple", solver, "weak_coupling", elapsed,
-           len(sol.reaction_sd), sol.status)
-    assert len(sol.reaction_sd) == 2, (
-        f"[{solver}] Expected 2 OptCouple solutions, got {len(sol.reaction_sd)}")
+    record("optcouple", solver, "weak_coupling", elapsed, len(sol.reaction_sd), sol.status)
+    assert len(sol.reaction_sd) == 2, (f"[{solver}] Expected 2 OptCouple solutions, got {len(sol.reaction_sd)}")
 
 
 # ===========================================================================
 # Standard suite (--medium): iMLcore gene-level MCS  (~47 s / solver each)
 # ===========================================================================
 
+
 @pytest.mark.medium
 @pytest.mark.parametrize("solver", STRONG_SOLVERS)
 @pytest.mark.timeout(90)
@@ -350,11 +341,8 @@ def test_imlcore_mcs_ethanol(solver, model_imlcore):
     """
     m = model_imlcore.copy()
     modules = [
-        sd.SDModule(m, SUPPRESS,
-                    constraints=["EX_etoh_e <= 1.0",
-                                 "BIOMASS_Ec_iML1515_core_75p37M >= 0.14"]),
-        sd.SDModule(m, PROTECT,
-                    constraints=["BIOMASS_Ec_iML1515_core_75p37M >= 0.15"]),
+        sd.SDModule(m, SUPPRESS, constraints=["EX_etoh_e <= 1.0", "BIOMASS_Ec_iML1515_core_75p37M >= 0.14"]),
+        sd.SDModule(m, PROTECT, constraints=["BIOMASS_Ec_iML1515_core_75p37M >= 0.15"]),
     ]
     t0 = time.perf_counter()
     sol = sd.compute_strain_designs(
@@ -366,10 +354,8 @@ def test_imlcore_mcs_ethanol(solver, model_imlcore):
         solver=solver,
     )
     elapsed = time.perf_counter() - t0
-    record("imlcore_ethanol", solver, "iMLcore", elapsed,
-           len(sol.reaction_sd), sol.status)
-    assert len(sol.reaction_sd) > 0, (
-        f"[{solver}] Expected ≥1 MCS for iMLcore ethanol scenario, got 0")
+    record("imlcore_ethanol", solver, "iMLcore", elapsed, len(sol.reaction_sd), sol.status)
+    assert len(sol.reaction_sd) > 0, (f"[{solver}] Expected ≥1 MCS for iMLcore ethanol scenario, got 0")
 
 
 @pytest.mark.medium
@@ -383,8 +369,7 @@ def test_imlcore_mcs_growth(solver, model_imlcore):
     """
     m = model_imlcore.copy()
     modules = [
-        sd.SDModule(m, SUPPRESS,
-                    constraints="BIOMASS_Ec_iML1515_core_75p37M >= 0.001"),
+        sd.SDModule(m, SUPPRESS, constraints="BIOMASS_Ec_iML1515_core_75p37M >= 0.001"),
     ]
     t0 = time.perf_counter()
     sol = sd.compute_strain_designs(
@@ -396,16 +381,15 @@ def test_imlcore_mcs_growth(solver, model_imlcore):
         solver=solver,
     )
     elapsed = time.perf_counter() - t0
-    record("imlcore_growth", solver, "iMLcore", elapsed,
-           len(sol.reaction_sd), sol.status)
-    assert len(sol.reaction_sd) > 0, (
-        f"[{solver}] Expected ≥1 MCS for iMLcore growth scenario, got 0")
+    record("imlcore_growth", solver, "iMLcore", elapsed, len(sol.reaction_sd), sol.status)
+    assert len(sol.reaction_sd) > 0, (f"[{solver}] Expected ≥1 MCS for iMLcore growth scenario, got 0")
 
 
 # ===========================================================================
 # Large suite (--large): iML1515 — known answer 393
 # ===========================================================================
 
+
 @pytest.mark.large
 @pytest.mark.parametrize("solver", [GUROBI])  # CPLEX is very slow on this one; skip to save time
 def test_iml1515_mcs_393(solver):
@@ -421,24 +405,22 @@ def test_iml1515_mcs_393(solver):
     t0 = time.perf_counter()
     sol = sd.compute_strain_designs(
         m,
-        sd_modules=[sd.SDModule(m, SUPPRESS,
-                                constraints="BIOMASS_Ec_iML1515_core_75p37M >= 0.001")],
+        sd_modules=[sd.SDModule(m, SUPPRESS, constraints="BIOMASS_Ec_iML1515_core_75p37M >= 0.001")],
         solution_approach=POPULATE,
         max_cost=3,
         gene_kos=True,
         solver=solver,
     )
     elapsed = time.perf_counter() - t0
-    record("iml1515_393", solver, "iML1515", elapsed,
-           len(sol.reaction_sd), sol.status)
-    assert len(sol.reaction_sd) == 393, (
-        f"[{solver}] Expected 393 MCS for iML1515, got {len(sol.reaction_sd)}")
+    record("iml1515_393", solver, "iML1515", elapsed, len(sol.reaction_sd), sol.status)
+    assert len(sol.reaction_sd) == 393, (f"[{solver}] Expected 393 MCS for iML1515, got {len(sol.reaction_sd)}")
 
 
 # ===========================================================================
 # Session teardown: write JSON + print comparison table
 # ===========================================================================
 
+
 @pytest.fixture(scope="session", autouse=True)
 def _write_results():
     yield
@@ -449,7 +431,9 @@ def _write_results():
         "git_sha": _git_sha(),
         "platform": platform.platform(),
         "python": platform.python_version(),
-        "solver_versions": {s: _solver_ver(s) for s in STRONG_SOLVERS},
+        "solver_versions": {
+            s: _solver_ver(s) for s in STRONG_SOLVERS
+        },
         "results": _RESULTS,
     }
     out_file = RESULTS_DIR / f"{_SESSION_TS}.json"
@@ -460,14 +444,12 @@ def _write_results():
     tests_ran = sorted({(r["test"], r["model"]) for r in _RESULTS})
     col = 16
 
-    header = (f"{'Test':<22} {'Model':<16}"
-              + "".join(f"  {s:>{col}}" for s in solvers_ran))
+    header = (f"{'Test':<22} {'Model':<16}" + "".join(f"  {s:>{col}}" for s in solvers_ran))
     print("=== Solver Performance Summary ===")
     print(header)
     print("-" * len(header))
     for (test, model) in tests_ran:
-        row = {r["solver"]: r for r in _RESULTS
-               if r["test"] == test and r["model"] == model}
+        row = {r["solver"]: r for r in _RESULTS if r["test"] == test and r["model"] == model}
         line = f"{test:<22} {model:<16}"
         for s in solvers_ran:
             if s in row:
diff --git a/tests/test_10_gene_design_validity.py b/tests/test_10_gene_design_validity.py
index 8f144ca..25fd85e 100644
--- a/tests/test_10_gene_design_validity.py
+++ b/tests/test_10_gene_design_validity.py
@@ -9,8 +9,7 @@
 from math import inf
 from cobra.io import read_sbml_model
 import straindesign as sd
-from straindesign.names import (MODULES, MAX_COST, MAX_SOLUTIONS, SOLUTION_APPROACH,
-                                KOCOST, GKOCOST, SOLVER)
+from straindesign.names import (MODULES, MAX_COST, MAX_SOLUTIONS, SOLUTION_APPROACH, KOCOST, GKOCOST, SOLVER)
 from straindesign import SUPPRESS, PROTECT
 
 TOL = 1e-6
@@ -37,11 +36,25 @@ def _apply_gene_design(model, design):
 
 
 def _gpr_mcs_setup(model, gko_cost, solver, approach):
-    modules = [sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]),
-               sd.SDModule(model, PROTECT, constraints=[[{'r_bm': 1.0}, '>=', 1.0]])]
-    return {MODULES: modules, MAX_COST: 3, MAX_SOLUTIONS: inf, SOLUTION_APPROACH: approach,
-            KOCOST: {'rs_up': 1.0, 'rd_ex': 1.0, 'rp_ex': 1.1},
-            GKOCOST: gko_cost, SOLVER: solver}
+    modules = [
+        sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]),
+        sd.SDModule(model, PROTECT, constraints=[[{
+            'r_bm': 1.0
+        }, '>=', 1.0]])
+    ]
+    return {
+        MODULES: modules,
+        MAX_COST: 3,
+        MAX_SOLUTIONS: inf,
+        SOLUTION_APPROACH: approach,
+        KOCOST: {
+            'rs_up': 1.0,
+            'rd_ex': 1.0,
+            'rp_ex': 1.1
+        },
+        GKOCOST: gko_cost,
+        SOLVER: solver
+    }
 
 
 @pytest.fixture
diff --git a/tests/test_11_gene_regulatory_compression.py b/tests/test_11_gene_regulatory_compression.py
index 54fc60e..8049c45 100644
--- a/tests/test_11_gene_regulatory_compression.py
+++ b/tests/test_11_gene_regulatory_compression.py
@@ -34,12 +34,23 @@ def _coupled_chain():
     m = Model("chain")
     mets = {x: Metabolite(x) for x in "ABCD"}
     m.add_metabolites(list(mets.values()))
-    vin = Reaction("vin"); vin.add_metabolites({mets["A"]: 1}); vin.bounds = (0, 10)
-    r1 = Reaction("r1"); r1.add_metabolites({mets["A"]: -1, mets["B"]: 1}); r1.bounds = (0, 1000)
-    r2 = Reaction("r2"); r2.add_metabolites({mets["B"]: -1, mets["C"]: 1}); r2.bounds = (0, 1000)
-    r3 = Reaction("r3"); r3.add_metabolites({mets["C"]: -1, mets["D"]: 1}); r3.bounds = (0, 1000)
-    vout = Reaction("vout"); vout.add_metabolites({mets["D"]: -1}); vout.bounds = (0, 1000)
-    m.add_reactions([vin, r1, r2, r3, vout]); m.objective = "vout"
+    vin = Reaction("vin")
+    vin.add_metabolites({mets["A"]: 1})
+    vin.bounds = (0, 10)
+    r1 = Reaction("r1")
+    r1.add_metabolites({mets["A"]: -1, mets["B"]: 1})
+    r1.bounds = (0, 1000)
+    r2 = Reaction("r2")
+    r2.add_metabolites({mets["B"]: -1, mets["C"]: 1})
+    r2.bounds = (0, 1000)
+    r3 = Reaction("r3")
+    r3.add_metabolites({mets["C"]: -1, mets["D"]: 1})
+    r3.bounds = (0, 1000)
+    vout = Reaction("vout")
+    vout.add_metabolites({mets["D"]: -1})
+    vout.bounds = (0, 1000)
+    m.add_reactions([vin, r1, r2, r3, vout])
+    m.objective = "vout"
     return m
 
 
@@ -48,7 +59,8 @@ def test_coupled_exemption_keeps_reaction_and_merges_rest():
     and flux space is preserved."""
     base = _fba(_coupled_chain(), {"vout": 1})
 
-    m_full = _coupled_chain(); compress_model(m_full)
+    m_full = _coupled_chain()
+    compress_model(m_full)
     assert len(m_full.reactions) == 1, "unprotected coupled chain should fully collapse"
 
     m_prot = _coupled_chain()
@@ -62,20 +74,31 @@ def test_coupled_exemption_keeps_reaction_and_merges_rest():
 
 def _gene_model():
     """g1 controls r1 (A->2B) and r2 (B->C), which are coupled 1:2."""
-    m = Model("g"); A, B, C = Metabolite("A"), Metabolite("B"), Metabolite("C")
+    m = Model("g")
+    A, B, C = Metabolite("A"), Metabolite("B"), Metabolite("C")
     m.add_metabolites([A, B, C])
-    vS = Reaction("vS"); vS.add_metabolites({A: 1}); vS.bounds = (0, 10)
-    r1 = Reaction("r1"); r1.add_metabolites({A: -1, B: 2}); r1.bounds = (0, 1000); r1.gene_reaction_rule = "g1"
-    r2 = Reaction("r2"); r2.add_metabolites({B: -1, C: 1}); r2.bounds = (0, 1000); r2.gene_reaction_rule = "g1"
-    vC = Reaction("vC"); vC.add_metabolites({C: -1}); vC.bounds = (0, 1000)
-    m.add_reactions([vS, r1, r2, vC]); m.objective = "vC"
+    vS = Reaction("vS")
+    vS.add_metabolites({A: 1})
+    vS.bounds = (0, 10)
+    r1 = Reaction("r1")
+    r1.add_metabolites({A: -1, B: 2})
+    r1.bounds = (0, 1000)
+    r1.gene_reaction_rule = "g1"
+    r2 = Reaction("r2")
+    r2.add_metabolites({B: -1, C: 1})
+    r2.bounds = (0, 1000)
+    r2.gene_reaction_rule = "g1"
+    vC = Reaction("vC")
+    vC.add_metabolites({C: -1})
+    vC.bounds = (0, 1000)
+    m.add_reactions([vS, r1, r2, vC])
+    m.objective = "vC"
     return m
 
 
 def _max_C_under_g1_le_1(protect):
     m = _gene_model()
-    cmap = compress_model(m, no_coupled_compress_reacs=({"r1", "r2"} if protect else set()),
-                          propagate_gpr=True)
+    cmap = compress_model(m, no_coupled_compress_reacs=({"r1", "r2"} if protect else set()), propagate_gpr=True)
     extend_model_gpr(m, use_names=False)
     m.reactions.g1.upper_bound = 1.0
     return _fba(m, _remap_obj({"vC": 1}, cmap))
@@ -84,7 +107,9 @@ def _max_C_under_g1_le_1(protect):
 def test_gene_regulatory_multiplicity_preserved_under_compression():
     """With g1 controlling coupled r1,r2, the regulatory bound g1<=1 must give the same max product
     compressed (with protection) as uncompressed; without protection it is wrong (3x)."""
-    mu = _gene_model(); extend_model_gpr(mu, use_names=False); mu.reactions.g1.upper_bound = 1.0
+    mu = _gene_model()
+    extend_model_gpr(mu, use_names=False)
+    mu.reactions.g1.upper_bound = 1.0
     ref = _fba(mu, {"vC": 1})
     assert abs(ref - 2.0 / 3.0) < 1e-4, "reference should be 2/3"
 
diff --git a/tests/test_12_save_load_embed_model.py b/tests/test_12_save_load_embed_model.py
index 176ace6..e3ef907 100644
--- a/tests/test_12_save_load_embed_model.py
+++ b/tests/test_12_save_load_embed_model.py
@@ -12,8 +12,7 @@
 # imported into the package), so reach the module via sys.modules.
 csd = sys.modules["straindesign.compute_strain_designs"]
 from straindesign import SUPPRESS, PROTECT
-from straindesign.names import (MODULES, MAX_COST, MAX_SOLUTIONS,
-                                SOLUTION_APPROACH, KOCOST, GKOCOST, SOLVER, SEED)
+from straindesign.names import (MODULES, MAX_COST, MAX_SOLUTIONS, SOLUTION_APPROACH, KOCOST, GKOCOST, SOLVER, SEED)
 
 GPR = os.path.join(os.path.dirname(__file__), "model_gpr.xml")
 TOL = 1e-6
@@ -36,12 +35,15 @@ def _compute(model, threshold=None):
     if threshold is not None:
         old, csd.LAZY_EXPANSION_THRESHOLD = csd.LAZY_EXPANSION_THRESHOLD, threshold
     try:
-        return sd.compute_strain_designs(
-            model,
-            sd_modules=[sd.SDModule(model, SUPPRESS,
-                                    constraints="Biomass_Ecoli_core >= 0.1")],
-            gene_kos=True, max_cost=3, max_solutions=8,
-            solution_approach="any", solver=_solver(), seed=1, compress=True)
+        return sd.compute_strain_designs(model,
+                                         sd_modules=[sd.SDModule(model, SUPPRESS, constraints="Biomass_Ecoli_core >= 0.1")],
+                                         gene_kos=True,
+                                         max_cost=3,
+                                         max_solutions=8,
+                                         solution_approach="any",
+                                         solver=_solver(),
+                                         seed=1,
+                                         compress=True)
     finally:
         if threshold is not None:
             csd.LAZY_EXPANSION_THRESHOLD = old
@@ -53,7 +55,7 @@ def test_embed_and_restore_roundtrip(model, tmp_path):
     ref_rsd = sols.get_reaction_sd()
     ref_gsd = sols.get_gene_sd()
     f = str(tmp_path / "sd.pkl")
-    sols.save(f)                                   # embed_model=True by default
+    sols.save(f)  # embed_model=True by default
 
     # default load: model NOT rebuilt, but solutions are there
     loaded = sd.SDSolutions.load(f)
@@ -67,14 +69,13 @@ def test_embed_and_restore_roundtrip(model, tmp_path):
     assert len(restored._model.reactions) == len(model.reactions)
     assert len(restored._model.genes) == len(model.genes)
     r = "AKGDH"
-    assert (restored._model.reactions.get_by_id(r).gene_reaction_rule
-            == model.reactions.get_by_id(r).gene_reaction_rule)
+    assert (restored._model.reactions.get_by_id(r).gene_reaction_rule == model.reactions.get_by_id(r).gene_reaction_rule)
 
 
 def test_explicit_model_takes_precedence(model, tmp_path):
     sols = _compute(model)
     f = str(tmp_path / "sd.pkl")
-    sols.save(f, embed_model=False)                # leaner file, no snapshot
+    sols.save(f, embed_model=False)  # leaner file, no snapshot
     loaded = sd.SDSolutions.load(f)
     assert loaded._embedded_model_dict is None
     # nothing to restore from, stays model-less
@@ -88,19 +89,19 @@ def test_live_model_and_solver_never_pickled(model, tmp_path):
     f = str(tmp_path / "sd.pkl")
     sols.save(f)
     with open(f, "rb") as fh:
-        raw = pickle.load(fh)               # must not raise (no live solver in pickle)
+        raw = pickle.load(fh)  # must not raise (no live solver in pickle)
     assert raw._model is None
     assert raw._embedded_model_dict is not None
 
 
 # ── lazy round-trip: save must not force-expand; restore then expand ─────
 def test_lazy_save_no_expand_then_restore_and_expand(model, tmp_path):
-    sols = _compute(model, threshold=1)            # force lazy expansion
+    sols = _compute(model, threshold=1)  # force lazy expansion
     assert sols.is_lazy
     materialized = sols.get_num_materialized()
     f = str(tmp_path / "lazy.pkl")
-    sols.save(f)                                   # must NOT expand_all (no hang)
-    assert sols.is_lazy                            # still lazy after save
+    sols.save(f)  # must NOT expand_all (no hang)
+    assert sols.is_lazy  # still lazy after save
 
     # default load: lazy, materialized reps available, expand errors clearly
     loaded = sd.SDSolutions.load(f)
@@ -123,6 +124,7 @@ def test_lazy_save_no_expand_then_restore_and_expand(model, tmp_path):
 #    on compressed representatives AND after expansion, then the whole
 #    computation is REPRODUCED from the artifact alone (embedded model + sd_setup).
 
+
 def _max_flux(model, rid):
     with model:
         model.objective = rid
@@ -157,18 +159,32 @@ def _validate_designs(reaction_designs):
 
 def _keyset(reaction_designs):
     """Canonical set of designs = set of frozensets of knocked-out reaction ids."""
-    return {frozenset(k for k, v in d.items() if v in (-1, -1.0, False))
-            for d in reaction_designs}
+    return {frozenset(k for k, v in d.items() if v in (-1, -1.0, False)) for d in reaction_designs}
 
 
 def _combined_setup(model, solver, approach):
-    modules = [sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]),
-               sd.SDModule(model, PROTECT, constraints=[[{"r_bm": 1.0}, ">=", 1.0]])]
-    return {MODULES: modules, MAX_COST: 3, MAX_SOLUTIONS: inf,
-            SOLUTION_APPROACH: approach,
-            KOCOST: {"rs_up": 1.0, "rd_ex": 1.0, "rp_ex": 1.1},         # reaction KOs
-            GKOCOST: {g.id: 1.0 for g in model.genes},                 # gene KOs
-            SOLVER: solver, SEED: 7}
+    modules = [
+        sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]),
+        sd.SDModule(model, PROTECT, constraints=[[{
+            "r_bm": 1.0
+        }, ">=", 1.0]])
+    ]
+    return {
+        MODULES: modules,
+        MAX_COST: 3,
+        MAX_SOLUTIONS: inf,
+        SOLUTION_APPROACH: approach,
+        KOCOST: {
+            "rs_up": 1.0,
+            "rd_ex": 1.0,
+            "rp_ex": 1.1
+        },  # reaction KOs
+        GKOCOST: {
+            g.id: 1.0 for g in model.genes
+        },  # gene KOs
+        SOLVER: solver,
+        SEED: 7
+    }
 
 
 @pytest.mark.parametrize("approach", ["any", "best", "populate"])
@@ -192,12 +208,12 @@ def test_full_reproducibility_roundtrip(tmp_path, approach, force_lazy):
     ref_gene_sd = sol.get_gene_sd()
     ref_keys = _keyset(sol.get_reaction_sd())
     assert ref_gene_sd
-    _validate_designs(sol.get_reaction_sd())            # sanity on the fresh result
+    _validate_designs(sol.get_reaction_sd())  # sanity on the fresh result
 
     # save a self-contained artifact (embed_model=True by default)
     f = str(tmp_path / "repro.pkl")
     sol.save(f)
-    assert sol.is_lazy == force_lazy                    # save never expanded a lazy result
+    assert sol.is_lazy == force_lazy  # save never expanded a lazy result
 
     # (A) reload WITHOUT a model: designs preserved; validate the COMPRESSED
     #     representatives directly via SUPPRESS/PROTECT FBA
@@ -210,11 +226,11 @@ def test_full_reproducibility_roundtrip(tmp_path, approach, force_lazy):
     #     expanded design set via FBA
     restored = sd.SDSolutions.load(f, model=True, cmp_model=True)
     assert restored._model is not None
-    restored.expand_all()                               # no-op if already non-lazy
+    restored.expand_all()  # no-op if already non-lazy
     assert not restored.is_lazy
     _validate_designs(restored.get_reaction_sd())
-    expanded_keys = _keyset(restored.get_reaction_sd())   # full design set
-    assert expanded_keys >= ref_keys                      # reps ⊆ full set
+    expanded_keys = _keyset(restored.get_reaction_sd())  # full design set
+    assert expanded_keys >= ref_keys  # reps ⊆ full set
 
     # (C) REPRODUCE the computation from the artifact ALONE: the embedded
     #     (uncompressed) model + the stored sd_setup must regenerate the exact
@@ -229,11 +245,27 @@ def test_full_reproducibility_roundtrip(tmp_path, approach, force_lazy):
 
 
 def _combined_gpr(model, solver, approach="any"):
-    modules = [sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]),
-               sd.SDModule(model, PROTECT, constraints=[[{"r_bm": 1.0}, ">=", 1.0]])]
-    return {MODULES: modules, MAX_COST: 3, SOLUTION_APPROACH: approach,
-            KOCOST: {"rs_up": 1.0, "rd_ex": 1.0, "rp_ex": 1.1},
-            GKOCOST: {g.id: 1.0 for g in model.genes}, SOLVER: solver, SEED: 7}
+    modules = [
+        sd.SDModule(model, SUPPRESS, constraints=["1.0 rd_ex >= 1.0 "]),
+        sd.SDModule(model, PROTECT, constraints=[[{
+            "r_bm": 1.0
+        }, ">=", 1.0]])
+    ]
+    return {
+        MODULES: modules,
+        MAX_COST: 3,
+        SOLUTION_APPROACH: approach,
+        KOCOST: {
+            "rs_up": 1.0,
+            "rd_ex": 1.0,
+            "rp_ex": 1.1
+        },
+        GKOCOST: {
+            g.id: 1.0 for g in model.genes
+        },
+        SOLVER: solver,
+        SEED: 7
+    }
 
 
 def test_networktools_model_dict_rational_exact():
@@ -243,18 +275,21 @@ def test_networktools_model_dict_rational_exact():
     from straindesign.networktools import model_to_dict, model_from_dict
     import json, math
     m = Model("rat")
-    A = Metabolite("A_c", compartment="c"); B = Metabolite("B_c", compartment="c")
-    r = Reaction("r1"); m.add_reactions([r])
+    A = Metabolite("A_c", compartment="c")
+    B = Metabolite("B_c", compartment="c")
+    r = Reaction("r1")
+    m.add_reactions([r])
     r.add_metabolites({A: Fraction(1, 3), B: Fraction(-7, 3)})
-    r.lower_bound = Fraction(1, 3); r.upper_bound = float("inf")
+    r.lower_bound = Fraction(1, 3)
+    r.upper_bound = float("inf")
     d = model_to_dict(m)
-    rt = model_from_dict(json.loads(json.dumps(d)))     # force through JSON primitives
+    rt = model_from_dict(json.loads(json.dumps(d)))  # force through JSON primitives
     lb = rt.reactions.r1.lower_bound
     coeff = rt.reactions.r1.metabolites[rt.metabolites.A_c]
-    assert isinstance(lb, Fraction) and lb == Fraction(1, 3)          # exact, still rational
+    assert isinstance(lb, Fraction) and lb == Fraction(1, 3)  # exact, still rational
     assert isinstance(coeff, Fraction) and coeff == Fraction(1, 3)
-    assert float(Fraction(1, 3)) != Fraction(1, 3)                    # sanity: float would differ
-    assert math.isinf(rt.reactions.r1.upper_bound)                   # inf round-trips
+    assert float(Fraction(1, 3)) != Fraction(1, 3)  # sanity: float would differ
+    assert math.isinf(rt.reactions.r1.upper_bound)  # inf round-trips
 
 
 def test_networktools_model_dict_float_matches_cobra(model):
@@ -263,7 +298,7 @@ def test_networktools_model_dict_float_matches_cobra(model):
     from straindesign.networktools import model_to_dict, model_from_dict
     import json
     d = model_to_dict(model)
-    json.dumps(d)                                        # must be JSON-clean
+    json.dumps(d)  # must be JSON-clean
     rt = model_from_dict(d)
     assert len(rt.reactions) == len(model.reactions)
     assert len(rt.genes) == len(model.genes)
@@ -278,12 +313,12 @@ def test_compressed_model_embedded_by_default_restore_optin(tmp_path):
     assert sol._cmp_model is not None and len(sol._cmp_model.reactions) < len(m.reactions)
     cm0 = sol._cmp_model
     f = str(tmp_path / "sd.pkl")
-    sol.save(f)                                          # embeds BOTH by default
+    sol.save(f)  # embeds BOTH by default
 
     # default load restores neither model
     o0 = sd.SDSolutions.load(f)
     assert o0.get_model() is None and o0.get_compressed_model() is None
-    assert o0._embedded_cmp_model_dict is not None       # ...but it IS embedded
+    assert o0._embedded_cmp_model_dict is not None  # ...but it IS embedded
 
     # model=True, cmp_model=True rebuild both; compressed model is exact + usable
     o = sd.SDSolutions.load(f, model=True, cmp_model=True)
@@ -296,18 +331,19 @@ def test_compressed_model_embedded_by_default_restore_optin(tmp_path):
         assert r0.lower_bound == r1.lower_bound and r0.upper_bound == r1.upper_bound
         for a in ("lower_bound", "upper_bound"):
             if isinstance(getattr(r0, a), Fraction):
-                assert isinstance(getattr(r1, a), Fraction)   # rational not float-ified
+                assert isinstance(getattr(r1, a), Fraction)  # rational not float-ified
         for mt, c0 in r0.metabolites.items():
             c1 = r1.metabolites[cm.metabolites.get_by_id(mt.id)]
             assert c0 == c1 and (not isinstance(c0, Fraction) or isinstance(c1, Fraction))
     assert any(isinstance(b, Fraction) for r in cm.reactions for b in (r.lower_bound, r.upper_bound))
     cm.objective = cm.reactions[0].id
-    assert cm.optimize().status == "optimal"             # a fully working model
+    assert cm.optimize().status == "optimal"  # a fully working model
 
     # explicit overrides take precedence
     assert sd.SDSolutions.load(f, cmp_model=cm0).get_compressed_model() is cm0
     # embed_model=False -> no compressed snapshot
-    f2 = str(tmp_path / "lean.pkl"); sol.save(f2, embed_model=False)
+    f2 = str(tmp_path / "lean.pkl")
+    sol.save(f2, embed_model=False)
     assert sd.SDSolutions.load(f2, model=True, cmp_model=True).get_compressed_model() is None
 
 
@@ -317,14 +353,15 @@ def test_compressed_solutions_analyzable_in_restored_cmp_model(tmp_path):
     solver = _solver()
     m = read_sbml_model(GPR)
     sol = sd.compute_strain_designs(m, sd_setup=_combined_gpr(m, solver))
-    f = str(tmp_path / "sd.pkl"); sol.save(f)
+    f = str(tmp_path / "sd.pkl")
+    sol.save(f)
     o = sd.SDSolutions.load(f, model=True, cmp_model=True)
     cm = o.get_compressed_model()
     cm_rxn_ids = {r.id for r in cm.reactions}
     non_empty = [cs for cs in o.compressed_sd if cs]
     assert non_empty, "expected at least one non-empty compressed solution"
     for cs in non_empty:
-        assert set(cs).issubset(cm_rxn_ids)              # analysable in the compressed model
+        assert set(cs).issubset(cm_rxn_ids)  # analysable in the compressed model
     # apply one compressed KO set and confirm FBA still runs in the small model
     cs = non_empty[0]
     with cm:
@@ -339,8 +376,12 @@ def test_networktools_preserves_objective_direction():
     """networktools serializers keep the objective SENSE, which cobra drops."""
     from cobra import Model, Metabolite, Reaction
     from straindesign.networktools import model_to_dict, model_from_dict
-    m = Model("obj"); A = Metabolite("A_c", compartment="c")
-    r = Reaction("r1", lower_bound=0, upper_bound=10); m.add_reactions([r])
-    r.add_metabolites({A: 1.0}); m.objective = "r1"; m.objective.direction = "min"
+    m = Model("obj")
+    A = Metabolite("A_c", compartment="c")
+    r = Reaction("r1", lower_bound=0, upper_bound=10)
+    m.add_reactions([r])
+    r.add_metabolites({A: 1.0})
+    m.objective = "r1"
+    m.objective.direction = "min"
     rt = model_from_dict(model_to_dict(m))
     assert rt.objective.direction == "min"