OptKnock: Bilevel Optimization for Metabolic Engineering

Table of Contents

1. Biological Background
2. The Bilevel Problem Structure
3. From Bilevel to Single-Level via LP Duality
4. The Complete MILP Formulation
5. Objective Function Variants
6. Greedy Baseline
7. Batch Production Metric
8. Implementation Overview
9. Experimental Results

---

1. Biological Background

OptKnock is an algorithm for metabolic engineering: given a genome-scale metabolic model of a microorganism, find a small set of gene knockouts that force the cell to overproduce a desired chemical — such as ethanol or succinate — as a consequence of its own growth optimization.

The key biological insight is growth-coupling. A cell's primary goal is to maximize its own growth rate. If you eliminate competing pathways, the cell may have no viable route to regenerate the energy carriers (ATP, NADH) it needs to grow other than routing flux through the target product pathway. The target becomes a growth-coupled byproduct: it is produced automatically whenever the cell divides.

This is fundamentally different from expressing a heterologous pathway or overexpressing an enzyme. A growth-coupled strain is evolutionarily stable — any mutant that stops producing the target would also grow more slowly and be outcompeted. The engineering problem is: which genes to delete?

The E. coli core model (e_coli_core.json) used in this implementation is the canonical COBRA metabolic reconstruction of Escherichia coli's central carbon metabolism — 95 reactions and 72 metabolites. The stoichiometric matrix $S$ encodes every reaction's mass balance.

---

2. The Bilevel Problem Structure

OptKnock has an inherently nested structure: two levels of optimization with opposing players.

The Inner Problem — What the Cell Does

Given a fixed set of gene knockouts (some reactions disabled), the cell distributes metabolic flux to maximize its own biomass production rate. This is standard Flux Balance Analysis (FBA), a well-established LP:

$$ \begin{aligned} \max_{v} \quad & v_{\text{biomass}} \\ \text{s.t.} \quad & S v = 0 && \text{(steady-state mass balance, for all metabolites)} \\ & \ell_j, r_j \leq v_j \leq u_j, r_j && \text{(flux bounds, scaled by knockout status)} \end{aligned} $$

where:

• $v_j$ is the flux through reaction $j$ (units: mmol/gDW/h)
• $S \in \mathbb{R}^{m \times n}$ is the stoichiometric matrix ($S_{ij}$ = net stoichiometry of metabolite $i$ in reaction $j$; positive = produced, negative = consumed)
• $\ell_j,, u_j$ are the lower and upper flux bounds on reaction $j$
• $r_j \in {0, 1}$ is the knockout indicator: $r_j = 0$ forces $v_j = 0$ (reaction silenced)

The mass balance constraint $Sv = 0$ enforces steady state: at every metabolite node, the net rate of production equals the net rate of consumption. This is valid because metabolic fluxes operate on a timescale much faster than cell growth — the metabolite concentrations are essentially at steady state during exponential growth.

The Outer Problem — What the Engineer Chooses

The engineer controls which genes to delete. The objective is to find at most $k$ gene knockouts such that the inner-optimal flux through the target reaction is maximized:

$$ \begin{aligned} \max_{r,,v} \quad & v_{\text{target}} \\ \text{s.t.} \quad & v \text{ solves the inner FBA given } r \\ & \sum_i (1 - y_i) \leq k && \text{(at most } k \text{ gene knockouts)} \\ & \text{GPR rules link } y \text{ to } r && \text{(gene-reaction logic)} \end{aligned} $$

The constraint "$v$ solves the inner FBA given $r$" is what makes this a bilevel program. The outer decision variable $r$ (or equivalently, the gene presence vector $y$) appears in the inner problem's constraint set. Solving one problem as a constraint inside another is not something a standard LP or MILP solver can do directly.

---

3. From Bilevel to Single-Level via LP Duality

The fundamental insight of OptKnock is that the inner FBA is a linear program, and every LP has an equivalent dual LP. Replacing the inner optimization by its optimality conditions — using LP duality theory — collapses the bilevel problem into a single-level MILP.

Step 1: Reversible Reaction Splitting

Before taking the dual, all reactions are made non-negative by splitting each reversible reaction ($\ell_j < 0$) into a forward part and a reverse part:

• Forward: $v^+j \in [0, u_j]$, stoichiometry column $+S{:,j}$
• Reverse: $v^-j \in [0, -\ell_j]$, stoichiometry column $-S{:,j}$

After splitting, every flux variable is non-negative ($\ell_j = 0$). This is essential because the dual variables $\mu^-_j$ (lower-bound duals) vanish when $\ell_j = 0$, greatly simplifying the formulation.

Step 2: Writing the Dual of the Inner LP

The inner FBA (after splitting) has the standard primal form $\max{c^\top v : Sv = 0,, 0 \leq v \leq u \cdot r}$ where $c_j = 1$ for the biomass reaction and $c_j = 0$ elsewhere. Taking the LP dual introduces one dual variable per constraint:

Primal constraint	Dual variable	Interpretation
$Sv = 0$ (equality, $m$ constraints)	$\lambda_i \in \mathbb{R}$ (free)	Shadow price of metabolite $i$: marginal increase in biomass if one unit of metabolite $i$ were added to the system
$v_j \leq u_j r_j$ ($n$ upper-bound constraints)	$\mu^+_j \geq 0$	Reduced cost: positive when reaction $j$ is operating at full capacity
$v_j \geq 0$ ($n$ lower-bound constraints)	$\mu^-_j \geq 0$	Reduced cost: positive when reaction $j$ carries zero flux

The dual LP is:

$$ \min_{\lambda,, \mu^+,, \mu^-} \quad \sum_j u_j r_j, \mu^+_j \quad \text{s.t.} \quad S^\top \lambda + \mu^+ - \mu^- = c,\quad \mu^+, \mu^- \geq 0 $$

Step 3: KKT Stationarity

The KKT stationarity condition for the primal maximization states that at any optimal primal-dual pair, the gradient of the Lagrangian with respect to $v$ must be zero. For each reaction $j$:

$$ \underbrace{\sum_i S_{ij}, \lambda_i}_{\text{metabolite contribution}} ;+; \mu^+_j ;-; \mu^-_j ;=; c_j \quad \forall j $$

This is the dual feasibility constraint. It holds as a strict equality for every reaction — both active ($r_j = 1$) and knocked-out ($r_j = 0$). Geometrically, the metabolite shadow prices (weighted by each reaction's stoichiometry) must exactly balance the objective coefficient $c_j$, with any slack absorbed by the bound dual variables.

Step 4: Strong Duality as the Optimality Certificate

For any LP, strong duality holds at the optimum: the primal and dual objectives are equal. For the inner FBA (after splitting, so $\ell_j = 0$):

$$ \underbrace{v_{\text{biomass}}}_{\text{primal objective}} ;=; \underbrace{\sum_j u_j, \mu^+_j}_{\text{dual objective at } \ell_j=0} $$

This single equality, combined with primal feasibility ($Sv = 0$, $0 \leq v \leq u \cdot r$) and dual feasibility (stationarity above), is a complete certificate that the inner LP is solved optimally. No separate inner LP solve is needed — we just enforce these constraints in the outer problem.

Why does this work? The Karush-Kuhn-Tucker theorem states that for a linear program, a primal-dual pair $(v^, \lambda^, \mu^{+}, \mu^{-})$ is optimal if and only if:

1. Primal feasibility: $Sv = 0$, $0 \leq v \leq u \cdot r$
2. Dual feasibility: $S^\top \lambda + \mu^+ - \mu^- = c$, $\mu^+, \mu^- \geq 0$
3. Complementary slackness: $\mu^+_j (u_j r_j - v_j) = 0$ and $\mu^-_j v_j = 0$ for all $j$

Strong duality (primal obj = dual obj) is equivalent to conditions 1–3 holding simultaneously. Enforcing conditions 1 and 2 plus the strong duality equality in the outer MILP is sufficient to guarantee that any feasible $(v, \lambda, \mu^+, \mu^-)$ is a primal-optimal FBA solution.

Step 5: Linearizing the Bilinear Products $r_j \cdot \mu^+_j$

The strong duality constraint nominally reads $v_{\text{biomass}} = \sum_j u_j r_j \mu^+_j$, which contains bilinear products: $r_j \in {0,1}$ (binary) times $\mu^+_j \geq 0$ (continuous). These are nonlinear and cannot appear in an LP or MILP directly.

Since $r_j$ is binary, we can exactly linearize each product using the McCormick substitution. Define auxiliary continuous variables:

$$ p_j ;=; r_j \cdot \mu^+_j, \qquad q_j ;=; r_j \cdot \mu^-_j $$

and replace each bilinear term with four linear inequalities (assuming $\mu^+_j,, \mu^-_j \in [0, M]$):

$$ p_j \leq M, r_j \qquad p_j \leq \mu^+_j \qquad p_j \geq \mu^+_j - M(1 - r_j) \qquad p_j \geq 0 $$

The logic is exact:

• When $r_j = 1$: the four constraints force $p_j = \mu^+_j$ exactly.
• When $r_j = 0$: the first constraint forces $p_j \leq 0$, so $p_j = 0$.

The strong duality constraint becomes the linear equality:

$$ v_{\text{biomass}} ;=; \sum_j \bigl(u_j, p_j - \ell_j, q_j\bigr) $$

which, since $\ell_j = 0$ after splitting, simplifies to $v_{\text{biomass}} = \sum_j u_j, p_j$.

Step 6: Adding Complementary Slackness Constraints

Complementary slackness (CS) is implied by primal/dual feasibility plus strong duality, so it is mathematically redundant. However, the LP relaxation of the MILP (obtained by relaxing the binary constraints) is much weaker without it — the branch-and-bound solver takes many more nodes to close the optimality gap. We add explicit CS constraints using additional binary indicator variables $w^+_j,, w^-_j \in {0,1}$:

$$ \mu^+_j ;\leq; M_{\text{dual}}, w^+_j \qquad\quad u_j r_j - v_j ;\leq; M_{\text{flux}},(1 - w^+_j) $$

$$ \mu^-_j ;\leq; M_{\text{dual}}, w^-_j \qquad\quad v_j - \ell_j r_j ;\leq; M_{\text{flux}},(1 - w^-_j) $$

The first pair enforces: either $\mu^+_j = 0$ (reaction not at its upper bound), or the upper bound is tight ($u_j r_j - v_j = 0$). The second pair enforces: either $\mu^-_j = 0$ (reaction not at its lower bound), or the lower bound is tight ($v_j = 0$). This dramatically tightens the LP relaxation and reduces branch-and-bound nodes.

---

4. The Complete MILP Formulation

Combining everything, the bilevel OptKnock reduces to the following single-level MILP over variables $v, y, r, \lambda, \mu^+, \mu^-, p, q, w^+, w^-, z$:

$$\max \quad v_{\text{target}} \tag{outer objective}$$

Primal feasibility (inner LP):

$$S,v = 0 \tag{mass balance}$$

$$0 ;\leq; v_j ;\leq; u_j, r_j \quad \forall j \tag{flux bounds with knockout}$$

Dual feasibility (KKT stationarity):

$$\sum_i S_{ij},\lambda_i ;+; \mu^+_j ;-; \mu^-_j ;=; c_j \quad \forall j \tag{stationarity}$$

$$\mu^+_j,; \mu^-_j ;\geq; 0 \quad \forall j$$

Optimality certificate (strong duality):

$$v_{\text{biomass}} ;=; \sum_j u_j, p_j \tag{strong duality, after splitting}$$

McCormick linearization of $r_j \cdot \mu^+_j$ and $r_j \cdot \mu^-_j$:

$$p_j \leq M r_j,\quad p_j \leq \mu^+_j,\quad p_j \geq \mu^+_j - M(1-r_j),\quad p_j \geq 0 \quad \forall j$$

$$q_j \leq M r_j,\quad q_j \leq \mu^-_j,\quad q_j \geq \mu^-_j - M(1-r_j),\quad q_j \geq 0 \quad \forall j$$

Complementary slackness (tightens LP relaxation):

$$\mu^+_j \leq M_{\text{dual}}, w^+_j, \qquad u_j r_j - v_j \leq M_{\text{flux}},(1 - w^+_j) \quad \forall j$$

$$\mu^-_j \leq M_{\text{dual}}, w^-_j, \qquad v_j \leq M_{\text{flux}},(1 - w^-_j) \quad \forall j$$

Gene-Protein-Reaction (GPR) logic (links gene presence $y$ to reaction activity $r$):

$$z_{j,c} \leq y_g \quad \forall g \in \text{clause}~c \text{ of reaction } j$$

$$z_{j,c} \geq \textstyle\sum_{g \in \text{clause}} y_g - (|\text{clause}| - 1)$$

$$r_j \leq \textstyle\sum_c z_{j,c} \qquad r_j \geq z_{j,c} \quad \forall c$$

Knockout budget:

$$\sum_i (1 - y_i) ;\leq; k \tag{at most $k$ gene knockouts}$$

Viability floor (optional):

$$v_{\text{biomass}} ;\geq; 0.1 \cdot v^{\text{WT}}_{\text{biomass}} \tag{cell must maintain $\geq$10% WT growth}$$

This is a Mixed-Integer Linear Program (MILP). The binary variables are $y$ (gene presence), $r$ (reaction activity), $w^+$, $w^-$ (CS indicators), and $z$ (GPR clause indicators). The continuous variables are $v$, $\lambda$, $\mu^+$, $\mu^-$, $p$, $q$. Gurobi solves it via branch-and-bound.

---

5. Objective Function Variants

The outer objective — what the engineer is optimizing — affects both what knockout sets are found and how the MILP behaves. Two variants are implemented.

Linear Objective

$$\max \quad v_{\text{target}} + t \cdot v_{\text{biomass}}$$

This maximizes a weighted sum of the target production rate and growth rate, with $t$ acting as a tradeoff parameter (units: hours). A larger $t$ favors fast-growing strains; a smaller $t$ favors high-producers. At $t = 0$, this is pure target maximization. This objective is a linear function of the variables, so the MILP structure is standard.

Log Objective

$$\max \quad \ln(v_{\text{target}}) + t \cdot v_{\text{biomass}} - \ln(v_{\text{biomass}})$$

This is derived from the batch production formula (see Section 7): maximizing the log of total batch production encourages the solver to balance target flux against growth rate rather than maximizing either to the exclusion of the other. This is a nonlinear objective that Gurobi handles via its general constraint feature (addGenConstrLog), which introduces piecewise-linear approximations internally. The log objective requires $v_{\text{target}} \geq \epsilon$ and $v_{\text{biomass}} \geq \epsilon$ constraints to keep the log well-defined.

The log objective typically takes longer to solve and has higher MILP objective values that do not directly translate to batch production — the relationship is nonlinear. This explains the large MILP error seen for the log objective in the results.

---

6. Greedy Baseline

The greedy algorithm is a sequential one-at-a-time heuristic. At each step, it evaluates every candidate gene for knockout, picks the one that maximizes the batch production score (see Section 7), adds it to the knockout set, and repeats for up to $k$ steps.

knocked_out = []
for step in 1..k:
    best_gene = argmax over genes not yet knocked out of:
        batch_production(q_p, μ, t)  subject to min_biomass ≥ 10% WT
    knocked_out.append(best_gene)

After selecting knockouts, a post-processing LP re-solves FBA under the knockout set and minimizes target flux at fixed optimal biomass — this reveals the forced-coupling floor: the minimum target production that must occur at maximum growth.

The greedy algorithm is fast (each step is just a COBRApy FBA call), but it is myopic: it cannot identify knockout combinations whose synergistic effect only appears when multiple genes are deleted simultaneously. OptKnock's MILP formulation considers all $k$ knockouts jointly.

---

7. Batch Production Metric

Both methods are ultimately evaluated on a common biological objective: total product accumulated during a batch fermentation of length $t$ hours, starting from a unit inoculum ($X_0 = 1$ g/L dry weight).

Assuming exponential growth at specific growth rate $\mu$ (= $v_{\text{biomass}}$) and constant specific productivity $q_p$ (= $v_{\text{target}}$):

$$P(t) ;=; q_p \cdot \frac{e^{\mu t} - 1}{\mu} \qquad (\mu > 0)$$

$$P(t) ;=; q_p \cdot t \qquad (\mu = 0, \text{ limiting case})$$

Units: mmol/gDW. This metric simultaneously rewards a high target flux and a high growth rate — the exponential amplifies the specific productivity by the total biomass produced during the fermentation. It is the standard objective used in the greedy's gene selection score.

---

8. Implementation Overview

Reversible Splitting (split_reversible)

Each reaction with $\ell_j < 0$ is decomposed into a forward part (non-negative, stoichiometry $+S_{:,j}$) and a reverse part (non-negative, stoichiometry $-S_{:,j}$). Both parts share the same knockout indicator $r_j$ via GPR rules. Splitting ensures all fluxes are non-negative, which is required for the McCormick linearization and simplifies the dual (the lower-bound dual term $\ell_j \mu^-_j$ vanishes).

GPR Encoding (reaction_gene_sets)

Gene-Protein-Reaction rules in the COBRA model are Boolean expressions (e.g., (b0351 and b0352) or b1702). The implementation converts each rule to Disjunctive Normal Form (DNF) using SymPy, then encodes each AND-clause as a binary auxiliary $z_{j,c}$:

• AND-clause: $z_{j,c} = 1$ iff all genes in the clause are present
• OR between clauses: reaction is active iff any clause is satisfied

Reactions with no GPR annotation are always active ($r_j = 1$, representing spontaneous reactions or transport).

Inner Dual Constraints (add_inner_dual_constraints)

This shared helper builds the dual variables ($\lambda$, $\mu^+$, $\mu^-$), mass balance constraints, stationarity constraints, and strong duality equality. It operates in two modes:

• FBA mode (r=None): for verification — strong duality uses $u^\top \mu^+$ directly (no McCormick)
• OptKnock mode (r provided): introduces McCormick variables $p, q$ and uses the linearized strong duality

Post-processing (postprocess_lp)

After the MILP finds a knockout set, postprocess_lp verifies it with a two-step LP:

1. Maximize biomass under the knockouts → $\mu^*$
2. Fix biomass at $\mu^*$, minimize target flux → production floor

The minimum target flux at maximum biomass is the true forced-coupling floor. This verification step catches cases where the MILP objective value does not correspond to actual forced production (particularly with the log objective).

---

9. Experimental Results

Experiments were run on the E. coli core model sweeping:

• Targets: ethanol (EX_etoh_e), succinate (EX_succ_e)
• Knockout budgets $k$: ${0, 2, 4, 7, 10}$
• Fermentation times $t$: ${1, 4, 8, 16}$ hours
• Methods: greedy, OptKnock (linear), OptKnock (log)
• Solver: Gurobi with a 300-second time limit per solve

Summary Statistics

Production (batch P(t), mmol/gDW)

                            mean     median
method   objective
greedy   batch      10407.647844   0.000000
optknock linear      5008.773718  37.935391
         log        46589.748987  17.421627

The greedy median of zero is striking. It reflects that for many parameter combinations — particularly $k=0$ (no knockouts allowed) and short fermentation times — the greedy algorithm finds no productive knockout set, leaving the cell producing negligible target flux. When $k=0$, the WT cell maximizes biomass and routes essentially no flux to ethanol or succinate. The high greedy mean (10407) is driven by a handful of large-$k$, large-$t$ cases where the greedy fortuitously finds a good sequential knockout.

The OptKnock linear objective has a non-zero median (37.9), indicating it consistently finds at least some forced production even at modest knockout budgets. The log objective's very high mean (46589) with a low median (17.4) reflects heavy skew from a few outlier cases where the nonlinear objective found unusual solutions — often not matching actual batch production (see MILP error below).

% Improvement of OptKnock over Greedy

              mean     median
objective
linear      260.06     49.91
log        2589.87     58.38

The median improvements — 50% for linear, 58% for log — tell the more reliable story: OptKnock consistently outperforms greedy by about half again in the typical case. The extreme means (260%, 2590%) are dominated by cases where greedy produced near-zero output while OptKnock found a genuinely growth-coupled solution, causing the ratio to blow up. Dividing by a near-zero greedy value amplifies any OptKnock output enormously.

MILP Error: |MILP Production − Actual Production|

              mean        median
objective
linear       17.12  3.31 × 10⁻⁶
log        4451.72  4.11 × 10⁻⁸

The near-zero medians for both objectives confirm that the MILP formulation is tight: in the vast majority of cases the MILP's internal production estimate exactly matches the post-verification LP. The large means reflect a small number of outlier cases:

• For the linear objective, errors are small overall (mean 17 mmol/gDW), with isolated large deviations likely caused by the solver hitting the time limit and returning a suboptimal integer solution.
• For the log objective, the large mean error (4452 mmol/gDW) arises because the log MILP objective is not directly comparable to batch production: the solver maximizes $\ln(v_{\text{target}}) + t v_{\text{biomass}} - \ln(v_{\text{biomass}})$, which can be large even when the actual batch production is moderate. The post-verification LP then reveals the true forced production, which can differ substantially.

MILP Solve Time

           mean    median
objective
linear     1.21     1.01
log        4.37     1.19

The linear objective solves in about one second at median, very fast for a MILP. The log objective takes slightly longer on average (4.4 s mean) due to the nonlinear constraint handling, but the median (1.2 s) is still near one second — a few hard instances pull the mean up. Both are well within the 300-second time limit.

Interpreting Knockout Strategies

The knockouts returned by OptKnock are the globally optimal set (within the budget $k$) for the specified objective. In the lab, these are candidate genes to delete via CRISPR or homologous recombination. Because the cell must produce the target to grow at all, the knockout strain is evolutionarily stable.

The greedy algorithm finds locally optimal knockouts sequentially, which is faster but can miss synergistic combinations. The MILP explores the full space of $\binom{n_{\text{genes}}}{k}$ combinations implicitly through branch-and-bound.