はじめに
初投稿 & 初アドベントカレンダー参加になります
普段はメーカーでバイオ系研究職として働いています。
ニッチな記事ですが、お手柔らかにお願いします...！🙌
今回は、微生物によるモノづくりで、ターゲット化合物の生産を向上させる代謝反応のノックアウト組み合わせを探索したい時に利用する、OptKnock (Burgard et al, 2003)1 について、原論文をもとに理論部分の解説をしていこうと思います。
20年以上前の古典的な手法ではありますが、以降の RobustKnock や OptGene などの重要な手法に大きく影響を与えていますので、じっくり紐解く価値は高いと思います。
次回の記事 ～②実装編～ では、数理最適化モデリングライブラリ optlang を使った実装（Python）をやる予定です。
FBA (Flux Balance Analysis)
まずは、通常の FBA (Flux Balance Analysis) を復習しましょう。
$\mathscr{M} = \lbrace 1, \ldots, M \rbrace$ 個の代謝物と, $\mathscr{N} = \lbrace 1, \ldots, N \rbrace$ 個の反応で構成される定常状態の代謝ネットワークにおいて、細胞成長を目的関数としたときのフラックス分布の一つは以下のような数理最適化問題として表現されます。
\begin{align}
& \textbf{maximize} \quad \boldsymbol v_{biomass} \\[0.5em]
& \textbf{subject to} \quad \sum_{j \in \mathscr{N}} S_{ij}v_{j} = 0, \quad \forall i \in \mathscr{M} \\
& \quad\quad\quad\quad\quad\quad v^{lb}_{j} \le v_j \le v^{ub}_{j}, \quad \forall j \in \mathscr{N} \\
& \quad\quad\quad\quad\quad\quad v_{glc} = v_{glc\_uptake} \\
& \quad\quad\quad\quad\quad\quad v_{atp} \ge v_{atp\_main} \\
& \quad\quad\quad\quad\quad\quad v_{biomass} \ge v^{target}_{biomass} \\
& \quad\quad\quad\quad\quad\quad v_j \ge 0, \quad \forall j \in \mathscr{N}_{irrev} \\
& \quad\quad\quad\quad\quad\quad v_j \ge 0, \quad \forall j \in \mathscr{N}_{secr\_only} \\
& \quad\quad\quad\quad\quad\quad v_j \in \mathcal{R}, \quad \forall j \in \mathscr{N}_{rev}
\end{align}
基本なのであまり説明は不要かと思いますが、FBAではおなじみの「定常状態」・「細胞の成長最大化」を数理最適化問題に反映させたものとなっています。
OptKnock の構造
ではでは、FBAの復習を終えたところで、本題である OptKnock の概要を見ていきましょう。
OptKnockは、細胞増殖の最大化を制約として、細胞につくらせたいターゲット化合物の生産を最大化するような、代謝反応（遺伝子）のノックアウト組み合わせを大域的に探索する手法です。
OptKnockの定式化は下図のように、2つの最適化問題がネストした、いわゆる bilevel-optimization (2レベル最適化問題) の形式をしています。
（原論文 Figure 1. より引用）
上位レベルの問題 (ターゲット化合物の生産最大化) の中に下位レベルの問題 (細胞の成長最大化) が組み込まれおり、下位レベルの最適解集合を上位レベルの制約として課しています。
OptKnock の定式化
それでは、具体的に数式を見ていきましょう！
まず、ノックアウト対象となる代謝反応（遺伝子）をバイナリ変数 $y_j$ で定義します。
y_j =
\begin{cases}
1 & (v_j \text{ is active}) \\
0 & (v_j \text{ is not active}, \quad \forall j \in \mathscr{N})
\end{cases}
そして、以下のように $v_j$ の上下限制約に組み込んで半連続変数として定義します。
v^{lb}_{j} \cdot y_j \le v_j \le v^{ub}_{j} \cdot y_j, \quad \forall j \in \mathscr{N}
$y_j = 0 \enspace \text{(knockout)}$ の時には $0 \le v_j \le 0$, つまり $v_j = 0$ となって破壊されており, $y_j = 1 \enspace \text{(active)}$ の時には通常の上下限制約 $v^{lb}_j \le v_j \le v^{ub}_j$ になっており、目的の挙動がちゃんとできていることが確認できますね！
定義したバイナリ変数 $y_j$ を使うと、OptKnock全体は以下のように書けます。
\begin{align}
& \underset{\boldsymbol{y}}{\textbf{maximize}} \quad v_{chemical} \\[0.5em]
& \textbf{subject to} \quad \underset{\boldsymbol{v}}{\textbf{maximize}} \quad\;\; v_{biomass} \quad \textbf{(Primal)}\\[0.5em]
& \quad\quad\quad\quad\quad\quad v_{biomass} \ge v^{target}_{biomass} \\[0.5em]
& \quad\quad\quad\quad\quad\quad y_j \in \lbrace0, 1\rbrace \\[0.5em]
& \quad\quad\quad\quad\quad\quad \sum_{j \in \mathscr{N}}(1-y_j) \le K \\[0.5em]
& \quad\quad\quad\quad\quad\quad \textbf{subject to} \quad \sum_{j \in \mathscr{N}} S_{ij}v_{j} = 0 \\[0.5em]
& \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad v^{lb}_{j} \cdot y_j \le v_j \le v^{ub}_{j} \cdot y_j \\[0.5em]
& \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad v_{glc} = v_{glc\_uptake} \\[0.5em]
& \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad v_{atp} \ge v_{atp\_main} \\[0.5em]
&
\end{align}
ここで、 $v_{chemical}$ はターゲット化合物の生産フラックス、$K$ はノックアウトする反応数の上限です。
このように定式化することで、「細胞が成長するためにはターゲット化合物を生産せざるを得ない」、という状況を作り出すノックアウト組み合わせ $\boldsymbol{y}$ を求めることができます。
どうやって解くの？
この2レベル最適化問題をそのままソルバーに投げたいところですが、残念ながらそれはできません❌
ソルバーが扱えるのは 数式（静的な条件） だけであり、「下位レベル問題（FBA）を解け！（$\text{maximize} \ v_{biomass}$）」 という 命令（動的な条件） を制約の中に含めることはできないからです。
そのため、何かしらの工夫を考える必要があります。
FBA の命令 "$\text{maximize} \ v_{biomass}$" を、何とかして数式だけの条件に書き換えたいですよね...🤔
「FBAを解くためにソルバーを回した結果」と、「ある数式を満たす解」が 完全に等価（必要十分） であれば、FBA を実行する代わりにその数式を上位レベル問題の制約として与えれば解決できそうです。
では、「細胞がバイオマスを最大化しきった状態（下位レベルの最適解）」であることを証明する数式とは一体なんでしょうか？
例えば、非線形凸関数であれば、微分がゼロという条件が最大値の証明として使えそうです。ただ、今回は線形計画問題（LP）なので、微分は使えません。
そこで登場するのが、強双対定理 です。
強双対定理の登場
LPでは、主問題（Primal）と双対問題（Dual）の間に、以下の強双対定理が成り立ちます。
強双対定理
主問題（P）に最適解 $x^*$ が存在すれば、双対問題（D）にも最適解 $y^*$ が存在し、以下が成り立つ。
\boldsymbol{c}^{T}\boldsymbol{x^*} = \boldsymbol{b}^{T}\boldsymbol{y^*}
つまり、
$\textbf{Pの目的関数値} = \textbf{Dの目的関数値}$
が成り立つということです。
下位レベル問題の双対問題を考えて、強双対定理を適用すれば、「下位レベル問題を解け（FBA）」という命令の代わりに、
$\textbf{FBAの最適値} = \textbf{双対問題の最適値}$
という等式制約を上位レベル問題に組み込めば、通常の MILP（混合整数線形計画問題）としてソルバーで処理できそうです！
双対問題を考えよう
方針が立ったところで、まずは下位レベル問題（細胞の挙動）の双対問題を考えていきましょう。
$v_{chemical}$ を最大化する $\boldsymbol{y}$ は上位レベルで最適化するので、下位レベルを考える段階ではひとまず $\boldsymbol{y}$ は定数ベクトルとして扱います。したがって、下位レベルで求めるのは、ある反応ノックアウト組み合わせ $\boldsymbol{y}$ が適用された条件下で、細胞成長 $v_{biomass}$ を最大化するようなフラックス分布 $\boldsymbol{v}$ です。
まずは、下位レベル問題（主問題）の最後の2つの式（$v_{glc}, v_{atp}$）について、それぞれ以下のように変形し、その他の反応の上下限制約の中に含めてしまいましょう。
\begin{align}
& v_{glc\_uptake} \le v_{glc} \le v_{glc\_uptake} \\[0.5em]
& v_{atp\_main} \le v_{atp} \le v_{atp}^{ub} \\[0.5em]
\end{align}
まとめると以下のようにコンパクトになりました。
\begin{align}
& \textbf{＜主問題＞} \\[0.5em]
& \underset{\boldsymbol{v}}{\textbf{maximize}} \quad\;\; v_{biomass} \\[0.5em]
& \textbf{subject to} \quad \sum_{j \in \mathscr{N}} S_{ij}v_{j} = 0 \\[0.5em]
& \quad\quad\quad\quad\quad\quad v^{lb}_{j} \cdot y_j \le v_j \le v^{ub}_{j} \cdot y_j \\[0.5em]
\end{align}
次は、主問題を標準形に直しましょう（目的関数の符号反転に注意⚠️）。
\begin{align}
& \textbf{＜主問題＞} \\[0.5em]
& \textbf{minimize} \quad -v_{biomass} \\[0.5em]
& \textbf{subject to} \quad \sum_{j \in \mathscr{N}} S_{ij}v_j = 0 \\[0.5em]
& \quad\quad\quad\quad\quad\quad v_j^{lb} \cdot y_j -v_j \le 0 \\[0.5em]
& \quad\quad\quad\quad\quad\quad v_j - v_j^{ub} \cdot y_j \le 0\\[0.5em]
\end{align}
そして、双対変数（ラグランジュ乗数）を以下のように導入し、
\begin{align}
& 等式制約 \left(S_{ij}v_{j}=0\right) \;: \; \mu_i \left(\text{自由変数}\right)\\[0.5em]
& 不等式制約 \left(v_j^{lb} \cdot y_j - v_j \le 0 \right) \; : \; \lambda_j^{lb} \ge 0\\[0.5em]
& 不等式制約 \left(v_j - v_j^{ub} \cdot y_j \le 0 \right) \; : \; \lambda_j^{ub} \ge 0\\[0.5em]
\end{align}
ラグランジュ関数 $L$ を作ったら、$\boldsymbol{v}$ について項を整理します（ここで、$v_j^{lb}, v_j^{ub}$ は定数であることに注意⚠️）。
\begin{align}
L(\boldsymbol{v}, \mu, \lambda^{lb}, \lambda^{ub}) &= -v_{biomass} + \sum_{i \in \mathscr{M}}\mu_i \left( \sum_{j \in \mathscr{N}}S_{ij}v_j \right) + \sum_{j \in \mathscr{N}}\lambda_j^{lb}(v_j^{lb} \cdot y_j - v_j) + \sum_{j \in \mathscr{N}}\lambda_j^{ub}(v_j - v_j^{ub} \cdot y_j) \\
&= \sum_{j \in \mathscr{N}} v_j \left( \sum_{i \in \mathscr{M}} \mu_i S_{ij} - \lambda_j^{lb} + \lambda_j^{ub} - \delta_{j, \;biomass} \right) + \sum_{j \in \mathscr{N}}\left(\lambda_j^{lb}v_j^{lb}y_j - \lambda_j^{ub}v_j^{ub}y_j\right) \quad (\delta \text{ : クロネッカーデルタ})\\[0.5em]
\end{align}
ここで、$\delta_{j, \ biomass}$ はクロネッカーデルタで、$j = \text{biomass}$ の時に $1$ となり、$j \not= \text{biomass}$ の時に $0$ となります。
このラグランジュ関数は、$L(\boldsymbol{v}) = A\boldsymbol{v}+B$ の形をしているため、$\boldsymbol{v}$ に対しては単なる一次関数です。主問題は「$\boldsymbol{v}$を動かして $L$ を最大化」しようとします。
いま、各 $v_j$ は両方向の向きに反応が進行することが許されているため、$v_j$ は正負両方の値をとる自由変数です。そのため、$L(\boldsymbol{v})$ の傾き $A$ がゼロ以外の場合は、$L$ の最大値が $\infty$ になってしまい有限な値をとれません。双対問題は、主問題の良い上界を見つける問題ですので、有限の値をとるためには「傾き $A$ がゼロ」の場合しかありません。つまり、
\frac{\partial L}{\partial v_j} = \sum_{i \in \mathscr{M}} \mu_i S_{ij} - \lambda_j^{lb} + \lambda_j^{ub} - \delta_{j, \;biomass} = 0
という等式が導かれます。これこそが双対問題の制約になります。
そして、$L$ の定数項である $\sum_{j \in \mathscr{N}}\left(\lambda_j^{lb}v_j^{lb}y_j - \lambda_j^{ub}v_j^{ub}y_j\right)$ が双対問題の目的関数になります。
\begin{align}
& \textbf{＜双対問題＞} \\[0.5em]
& \textbf{maximize} \quad \sum_{j \in \mathscr{N}}\left(\lambda_j^{lb}v_j^{lb}y_j - \lambda_j^{ub}v_j^{ub}y_j\right) \\[0.5em]
& \textbf{subject to} \quad \sum_{i \in \mathscr{M}} \mu_i S_{ij} - \lambda_j^{lb} + \lambda_j^{ub} - \delta_{j, \;biomass} = 0 \\[0.5em]
\end{align}
強双対定理の再登場
さて、ここで強双対定理を思い出しましょう。
LPにおいて最適解が存在する場合、「主問題の目的関数値」と「双対問題の目的関数値」は一致する、でしたね。
ということで、ようやく OptKnock を解くキーとなる以下の式にたどりつきました！（ここでも符号に注意⚠️）
v_{biomass} = \sum_{j \in \mathscr{N}}\left(\lambda_j^{ub}v_j^{ub}y_j - \lambda_j^{lb}v_j^{lb}y_j\right)
OptKnockの下位レベル問題でFBAを解く代わりに、上記の等式を導入すれば、もともと bilevel-optimization だった問題が見事に single-level optimization へと変身しました！
\begin{align}
& \textbf{＜OptKnock (single-level)＞} \\[0.5em]
& \underset{\boldsymbol{y}}{\textbf{maximize}} \quad v_{chemical} \\[0.5em]
& \textbf{subject to} \quad v_{biomass} = \sum_{j \in \mathscr{N}}\left(\lambda_j^{ub}v_j^{ub}y_j - \lambda_j^{lb}v_j^{lb}y_j\right)\\[0.5em]
& \quad\quad\quad\quad\quad\quad \sum_{i \in \mathscr{M}} \mu_i S_{ij} - \lambda_j^{lb} + \lambda_j^{ub} - \delta_{j, \;biomass} = 0 \\[0.5em]
& \quad\quad\quad\quad\quad\quad y_j \in \lbrace0, 1\rbrace \\[0.5em]
& \quad\quad\quad\quad\quad\quad \sum_{j \in \mathscr{N}}(1-y_j) \le K \\[0.5em]
& \quad\quad\quad\quad\quad\quad \sum_{j \in \mathscr{N}} S_{ij}v_{j} = 0 \\[0.5em]
& \quad\quad\quad\quad\quad\quad v^{lb}_{j} \cdot y_j \le v_j \le v^{ub}_{j} \cdot y_j \\[0.5em]
& \quad\quad\quad\quad\quad\quad v_{biomass} \ge v^{target}_{biomass} \\
&
\end{align}
ここまで来ればただの MILP なので、ソルバーに投げて最適解を求めることができますね😎
めでたしめでたし！...と思いきや、最後の困難が待ち構えていました😵
最後の仕上げ Big-M
先ほどまでは定数として扱っていたバイナリ変数 $y$ ですが、single-level にした後はちゃんと変数として扱わないといけません。すると、強双対定理で導いた等式制約の中に、$\lambda_j \cdot y_j$ という非線形項が含まれていることになります。こいつを何とかしましょう！
幸いなことに、この非線形項は Big-M法 による線形化で解決できます。
$z_j = \lambda_j \cdot y_j$ とおくと、
\begin{align}
& y = 0 \; のとき \quad z = 0 \\
& y = 1 \; のとき \quad z = \lambda \\
\end{align}
という挙動を、掛け算を使わずに、以下のような不等式制約で表現することを考えます（$M$は十分大きな値）。
\begin{align}
& 0 \le z_j \le \lambda_j \\[0.5em]
& \lambda_j - M\left(1-y_j\right) \le z_j \le M \cdot y_j \\[0.5em]
\end{align}
こうして無事にソルバーで求解可能な形式の OptKnock を得ることができました！
\begin{align}
& \textbf{＜OptKnock (final)＞} \\[0.5em]
& \underset{\boldsymbol{y}}{\textbf{maximize}} \quad v_{chemical} \\[0.5em]
& \textbf{subject to} \quad v_{biomass} = \sum_{j \in \mathscr{N}}\left(v_j^{ub}z_j^{ub} - v_j^{lb}z_j^{lb}\right) \\[0.5em]
& \quad\quad\quad\quad\quad\quad 0 \le z_j^{lb} \le \lambda_j^{lb} \\[0.5em]
& \quad\quad\quad\quad\quad\quad \lambda_j^{lb} - M\left(1-y_j\right) \le z_j^{lb} \le M \cdot y_j \\[0.5em]
& \quad\quad\quad\quad\quad\quad 0 \le z_j^{ub} \le \lambda_j^{ub} \\[0.5em]
& \quad\quad\quad\quad\quad\quad \lambda_j^{ub} - M\left(1-y_j\right) \le z_j^{ub} \le M \cdot y_j \\[0.5em]
& \quad\quad\quad\quad\quad\quad \sum_{i \in \mathscr{M}} \mu_i S_{ij} - \lambda_j^{lb} + \lambda_j^{ub} - \delta_{j, \;biomass} = 0 \\[0.5em]
& \quad\quad\quad\quad\quad\quad \lambda_j^{lb} \ge 0, \quad \lambda_j^{ub} \ge 0 \\[0.5em]
& \quad\quad\quad\quad\quad\quad y_j \in \lbrace0, 1\rbrace \\[0.5em]
& \quad\quad\quad\quad\quad\quad \sum_{j \in \mathscr{N}}(1-y_j) \le K \\[0.5em]
& \quad\quad\quad\quad\quad\quad \sum_{j \in \mathscr{N}} S_{ij}v_{j} = 0 \\[0.5em]
& \quad\quad\quad\quad\quad\quad v^{lb}_{j} \cdot y_j \le v_j \le v^{ub}_{j} \cdot y_j \\[0.5em]
& \quad\quad\quad\quad\quad\quad v_{biomass} \ge v^{target}_{biomass} \\[0.5em]
\end{align}
今度こそこれで終了です！🥳
あとは、この single-level OptKnock を optlang などの数理最適化モデリングライブラリを使ってゴリゴリ書いていってソルバーに投げれば、欲しかった反応ノックアウト組み合わせを得ることができます（実際には色々留意点がありますが、それはまた今度...）。
以上になります。最後までお読みいただき、ありがとうございました。
次回は ②実装編 に続く予定です！（たぶん）