生体分子のシミュレーションにおけるImplicit溶媒モデルの適用と開発(2014)
以下に書いたものは、論文を基にした私の理解なので、正確性を欠きます。正確な情報を知りたい場合は、下記の論文を参照してください。また、下記の論文の図を参照することで、理解が深まると思います。
Design and application of implict solvent models in biomlecular simulations
前置き
マルチスケールの分子動力学法の計算は、様々な有益な情報を生み出すことが出来る。分子動力学の実際の計算の90％は水であると言われている。そのために水のスケーリングを変えることは、分市道力学法で扱える空間スケール、時間スケール共に大きく広げる可能性を秘めている。また、2013年のNovel prizeは、分子動力学法のマルチスケール化のパイオニアであるMartin Karplus,Michael Levitt,Ariel Warshelらに与えられた。得に、このREVIEWで紹介するImplicit溶媒モデルのパイオニアはWarshelである。
一番、簡単な方法として、溶質の第一水和殻の影響はSASA(Solvent-Accessible Surface Area)に比例することから、相互作用にSASAを含めるモデルがある。
Implicit溶媒のモデル
Free energy of solvation
溶媒和自由エネルギー\(\Delta G_{sol}\)は以下のように分解することが出来る。溶質の穴に溶媒が入るときのエネルギー\(\Delta G_{cav}\)、その穴に溶質が入るときのエネルギー\(\Delta G_{vdW}\)と溶質と溶媒間の静電相互作用\(\Delta G_{ele}\)。つまり、
\(
\begin{eqnarray*}
\Delta G_{sol}=\Delta G_{cav}+\Delta G_{vdW}+\Delta G_{ele}
\end{eqnarray*}
\)
となる。SASA法は、\(\Delta G_{cav}+\Delta G_{vdW}\)や全体のモデル化をすることが可能である。一方で、Poisson-BoltzmannやGeneralized-Bornは、\(\Delta G_{ele}\)をモデル化することが出来る。
\begin{eqnarray*}
\Delta G_{sol}=\Delta G_{cav}+\Delta G_{vdW}+\Delta G_{ele}
\end{eqnarray*}
\)
Solvent-accessible surface area(SASA) models
この方法はオクタノールと水間の移動による自由エネルギー差に基づいている。SASA法は溶媒と溶質の表面の面積に溶媒和自由エネルギーが比例することを想定しており、\(\Delta G_{solv}\)は、ポテンシャル\(V^{SASA}_{solv}\)で導かれる。微小面積当たりの溶媒和自由エネルギーを\(σ^{SASA}_i\)とし、原子iのSASAを\(SASA_i\)とすると
\(
\begin{eqnarray*}
V^{SASA}_{solv}({\bf r})=\sum_i σ^{SASA}_i\cdot SASA_i({\bf r})
\end{eqnarray*}
\)
となる。さらに厳密にモデル化するために、GB/SASAやPB/SASAがある。\begin{eqnarray*}
V^{SASA}_{solv}({\bf r})=\sum_i σ^{SASA}_i\cdot SASA_i({\bf r})
\end{eqnarray*}
\)
また溶質内の相互作用を長距離まで考えることが出来る。以下のようなSASA/VOLモデルがある。以下の式の\(g({\bf r}_i\)はスイッチング関数である。
\(
\begin{eqnarray*}
V^{VOL}_{solv}({\bf r})=\sum_i σ^{VOL}_i\cdot g({\bf r}_i) \cdot \frac{4}{3}\pi R^3_i
\end{eqnarray*}
\)
実際に計算する方法として、POPSやPOPSCOMPなどがある。またSASAの代わりに単に溶媒和エネルギーを部分体積に分割する方法としてEEF1やEEF1-SBなどがある。\begin{eqnarray*}
V^{VOL}_{solv}({\bf r})=\sum_i σ^{VOL}_i\cdot g({\bf r}_i) \cdot \frac{4}{3}\pi R^3_i
\end{eqnarray*}
\)
\(
\begin{eqnarray*}
\Delta G_{i}^{sol}=\Delta G_{i}^{ref}-\sum_{j\neq i}f_i({\bf r}_{ij})\cdot VOL_j
\end{eqnarray*}
\)
\begin{eqnarray*}
\Delta G_{i}^{sol}=\Delta G_{i}^{ref}-\sum_{j\neq i}f_i({\bf r}_{ij})\cdot VOL_j
\end{eqnarray*}
\)
Poisson-Boltzmann
静電ポテンシャルは以下のようになり、
\(
\begin{eqnarray*}
&&\nabla\cdot[\varepsilon({\bf r})\nabla\Phi({\bf r})]=-\frac{\rho({\bf r})}{\varepsilon_0}\\
&&\Phi({\bf r})=\frac{1}{4\pi\varepsilon_0\varepsilon}\int \frac{\rho({\bf r})}{|{\bf r}|}dV({\bf r})
\end{eqnarray*}
\)
Poisson-Boltzmannは特殊な状況を考えており、平衡状態において電荷はボルツマンエネルギー分布に従うとすると、\begin{eqnarray*}
&&\nabla\cdot[\varepsilon({\bf r})\nabla\Phi({\bf r})]=-\frac{\rho({\bf r})}{\varepsilon_0}\\
&&\Phi({\bf r})=\frac{1}{4\pi\varepsilon_0\varepsilon}\int \frac{\rho({\bf r})}{|{\bf r}|}dV({\bf r})
\end{eqnarray*}
\)
\(
\begin{eqnarray*}
\rho_B({\bf r})=\sum_ic^b_iz_iqe^{-z_iq\Phi({\bf r})/k_BT}
\end{eqnarray*}
\)
これを上の式に代入すると、\begin{eqnarray*}
\rho_B({\bf r})=\sum_ic^b_iz_iqe^{-z_iq\Phi({\bf r})/k_BT}
\end{eqnarray*}
\)
\(
\begin{eqnarray*}
&&\nabla\cdot[\varepsilon({\bf r})\nabla\Phi({\bf r})]=-\rho_B({\bf r})-\sum_ic^b_iz_iqe^{-z_iq\Phi({\bf r})/k_BT}
\end{eqnarray*}
\)
\((e^{-z_iq\Phi({\bf r})/k_BT}\ll1)\)という条件で上記の式は\((1-z_jq\Phi({\bf r})/k_BT)\)のように線形化することできる。よって、\begin{eqnarray*}
&&\nabla\cdot[\varepsilon({\bf r})\nabla\Phi({\bf r})]=-\rho_B({\bf r})-\sum_ic^b_iz_iqe^{-z_iq\Phi({\bf r})/k_BT}
\end{eqnarray*}
\)
\(
\begin{eqnarray*}
&&\nabla\cdot[\varepsilon({\bf r})\nabla\Phi({\bf r})]=-\rho_B({\bf r})+\sum_ic^b_iz_i^2\frac{q^2\Phi({\bf r})}{k_BT}
\end{eqnarray*}
\)
\begin{eqnarray*}
&&\nabla\cdot[\varepsilon({\bf r})\nabla\Phi({\bf r})]=-\rho_B({\bf r})+\sum_ic^b_iz_i^2\frac{q^2\Phi({\bf r})}{k_BT}
\end{eqnarray*}
\)
Generalized Born
溶質が半径αの球であるとすると、内部の誘電率を\(\varepsilon_{int}\)と外部の誘電率を\(\varepsilon_{ext}\)とすると、
\(
\begin{eqnarray*}
\Delta G_{solv}=-\frac{1}{2}\left(\frac{1}{\varepsilon_{int}}-\frac{1}{\varepsilon_{ext}}\right)\frac{q^2}{\alpha}
\end{eqnarray*}
\)
α(Born effective radius)は\(\Delta G_{solv}\)と合うように決定されるパラメータである。電荷を持つi,j間の溶媒和エネルギーは\begin{eqnarray*}
\Delta G_{solv}=-\frac{1}{2}\left(\frac{1}{\varepsilon_{int}}-\frac{1}{\varepsilon_{ext}}\right)\frac{q^2}{\alpha}
\end{eqnarray*}
\)
\(
\begin{eqnarray*}
\Delta G_{solv}=-\frac{1}{2}\left(\frac{1}{\varepsilon_{int}}-\frac{1}{\varepsilon_{ext}}\right)\frac{q_iq_j}{\sqrt{r^2_{ij}+\alpha_i\alpha_je^{-r^2_{ij}/4\alpha_i\alpha_j}}}
\end{eqnarray*}
\)
GROMACS,AMBER,NAMDにそれぞれパッケージとして入っている。NAMDのパッケージはGPU/CPUのハイブリッドで計算できる。
\begin{eqnarray*}
\Delta G_{solv}=-\frac{1}{2}\left(\frac{1}{\varepsilon_{int}}-\frac{1}{\varepsilon_{ext}}\right)\frac{q_iq_j}{\sqrt{r^2_{ij}+\alpha_i\alpha_je^{-r^2_{ij}/4\alpha_i\alpha_j}}}
\end{eqnarray*}
\)
Combined PB/SASA,GB/SASA and related methods
AMBERにはHTC、OBCが組み込まれ、CHARMMにはGBMV、GBSWやFACTSが組み込まれている。
Force matchingによるパラメータ化
Implicit solvent modelのパラメータの決定にはいくつかの方法がある。一つには、PBを参照モデルにしてGBのパラメータを決定する方法。また、参照モデルとして実験から得られる溶媒和自由エネルギー差を再現するパラメータを採用する方法。さらに、溶媒と溶質の系において、アンサンブルを再現する方法。最後に、explicitで溶媒を扱った時の溶媒和力？を合わせる方法。
最後のforce-matchingによる方法として
\(
\begin{eqnarray*}
σ^{SASA}_i=-\frac{\partial A_i/\partial r_i}{|\partial A_i /\partial r_i|^2}\cdot \langle f_i^{expl}\rangle
\end{eqnarray*}
\)
\( \langle f_i^{expl}\rangle\)はexplicitの溶媒和力を表しており、\(\partial A_i/\partial r_i\)は距離依存的に変化する溶媒に晒されている表面を表している。
\begin{eqnarray*}
σ^{SASA}_i=-\frac{\partial A_i/\partial r_i}{|\partial A_i /\partial r_i|^2}\cdot \langle f_i^{expl}\rangle
\end{eqnarray*}
\)
最後に
Implicit溶媒モデルの大まかな分類は図ようになる。
