Constant pH MD
Constant pH MDは，あるpH条件（例えばpH=7.0）において脱プロトンorプロトン状態の存在比をMD計算で予測するための手法です。複数のpH条件下で脱プロトンorプロトン状態の存在比を計算（=滴定曲線を計算）することにより，pKaを予測することに繋がります。本記事では，参考文献１で採用されているConstant pH MDの背景理論を概説したいと思います。
pH下を表現するポテンシャル関数
Constant pH MDを実現するためには，脱プロトンorプロトン状態の両方を表現することができ，かつpHという条件を盛り込んだポテンシャル関数を用意する必要があります。そのため，以下のような関数が用いられます。
$$
\begin{array}{}U(\mathbf{X},\lambda,{\rm pH})&=&\lambda U_0(\mathbf{X})+(1-\lambda)U_1(\mathbf{X})+C(\lambda,{\rm pH})\\U_0(\mathbf{X})&:&脱プロトン状態を表すポテンシャル関数\\U_1(\mathbf{X})&:&プロトン状態を表すポテンシャル関数\\C(\lambda,{\rm pH})&:&ある{\rm pH}下でのサンプリングを表現するためのバイアスポテンシャル\end{array}
$$
$${\lambda}$$は脱プロトンorプロトン状態を一つのポテンシャル関数で表す役割を担っており，０から１までの値を取ります。
$${C(\lambda,{\rm pH})}$$は以下のような内容で構成されます。
$$
\begin{array}{}C(\lambda,{\rm pH})&=&-G_{\rm calc}^{\rm model}(\lambda)+\lambda\Delta G_{\rm exp}^{\rm model}({\rm pH})+U^{\rm barr}(\lambda)\\\Delta G_{\rm exp}^{\rm model}({\rm pH})&=&RT({\rm pKa}_{\rm exp}^{\rm model}-{\rm pH}\ln 10)\\U^{\rm barr}(\lambda)&=&-4\beta(\lambda-0.5)^2\end{array}
$$
pH条件を記述する役割を担う項は，$${G_{\rm calc}^{\rm model}(\lambda)}$$と$${\Delta G_{\rm exp}^{\rm model}({\rm pH})}$$の２つになります。$${U^{\rm barr}(\lambda)}$$は物理的に意味のある$${\lambda=0\ {\rm or}\ 1}$$近傍の状態を効率良くサンプリングする役割を担います。$${C(\lambda,{\rm pH})}$$としてなぜこのような関数が採用されるかを直感的に理解しがたいと感じる方もいると思いますが，それについては後ほど説明します。
$${G_{\rm calc}^{\rm model}(\lambda)}$$は，モデルケースの$${\lambda}$$に対する平均力ポテンシャルであり，Constant pH MD実行前に予め計算しておく必要があります。ここで言及しているモデルケースとは，計算対象と同一の化学種で既に実験値が分かっているものになります。例えば，あるタンパクのグルタミン酸のpKaが予測対象であった場合，グルタミン酸の単分子がモデルケースに相当します。同様にして，$${\Delta G_{\rm exp}^{\rm model}({\rm pH})}$$はモデルケースの脱プロトン状態とプロトン状態の自由エネルギー差を$${{\rm pH}}$$の関数として表したものになります。
λ Dynamics
物理的に意味ある$${\lambda=0}$$と$${\lambda=1}$$の状態の自由エネルギー差をMD計算で求めたいわけですが，$${\lambda}$$に関するサンプリングを$${\lambda}$$ Dynamicsで実現することを考えます。
$${\lambda}$$ Dynamicsは，$${\lambda}$$を力学変数として取り扱うことにより，１つのMD計算で$${0\leq\lambda\leq 1}$$の状態のサンプリングを実現することができます。
取る値を0～1に制限し，かつ連続変数であることが望ましいことを考えると，$${\lambda}$$そのものではなく$${\lambda=\sin\theta}$$とし，$${\theta}$$を力学変数とした方が便利です。
$${\lambda}$$ Dynamicsのメリットの一つとして，$${\lambda(=\sin\theta)}$$を導入する前後でMD計算として計算コストがさほど変わらないことが挙げられます。その結果，タンパクのようなプロトン状態を決めたい箇所が複数あるような分子を扱う場合においても，通常のMD計算と同等の計算コストで$${\lambda}$$ Dynamicsの計算を実行できます。
脱プロトンorプロトン状態の自由エネルギー差
$${\lambda}$$ Dynamicsで得られたサンプリングデータから，$${\lambda=0}$$と$${\lambda=1}$$の頻度比を計算することを考えます。
$$
\begin{array}{}\frac{\rho(\lambda=0, {\rm pH})}{\rho(\lambda=1, {\rm pH})}&=&\frac{\int{\rm d}\mathbf{X}\exp\left(-\beta U(\mathbf{X},\lambda=0,{\rm pH})\right)}{\int{\rm d}\mathbf{X}\exp(-\beta U(\mathbf{X},\lambda=1,{\rm pH}))}\\&=&\frac{\exp\left(\beta G_{\rm calc}^{\rm model}(\lambda=0)\right)\int{\rm d}\mathbf{X}\exp(-\beta U_1(\mathbf{X}))}{\exp\left(\beta\left(G_{\rm calc}^{\rm model}(\lambda=1)-\Delta G_{\rm exp}^{\rm model}({\rm pH})\right)\right)\int{\rm d}\mathbf{X}\exp(-\beta U_0(\mathbf{X}))}\\&=&\exp\left(-\beta\left(\Delta G_{\rm calc}-\Delta G_{\rm calc}^{\rm model}+\Delta G_{\rm exp}^{\rm model}({\rm pH})\right)\right)\\\Delta G_{\rm calc}&:=&\frac{\int{\rm d}\mathbf{X}\exp(-\beta U_1(\mathbf{X}))}{\int{\rm d}\mathbf{X}\exp(-\beta U_0(\mathbf{X}))}\\\Delta G_{\rm calc}^{\rm model}&:=&G_{\rm calc}^{\rm model}(\lambda=0)-G_{\rm calc}^{\rm model}(\lambda=1)\end{array}
$$
ここで，$${\beta}$$は逆温度になります。
さらに以下の等式が成立すると仮定します。
$$
\Delta G_{\rm exp}({\rm pH})-\Delta G_{\rm exp}^{\rm model}({\rm pH})=\Delta G_{\rm calc}-\Delta G_{\rm calc}^{\rm model}
$$
この式は，計算で使用したポテンシャルが，計算対象とした化合物とモデルケースの脱プロトン→プロトン状態の自由エネルギー差を正確に記述できると仮定したことと等価です。
上記の仮定の下，$${\lambda}$$ Dynamicsで得られた$${\lambda=0}$$と$${\lambda=1}$$の頻度比は以下の情報を計算していることになります。
$$
\frac{\rho(\lambda=0, {\rm pH})}{\rho(\lambda=1, {\rm pH})}=\exp(-\beta\Delta G_{\exp}(\rm pH))
$$
以上により，Constant pH MD計算はあるpH条件下での脱プロトン→プロトン状態の自由エネルギー差の実験値を計算で予測していることになります。
モデルケースのλに対する平均力ポテンシャル
$${\lambda}$$ DynamicsによるConstant pH MD計算を実現するためには，$${0\leq\lambda\leq 1}$$の範囲で$${G_{\rm calc}^{\rm model}(\lambda)}$$の値を決めなければなりません。幸いにも，$${G_{\rm calc}^{\rm model}(\lambda)}$$の形状は比較的に単純な場合が多く，FEP計算等で離散的に$${G_{\rm calc}^{\rm model}(\lambda)}$$の値を求めた後，2～4次の多項式フィッティングで事足りるみたいです。あるモデルケースにおいて，多項式近似の$${G_{\rm calc}^{\rm model}(\lambda)}$$を求めてしまえば，同種の化学種の計算に対しては使いまわすことができます。
ヒルの式による滴定曲線のフィッティング
滴定曲線の形状からpKaの値を求めることを考えると，滴定曲線をpHを連続変数とした関数で表現すると便利です。
ヒルの式で各pH下における脱プロトン状態とプロトン状態の存在比をフィッティングすることを考えます。
$$
S({\rm pH}):=\frac{\rho(\lambda=0, {\rm pH})}{\rho(\lambda=1, {\rm pH})}=\frac{1}{1+10^{n({\rm pKa-pH})}}
$$
$${n, {\rm pKa}}$$がフィッティングパラメータになります。フィッティングによって決まった$${\rm pKa}$$の値がConstant pH MD計算での最終的なアウトプットである$${\rm pKa}$$予測値になります。
まとめ
本記事の要点を以下にまとめます。
Constant pH MDは，あるpH条件において脱プロトンorプロトン状態の存在比をMD計算で予測するための手法である
脱プロトンorプロトン状態は，$${\lambda}$$変数で表現され，また$${\lambda}$$の値については$${\lambda}$$ Dynamicsでサンプリングされる
モデルケースのpKaの実験値，及びλに対する平均力ポテンシャルが事前に必要となる