MDで標準結合自由エネルギーを求める際の補正項について

これはなに

分子動力学シミュレーション（Molecular Dynamics simulation, 以下MD）を行うことにより、タンパク質とリガンドの標準結合自由エネルギーを見積もることができます。その際に、標準結合自由エネルギーを求める際に考える必要のある体積補正(volume correction)をする必要がある場合があります。しかしこの補正について、chatGPTをはじめとするAIも良い答えを出してくれず、（特にMDを分散、並列して実行した際の）理論と実装の両方で日本語で容易に読めるものもすぐには見当たりませんでした。そこで、整理のために書くことにした文章がこれです。

背景

結合親和性の大きな分子が見つかると嬉しいです。（なにが嬉しいか分からない人はそこら辺の人をつかまえて聞いてください。）結合親和性、つまり結合に伴う自由エネルギー差は、MDを用いて見積もることができる物理量です。またMDでは詳細な結合ポーズが分かることもあり、理論、実験と異なる第3の方法として創薬に役に立つと（おそらく）考えらえています（知りませんが）。そのため、結合自由エネルギー、また標準結合自由エネルギー計算について抑えておく必要があります。

標準結合自由エネルギーの計算について

実験の観点から

タンパク質AとリガンドBが存在して、溶媒中でくっついたり離れたりしている状況を考えます。このとき実験で得られる値として、学生実験でやるように、解離定数が挙げられます。解離定数$K_d$は、化学平衡$A+B⇌AB$に対して次のように定義されます。

$$K_d=\dfrac{[A][B]}{[AB]}$$

ここで、$[A], [B], [AB]$はそれぞれタンパク質A、リガンドB、複合体ABの濃度を表します。

また、平衡定数$K$との間には次の関係が成り立ちます。

$$K = \dfrac{1}{C^o} K_d$$

ここで、$C$は基準濃度です。（活量の代わりに濃度を使う妥当性については、適当な教科書や学生実験のテキストを参照してください。）このとき標準結合自由エネルギーとの間には次の関係が成り立ちます。

$$ΔG^o = -k_B T \log \dfrac{K_d}{1M}$$

ここで、$k_B$はボルツマン定数、$T$は温度です。また、基準状態として標準状態$C^o=1M≒1/1661 Å^{-3}$です。

シミュレーションの観点から

MDを行うことにより計算できる値は、結合自由エネルギー差$ΔG$です。実験値$K$と比較するためには、標準状態における結合自由エネルギー$ΔG^o$を求める必要があります。$ΔG$に対して、$ΔG^o$を計算するための補正を、体積補正(volume correction)と呼びます。（この文章が分からない人は例えば次を考えてください。化学熱力学の教科書では$ΔG=ΔH-TΔS$と書かれます。ここでエントロピーSはシミュレーションを行う系の大きさ（箱の大きさ）によって変わります。任意の条件でのMDを適切に比較するためには補正をする必要がある気がしてきませんか。）

MDでのvolume correctionは、おそらくdoudouの論文[1]で初めて扱われました。論文を読めそうな人は論文を読んでください。umbrella samplingなどにより拘束をかけた系に対してvolume correctionを行っています。拘束のためのポテンシャルも追加で考えており、式の見た目が結構難しくなっています。より簡潔になっている導出として、buchらの論文[2]および、そのSupporting Informationが挙げられます。次節ではこの導出を追います。

導出

buchの論文のSupporting Informationにある導出を追っていきます。なお、これはWe follow the methodology of Doudou et al. (1)から始まるもので、Doudouらの論文の導出を整理したものです。

まず、自由エネルギーは次のように計算されます。

$$ΔG_{PMF} = -k_B T \log \dfrac{Q_{b}}{Q_{u}} \tag{4}$$

ここで、$Q_{b}, Q_{u}$はそれぞ結合状態、非結合状態における分配関数です。正規化の項が共通であることを踏まえると、その比は次のようになります。

$$\dfrac{Q_b}{Q_u}=\dfrac{\int_b \exp(-βW(r))dr}{\int_u \exp(-βW(r))dr} \tag{5}$$

ここで、$b, u$はそれぞれ結合状態、非結合状態を、また$β=1/k_B T$、$W(r)$は3次元空間でのPMF(potential of mean force)です。簡単のためにこれ以降$W(r)$の最小値を0とします。結合状態、非結合状態の定義については、PMFが平坦になる領域を領域を非結合状態とし、それ以外を結合状態とします。このとき、$Q_u$について、次が言えます。

$$\int_u \exp (-\beta W(r)) dr = \exp (-\beta ΔW) V_u \tag{6}$$

ここで、$ΔW$は結合状態、非結合状態のPMFの差であり（$W$の最小値が0になるように移動させていることから、$ΔW$はPMFが平坦になった時の値であり、非結合状態においては$W(r)=ΔW$で固定される）、$V_u$は非結合状態の体積、つまり、$V_u = \int _u dr$です。（これが$V_u$の定義です)。以上から、${Q_b}/{Q_u}$は次の形で書けます。

$$\begin{align} \dfrac{Q_b}{Q_u} &= \dfrac{\int_b \exp(-βW(r))dr}{\exp (-\beta ΔW) V_u} \\ &= \dfrac{V_b}{V_u} \exp (\beta ΔW) \end{align} \tag{7}$$

ここで、$V_b=\int_b \exp(-βW(r))dr$です。($V_u$のように$V_b = \int_b dr$ではありません。単純な体積ではなく重みづけされた値です）。

さて、標準濃度における結合自由エネルギー$ΔG^0$は次のように書けます。

$$ΔG^o = ΔG_{PMF} + ΔG_V \tag{8}$$

ここで、

$$ΔG_v = -k_B T \log \dfrac{V_u}{V^{o}} \tag{9}$$

です。(4)、(7)、(8)、(9)から$ΔG^o$は次のように書けます。

$$\begin{align} ΔG^o &= ΔG_{PMF} + ΔG_v \\ &= -k_B T \log (\dfrac{V_b}{V_u} \exp (\beta ΔW)) -k_B T \log \dfrac{V_u}{V^{o}} \\ &= -k_B T \log \dfrac{V_b}{V^{u}} -k_B T \log \exp (\beta ΔW) - k_B T\log \dfrac{V_u}{V^{o}} \\ &= -k_B T \log \dfrac{V_b}{V_u} \dfrac{V_u}{V^o} - ΔW \\ &= -ΔW-k_B T \log \dfrac{V_b}{V^{o}} \\ \end{align} \tag{11}$$

これにより、系の大きさ（箱の大きさ）により値が変わる$V_u$は打ち消され、結合状態の体積$V_b=\int_b \exp(-βW(r))dr$を求めることで、標準濃度における結合自由エネルギーを求められていることが分かります。

導出を受けての整理

改めて$ΔG^o$の表記を整理します。式(11)から分かるように、$ΔG^o$の計算は大きく1,2の2つに分けられ、補正項は$V_u$、$V_b$という異なる2つの値のどちらかを計算する必要があります。

1. $ΔG^o = ΔG_{PMF}-k_B T \log \dfrac{V_u}{V^{o}}$
2. $ΔG^o = -ΔW - k_B T \log \dfrac{V_b}{V^{o}}$

ここで、$ΔW$と$ΔG_{PMF}$の関係について一応整理しておきます。まず重要なこととして、上の式から明らかですがこの2つは同一の値ではありません。式(4)、(7)から

$$\begin{align} ΔG_{PMF} &= -k_B T \log (\dfrac{V_b}{V_u} \exp (\beta ΔW)) \\ &= -k_B T \log \dfrac{V_b}{V_u} - ΔW \end{align}\tag{12}$$

であるように、(符号の違いを無視しても)$V_b=V_u$の場合を除き、この2つが一致することはありません。

さて、ここで$W$の計算について考えてみます。doudouやbuchの論文では触れられていませんが、Pandeのfolding@homeに代表されるMDの分散実行、またそれに続く、F.Noeに代表されるマルコフ状態(MSM)を用いたMDの解析では$W$はトラジェクトリの確率を用いて計算されます。
具体的には、MDにより得られたトラジェクトリを適当にクラスタリングして、各クラスターを$s_i$とします。また、MSMにおける定常分布での存在確率を$π_{s_i}$とします。このとき状態$s_i$のPMF $W(s_i)$は次のように計算されます。

$$W(s_i) = -k_B T \log π_{s_i} \tag{13}$$

特に結合、解離を行うMDにおいて、$ΔW$は1次元、および3次元での計算について分かっていれば何とかなることが多いです。1次元においては、$ΔW$は$W(s_i)$の最小値と最大値の差を求められます。つまり、

$$ΔW= -k_B T \log \dfrac{\min (π_{s_i})}{\max (π_{s_i})} \\$$

です。一般には結合状態で確率最大、非結合状態で確率最小となり、この2点（2クラスター）の差を考えればいいです。。3次元空間での$ΔW$の計算は、クラスター単位で最大値と最小値の差を取るのではなく、結合サイトを中心とする球面で周辺化する必要があります。つまり

$$W(d_i)= -k_B T \log \sum_{r=d_i} ^ {d_i+Δd} π_{s_i} \\$$

としたとき、$ΔW=\max(W(d_i))-\min(W(d_i))$として計算できます。なお、2つ目の式(3次元での$W$)は、$Δd$を適切に取ると最初の式(1次元でのW)と一致することは直感的に分かります。
続いて、$ΔG_{PMF}$について考えます。$ΔW$と同様に、$ΔG_{PMF}$も$π$を用いて書くと次のようになります。

$$\begin{align} ΔG_{PMF} &= -k_B T \log \dfrac{Q_{b}}{Q_{u}} \\ &= -k_B T \log \dfrac{\int_b \exp(-βW(r))dr}{\int_u \exp(-βW(r))dr} \\ &= -k_B T \log \dfrac{\int_b \exp(-β×(-k_BT\log π_{s_i}))dr}{\int_u \exp(-β×(-k_BT\log π_{s_i}))dr} \\ &= -k_B T \log \dfrac{\Sigma_b π_{s_i}}{\Sigma_u π_{s_i}}\\ \end{align} \tag{14}$$

最後の式変形においては、計算実行時は離散化されることを踏まえて、積分を和に置き換えています。$ΔW$の計算は単純な2クラスターの差分を計算していましたが、$ΔG_{PMF}$は結合領域、非結合領域それぞれの確率を足し合わせたものの比を取っています。

実装

ここまでを踏まえると、改めて$ΔG^o$は次のように書くことができます。

1. $ΔG^o = -k_B T \log \dfrac{\Sigma_b π_{s_i}}{\Sigma_u π_{s_i}}-k_B T \log \dfrac{\Sigma_u v_i}{V^{o}}$

2. $ΔG^o = -k_B T \log \dfrac{\min \sum_{r=d_i} ^ {d_i+Δd} π_{s_i}}{\max\sum_{r=d_i} ^ {d_i+Δd} π_{s_i}} - k_B T \log \dfrac{\Sigma_b \pi_{s_i} v_i}{V^{o}}$

計算機を用いて計算する都合上、離散化した形で書きました。それぞれについて、実装の詳細を考えます。

1について

$ΔG^o = -k_B T \log \dfrac{\Sigma_b π_{s_i}}{\Sigma_u π_{s_i}}-k_B T \log \dfrac{\Sigma_u v_i}{V^{o}}$
初項の計算は面倒です。まずハイパーパラメータとして、結合領域、非結合領域の定義が必要です。例えば結合領域の定義を、結合距離が2nm以下といった距離を基準に決めたり、$π_i<\max(π)±\epsilon$や、$π_i!=\min(π)±\epsilon$といったPMFを基準に決めるものが考えられます。この定義ができた後は、マルコフ状態モデルの定常分布における確率を足し合わせ、その比を計算することで初項が計算できます。
補正項の計算は簡単です。$V_u=\Sigma_u v_i$は重みが存在しません。そのため、非結合領域の体積を計算しさえすればいいです。結合領域の体積があまり大きくないことと、対数を取ることから、$\Sigma_u v_i$はサンプリング領域について、その凸包を求め、その体積を計算するといった近似でそれなりの解を得ることができると思われます（初項のハイパラのほうが最終的な値に与える影響が大きいため、結合領域も$V_u$として含めるこの近似でも問題ありません。）

2について

$ΔG^o = -k_B T \log \dfrac{\min \sum_{r=d_i} ^ {d_i+Δd} π_{s_i}}{\max\sum_{r=d_i} ^ {d_i+Δd} π_{s_i}} - k_B T \log \dfrac{\Sigma_b \pi_{s_i} v_i}{V^{o}}$

初項の計算は1次元においては簡単です。$π_i$の最大値、最小値を求め、その比を取ればいいです(Wの形にして差をとってもいいです）。3次元の場合は球面の確率を足し合わせて、それが最大となるような半径と、最小となるような半径について、その確率の差を取ればいいです。
補正項の計算は面倒です。結合領域の重み付き体積を計算する必要があります。こちらもどうせ雑な近似でもほとんど影響がないと考えられますが、例えば次のようにして計算することができます。

• サンプリング領域についてその凸包を求める。
• 各クラスターについて、その存在確率を計算する。
• クラスター中心でボロノイ分割を行い、その体積を計算する。このとき、分割の端については、その体積が発散してしまうため、凸包の体積を境界に用いる。
• 各クラスターについて、その体積を存在確率で重みづけし、足し合わせる。

ボロノイ分割はクラスターに対して行うため、計算量は（それ以前にトラジェクトリ全体を見ていることを考えると）無視できるほど小さいです。また、非結合領域の存在確率は十分に小さいことから（重みづけした際に非結合領域の体積は無視できるほど小さくなることから）、結合/非結合領域の定義をする必要がなくトラジェクトリ全体に対して凸包やボロノイ分割を考えて問題ありません。

まとめ

標準結合自由エネルギーの計算について、その導出と実装について整理しました。計算に$W$を使うのか、$ΔG$を使うのかで補正項として計算する値も変わることが分かると思います。この文章を読んで、標準結合自由エネルギーの計算について理解していただけたら幸いです。なにか間違いがあれば、ご指摘いただけると幸いです。

参考文献

1. Doudou, S., Burton, N. A., Henchman, R. H., & Oostenbrink, C. (2009). Standard free energy of binding from a one-dimensional potential of mean force. Journal of Chemical Theory and Computation, 5(10), 2764-2770.
2. Buch, I., Harvey, M. J., Giorgino, T., Anderson, D. P., De Fabritiis, G., & Chipot, C. (2011). Complete reconstruction of an enzyme-inhibitor binding process by molecular dynamics simulations. Proceedings of the National Academy of Sciences, 108(25), 10184-10189.