はじめに
初投稿 & 初アドベントカレンダー参加になります。
ニッチな記事ですが、お手柔らかにお願いします...！
今回は、微生物によるモノづくりで、ターゲット化合物の生産を向上させる代謝反応のノックアウト組み合わせを探索したい時に利用する、OptKnock (Burgard et al, 2003)1 について、原論文をもとに理論部分の解説をしていこうと思います。
次回の記事 ～②実装編～ では、数理最適化モデリングライブラリ optlang を使った実装（Python）をやる予定です。
FBA (Flux Balance Analysis)
まずは、通常の FBA (フラックスバランス解析) を復習しましょう。
代謝物集合 $\mathcal{I}$ と反応集合 $\mathcal{J}$ で構成される定常状態の代謝ネットワークを考えます。細胞成長を目的関数としたときの最適なフラックス分布は、以下のような数理最適化問題を解くことで求めることができます。
\begin{align}
& \textbf{maximize} \quad v_{biomass} \\[0.5em]
& \textbf{subject to} \quad \sum_{j \in \mathcal{J}} S_{ij}v_{j} = 0, \quad \forall i \in \mathcal{I} \\
& \quad\quad\quad\quad\quad\quad v^{lb}_{j} \le v_j \le v^{ub}_{j}, \quad \forall j \in \mathcal{J} \\
& \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} \\
\end{align}
FBAではおなじみの「定常状態」・「細胞成長の最大化」を数理最適化問題に反映させたものとなっています。
OptKnock の構造
FBAの復習を終えたところで、本題である OptKnock の概要を見ていきます。
OptKnockは、細胞増殖の最大化を制約として、細胞につくらせたいターゲット化合物の生産を最大化するような、反応のノックアウト組み合わせを大域的に探索する手法です。
OptKnockの定式化は下図のように、2つの最適化問題がネストした、いわゆる bilevel-optimization (2レベル最適化問題) の形式をしています。
（原論文 Figure 1. より引用）
上位問題 (ターゲット化合物の生産最大化) の中に下位問題 (細胞成長の最大化) が組み込まれおり、下位問題の最適解集合を上位問題の制約として課しています。
OptKnock の定式化
それでは、具体的に数式を見ていきましょう。
まず、ノックアウト対象となる反応をバイナリ変数 $y_j$ で定義します。
y_j =
\begin{cases}
1 & (\text{reaction} ~ j ~ \text{is active}) \\
0 & (\text{reaction} ~ j ~ \text{is not active})
\end{cases}
そして、以下のように $v_j$ の上下限制約に組み込む形で定義します。
v^{lb}_{j} \cdot y_j \le v_j \le v^{ub}_{j} \cdot y_j, \quad \forall j \in \mathcal{J}
$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 \\[0.5em]
& \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad y_j \in \lbrace0, 1\rbrace, \quad \forall j \in \mathcal{J} \\[0.5em]
& \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \sum_{j \in \mathcal{J}}(1-y_j) \le K \\[0.5em]
& \quad\quad\quad\quad\quad\quad \textbf{subject to} \quad \sum_{j \in \mathcal{J}} 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_{biomass} \ge v^{target}_{biomass} \\[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}$ を求めることができます。
強双対定理
定式化したものの、ソルバーは最適化問題が入れ子になった状態では処理できません。そこで、以下の強双対定理 を用います。
強双対定理
主問題（P）に最適解 $x^*$ が存在すれば、双対問題（D）にも最適解 $y^*$ が存在し、以下が成り立つ。
\boldsymbol{c}^{T}\boldsymbol{x^*} = \boldsymbol{b}^{T}\boldsymbol{y^*}
下位問題の双対問題を考えて強双対定理を適用し、
$\textbf{下位問題の最適値} = \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 \mathcal{J}} 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 \mathcal{J}} 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 \mathcal{I}}\mu_i \left( \sum_{j \in \mathcal{J}}S_{ij}v_j \right) + \sum_{j \in \mathcal{J}}\lambda_j^{lb}(v_j^{lb} \cdot y_j - v_j) + \sum_{j \in \mathcal{J}}\lambda_j^{ub}(v_j - v_j^{ub} \cdot y_j) \\
&= \sum_{j \in \mathcal{J}} v_j \left( \sum_{i \in \mathcal{I}} \mu_i S_{ij} - \lambda_j^{lb} + \lambda_j^{ub} - \delta_{j, \;biomass} \right) + \sum_{j \in \mathcal{J}}\left(\lambda_j^{lb}v_j^{lb}y_j - \lambda_j^{ub}v_j^{ub}y_j\right) \\[0.5em]
\end{align}
ここで、$\delta_{j, \ biomass}$ はクロネッカーデルタで、$j = \text{biomass}$ の時に $1$ となり、$j \not= \text{biomass}$ の時に $0$ となります。
このラグランジュ関数は、$\boldsymbol{v}$ に関して $L(\boldsymbol{v}) = A\boldsymbol{v} + B$ という一次式になっています。
ここで、主問題はすでに
\textbf{minimize} \quad -v_{biomass}
という最小化問題の形に直してあるので、双対関数は
g(\mu,\lambda^{lb},\lambda^{ub})
=
\inf_{\boldsymbol v}
L(\boldsymbol{v}, \mu, \lambda^{lb}, \lambda^{ub})
で定義されます。
また、フラックスの上下限制約はすでに不等式制約として外に出しているため、ラグランジュ関数の中では各 $v_j$ を自由変数として扱えます。したがって、もし $v_j$ の係数
\sum_{i \in \mathcal{I}} \mu_i S_{ij} - \lambda_j^{lb} + \lambda_j^{ub} - \delta_{j,\;biomass}
が 0 でなければ、$v_j$ を正または負の無限大方向へ動かすことで $L(\boldsymbol v)$ はいくらでも小さくできてしまい、
\inf_{\boldsymbol v} L(\boldsymbol v,\mu,\lambda^{lb},\lambda^{ub}) = -\infty
となります。
双対問題では有限な双対関数値をもつ必要があるため、各 $v_j$ の係数はすべて 0 でなければなりません。したがって、
\sum_{i \in \mathcal{I}} \mu_i S_{ij}
- \lambda_j^{lb}
+ \lambda_j^{ub}
- \delta_{j,\;biomass}
=0
\qquad \forall j \in \mathcal{J}
という双対問題の制約が導かれます。
そして、$L$ の定数項である $\sum_{j \in \mathcal{J}}\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 \mathcal{J}}\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 \mathcal{I}} \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 \mathcal{J}}\left(\lambda_j^{ub}v_j^{ub}y_j - \lambda_j^{lb}v_j^{lb}y_j\right)
OptKnockの下位問題を解く代わりに上記の等式を導入すれば、もともと 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 \mathcal{J}}\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 \mathcal{I}} \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 \mathcal{J}}(1-y_j) \le K \\[0.5em]
& \quad\quad\quad\quad\quad\quad \sum_{j \in \mathcal{J}} 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}
ただ、この時点では $\lambda_j \cdot y_j$ という非線形項が含まれているため、これを線形化する必要があります。
最後の仕上げ 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＞} \\[0.5em]
& \underset{\boldsymbol{y}}{\textbf{maximize}} \quad v_{chemical} \\[0.5em]
& \textbf{subject to} \quad v_{biomass} = \sum_{j \in \mathcal{J}}\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 \mathcal{I}} \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 \mathcal{J}}(1-y_j) \le K \\[0.5em]
& \quad\quad\quad\quad\quad\quad \sum_{j \in \mathcal{J}} 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}
あとはこの定式化に沿ってモデリングしていけば、最適な反応ノックアウト組み合わせを得ることができます。
以上になります。最後までお読みいただき、ありがとうございました！
次回は ②実装編 に続く予定です（たぶん）。