粗視化力場を最適化するための遺伝的アルゴリズム(2013)

以下に書いたものは、論文を基にした私の理解なので、正確性を欠きます。正確な情報を知りたい場合は、下記の論文を参照してください。また、下記の論文の図を参照することで、理解が深まると思います。

Evolutionary Algorithm in the Optimization of a Coarse-Grained Force Field

前置き

　計算機科学の発展に伴い、分子動力学法で扱える空間スケール、時間スケールは大きく広がった。しかし、それには限界が有る。なぜなら、対象とする系は如何様にも大きくできるからである。そのような限界を突破するための一つの方法として、粗視化分子動力学法がある。この方法は分子動力学法において、全原子の相互作用を計算するのをやめ、原子のグループを作り、それを一つの粒子として扱う方法である。このような自由度を削減する方法により、計算コストを大幅に減らすことが出来る
　しかし、上記のような方法は、粗視化力場の新たな構築が必須となる。その作業は大変でつまらないものであり、最適な力場を決定できているかわかない。そのため、本論文では、そのパラメータを自動的に決定するために遺伝的アルゴリズムを導入した。遺伝的アルゴリズムとは、生物の遺伝システムを模したアルゴリズムであり、その過程は選択、交差、突然変異を含んでいる。これらの過程により、次世代では、前世代の結果より良い結果を生み出すことができるアルゴリズムである。

前置きの前置き

　粗視化力場を最適化する方法として、遺伝的アルゴリズムがメタヒューリスティックな(メタヒューリスティクスとは特定の問題に限定されず、どのような問題に対しても汎用的に対応できるように設計された、アルゴリズムの基本的な枠組みのことである。by wikipedia)アルゴリズムに対して、他に解析的な方法とシステマティックな局所探索法がある。

解析的な方法

• Boltzmann iversion(BI)
• Force-matching method

システマティックな局所探索法

• CG-OPT
• GROW
• Renormalization group method

方法

本論文では、1ヌクレオチド=1粒子として粗視化したRNAを対象とし、パラメータの最適化を試みた。粗視化粒子の位置は、リン原子に配置した。またその際のRNAのポテンシャルは以下のようにした。

\(\begin{eqnarray*}
&&U=U_b+U_{\theta}+U_{\phi}+U_{bp}+U_{nb}\\
&&U_{bp}=U_{i:j}+U_{i:j+1}+U_{i+1:j+1}
\end{eqnarray*}\)

これらのポテンシャルの詳細は、まず分子内のポテンシャルから

\(\begin{eqnarray*}
&&U_b(r)=\frac{1}{2}k_r(r-r_0)^2\\
&&U_{\theta}(r)=\frac{1}{2}k_{\theta}(\theta-\theta_0)^2\\
&&U_{\phi}(\phi)=k_{\phi}(1-cos(\phi-\phi_0))
\end{eqnarray*}\)

次に、ベースペア間のポテンシャルは2つ存在し、以下のように定義する。

\(\begin{eqnarray*}
&&U_{i:j}(r)=k^{i:j}\frac{1}{2}(r-r^{i:j}_0)^2\\
&&U_{i:j}(r)=(U^{i:j}_{Morse}(r)-U^{i:j}_0 c^{i:j})sw^{i:j}(r)\\
&&U^{i:j}_{Morse}(r)=U^{i:j}_0[1-exp(-\alpha^{i:j}(r-r^{i:j}_0))]^2\\
&&sw^{i:j}(r)=\frac{1}{2}[1-ｔａｎｈ (\lambda^{i:j}(r-r^{i:j}_1))]
\end{eqnarray*}\)

\(U^{i:j}_0 c^{i:j}\)は、相互作用している時としていない時のエネルギー差で、\(\lambda_{i:j}\)はスイッチ関数をコントロールするパラメータである。
次に、ノンローカルな相互作用として、静電相互作用と疎水性相互作用と反発力を表現したポテンシャルは

\(\begin{eqnarray*}
&&U_{nb}(r)=\frac{q_1q_2}{4\pi\varepsilon_0\varepsilon_rr}
\end{eqnarray*}\)

\(
\begin{eqnarray*}
U_{nb}(r) = \begin{cases}
k(r-r_0)^2+U_{nb}^2 & if\ r＜r_0\\
(U_{nb}^{bar}-U_{nb}^0)\exp(-σ_1(r-r_1)^2)+U_{nb}^2 & if\ r_0＜r＜r_1\\
U_{nb}^{bar} & if\ r_1＜r＜r_2\\
U_{nb}^{bar}\exp(-σ_2(r-r_2)^2) & otherwise
\end{cases}
\end{eqnarray*}
\)

Fitnes Function.この関数を使って、CG-FFと実験結果やall-atomの結果の一致度を調べる。

\(\begin{eqnarray*}
&&f_{tot}=\sum^n_{i=1}w_if'_i
\end{eqnarray*}\)

\(
\begin{eqnarray*}
f'_i = \begin{cases}
0 & if\ f_i＜l_i\\
\frac{u_i-f_i}{u_i-l_i} & if\ l_i＜f_i＜u_i\\
1 & if\ f_i＞u_i
\end{cases}
\end{eqnarray*}
\)

以下、具体的なFitnes Functionの具体例を見ていく。Distance Distribution Functions.距離の分布間の距離の比較。

\(
\begin{eqnarray*}
f_{HI} =\frac{2\sum^N_{i=1}d_1(i)d_2(i)}{\sum^N_{i=1}d_1^2(i)+\sum^N_{i=1}d_2^2(i)}
\end{eqnarray*}
\)

RMSD(Root Mean Square Deviation)の比較

\(
\begin{eqnarray*}
RMSD(t)=\sqrt{\frac{1}{N}\sum^N_{n=1}({\bf r}(n,t)-{\bf r}_{ref}(n))^2}
f_{RMSD}=\frac{1}{T}\sum_{t=0}^{T}RMSD(t)
\end{eqnarray*}
\)

RMSF(Root Mean Square Fluctuations)の比較

\(
\begin{eqnarray*}
RMSF(n)=\sqrt{\frac{1}{N}\sum^N_{n=1}({\bf r}(n,t)-{\bf r}_{ref}(n))^2}
\beta(n)=\frac{8\pi^2}{3}RMSF(n)^2\\
f_{RMSF}=\frac{1}{N}\sum^N_{n=1}| RMSF_{simulation}(n)-\sqrt{\frac{3}{8\pi^2}\beta(n)} |
\end{eqnarray*}
\)

Distance Constraints(CG-FFでは、あるボンドに拘束をかけることがある。)

\(
\begin{eqnarray*}
f_{dist}=|\overline{x}-x_{ref}|
\end{eqnarray*}
\)

遺伝的アルゴリズム

1. ランダムにCG-FF作成する。CG-FFの数はポピューレションサイズとなる。
2. 上のCG-FFを使って、シミューレションを実行する。
3. CG-MDの上記のfitness functionを計算する。
4. 2-3のステップを繰り返し、結果を平均化する。
5. fitnes functionを使い、CG-FFをランキングし、1番だけを残し、変異と交叉を行い、新しい世代を作成する。
6. 2からのステップを繰り返す。

思うところ

• 遺伝的アルゴリズムをMDのパラメータ決定に用いた面白い例。2008年にDNAのパラメータ決定をしている。何故。タンパク質のパラメータ決定をしないのか気になるところ。
• 結局、出てきたパラメータがいいのか悪いのかわからない。遺伝的アルゴリズムを繰り返し得られたパラメータが過学習しているかもしれない。
• どうせなら関数の選択にも遺伝的アルゴリズムを適用してみたらいいと思う。
• Go-modelと違いフォールディング中間体を再現しているか不明。Φ値の実験と比較してみると面白いかも。