前回の記事でGoogle Colab上で分子動力学（MD）計算を行うMaking-it-rainで遊んでみました。計算したい系のパラメータをいれるだけなのでとても便利です。
・・・「でもそのパラメータが何かよくわからない！」「なんならMD計算が何かもよくわかってない！」
というわけで、パラメータの理解に必要そうなMD計算の基本的な用語をお勉強したメモです。せめてパラメータが何を言っているかわかりたい！
参考資料
参考にさせていただいた書籍はこちらです。
いろいろとググった結果、本棚にあったこちらの書籍に戻ってきました。2009年出版の書籍ですが、知識ゼロから読めて「なんかこの単語見たことあるぞ！」って状態になれるので個人的にとてもオススメです。*1
また、以下のフリーでアクセスできる資料も参考にさせていただきました。*2
- 古明地 勇人, 上林 正巳, 長嶋 雲兵 生体分子の分子動力学シミュレーション(1)方法
Journal of Chemical Software 2000(6)1-36 - 古明地 勇人, 田島 澄恵, 原口 誠, 高橋 伸幸, 上林 正巳, 長嶋 雲兵 生体分子の分子動力学シミュレーション (2) 応用
Journal of Chemical Software 2001(7)1-28
1. MD計算の基本
1-1. MD計算って？
とりあえずWikipedia！
分子動力学法（ぶんしどうりきがくほう、英: molecular dynamics、MD法）は、原子ならびに分子の物理的な動きのコンピューターシミュレーション手法である。（分子動力学法-Wikipedia）
というわけで、分子の動的な性質がわかります。
タンパク質を対象としたMD計算は、1977年カープラス先生らの論文に端を発するそうです。 *3
ウシ膵臓トリプシン阻害剤（Bovine Pancreatic Trypsin Inhibitor, BPTI）の9.2ps（ピコ秒）のシミュレーションで、これにより「タンパク質は動的な系である」という概念を提唱し、それまでの「タンパク質は固い」という概念を転換したそうです。新しい手法を適応してものの見方を変えてしまう研究、2013年のノーベル化学賞受賞につながっているのも納得ですね。*4
さて、この記事ではMD計算の用語について大体下図みたいに追っていきます。
1-2. 動的な性質がなぜわかるの？
まず「どうやって動的な性質がわかるか？」という点ですが、「ある条件において、立体構造が時間に依存してどのように変化するか？」シミュレーションを行います。この構造変化の時系列のことをトラジェクトリ（軌跡, trajectory）とよびます。
より具体的には、ある時点での「立体構造の状態（座標、速度）」をもとに、ニュートンの運動方程式を解いて次の時点で立体構造の状態を導きます。この作業を繰り返していくことで各ステップの構造のスナップショット（"静止画"）が得られます。スナップショットを繋げば動画になるので、動的な状態を観察することができます。なお、ステップ間の時間間隔のことを時間刻み（time step）とよぶそうです。
1-3. 運動方程式で位置と速度を求める
MD計算ではニュートンの運動方程式を解いて、次のステップの構造（座標、速度）を求めます。図のように「段階的に時間が進む」ことを時間発展（time evolution）とよび、「位置や速度を１ステップごとに求めて時間発展する」ことを時間積分（time integration）とよびます。
上図のように、ある時刻の位置と速度は一つ前のステップの情報をもとに求まります。新しい位置と速度で構造を更新していくことで構造変化の時系列（トラジェクトリ）が得られます。
ただし、上記は２次の微少量（dvdt）を無視しているので、コンピュータによる数値計算では時間発展の間に誤差が蓄積して計算が破綻するそうです。誤差の蓄積を防ぐ時間積分法として「差分方程式から質点の位置と速度を求める方法」があり、Verlet法や速度Verlet法、蛙飛び（leap-frog）法といった方法があります。
1-4. 分子力場
時間積分法で速度を求める際に各原子に働く力（F）を利用していますが、これは分子力場（Force field）から求めることができます。
力のベクトルはポテンシャルエネルギー（U)の勾配（偏微分してマイナス（-∇U））で求まります。古典論に基づくMD計算では、原子の座標を関数に代入することでポテンシャルエネルギーもとめます。分子のポテンシャルエネルギーを与える関数とその係数が分子力場です。
分子力場の例としてはAMBERやCHARMM、GROMOSといったものがあります。たとえばAMBER力場の関数形式は以下のような感じです（AMBER (分子動力学)-Wikipedia）。
５つの項のうち、最初の３つは原子間の結合に関係する項（結合項）で、残り２つは結合の有無に関わらず関係する項（非結合項）です。
MD計算では、基本的に各原子が剛体球とばねで繋がれているとするモデルを前提としています。高校物理でよくみたやつですね！（格好良くいうと調和振動子？）
結合項のうち距離(r)や角度(θ)は調和ポテンシャルとなっていて、理想的な状態（平衡結合距離(r_0)、平衡結合角(θ_0))からズレると２乗に比例してエネルギーが大きくなります。また、二面角は回転ポテンシャルに基づいていて、「n : 二面角を作る４原子の対称性」「γ : 位相」です。
非結合項はいずれも距離に依存しますが、減衰の速いファンデルワールス・エネルギーと比較して、静電エネルギーは距離の逆１乗に比例するため減衰が遅く、長距離の相互作用にも寄与します。
非結合項では（結合項と異なり）計算する原子のペア数が系の大きさに応じて膨大に増えるため、計算時間に大きく影響するという課題があります。そこで、計算を簡単にするため「一定の距離（カットオフ, cut-off）以内にある場合のみを考慮する」（カットオフ法）という方法が以前はよく用いられたそうです。ところが、上記の通り「静電エネルギーは長距離にも影響する」ので、一定の距離でバサッと切ってしまうと相互作用を無視してしまい、計算精度が低下してしまいます。現在では、これを改善するための効率的で精度の高い手法が色々と実装されているそうです。
1-5. 境界条件と(長距離)静電相互作用の計算
静電相互作用の計算の課題について述べたので、もう少し解決方法を眺めます。まず、前提知識として境界条件にふれておきます。
1-5-1. 境界条件
MD計算の計算系には孤立系と周期系があります。孤立系では「計算系の外は真空」となっており、周期系ではずらっとくっついて並んだ「周期的な容器（セル）」を考えます。
孤立系で溶媒を明示的に扱う（explicit solvent）モデルでは、球境界条件という「球状の容器の中にタンパク質やリガンド、水分子を閉じ込める方法」がよく用いられるそうです。この系では境界付近で水の挙動が不自然になってしまうという課題があります。
周期系では「周期的な箱」のなかに分子を詰めます。中心を基本セル、それ以外をイメージセルとよび、セルの形状は立方体や直方体が一般的だそうです。
上図のようなイメージで、セルから外に飛び出た分子は反対側からセルに入ってきます*5。周期性を取りいれることで、「実験系よりも圧倒的に少ない原子数」を取り扱うだけで、「マクロな系の性質を求められる」というのが基本的な考え方のようです。周期的境界条件では、球境界条件の境界付近における不自然な挙動が解消できるという利点があります。
1-5-2. （長距離）静電相互作用の計算
先に見た通り、静電相互作用は考えるペアの組み合わせが多く、MD計算で計算時間がかかる部分で、なんとかして効率的に解きたい問題です。
基本的な解決方法の考え方は「近いところ」と「遠いところ」に領域を分けて、寄与の大きい前者は真面目に厳密計算し、少ない後者は近似計算することです。極端なものは「遠いところ＝寄与ゼロ」とするもので、先のカットオフ法です。
もう少し真面目に計算しよう！ということで出てきたEwald法です。Ewald法は周期的境界条件の周期性をうまく利用した高速計算方法です。うまいこと式を分割して収束を速くする感じです。具体的な式は書籍やWikipedia（エバルトの方法-Wikipedia）をご参照ください。*6
Ewald法をさらに高速化したものとして、Particle Mesh Ewald法（PME法）があります。こちらは波数空間の式（E_wave）を「電荷分布をグリッド上に内挿する形」に書き換え、高速フーリエ変換(FFT)することで高速化しています。
Ewald法は周期境界条件を利用したものでしたが、球境界条件でも適用できる計算方法として多重極展開法があります。「遠いところ」にある電荷を持つ点群を多重極展開により近似する方法で、代表点との相互作用のみの計算に単純化することで、考慮する相互作用ペアを大幅に削減することができます。この手法を使った例に、セルマルチポール法（Cell Multipole Method, CMM）やファストマルチポール法（Fast Multipole Method, FMM）といったものがあるそうです。*7
以上が計算コストの高い静電相互作用についての工夫でした。
1-6. マクロな視点で（統計アンサンブル）
ここまでで、MD計算では運動方程式を解いて系の逐次変化をたどっていくこと、運動方程式の力はポテンシャルエネルギーから決まり、エネルギーを求めるために分子力場が使われることがわかりました。最後にもう少しマクロな視点でのお話にふれておきます。
多くの実験は定温、定圧条件でおこなわれ、また、生体内で起こる現象も近似的には定温・定圧条件で起こります。これらと比較し生物学的な現象を解き明かす上で、MD計算結果をより自然（リアル）な系のシミュレーションとするために、さまざまな系の制御が行われます。
1-6-1. アンサンブル
「運動量や座標のある値の組」で指定される系の状態を微視的状態とよび、「多数の微視的状態の集まり」をアンサンブルとよびます。MD計算ではトラジェクトリがアンサンブルに対応します。
アンサンブルとしては以下があります。
先に見た通り、定温や定圧といった条件が意味のあるシミュレーションとして大事になってくるので、さまざまな温度制御や圧力制御の工夫によりカノニカルアンサンブル、NPTアンサンブルを得ることが目指されます。
1-6-2. 温度制御（カノニカルアンサンブル）
系の温度は「系を構成する原子の運動エネルギーの総和」で決まるので、原子の速度が制御できれば温度を制御できることになります。なかでも自然な速度分布に制御するには、統計力学的に妥当な条件を満たす必要があります。
いろいろな温度制御法が提案されていますが、共通して「仮想的な巨大熱浴と系との相互作用を想定して系の温度を制御している」という特徴があるそうです。
温度制御法の例はこんな感じです。*10
難しくてよくわからなかったので正しい情報は書籍や、自然科学研究機構 分子科学研究所の奥村先生が研究室HPで公開してくださっている資料（分子シミュレーション研究会会誌 アンサンブル 連載解説「分子動力学シミュレーションにおける温度・圧力制御」）などをご参照ください。*11
1-6-3. 圧力制御（NPTアンサンブル）
温度制御では速度が大事でしたが、圧力制御では系の体積を増減させることで系の圧力がコントロールされます。代表的な圧力制御の方法としてアンダーソン法があります。温度制御の能勢法（能勢ーフーバー法）と組み合わせた能勢-アンダーソン法では、NPTアンサンブルが実現できます。
・・・詳細はよくわかりませんでした。
以上、アンサンブルについてでした。細かな理論については理解できませんでしたが、MD計算で得られるトラジェクトリ（アンサンブル）が意味のあるもの（熱力学的な平衡状態からのサンプリング）とするために、一定の条件を保つよう制御する工夫が行われている、ということがわかりました。
1-7. 古典MDと第一原理MD（と機械学習MD）
以上でMDの基本的な用語はおしまいです。ところで、ここまで見てきたものは古典分子動力学とよばれるもので、ニュートンの運動方程式（古典力学）に従い、経験的なポテンシャル（分子力場）を利用して原子間の相互作用を求めています。
これに対して、原子間相互作用を量子力学に基づいて計算する第一原理分子動力学というものもあるそうです。量子力学で電子の状態を求め、原子核、電子のすべての相互作用を考慮して原子どうしの相互作用を計算します。このような量子論に基づくMDは1985年に提案されたカー - パリネロ（Car - Parrinello, CP）法に遡るそうです。
古典MDでも「スパコンを使って計算！」みたいな大変そうな感じなのに、時間ステップごとに第一原理計算をして力を求めるなんてめちゃくちゃ大変そうですね。
「計算が大変なところには機械学習がやってくる！」というわけで、この力の計算部分を機械学習に基づくポテンシャル（ニューラルネットワークポテンシャル）に置き換えて計算を高速化した機械学習分子動力学というのもあるそうです。*12
また、第一原理MDと厳密に等しい結果が得られる自己学習ハイブリッドモンテカルロ（SLMC）法というのもあり、機械学習MDほど計算コストは下がりませんが、学習結果に依存せずに高い計算精度が得られるという利点があるそうです。*13
2. MD計算の流れと設定
2-1. MD計算の流れ
MDの基本的な内容がだいたいわかったので、具体的な計算の流れをたどりたいと思います。だいたいこんな感じです。*14
シミュレーションボックスの準備で「系の中和」と書いていますが、これは周期的境界条件でシミュレーションを行う際に扱う時に、電荷が中和されていないと静電相互作用の計算がうまくいかなくなることに関係しているそうです。
2-2. 力場の選択
分子力場の例としてAMBERやCHARMM、GROMOSといったものがあることをみました。*15
ただし、それぞれの力場の中でも、どのような種類の分子を対象に計算を行うかによって選択肢が変わってくるようです。例えば、Amberの Webサイトの説明では以下のようになっていました。
2-3. 明示的な水分子モデル
水分子の配置で明示的、非明示的とかきましたが、これは水も分子として扱うか、それとも連続体モデルのように扱うか、という選択です。
明示的な水分子としてはSPC*16、TIP系（TIP3P、TIP4P、TIP5P）、OPC*17といったモデルがあり、Wikipediaの水モデル項にも詳しく説明されています。
TIP系の数字（3、４、５）は、それぞれ「いくつの点をつかって水分子をモデル化するか？」と関係します。単純に考えると、酸素原子１つと水素原子2つの3点モデル（TIP3P、SPC）となりそうですが、4点モデルではH-O-H角の垂直二等分線上の点Mもさらに考慮し（TIP4P）、5点モデルでは酸素原子ローンペア方向の2点Lを考慮する（TIP5P）モデルとなっています。
以上が、MDの基本でふれていなかった設定項目についてです。
3. Making-it-rainのパラメータ振り返り
MDでやりたいことと、計算の流れがだいたいつかめたので、前回の記事で遊んだMaking-it-rainチュートリアルのパラメータを振り返ってみましょう！
3-1. シミュレーションボックス（系）の準備
まずは「①シミュレーションボックス（系）の準備」に相当する段階で、Amberトポロジーファイルを作成している部分です。
タンパク質対象なので力場（Force_field）として「ff19SB」を選択しているのは先にAmberのマニュアルで見た通りですね。また、水分子のモデル（Water_type）としては明示的な「OPC」を選択しています。
選択項目の中にありませんが、コードの中身を見るとMaking-it-rainは静電相互作用の計算として「PME（Particle Mesh Ewald）」を採用しているようです。ですので、境界条件は周期的境界条件となっていて、「Size_box」でセルの大きさを設定することになります。またEwald法と関連するので系を中和するための設定（Ions, Concentration）も入っています。
3-2. エネルギー最小化と平衡化
続いて「②エネルギー最小化」と「③平衡化」に相当する段階です。パラメータ設定はEquilibration protocolでまとめられています。
平衡化シミュレーションの時間を設定する項目で、「Time」の他に「Integration_timestep」があります。MD計算で運動方程式を解くときは、微少時間の時間積分を考えるという話だったので、微少時間に相当するパラメータとなりそうです。（上図では0.2nsのシミュレーション中、2fs間隔で積分計算を行う感じ？）
また、計算で得られるアンサンブルに関係するパラメータとして温度と圧力があります。Making-it-rainはNPTアンサンブルのみとなっています。
MD計算で欲しいのは「一定の時間間隔ごとの構造のスナップショット（トラジェクトリ）」ですが、「Write_the_trajectory」がその時間間隔のパラメータになっています。運動方程式の計算はフェムト秒間隔ですが、すべての計算結果を保存したらデータ量が膨大になるので、ピコ秒オーダーで出力する感じでしょうか？
3-3. シミュレーション本番
最後に「④シミュレーション本番」に相当する部分です。ここのパラメータは平衡化の部分とほぼ一緒です（Making-it-rain特有の時間設定パラメータ（Stride）が増えているくらい）。
ところで、Making-it-rainはNPTアンサンブルと書きましたが、「どのように温度・圧力制御を行うか？」の設定の選択肢はありませんでした。
コードの中身を見ると、MonteCarloBarostat(pressure, temperature)やLangevinIntegrator(temperature, friction, dt)というそれっぽい箇所がみられます。Making-it-rainがMD計算に使用しているOpenMMのAPIを参照してみましょう。
MonteCarloBarostatは「モンテカルロアルゴリズムで周期的なセルのサイズを調整し、一定の圧力の効果をシミュレーションする」とあります。これが圧調整器（バロスタット）みたいですね。ただし「温度を維持する機能はないので、LangevinIntegratorやAndersenThermostatといったものと合わせて使うように」とあります。
また、LangevinIntegratorは、「ランジュバン動力学を使って系のシミュレーションを行う積分器」とあります。「ランジュバン動力学 - Wikipedia」によると、「サーモスタットを用いた時のように温度を制御できるため、正準集団（カノニカルアンサンブル）を近似できる。」とのことですので、こちらが一定の温度に保つための温度調整器（サーモスタット）の役割を果たしているようです。
ところで、先に見た温度制御の例で、「サーモスタットはだいたい速度をコントロールしている」「ランジュバン方程式類似で摩擦力を使う」といった内容を見ました。LangevinIntegratorの引数に摩擦係数（friction coefficient）があるのを見ると、確かに摩擦力で速度制御して温度をコントロールしてそうだな、という印象です。*18
以上、Making-it-rainのパラメータ振り返りでした。MD計算の考え方や手順をぼんやり眺めた後だと、パラメータの意味が少しわかった気がします！
4. Making-it-rainトラジェクトリ解析の振り返り
Making-it-rainには基本的なトラジェクトリ解析もありましたが、これらもさっぱり意味がわかりませんでした。ついでに振り返っておきましょう。
一般的な解析手法について素晴らしい資料がありました！理化学研究所 計算科学研究センターのeラーニング「計算生命科学の基礎６ 生体系分子シミュレーションの新展開」です！
横浜市立大学 池口満徳先生（研究室HP）の講義動画（YouTube）とスライド資料（SlideShare）を公開してくださっています。うれしい！*19
4-1. RMSDとRMSF
まずは、根平均二乗変位(Root Mean Square Deviation, RMSD)と根平均二乗揺らぎ(Root Mean Square Fluctuation, RMSF）です。どっちもズレを評価している感じ。
RMSDについてMaking-it-rainの結果と説明を並べてみるとこんな感じ。
上の例だと、PDBから取得した初期構造（参照構造）からのズレを可視化できていることになりそうです。また、トラジェクトリ前半でRMSDは増加傾向にあるものの、後半では一定の範囲に落ち着いているので、確かに収束具合がわかります。
RMSFについてはこんな感じ。
RMSDとRMSFの用語が似すぎていて違いがよくわかっていませんでしたが、参照構造との比較か平均構造との比較かという違いだったんですね。横軸が違う（time, residue）意味もやっとわかりました。
4-2. 回転半径
回転半径はこんな感じ。
実験で得られるデータと関係しているというのは知りませんでした。上記の例は時間も短く、フォールディングの検証でもないのであまりこの指標には合わないのかもしれません。
4-3. 主成分解析（PCA）
PCAはこんな感じ。
タンパク質を相手にPCAを行う意味がよくわかっていなかったですが、揺らぎの解析だったんですね。平面プロット（PC1, PC2）で近くにあるCαは、一緒に動いている残基の可能性が高いということでしょうか？？
4-4. 相関行列
Making-it-rainではPearson's Cross Correlation(CC)となっていましたが、出力結果を見た感じ相関行列(Correlation Matrix)や動的相関図（Dynamical Cross Correlation Map, DCCM）と呼ばれているものと同じでしょうか？
この解析にはいくつか問題点も指摘されているそうです。式の通り２つの残基の平行な協同運動しか観測できず、また中心付近のコアドメイン内の協同運動を観測することもできません。また、揺らぎの相関なので収束が遅く、再現性が低いそうです。
詳しい分子内協同運動を調べるには、相関行列よりも先の主成分解析やEssential Dynamics(ED)といったものの方が良いそうです。
以上、トラジェクトリ解析の振り返りでした。解釈は難しいですが、とりあえず何を目的とした解析なのかを把握できたのでヨシ！
5. おわりに
以上、「今回はMD計算のパラメータがわからない！っていうか、そもそもMD計算て何をやっているんだ？？」ということで、MD計算の基本的な内容とパラメータ等の用語のお勉強をしてみました。
ついでに、前回の記事で遊んだMaking-it-rainのチュートリアルを振り返って、パラメータや操作の意味・解析方法について復習してみました。
「計算がうまくいっているのか？」「妥当な結果なのか？」といったところまではわかりませんでしたが、とりあえず何をやろうとしているかぼんやりわかったのでヨシ！
書籍の参考部分が多くなってしまったので、最後にもう一度おすすめしときます。
タンパク質計算科学 ―基礎と創薬への応用―（CD-ROM付き） / 神谷 成敏 肥後 順一 福西 快文 中村 春木 著 | 共立出版
今回も色々と間違いが多そうなのでご指摘いただけると嬉しいです。（初心者向けの本もお薦めいただけると嬉しいです。）
ではでは！
*1:いろいろな資料を参考にさせていただいていますが、ベースは書籍を参考にしているので先に宣伝しときます。PDBファイルの読み方を含めてタンパク質をコンピューターで扱う方法がざっとさらえて素敵な書籍です。基本的な内容は変わっていないと思うので、まだまだ現役の書籍とは思いつつ、アップデートを期待している書籍でもあります。（もちろん私は全くの部外者で、アフィリエイトでもないので（安心して？）リンクをクリックしてしてください)。なお本記事中の図表の作成にあたって参考にさせていただいた書籍の対応箇所はキャプションに記載しています。
*2:こっちはさらに前の2000年、2001年の資料です。日本語の素晴らしい資料が20年経ってもフリーでアクセスできるのは最高ですね！
*3:McCammon, J. A., Gelin, B. R. & Karplus, M., Dynamics of folded proteins Nature 267, 585-590 (1977).
*4:なお、この記念すべきシミュレーションの動画はカープラス研究室HP（Karplus group HP）でみることができます。
*5:SF映画でドアを開けたら元の部屋につながっている感じ？？
*6:フーリエ級数という言葉が出てきた瞬間に理解を諦めました。
*7:湯どうふさんのYouTube動画が簡潔でわかりやすいです。多重極展開法のYouTube動画はこちらで、セルマルチポール法のYoutTube動画はこちらです。
*8:東京大学大学院農学生命科学研究科 アグリバイオインフォマティクス教育研究ユニット 講義「分子モデリングと分子シミュレーション」 2019年5月23日の資料「分子動力学法とモンテカルロ法」の中の図表と記述を参考にしました
*9:イメージはいらすとやさんの図から
*10:下図の特徴に書いている部分がわかりにくいのは私が全く理解していないからです。ごめんなさい。
*12:たとえばAENET（The Atomic Energy Network）というオープンソースソフトウェアがあるそう（生体分子というよりもマテリアル寄り？）。
*13:原子力機構システム計算科学センター 永井佑紀先生が講演資料をGitHubで公開してくださっていて面白かったのでご紹介しておきます（スライドPDFのURL(GitHub) 。
*14:参考：@Ag_smith先生 Qiita記事 MDシミュレーションのチュートリアル〜PDB: 1LKEの場合〜 より
*15:Wikipediaの「力場（化学）)」項に他にもいろいろな力場の例が載っています。
*16:単純点電荷, simple point-chargeの略？
*17:Oはoptimal？参考: J. Phys. Chem. Lett. 2014(5)3863–3871
*18:ランジュバン方程式・動力学が全くわかっていないので間違っていたらごめんなさい。
*19:以下の説明では池口先生の資料と、冒頭にあげた古明地先生らの資料の記述を参考に引用させていただいています。
