MDシミュレーションのチュートリアル記事を書いてみました。
今回題材として選んだタンパク質はPDB: 1LKEです。
この構造のポイントは以下の3つです。
- リガンドとしてDigoxigeninが結合している
- TYR-117とGLY-123の間に構造データが存在しないミッシングループが存在し、SER-166の後ろにもいくつか残基が続く（参考：PDB: 1LNM）
- Cys-8とCys-115でジスルフィド結合が形成されている。さらに、Cys-42とミッシング残基のCys-170でもジスルフィド結合が形成されている（参考：PDB: 1LNM）
これらにどう対処すればよいかを知っていればほとんどのタンパク質MDシミュレーションに対応できるからです。
（2018年9月13日） GitHub（ https://github.com/YoshitakaMo/MDsimulations ）に使用した全ファイルと完成済みファイル例を置きました（minimizationまで）。参考にしてください。
では行ってみましょう。
環境
- macOSまたはLinux環境推奨。Windows 10の場合はWindows Subsystem for Linux（通称WSL）をインストールすれば、Ubuntuのような環境を構築することができる。GPUがあると計算が最大10倍速くらいになる。
- AmberTools 22（AMBER20があればそれも）がインストールされている
- Gromacs 2023.2がインストールされている
- Modellerがインストールされている（ミッシング残基が存在する場合）
環境設定についてはそれぞれのリンクを読んで下さい。Modellerについてはすぐ下で述べます。
Modellerのインストール(macOSの場合)
macOSの場合はターミナルからHomebrewをインストールしてある状態で、以下のコマンドを打ちます。
brew tap salilab/salilab
brew install modeller
これだけでインストールできます。
ただし、使用開始する前にファイル /usr/local/opt/modeller/modlib/modeller/config.py にLicense Keyを入れておく必要がある(XXはバージョン。現在は10.1)。vimやVS Codeなどでこのファイルを開いて中身を以下のように書き換える。
install_dir = r'/usr/local/Cellar/modeller/10.1'
license = r'MODELIRANJE'
これで完了です。
Windows 10、WSLまたはUbuntuを使っている方は、ターミナルを開いてModeller公式ページのUbuntu用のインストール方法を試してみてください。
PDBファイルのダウンロード
適当な作業ディレクトリを作成します。ここでは~/Desktop/1LKEとします。まずはそのディレクトリを作り、そこにRCSB PDBのウェブサイトからPDBファイルのPDB: 1LKEをダウンロードしてきます。
# ディレクトリの作成
mkdir -p ~/Desktop/1LKE/PDB
# ディレクトリへの移動
cd ~/Desktop/1LKE/PDB
# PDB:1LKEのファイルをcurlコマンドでダウンロード
curl -o 1LKE.pdb "https://files.rcsb.org/download/1LKE.pdb"
このPyMOLなどの可視化ソフトで構造を見ておきましょう。この構造で確認しておくポイントは以下の3つです。
- リガンドとしてDigoxigeninが結合している
- TYR-117とGLY-123の間に構造データが存在しないミッシングループが存在し、SER-166の後ろにもいくつか残基が続く（参考：PDB: 1LNM）
- Cys-8とCys-115でジスルフィド結合が形成されている。さらに、Cys-42とミッシング残基のCys-170でもジスルフィド結合が形成されている（参考：PDB: 1LNM）
このMDシミュレーションのチュートリアルではこれらの部分の設定方法も含めながら進んでいきます。
AmberToolsを使った1LKEのトポロジーファイル作成
トポロジーファイルとは
MDシミュレーションの基本的な考えは、タンパク質などの生体分子を構成する各原子が剛体球とばねで繋がれているとするモデルを大前提としています。これに加え、3原子間の角度についても同様な球とばねモデルを適用し（例として水分子のH-O-H角度が約104.5°となるように制御するなど）、さらに4原子間の二面角についても設定し、分子の平面性を制御したりしています。また、例えば同じ炭素原子であってもsp3炭素とsp2炭素の振る舞いが違うように、各原子について個別に適切なパラメータを設定する必要があります。
幸いなことに、タンパク質のシミュレーションではタンパク質を構成するアミノ酸はおよそ20種類に限られ、かつパラメータもアミノ酸ごとに同じものを使いまわしても大きな問題にはならないことが判明していますので、これらのパラメータセットは**力場（forcefield）**としてライブラリ化されており、いくつかの信頼と実績のある力場が世界中で広く使われています。この力場を使ってタンパク質のMDシミュレーションの系に対して作成したファイルのことをトポロジーファイルと言います。
トポロジーファイルはタンパク質のMDシミュレーションをする上で必須のファイルであり、パラメータは動きと精度に大きく関わるので非常に重要なものとなります。
Step 1-1. PDBファイルの原子の座標部分だけ取り出す
AMBERとAmberToolsに付属のtleapというモジュールを使うと、PDBファイルからトポロジーファイルを作成することができます。さっそくこれを使ってみましょう。
まずcp 1LKE.pdb 1LKE_truncated.pdbで先程ダウンロードしたファイルのバックアップを確保しておきます。続いてこの1LKE_truncated.pdbファイルをテキストエディタで見てみます。すると、462行目から1793行目まで、ATOM、TER、HETATMで始まる文字列があるのがわかります。この部分だけを確保し、残りは消してしまいましょう。
ATOM 1 N ASP A 5 2.533 -40.650 52.271 1.00 69.29 N
ATOM 2 CA ASP A 5 1.130 -40.378 52.680 1.00 70.86 C
ATOM 3 C ASP A 5 1.018 -39.026 53.360 1.00 66.11 C
ATOM 4 O ASP A 5 1.950 -38.566 54.021 1.00 68.69 O
ATOM 5 CB ASP A 5 0.598 -41.480 53.609 1.00 74.53 C
（中略）
HETATM 1330 O HOH A 556 28.620 -15.856 57.423 1.00 53.96 O
HETATM 1331 O HOH A 557 4.205 -24.295 47.080 1.00 39.36 O
HETATM 1332 O HOH A 558 27.529 -13.304 59.699 1.00 58.19 O
END
として切り出しておきます。これを1LKE_truncated.pdbとして保存します。
Step 1-2. Modellerを使ったミッシング残基の補完（存在する場合のみ）
このステップの動作は、上に引き続き~/Desktop/1LKE/PDBを想定しています。
今回のタンパク質1LKEには一部の残基が見えていないミッシング領域が存在します。
RCSB PDBのウェブサイトでタンパク質構造1LKEのsequenceタブを見ると、TYR-117とGLY-123の間にはDEDKKのアミノ酸が存在することがわかります。
この画像でUNMODELEDの行に注目します。この行は座標情報が存在していないアミノ酸残基の場所を示しています。UNMODELEDの行の118-122の灰色の四角にマウスカーソルを載せると情報が表示されます。
118番から122番までのアミノ酸残基、DEDKKがミッシング（UNMODELED）であることがわかります。この箇所は結晶中でさえもゆらぎが大きいためか、この相当する部位の電子雲がほとんど見えず、その座標部分がPDBファイル内に記載されていません。しかし、このDEDKK配列は必ず存在するはずであり、またMDシミュレーションをする場合には当然ながらこれらを含んだ構造にしている必要があります。
またジスルフィド結合（DISULFIDE BRIDGE）についても注目します。上と同じSequenceタブでマウスをDISULFIDE BRIDGEの部分に載せると、CYS-8とCYS-115の間でジスルフィド結合が形成されていることがわかります。
ジスルフィド結合はタンパク質構造を保つ上でとても重要な働きをするので、必ず入れる必要があります。
ここで、やや発展的な内容ですが、ジスルフィド結合が他にも存在していないか調べてみましょう。RCSB PDBのStructureタブに戻って下の方にはLiteratureという情報欄があります。そこには、PDB: 1LKEに関係するタンパク質がいくつか示されています。
ここにPDB: 1KXO, 1LNMのことが書かれています。こちらのいずれかのPDBデータについてSequenceタブを調べてみると、そちらの構造にはCYS-42, CYS-170間にもう1つのジスルフィド結合が存在することがわかります。（参考：http://www.rcsb.org/sequence/1LNM）
アミノ酸配列が同じであればタンパク質構造の概形は同じなので、こちらでもう1つのジスルフィド結合が見えたのであれば、1LKEでもこのジスルフィド結合が存在している可能性は極めて高いと言えます。そしてジスルフィド結合はタンパク質構造を保つ上でとても重要なので（2回め）、この情報を入れた上でMDシミュレーションをしたいということになります。ただし、このような座標情報がないアミノ酸残基の座標情報を人間が手で適切に補完してあげるのは困難です。そこで機械的に行いたいわけですが、こうしたミッシング領域の座標の補完を行うとき、よく用いられるのがModellerというソフトです。ModellerはSali Labが開発しているPythonベースのタンパク質構造編集ソフトウェアです。使いこなせれば非常に強力にタンパク質構造を編集できます。
今回はModellerで以下の3つの設定を踏まえたモデリングを行います。
- TYR-117とGLY-123の間のミッシング領域についてDEDKKのアミノ酸ループ構造を補完する。
- さらに、SER-166の後ろにもいくつか残基が続くので、EAACKVNのタンパク質構造を補完する（今回はとりあえずCYS-170で形成されるジスルフィド結合周辺の構造を補完するのが狙いです）。
- PDBファイル上のCys-8とCys-115でジスルフィド結合が形成されている。さらに、Cys-42とミッシング残基のCys-170（EAACKVN配列中のC）でもジスルフィド結合が形成されている。
Modellerの実行には実行バイナリmod10.1（10.1はバージョン番号）にPATHが通っていることの他に（which mod10.1とコマンドを入力してみてPATHが表示されれば成功）、Modellerに読み込ませるPythonスクリプトのalign.pyと配列アラインメントファイルのalignment.aliの3つが必要です。
では、まずは以下のようなalign.pyファイルを用意します。
#!/usr/bin/env python3
import modeller
import modeller.environ
import modeller.util.logger
import modeller.automodel
import modeller.parallel
# 並列計算用のコア数を設定(2021年3月17日追加)
# https://salilab.org/modeller/9.15/manual/node457.html
ncpus = 8
j = modeller.parallel.job(host="localhost")
for i in range(ncpus):
j.append(modeller.parallel.LocalWorker())
modeller.util.logger.log.verbose()
env = modeller.environ()
# HETATMの座標を保存するようにフラグを追加
env.io.hetatm = True
# 水分子の座標を保存するようにフラグを追加
env.io.water = True
# directories for input atom files
env.io.atom_files_directory = ['.', '1LKE_truncated.pdb']
a = modeller.automodel.AutoModel(env, alnfile = 'alignment.ali', knowns = '1LKE', sequence = '1LKE_fill')
a.starting_model= 1
a.ending_model = 8
a.md_level = modeller.automodel.refine.fast
# 並列計算の実行（2021年3月17日追加）
a.use_parallel_job(j)
a.make()
現行のmodellerでは、マルチコアCPUを利用して並列に構造モデリングを行う機能が実行されています。これを利用することでたくさんのモデリングを行うときに大幅に計算時間を短縮することができます。
残基番号が4つずれていることに注意してください。これはModellerでは以下のalignment.aliで定義した構造情報を持つアミノ酸残基を必ず1番めに再設定するという性質があるためです（言い換えれば、アウトプットされるPDBの残基番号を基準としてここの数字を設定します）。PDB: 1LKEの構造情報はChain Aの残基番号5から始まっています。
ミッシング領域に含まれるCys-170（完成後のモデル構造上では166）についてもSS結合を定義しています。
次にアラインメントファイルのalignment.aliです。
>P1;1LKE
structureX:1LKE_truncated:5:A::A:undefined:undefined:-1.00:-1.00
DGACPEVKPVDNFDWSQYHGKWWQVAAYPDHITKYGKCGWAEYTPEGKSVKVSRYSVIHGKEYFSEGTAYPVGDSK
IGKIYHSYTIGGVTQEGVFNVLSTDNKNYIIGYFCSY-----GHMDLVWVLSRSMVLTGEAKTAVENYLIGSPVVDSQKL
VYSDFS-------...........................................................*
>P1;1LKE_fill
sequence:::::::::
DGACPEVKPVDNFDWSQYHGKWWQVAAYPDHITKYGKCGWAEYTPEGKSVKVSRYSVIHGKEYFSEGTAYPVGDSK
IGKIYHSYTIGGVTQEGVFNVLSTDNKNYIIGYFCSYDEDKKGHMDLVWVLSRSMVLTGEAKTAVENYLIGSPVVDSQKL
VYSDFSEAACKVN...........................................................*
上の方が構造情報を持つテンプレート配列（PDBの構造と同じに揃える）、下のほうが補完した配列です。必ず配列を対応させておく必要があり、ミッシングに相当する配列は-記号で埋めます。また、align.pyでenv.io.hetatm = True、env.io.water = Trueを定義していた場合、この配列情報中の.記号でそれぞれの残基を代用することができます（この後の計算中のModeller中間ファイル*.iniではBLKとして処理しているようです。）。上の配列ファイルでの.の数は、Digoxigeninが1つ分、水分子が58つで計59個あります。
これの2つのファイルを用意したら、cd ~/Desktop/1LKE/PDBでそのディレクトリを移動し、1LKE_truncated.pdbとalign.py, alignment.aliの3つのファイルが同じディレクトリに存在することを確認した上で、ターミナル上で
$ python3 align.py > align.log
を実行します。
並列計算の数値によりますが、およそ30秒くらいで計算は終了し、8つのモデリング構造がアウトプットされます。
align.logの最後の方に書かれている構造スコアを見てみます。この値が低いほどよい構造とされます。
>> Summary of successfully produced models:
Filename molpdf
----------------------------------------
1LKE_fill.B99990001.pdb 1138.31165
1LKE_fill.B99990002.pdb 1333.03198
1LKE_fill.B99990003.pdb 1128.88806
1LKE_fill.B99990004.pdb 1168.13550
1LKE_fill.B99990005.pdb 1195.39551
1LKE_fill.B99990006.pdb 1367.64893
1LKE_fill.B99990007.pdb 1110.27405
1LKE_fill.B99990008.pdb 1391.83472
モデリングにはある程度の乱数を利用しているため、この中で最も良かった構造の番号は毎回異なります。
値が最もよかったものを選んでPyMOLなどで構造を見てみましょう。またここで、構造を見てみたときに注意してほしい点があります。
このようにモデリングしたループがSS結合の中に入ってしまっている構造がいくつかあります。この構造は明らかに良くないので選ばないようにしましょう。
また、水分子とDigoxigeninがモデル後の構造にも保存されているでしょうか？大丈夫であれば、ここでは最もスコアの小さかった1LKE_fill.B99990007.pdbを選んで話を続けます。
Modellerで作られたこのファイルを開いてみると、残基番号が1番からに振り直されています。ここで（必要というわけではないですが、話をわかりやすくするため）、Gromacsのgmx editconfコマンドを使って残基番号を改めて5番からに振り直します。
# gmx editconfコマンドで、残基番号を5から振り直したものを1LKEmodel.pdbに出力する
gmx editconf -f 1LKE_fill.B99990007.pdb -o 1LKEmodel.pdb -resnr 5
入力-fを先程のモデル構造ファイルに、-oに出力ファイル名を指定し、renumberingを5番からに指定します。こうして得られた1LKEmodel.pdbファイルを使って次のStepに進みます。
Step 1-3. His残基の水素の位置の設定
ヒスチジン（His）残基は図のようなpKa曲線を持っており、pH = 6.0以下で側鎖のイミダゾール基に水素を2つもつ状態になり、それ以上になると水素を1つ、δ位かε位のどちらか1つについた状態になります。MDシミュレーションは通常pH 6.0〜9.0程度の溶液を想定して行うことが多いため、水素原子が1つついた状態のヒスチジンを取り扱うことがほとんどになります。
さて、イミダゾール基のδ位かε位のどちらか1つに水素原子がつくわけですが、水素原子はとても小さいのでタンパク質結晶構造解析であってもなかなか水素原子の位置を捉えることはできません。しかしどちらに設定するべきかというのは、場合によっては重要な要素になります。そこで通常は、タンパク質結晶構造中で、水素結合が形成されていそうな原子の方に水素原子を付加し、水素結合が形成されてそうでなければ適当にどちらか一方（通常ε位）に水素原子を設定するという手法を取ります。
今回使用するAMBERのtleapモジュールの場合、読み込むPDBファイルのHIS残基名をHIDに変更するとδ-protonatedのヒスチジン、HIEにするとε-protonatedのヒスチジンとして認識してくれるようになります。
では1LKE_model.pdbについてやってみます。このタンパク質内には5つのHis残基があります。このうち注目するべきはHis-86です。
His-86をみると、δ位に水素原子があったほうがdigoxigeninと水素結合を形成していい感じになりそうなのがわかります。His-35はε位に水素があったほうが水素結合を形成しやすそうなのがわかります（δ位ってのはイミダゾールの五角形の根元から近い方の窒素原子がδ、遠い方がεです）。他のヒスチジン残基はどっちでもよさそうです。なので、His-86だけは残基名をHIDにしてδ-protonatedということにしておきましょう。
先程の1LKEmodel.pdbファイルのHIS 86あたりの残基名をHIDにします。
ATOM 654 N HID 86 4.953 -19.957 36.923 1.00 79.07 N
ATOM 655 CA HID 86 5.853 -20.098 35.812 1.00 79.07 C
ATOM 656 ND1 HID 86 9.513 -20.269 37.045 1.00 79.07 N
ATOM 657 CG HID 86 8.205 -19.890 36.843 1.00 79.07 C
ATOM 658 CB HID 86 7.209 -20.754 36.129 1.00 79.07 C
ATOM 659 NE2 HID 86 9.292 -18.249 37.946 1.00 79.07 N
ATOM 660 CD2 HID 86 8.086 -18.655 37.401 1.00 79.07 C
ATOM 661 CE1 HID 86 10.118 -110.12 37.707 1.00 79.07 C
ATOM 662 C HID 86 5.229 -21.009 34.812 1.00 79.07 C
ATOM 663 O HID 86 4.616 -22.013 35.174 1.00 79.07 O
ここのHIS残基をすべてHIDに変更して保存したファイルを1LKE_trun_hid.pdbとしておきましょう。あとは自動的にε-protonatedになります。
Step 1-4. SS結合の設定（存在する場合のみ）
今回のケースのように、SS結合が存在する場合はもういくつか設定する必要があります。AmberToolsの場合、SS結合を形成させたい場合は以下の2つを行う必要があります。
- SS結合を形成する2つのCYS残基の名前をCYXに変更する
- SS結合を形成させる2残基のCYXのSGを結合させる設定を書く
1つめは上のHIS→HIDと同じように名前をCYS→CYXに変えるだけです。
ATOM 15 C ALA A 7 -0.641 -33.155 52.030 1.00 49.97 C
ATOM 16 O ALA A 7 0.142 -33.209 52.966 1.00 42.23 O
ATOM 17 CB ALA A 7 -3.095 -33.395 52.296 1.00 58.78 C
ATOM 18 N CYX A 8 -0.532 -32.270 51.035 1.00 46.12 N
ATOM 19 CA CYX A 8 0.517 -31.252 51.058 1.00 43.44 C
ATOM 20 C CYX A 8 0.439 -30.449 52.312 1.00 46.59 C
ATOM 21 O CYX A 8 -0.638 -29.981 52.691 1.00 51.88 O
ATOM 22 CB CYX A 8 0.410 -30.329 49.847 1.00 38.64 C
ATOM 23 SG CYX A 8 0.925 -31.183 48.413 1.00 36.24 S
ATOM 24 N PRO A 9 1.577 -30.296 52.978 1.00 45.14 N
ATOM 25 CA PRO A 9 1.646 -29.492 54.192 1.00 47.94 C
ATOM 26 C PRO A 9 1.342 -28.044 53.866 1.00 53.73 C
（中略）
ATOM 880 CE1 PHE A 114 3.892 -29.444 41.508 1.00 29.36 C
ATOM 881 CE2 PHE A 114 3.175 -27.133 41.524 1.00 28.15 C
ATOM 882 CZ PHE A 114 2.961 -28.488 41.165 1.00 23.55 C
ATOM 883 N CYX A 115 4.434 -28.614 45.678 1.00 20.99 N
ATOM 884 CA CYX A 115 3.916 -29.898 46.202 1.00 20.88 C
ATOM 885 C CYX A 115 2.614 -30.206 45.489 1.00 31.70 C
ATOM 886 O CYX A 115 1.789 -29.310 45.288 1.00 24.52 O
ATOM 887 CB CYX A 115 3.632 -29.806 47.726 1.00 26.63 C
ATOM 888 SG CYX A 115 2.960 -31.364 48.464 1.00 28.36 S
ATOM 889 N SER A 116 2.467 -31.439 45.049 1.00 32.36 N
ATOM 890 CA SER A 116 1.156 -31.884 44.551 1.00 48.05 C
今回、このファイルにはCYSが8, 42, 115, 170にあるのですが、いずれもCYXに変更します。115, 170のCYSもCYXにします（大切なことなので2回言いました）。 そしてこのファイルを1LKE_trun_hid_cyx.pdbとして保存しておきます。
2つめのSS結合の設定ですが、これについては以下のStep 2. tleapで紹介します。
Step 1-5. リガンドのパラメータファイルの用意（存在する場合のみ）
タンパク質は化合物と結合して機能を発揮するタイプのものが非常にたくさんあります。1LKEにもリガンド（Digoxigenin）が結合しています。しかし、20種類の"カノニカル"な（正統な）アミノ酸と違って、一般の化合物の形状は非常に多くの形状を持つので、その都度それについて適切な力場や電荷を用意してあげなければなりません。
この化合物のパラメータファイルの用意自体は非常に長くなってしまうため、作成方法は別の記事を参照してください。ここではすでに完成されたファイルdog.mol2が与えられているものとして次へ進みます。
@<TRIPOS>MOLECULE
DOG
62 66 1 0 0
SMALL
resp
@<TRIPOS>ATOM
1 O23 7.2640 0.6350 0.3660 o 1 DOG -0.603790
2 C23 6.1540 0.2200 0.1420 c 1 DOG 0.913720
3 C22 4.9650 0.9130 -0.3940 ce 1 DOG -0.552430
4 O21 5.7960 -1.0930 0.3690 os 1 DOG -0.543340
5 C21 4.4240 -1.2860 0.0110 c3 1 DOG 0.253200
6 C20 3.9340 0.0590 -0.4810 c2 1 DOG 0.088750
7 C17 2.5540 0.3660 -0.9980 c3 1 DOG -0.002670
8 C16 2.0410 -0.6750 -2.0480 c3 1 DOG -0.263190
9 C13 1.4010 0.4740 0.0890 c3 1 DOG 0.635860
10 C18 1.9120 0.6990 1.5210 c3 1 DOG -0.416180
11 C12 0.4900 1.6560 -0.3400 c3 1 DOG 0.078140
12 O12 1.2290 2.8590 -0.1190 oh 1 DOG -0.703970
13 C14 0.5950 -0.8510 -0.0860 c3 1 DOG -0.051480
14 C15 0.6220 -1.0940 -1.6100 c3 1 DOG 0.108190
15 O14 1.4070 -1.8660 0.5410 oh 1 DOG -0.712600
16 C8 -0.8080 -0.8380 0.5810 c3 1 DOG 0.042480
17 C7 -1.5700 -2.1580 0.3480 c3 1 DOG -0.040690
18 C9 -1.6660 0.3920 0.1730 c3 1 DOG 0.040630
19 C11 -0.8550 1.6840 0.3920 c3 1 DOG -0.475450
20 C10 -3.0760 0.4040 0.8680 c3 1 DOG 0.163710
21 C19 -2.9500 0.6320 2.3920 c3 1 DOG -0.573120
22 C5 -3.8090 -0.9540 0.6040 c3 1 DOG 0.141680
23 C6 -2.9460 -2.1570 1.0260 c3 1 DOG -0.331830
24 C4 -4.2990 -1.1120 -0.8530 c3 1 DOG 0.042980
25 C3 -5.1640 0.0520 -1.3340 c3 1 DOG 0.272340
26 O32 -6.3830 0.0000 -0.5830 oh 1 DOG -0.710140
27 C2 -4.4320 1.3850 -1.1280 c3 1 DOG -0.171740
28 C1 -3.9800 1.5490 0.3280 c3 1 DOG -0.056550
29 H1 4.9820 1.9620 -0.6620 ha 1 DOG 0.218680
30 H2 3.8680 -1.6470 0.8800 h1 1 DOG 0.060400
31 H3 4.3750 -2.0640 -0.7600 h1 1 DOG 0.060400
32 H4 2.6320 1.3530 -1.4600 hc 1 DOG 0.057230
33 H5 2.6940 -1.5500 -2.0730 hc 1 DOG 0.070310
34 H6 2.0490 -0.2540 -3.0590 hc 1 DOG 0.070310
35 H7 1.0840 0.8150 2.2280 hc 1 DOG 0.103440
36 H8 2.5110 -0.1450 1.8640 hc 1 DOG 0.103440
37 H9 2.5160 1.6060 1.5680 hc 1 DOG 0.103440
38 H10 0.2860 1.5600 -1.4200 h1 1 DOG 0.036850
39 H11 0.6880 3.5980 -0.4370 ho 1 DOG 0.444350
40 H12 0.4090 -2.1410 -1.8550 hc 1 DOG -0.017340
41 H13 -0.1430 -0.4960 -2.1160 hc 1 DOG -0.017340
42 H14 1.0310 -2.7280 0.3040 ho 1 DOG 0.443700
43 H15 -0.6050 -0.7670 1.6580 hc 1 DOG 0.058100
44 H16 -0.9860 -3.0040 0.7360 hc 1 DOG 0.031780
45 H17 -1.6900 -2.3330 -0.7290 hc 1 DOG 0.031780
46 H18 -1.8490 0.3110 -0.9080 hc 1 DOG 0.078700
47 H19 -1.4190 2.5530 0.0280 hc 1 DOG 0.151810
48 H20 -0.6740 1.8570 1.4580 hc 1 DOG 0.151810
49 H21 -3.9280 0.5080 2.8710 hc 1 DOG 0.137930
50 H22 -2.2570 -0.0600 2.8790 hc 1 DOG 0.137930
51 H23 -2.6080 1.6490 2.6160 hc 1 DOG 0.137930
52 H24 -4.7110 -0.9490 1.2300 hc 1 DOG 0.041280
53 H25 -2.8100 -2.1520 2.1150 hc 1 DOG 0.095400
54 H26 -3.4760 -3.0890 0.7900 hc 1 DOG 0.095400
55 H27 -4.8810 -2.0380 -0.9330 hc 1 DOG -0.016620
56 H28 -3.4560 -1.2110 -1.5480 hc 1 DOG -0.016620
57 H29 -5.3880 -0.0840 -2.4050 h1 1 DOG 0.002960
58 H30 -6.9330 0.7500 -0.8610 ho 1 DOG 0.433910
59 H31 -5.0920 2.2190 -1.4090 hc 1 DOG 0.028770
60 H32 -3.5790 1.4310 -1.8190 hc 1 DOG 0.028770
61 H33 -3.4880 2.5190 0.4630 hc 1 DOG 0.039300
62 H34 -4.8800 1.5760 0.9550 hc 1 DOG 0.039300
@<TRIPOS>BOND
1 1 2 2
2 2 3 1
3 2 4 1
4 3 6 2
5 3 29 1
6 4 5 1
7 5 6 1
8 5 30 1
9 5 31 1
10 6 7 1
11 7 8 1
12 7 9 1
13 7 32 1
14 8 14 1
15 8 33 1
16 8 34 1
17 9 10 1
18 9 11 1
19 9 13 1
20 10 35 1
21 10 36 1
22 10 37 1
23 11 12 1
24 11 19 1
25 11 38 1
26 12 39 1
27 13 14 1
28 13 15 1
29 13 16 1
30 14 40 1
31 14 41 1
32 15 42 1
33 16 17 1
34 16 18 1
35 16 43 1
36 17 23 1
37 17 44 1
38 17 45 1
39 18 19 1
40 18 20 1
41 18 46 1
42 19 47 1
43 19 48 1
44 20 21 1
45 20 22 1
46 20 28 1
47 21 49 1
48 21 50 1
49 21 51 1
50 22 23 1
51 22 24 1
52 22 52 1
53 23 53 1
54 23 54 1
55 24 25 1
56 24 55 1
57 24 56 1
58 25 26 1
59 25 27 1
60 25 57 1
61 26 58 1
62 27 28 1
63 27 59 1
64 27 60 1
65 28 61 1
66 28 62 1
@<TRIPOS>SUBSTRUCTURE
1 DOG 1 TEMP 0 **** **** 0 ROOT
上では表現できませんが、**最後に空行をはさんでおきます。**これをdog.mol2としてPDBディレクトリに入れて保存します。
参考： antechamberを使ったMD用小分子化合物の電荷パラメータの作り方
Step 2. tleapへのインプットファイル
続いてはtleapを用いたトポロジーファイル作成をするために、topディレクトリを作成します。
# top（トポロジー）のディレクトリ作成と移動
mkdir -p ~/Desktop/1LKE/top
cd ~/Desktop/1LKE/top
ここに以下の内容でleap.inファイルを作成します。Step 1-4で紹介したジスルフィド結合の設定はこのファイルの中で行います（bond mol~ mol~の部分に対応します）。
(2023年2月8日、小分子向けの汎用力場をGAFF2に変更)
#AMBER の力場パラメータff14SBを読み込む
source leaprc.protein.ff14SB #for AmberTools22
source leaprc.gaff2 #小分子向けの汎用力場「GAFF2」を読み込む
source leaprc.water.tip3p #for AmberTools22
#追加のイオンの力場の導入
loadAmberParams frcmod.ionsjc_tip3p
#.mol2ファイルの読み込み
DOG = loadMol2 ../PDB/dog.mol2
#pdbを"mol"として読み込む
mol = loadPDB ../PDB/1LKE_trun_hid_cyx.pdb
#SS-bondを設定
bond mol.8.SG mol.115.SG
bond mol.42.SG mol.170.SG
#イオンの追加。末尾の数だけ対応するイオンを配置する。0を指定すると系全体の電荷を0にするようにイオンを入れるいう設定になる。
addIons2 mol Na+ 20
addIons2 mol Cl- 0
#直方体型の周期境界の溶媒和ボックスを形成する。入力PDBファイルのサイズに対し、最低10A余裕を持って形成する。
solvateBox mol TIP3PBOX 10.0
#最後に、"mol"という溶媒和ボックスの系の電荷情報を表示する。0.00000になっていることが理想。
charge mol
#溶媒和された系のトポロジー・初期座標をleap.parm7, leap.rst7としてそれぞれ保存
saveAmberParm mol leap.parm7 leap.rst7
#溶媒和された系のPDBファイルをleap.pdbとして保存
savePDB mol leap.pdb
quit
このファイルを作成したら、このディレクトリでtleap -f leap.inのコマンドを実行します。しかし、次のようなメッセージが出るはずです。
Total unperturbed charge: -0.000000
Total perturbed charge: -0.000000
Checking Unit.
FATAL: Atom .R<ASN 173>.A<OXT 15> does not have a type.
Failed to generate parameters
Parameter file was not saved.
Writing pdb file: leap.pdb
printing CRYST1 record to PDB file with box info
Converting N-terminal residue name to PDB format: NASP -> ASP
Quit
Parameter file was not saved.と表示されるとトポロジーファイルの作成に失敗しています。まずはここでTotal perturbed charge: -0.000000となっていることを確認することと、そしてFATALが表示されていないことを確認します（上の例では表示されています）。この例では読み込んだPDBファイルのASN 173に謎のOXT15という原子が存在することが原因で失敗しているようです。XXX does not have a type.ってのは、力場に対応する原子タイプが存在しないために発生するエラーメッセージです。
読み込んでいるPDBファイルを編集し、問題が起きている箇所の原子を削っておきましょう。削った箇所には区切れを意味するTERを書いておくと良いです。
ATOM 1335 ND2 ASN 173 30.509 -41.532 33.732 1.00 51.92 N
ATOM 1336 C ASN 173 27.959 -39.594 35.439 1.00 51.92 C
ATOM 1337 O ASN 173 27.577 -39.210 36.576 1.00 51.92 O
TER
HETATM 1339 O23 DOG 174 15.551 -21.450 37.436 1.00105.39 O
HETATM 1340 C23 DOG 174 14.900 -22.454 37.433 1.00105.39 C
HETATM 1341 C22 DOG 174 13.769 -22.807 36.606 1.00105.39 C
こうして再びtleap -f leap.inコマンドを実行します。
Total unperturbed charge: 0.000000
Total perturbed charge: 0.000000
Checking Unit.
Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
total 556 improper torsions applied
Building H-Bond parameters.
Incorporating Non-Bonded adjustments.
Not Marking per-residue atom chain types.
Marking per-residue atom chain types.
(Residues lacking connect0/connect1 -
these don't have chain types marked:
res total affected
CASN 1
DOG 1
NASP 1
WAT 8527
)
(no restraints)
Writing pdb file: leap.pdb
printing CRYST1 record to PDB file with box info
~/apps/amber18/bin/teLeap: Warning!
Converting N-terminal residue name to PDB format: NASP -> ASP
~/apps/amber18/bin/teLeap: Warning!
Converting C-terminal residue name to PDB format: CASN -> ASN
Quit
Exiting LEaP: Errors = 0; Warnings = 2; Notes = 0.
こうしてBuilding topology.と出力され、0バイトではないleap.parm7とleap.rst7ファイルが作成されていればトポロジーファイルの作成ステップは完了です。一応leap.logの中にWARNINGなどがないか確認しておき、あってもそれが大きな問題でなさそうかを確認しておきましょう。
さて、ここで作成されたleap.pdbがどのようになっているか、分子描画ソフトウェアPyMOLなどで開いてみましょう。するとこのようになっているはずです。
図のように、タンパク質の周りに水分子、Naイオン、Clイオンが散りばめられている様子が見えます。特に、タンパク質の周りにたくさん存在する水分子は、ある種の見えない箱みたいなものに閉じ込められているかのように映るはずです。これは、今のleap.inの中で行っていたsolvateBoxによって、周期境界条件の溶媒和ボックスを生成したためです。
なぜこんな水ボックスを生成する必要があるかというと……現実のタンパク質1分子に対して水分子の占める体積は非常に大きいものになっていてほしいはずなのですが、コンピュータの計算の制約上、それをすべて表現するというのは困難です。そこで、無限ループする箱の中にタンパク質と水分子を入れるということで、擬似的に無限を表現することで一応の解決を得ています。図にするとこんな感じ。
この図で水色のボックスを先程生成したタンパク質と水分子の入った基準ボックスとしたとき、$x, y, z$方向に基準ボックスのコピーが無限に広がっているというイメージです。そしてボックス間で原子・分子の移動は可能となっており、ボックスの境界を超えた原子は見かけ上ボックスの反対側からまた流入するというような動きを見せます。また、静電相互作用のような長距離の力も適切に作用します。したがって、あまり小さな周期境界ボックスにタンパク質を入れてしまうと、各ボックスに存在するタンパク質のコピー同士間の静電相互作用が強まってしまいます。そこで上のleap.inでは10 Åの余裕を持って周期境界条件ボックスを作るような設定にしていました（本当はもうちょっと大きな値にしておいた方がいいような気がしますがここではチュートリアルなので勘弁してください）。つまり、隣り合うタンパク質とは20 Å以上常時離れていることになります。この影響を小さくするためにはボックスを大きく取ればいいのですが、そうすると今度は水分子の個数が加速度的に増えていきますのでシミュレーションの計算速度が落ちてしまうことになり、そこはトレードオフとなってしまいます。
この周期境界ボックスのサイズはleap.pdbの1行目にかかれてあります。このファイルをテキストエディタで開いてみますと、
CRYST1 71.492 72.959 68.898 90.00 90.00 90.00 P 1 1
という記述があり、$x=71.492, y=72.959, z=68.898$（Å）のサイズになっていることがわかります。
Step 3. GROMACS形式への変換
Step 2でAMBER形式のトポロジーファイルが完成しました。このままAMBERでのMDシミュレーションを開始することができるのですが、AMBERでのMD計算はGROMACSに比べ計算が2〜3倍遅いので、私はGROMACS形式にトポロジーファイルを変換してからGROMACSでのシミュレーションを使うことを好んでいます（最近はともにGPU対応になり、あまり差がなくなりました）。
GROMACS/ACPYPE.pyを用いたAMBER→Gromacs形式へのtopologyファイル変換のページを参考にして、~/appsディレクトリにacpype.py, groconvert.sh, gen_posre.plの3つのファイルを置き、以下のコマンドで変換します。
# macOSのみ、gsed, gawkコマンドのインストールを行っておく (デフォルトのsed, awkは流派が異なり使えない機能があるため)
brew install gnu-sed gawk
~/apps/groconvert.sh -i leap -o 1lke -r
これで1lke.groと1lke.topファイルを作成します。また、posreXX.itpというファイルがいくつか作成されます。これは後のMDシミュレーションで使う、各重原子（＝水素以外の原子）に対しての位置拘束(position restraint)ファイルです。
1lke.groファイルを開いてみると、中はこのようになっています。
leap_GMX.gro created by acpype (Rev: 0) on Fri Apr 6 11:52:56 2018
28284
1 ASP N 1 2.444 1.727 4.364
1 ASP H1 2 2.350 1.760 4.349
1 ASP H2 3 2.442 1.627 4.377
1 ASP H3 4 2.504 1.768 4.294
（中略）
8732 WAT H228281 1.518 1.755 0.175
8733 WAT O28282 1.421 1.671 1.497
8733 WAT H128283 1.500 1.622 1.477
8733 WAT H228284 1.450 1.739 1.558
7.14920 7.29594 6.88980
GROMACSの座標ファイル形式である.groの記述はシンプルで、1行目はコメント行（何を書いてもよい）、2行目は系の中に存在する原子数（超重要）、3行目以降は原子の記述（残基番号、残基名、atom type、原子の通し番号、x, y, z座標（単位はÅではなくnm））です。2行目の原子数は必ず実際にファイル内に存在する原子の数と一致させておく必要があります。3行目以降の記述は文字の列認識が固定されており、一文字でもずれると正しく認識してくれないことに注意です。ちなみに残基番号と原子の通し番号については実は飾りで、どんな数字を付けていてもエラーにはならず、上から順に認識します。そのため100000を超える原子を含むgroファイルであっても問題なく作動します。最終行は周期境界ボックスのサイズ指定であり、この指定は必須です。先程leap.pdbで見た値から1/10となっていますが、これはGROMACS形式がnm表記で統一していることと、AMBER形式ではÅ表記で統一しているために1/10ずれていることに起因します。なので問題ありません。
Step 4. gromacs indexファイルの作成
後のMDシミュレーションのために、インデックスファイル（通常index.ndx）を作成します。MDシミュレーション内に含まれる原子は、タンパク質に属するものやリガンド、水もイオンも含まれます。また「水＋イオン」とか「（水＋イオン）以外」というグループを作りたいこともあります。このインデックスファイルを作っておくと、「リガンドだけの座標データを切り出して取得」「温度の制御を水イオングループとそれ以外に別々に設定」などということが可能になります。このチュートリアル方式では、MDシミュレーションを始めたら早速このインデックスファイルを要求されることになるので今のうちに作っておきます。
インデックスファイルを作るコマンドはgmx make_ndxです。
$ gmx make_ndx -f 1lke.gro
すると、続いて追加入力を求められます。
There are: 8527 Water residues
Analysing Protein...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
0 System : 28284 atoms
1 Protein : 2605 atoms
2 Protein-H : 1338 atoms
3 C-alpha : 169 atoms
4 Backbone : 507 atoms
5 MainChain : 677 atoms
6 MainChain+Cb : 830 atoms
7 MainChain+H : 842 atoms
8 SideChain : 1763 atoms
9 SideChain-H : 661 atoms
10 Prot-Masses : 2605 atoms
11 non-Protein : 25679 atoms
12 Other : 62 atoms
13 DOG : 62 atoms
14 Na+ : 20 atoms
15 Cl- : 16 atoms
16 Ion : 36 atoms
17 DOG : 62 atoms
18 Na+ : 20 atoms
19 Cl- : 16 atoms
20 Water : 25581 atoms
21 SOL : 25581 atoms
22 non-Water : 2703 atoms
23 Water_and_ions : 25617 atoms
nr : group '!': not 'name' nr name 'splitch' nr Enter: list groups
'a': atom '&': and 'del' nr 'splitres' nr 'l': list residues
't': atom type '|': or 'keep' nr 'splitat' nr 'h': help
'r': residue 'res' nr 'chain' char
"name": group 'case': case sensitive 'q': save and quit
'ri': residue index
>
ここでは>を求められるたびに、以下のように追加入力を行います。qを押すと終了です。
23 Water_and_ions : 25617 atoms
nr : group '!': not 'name' nr name 'splitch' nr Enter: list groups
'a': atom '&': and 'del' nr 'splitres' nr 'l': list residues
't': atom type '|': or 'keep' nr 'splitat' nr 'h': help
'r': residue 'res' nr 'chain' char
"name": group 'case': case sensitive 'q': save and quit
'ri': residue index
> !23
Copied index group 23 'Water_and_ions'
Complemented group: 2667 atoms
24 !Water_and_ions : 2667 atoms
> name 24 nowation
> q
一度目は!23を入力します。これは23: Water_and_ionsに含まれる原子の補集合を表す!をつけた設定です。これを指定すると上記のようなComplemented group 2667 atomsが24 !Water_and_ions : 2667 atomsとして出てきます（番号は系によって異なります）。この24番のグループを改名したい場合はname 24 nowationとします。最後にqを入力して終了です。
長かったですが、これでトポロジーファイルの生成は終了です。操作に慣れてくればここまでで10〜20分くらいで終わります。
この時点でのこの作業ディレクトリ~/Desktop/1LKE/topのファイル構成はこのようになっているはずです。ls -ltコマンドで表示してみますと
$ ls -lt ~/Desktop/1LKE/top
total 25512
-rw-r--r-- 1 Agsmith staff 833741 4 6 11:52 index.ndx
-rw-r--r-- 1 Agsmith staff 42407 4 6 11:52 posre0.itp
-rw-r--r-- 1 Agsmith staff 42407 4 6 11:52 posre10.itp
-rw-r--r-- 1 Agsmith staff 42407 4 6 11:52 posre100.itp
-rw-r--r-- 1 Agsmith staff 42407 4 6 11:52 posre1000.itp
-rw-r--r-- 1 Agsmith staff 42407 4 6 11:52 posre20.itp
-rw-r--r-- 1 Agsmith staff 42407 4 6 11:52 posre200.itp
-rw-r--r-- 1 Agsmith staff 42407 4 6 11:52 posre50.itp
-rw-r--r-- 1 Agsmith staff 42407 4 6 11:52 posre500.itp
-rw-r--r-- 1 Agsmith staff 2056681 4 6 11:52 1lke.top
-rw-r--r-- 1 Agsmith staff 1272891 4 6 11:52 1lke.gro
-rw-r--r-- 1 Agsmith staff 1032459 4 6 11:50 leap.rst7
-rw-r--r-- 1 Agsmith staff 257958 4 6 11:50 leap.log
-rw-r--r-- 1 Agsmith staff 1955061 4 6 11:50 leap.pdb
-rw-r--r-- 1 Agsmith staff 5274680 4 6 11:50 leap.parm7
-rw-r--r-- 1 Agsmith staff 1294 4 6 11:45 leap.in
GROMACSを使ったエネルギー最小化と平衡化
ここからは水溶液中に存在するタンパク質を物理演算で動かしていくことになります。しかし、これまでの操作で作り出した初期座標では水溶液の様子を再現しているとは言えません。その初期座標ファイル1lke.groをPyMOLなどで見てもらえばわかりますが、このファイル内の水分子は規則正しく配置されたままですし、何より原子に速度情報が与えられていません。そこで、まずは系のエネルギーをある程度小さくしてから（＝エネルギー最小化）、各原子に（温度に対するBoltzmann-distributionに従う範囲で）ランダムに速度を与え、徐々に系全体を設定温度まで引き上げていく（＝平衡化）という予備動作を行うことで水溶液の状態に近づけていくという操作を行います。流れは
となり、後のProduction Runというステップにつなげていきます。
エネルギー最小化と平衡化の操作はほとんどすべてのMDシミュレーションにおいて使いまわすことができます。
さて、ここからは計算に時間がとてもかかることになります。しかし近年のGPUを使った計算手法の発展（GPGPU: GPUの演算資源を画像処理以外の目的に応用する技術のこと）のおかげで、GROMACSやAMBER、NAMDなどのMDシミュレーションのプログラムはGPUを使って計算を大幅に加速させることができます。その結果、GPUの取り付けてある計算機とそうでないマシンでは計算速度に雲泥の差が出てくるようになりました。以降の操作はGPUの有無にかかわらず可能ですが、GPUのあるマシンで行うことを推奨します。
上の最小化の操作を ~/Desktop/1LKE/gromacs/minimize で行うものとし、ここにmin1.mdp, min2.mdpとエネルギー最小化実行スクリプトrun.shを置き、NVT, NPT平衡化の操作は ~/Desktop/1LKE/gromacs/heatで行うものとし、ここにmd1.mdp ~ md9.mdpと平衡化実行スクリプトrun.shを置いて実行します。
Production Runの実行
平衡化が終わったら、あとはもう時間とハードディスク容量が許す限りいくらでもMDシミュレーションを行うことができます。この平衡化以降のMDシミュレーションのことをProduction Runと呼びます。実行ディレクトリを~/Desktop/1LKE/gromacs/prに移します。
mkdir -p ~/Desktop/1LKE/gromacs/pr
cd ~/Desktop/1LKE/gromacs/pr
GROMACSの場合、Production Runの設定ファイルmdrunset.mdpは以下のようになります。
define =
; Run parameters
integrator = md ; leap-frog integrator
dt = 0.002 ; 2 fs
nsteps = 50000000 ; 2 fs * 50000000 = 100 ns
cutoff-scheme = Verlet
; Output control
nstxout = 5000 ; save coordinates every 10 ps
nstvout = 5000 ; save velocities every 10 ps
nstenergy = 2500 ; save energies every 5 ps
nstlog = 5000 ; update log file every 10 ps
; Bond parameters
continuation = yes ; Restarting after previous RUN
constraint_algorithm = LINCS ; holonomic constraints
constraints = hbonds ; Convert the bonds with H-atoms to constraints
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Langevin Dynamics
bd_fric = 2.0 ; Brownian dynamics friction coefficient(correspond to gamma_ln in AMBER)
ld_seed = 1732 ; random seed
; Neighborsearching
ns_type = grid ; search neighboring grid cells
nstlist = 20 ; 10 fs
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.12 ; grid spacing for FFT(nm)
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = nowation Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation
gen_temp = 300.0 ; temperature for Maxwell distribution
gen_seed = 1732 ; used to initialize random generator for random velocities
特に重要なのはdt = 0.002とnstepsです。この値は結合長を拘束するSHAKEアルゴリズムやLINCSアルゴリズム（GROMACSではconstraint_algorithmで設定を切り替えられる）を利用している場合に、1ステップあたり2フェムト秒（dt = 0.002）後まで時間をワープさせても計算値が発散せずにシミュレーションを進めることが可能になっています。この繰り返しを規定するのがnstepsの値であり、掛け算した値がシミュレーション時間ということになります。現在、AMBER, NAMD, GROMACS、どのソフトウェアの設定をみてもdt = 0.002としてあることが多いです（一応0.003までやってもいい方法が開発されているのですけど、これをやると計算の解析のときに微妙に切りの悪い数字になりがちなので結局使われることは少ないです）。
GROMACSは初めにnstepsで規定しておいたMDシミュレーション分の時間に達するまで計算し続けます。一応、後からでもgmx convert-tprコマンドでこの時間を変更することが可能なのですが、ふつうは最初の時点でどれくらいのシミュレーション時間をしようか、よく考えてから設定しましょう。今回は100 nsのシミュレーションを行うことにし、nsteps = 50000000（5000万回ステップ計算とdt = 0.002の設定にしておきましょう。
温度制御と圧力制御はそれぞれV-rescale法とParrinello-Rahman法を使っています。これの詳細はものすごく物理学チックな話なので詳細は省きますが、圧力制御のParrinello-Rahman法を使うときは、すでに圧力がある程度落ち着いた状態から始めなければならないという点に気をつけておきましょう。詳細はGROMACS/Production Runの方法 を読んで下さい。
mdrunset.mdpファイルを用意したら、この実行ファイルrun.shを用意してコンピュータ上（NVIDIA GPUがあるマシンを強く推奨）で実行します。
#!/bin/sh
input=1lke
# if md${input}.tpr doesn't exist, grompp will make the tpr file.
if [ ! -e md${input}.tpr ];then
gmx grompp \
-f mdrunset.mdp \
-c ../heat/md9.tpr \
-r ../heat/md9.tpr \
-t ../heat/md9.cpt \
-p ../../top/${input}.top \
-n ../../top/index.ndx \
-o md${input}.tpr \
-po md${input}.mdp
fi
# gmx mdrun
gmx mdrun \
-v -s md${input} -deffnm ${input} -cpi ${input}.cpt
例によってgmx gromppで先にMDシミュレーションの実行バイナリファイルmd1lke.tprを作り出します。今回は-c ../heat/md9.tprと-t ../heat/md9.cptで座標情報と速度情報をheatの最終ステップから読み込む、つまりheatの計算直後からデータを引き継いで計算を実行すると設定しています。またここではif文を使っており、このファイルが同じディレクトリ上にない場合にはgmx gromppを実行する方式です。
これができたら続けてgmx mdrunでProduction Runを実行します。
# 作成したrun.shに実行許可のパーミッションを付与する
chmod +x run.sh
# 上で実行許可してあれば、そのプログラムは ./hoge という記法でターミナルから実行可能となる。
./run.sh
です。
この系を構成する原子数はおよそ28000ほどで、計算時間はNVIDIA GPUの1つ、GTX 1080Tiを積んだマシンでおよそ8〜10時間程度かかります。このような計算時間が長い場合には、計算途中でもしマシンがなんらかの理由でシャットダウンさせられてしまったり、もしくはユーザーの都合で計算を打ち切るというような事がありえます。そんなとき、計算を再スタートさせるときにまたProduction Runを一からやり直すというのはとても大変なことです。そこで、たいていのMDソフトウェアには計算の再スタートを円滑にするためのチェックポイント機能が備わっています。Gromacsのデフォルト設定では15分に一度、$input.cptという名前のチェックポイントファイルを作り出してくれます。このファイルにはその計算ステップ時点での時間情報、座標情報と速度情報を保存してくれているので、これさえ定期的に作り出していればもし計算が異常終了してもここから再スタートすることができるという仕組みです。
以上でProduction Runの説明は終了です。気の済むまでProduction Runをやりましょう。
ちなみにtprファイルを作成してしまった後で、「やっぱりMDのProduction Runの時間を変更したいな……」ってときは、 gmx convert-tpr コマンドでその時間を変更することができます。
# 上で作成したmd1lke.tprファイルの設定時間を300,000 ps (=300 ns)に変更したものをhoge.tprに書き出す
gmx convert-tpr -s md1lke.tpr -o hoge.tpr -until 300000
# hoge.tprをmd1lke.tprにリネーム（上書き）。
mv hoge.tpr md1lke.tpr
こうして作り変えた md1lke.tpr を作って再び上記 run.sh スクリプトを実行すれば300 nsまでのシミュレーションが行われるようになります。
AmberToolsを使ったトラジェクトリファイルの解析
Production Runで得られたトラジェクトリファイル1lke.trrを解析するには再びAmberToolsを使います。特に、AmberToolsのcpptrajというモジュールが大活躍します。
cpptrajというモジュールはトラジェクトリの処理や各種座標情報の解析に用います。詳しくはAmberToolsのマニュアルのcpptrajの章を読んでください。amber22のマニュアルはamber22ディレクトリの下のamber22_src/doc/Amber22.pdfにあります。2023年5月以降はAmberTools23が登場したことでAmber23.pdfという名前になっています。
# /path/to/はamber22_srcがあるところに読みかえる。
# openコマンドはmacOSのみ
$ open /path/to/amber22_src/doc/Amber23.pdf
これでマニュアルが開きます。
上記のcpptrajコマンドはbashやpythonのようなインタプリタ言語式でスクリプトを記述します。基本的にはtrajinコマンドで
## 1フレーム目だけを取り出す準備
# trajinで処理させたいトラジェクトリ名を指定。その後の1 1 1は、[最初のフレーム] [最後フレーム] [サンプリング間隔]いう意味。ここでは最初の1枚めのフレームだけ取り出す。
trajin 1lke.trr 1 1 1
# reference構造の指定
reference ../../top/leap.rst7
# 残基番号で1から170に対応する残基を、同じBoundary box内に位置するようにする（ボックスの外に行ってしまっても途中でループさせない）
unwrap :1-170
# 残基番号1-170のCα原子の重心を原点上にセンタリングするよう、各フレームのトラジェクトリを平行移動させる
center :1-170@CA mass origin
# 1フレーム目（first）をリファレンスとして、rmsが最小になるよう回転操作を行う（上記4段階目の操作に対応）
rms first out rmsd.dat @CA
# stripコマンドで、SOD, WAT……などの残基に対応するトラジェクトリを出力対象から外す（水・イオンのないトラジェクトリを作成するときに有用）。
strip :SOD,CLA,WAT,TIP3,Na+,Cl-
#以上の処理を施したものを1lke_init.pdbという名前で出力する
trajout 1lke_init.pdb pdb nobox
# 1フレーム目の取り出し実行
go
# これまでの入力情報をリセット
clear trajin
## 2段階目のtrajin処理の開始準備
trajin 1lke.trr 1 last 1
# もし複数のトラジェクトリファイルに分割している場合はこれに続けてtrajinで追加する
# trajin 1lke_2.trr 1 last 1
# trajin 1lke_3.trr 1 last 1
unwrap :1-170
trajout 1lkenj.trr
## 2段階目のtrajin処理実行
go
# これまでの入力情報をリセット
clear trajin
## 3段階目のtrajin処理の開始準備。
trajin 1lkenj.trr 1 last 10
unwrap :1-170
center :1-170@CA mass origin
# 1フレーム目（first）をリファレンスとして、rmsが最小になるよう回転操作を行う（上記4段階目の操作に対応）
rms first out rmsd.dat @CA
# stripコマンドで、SOD, WAT……などの残基に対応するトラジェクトリを出力対象から外す（水・イオンのないトラジェクトリを作成するときに有用）。
strip :SOD,CLA,WAT,TIP3,Na+,Cl-
#以上の処理を施したものを1lkenj.trrという名前で出力する
trajout 1lkenj2.trr
## 3段階目のtrajin処理実行
go
このファイルtrajfix.inを~/Desktop/1LKE/gromacs/prディレクトリに用意した後、以下のコマンドで実行します。-pにはtleapのときに作成したparm7ファイルを指定します。
cpptraj -p ../../top/leap.parm7 -i trajfix.in
これで1lkenj.trrと1lkenj2.trrという名前のファイルが出力されます。1lkenj.trrには周期境界条件を修正した状態、1lkenj2.trrには周期境界条件の修正と水・イオン分子の削除を行い、それを10ステップ刻みで出力した状態のトラジェクトリデータが含まれます。また、1lke_init.pdbファイルにはMD開始時点の構造情報が入っています。以降の可視化では1lkenj2.trrの方を用います。
ちなみに1lkenj.trrの処理は必要でない場合もありますが、あらかじめ一度unwrapかつステップ刻みを1ずつに設定して出力しておいた方が、第2段階以降のトラジェクトリ変換に失敗しにくくなるので、あえてこのようにしています。
PyMOLでトラジェクトリを表示する
タンパク質構造可視化ソフトウェアであるPyMOLを使ってMDシミュレーションの結果を可視化します。インストール方法はmacOS/Ubuntu 22.04へのオープンソース版PyMOL 2.5のインストール方法この記事を参考にして下さい。
1lke_init.pdbと1lkenj2.trrというファイルを使って、手持ちのPCのPyMOL上でトラジェクトリファイルを表示させることで、タンパク質の動きを様々な角度から見ることができます。ここではPyMOLを用います。
まずPyMOLで1lke_init.pdbのファイルを開いた後、さらにその構造オブジェクトの上に1lkenj2.trrのトラジェクトリファイルをドラッグ＆ドロップで載せることでロードすることができます。コマンドで行う場合、PyMOLのコマンドウィンドウのところにload 1lkenj2.trr, 1lke_init, 1と入力することで可能です。
cpptrajで距離・角度・二面角と残基ゆらぎの値などを取得する
cpptrajモジュールはトラジェクトリについての情報を解析し、ファイルに出力させることができます。よく出力に用いられるのはトラジェクトリ中の2原子間の距離とか3原子間の角度、4原子間の二面角などです。
以下がそのサンプルコードです。HIE31のNE2原子とDOG170のO14原子の間の距離とDOG170のC15, C14, O14原子とHIE31のNE2原子が形成する二面角を、トラジェクトリのフレームごとに出力します。また、atomicfluctで各残基のゆらぎ(RMSF)を計算させることができます。
trajin 1lkenj2.trr 1 last 1
# HIE31のNE2原子とDOG170のO14原子の間の距離をfoo.txtに出力
distance dist1 :31@NE2 :170@O14 out foo.txt
# MDシミュレーション中の各残基のゆらぎ（RMSF）を計算するatomicfluct
atomicfluct out rmsf.txt @CA byres
# DOG170のC15, C14, O14原子とHIE31のNE2原子が形成する二面角の値をfoo.txtに出力
dihedral dih1 :170@C15 :170@C14 :170@O14 :31@NE2 out foo.txt
strip :SOD,CLA,WAT,TIP3,Na+,Cl-
go
ところで、この1lkenj2.trrはすでに水・イオンを抜いた状態のトラジェクトリになっているため、その状態に対応するトポロジーファイル(.parm7)を用意してあげる必要があります。これはAmberToolsのante-MMBSA.pyで以下のようにして用意することができます。
cd ~/Desktop/1LKE/top/ # leap.parm7があるtopディレクトリに移動
ante-MMPBSA.py -p leap.parm7 -s :WAT,Na+,Cl- -c leap_stripped.parm7
これを行うと
Stripping :WAT,Na+,Cl- (solvent) from original topology, output is leap_stripped.parm7
Done stripping solvent!
と表示され、同ディレクトリにleap_stripped.parm7というファイルが生成されます。これはante-MMPBSA.pyのときに-sで指定した残基名を抜いた状態のparm7ファイルとなっています。
これらを用いて、再度~/Desktop/1LKE/gromacs/prディレクトリで上記analysis.inファイルとともに
cpptraj -p ../../top/leap_stripped.parm7 -i analysis.in
を実行すればOKです。
ちなみに、先述のトラジェクトリを書き出すコードtrajfix.inの中に組み込んでおいてそれにまとめることもできます。この場合はleap_stripped.parm7ではなくleap.parm7を用います。
## 1フレーム目だけを取り出す準備
# trajinで処理させたいトラジェクトリ名を指定。その後の1 1 1は、[最初のフレーム] [最後フレーム] [サンプリング間隔]いう意味。ここでは最初の1枚めのフレームだけ取り出す。
trajin 1lke.trr 1 1 1
# reference構造の指定
reference ../../top/leap.rst7
# 残基番号で1から170に対応する残基を、同じBoundary box内に位置するようにする（ボックスの外に行ってしまっても途中でループさせない）
unwrap :1-170
# 残基番号1-170のCα原子の重心を原点上にセンタリングするよう、各フレームのトラジェクトリを平行移動させる
center :1-170@CA mass origin
# stripコマンドで、SOD, WAT……などの残基に対応するトラジェクトリを出力対象から外す（水・イオンのないトラジェクトリを作成するときに有用）。
strip :SOD,CLA,WAT,TIP3,Na+,Cl-
#以上の処理を施したものを1lke_init.pdbという名前で出力する
trajout 1lke_init.pdb pdb nobox
# 1フレーム目の取り出し実行
go
# これまでの入力情報をリセット
clear trajin
## 2段階目のtrajin処理の開始準備
trajin 1lke.trr 1 last 1
# もし複数のトラジェクトリファイルに分割している場合はこれに続けてtrajinで追加する
# trajin 1lke_2.trr 1 last 1
# trajin 1lke_3.trr 1 last 1
unwrap :1-170
trajout 1lkenj.trr
## 2段階目のtrajin処理実行
go
# これまでの入力情報をリセット
clear trajin
## 3段階目のtrajin処理の開始準備。
## Gromacsは1フレーム目に初期情報が入っているので2フレームめ抽出すると良い。
trajin 1lkenj.trr 2 last 10
unwrap :1-170
center :1-170@CA mass origin
# 1フレーム目（first）をリファレンスとして、rmsが最小になるよう回転操作を行う（上記4段階目の操作に対応）
rms first out rmsd.txt @CA
# stripコマンドで、SOD, WAT……などの残基に対応するトラジェクトリを出力対象から外す（水・イオンのないトラジェクトリを作成するときに有用）。
strip :SOD,CLA,WAT,TIP3,Na+,Cl-
# HIE31のNE2原子とDOG170のO14原子の間の距離をoutput.txtに出力
distance dist1 :31@NE2 :170@O14 out output.txt
# MDシミュレーション中の各残基のゆらぎ（RMSF）を計算するatomicfluct
atomicfluct out rmsf.txt @CA byres
# DOG170のC15, C14, O14原子とHIE31のNE2原子が形成する二面角の値をoutput.txtに出力
dihedral dih1 :170@C15 :170@C14 :170@O14 :31@NE2 out output.txt
#以上の処理を施したものを1lkenj.trrという名前で出力する
trajout 1lkenj2.trr
## 3段階目のtrajin処理実行
go
これを用いて同様にcpptraj -p ../../top/leap.parm7 -i trajfix.inコマンドを実行すれば、距離・二面角情報の書き出しとトラジェクトリの書き出しが同時に行われます。
出力されたoutput.txtを見てみましょう。
#Frame dist1 dih1
1 2.9622 0.4116
2 2.9025 -4.1638
3 2.8128 -30.7028
4 2.8926 -17.0340
5 2.9170 -20.7868
6 3.0112 -4.1119
...
このように2つの値が並べて1000フレーム分（100ナノ秒）出力されています。これとPyMOLで1lkenj2.trrを表示してみたときの結果を表示してみると正しく値を取り出せていることも確認できるでしょう。
また、rmsd.txtの方もこうなっています。
#Frame RMSD_00002
1 0.0000
2 1.2040
3 1.2095
4 1.2948
5 1.3648
6 1.6311
7 1.8407
8 1.8810
9 1.8093
10 1.7317
...
これは後ほどRMSDの時間発展として使うことができます。
RMSFの値はrmsf.txtに出力されます。
#Res AtomicFlx
1.000 3.9846
2.000 3.0131
3.000 2.1024
4.000 1.4635
5.000 2.1263
6.000 1.8951
7.000 1.4176
...
matplotlibで結果を図にする
cpptrajで得られる角度情報をmatplotlibでプロットしてみましょう。この技術をマスターすれば、論文に載せられる図も簡単に作成することができます。
例として以下のrmsd.txtファイルのRMSDの時間変化を図にすることを考えます。
#Frame RMSD_00002
1 0.0000
2 1.2040
3 1.2095
4 1.2948
5 1.3648
6 1.6311
7 1.8407
8 1.8810
9 1.8093
10 1.7317
...
以下のplot.pyでプロットすることができます。当然python3にmatplotlibとnumpyがインストールされていることが条件です。まだインストールしていない方はターミナルからpipでインストールしておきましょう。
$ python3.11 -m pip install numpy matplotlib
#!/usr/bin/env python3.11
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update(
{
"text.usetex": False,
"font.family": "sans-serif",
"font.sans-serif": "Arial",
"font.size": 18,
"lines.color": "b",
"lines.linewidth": 1.0,
}
)
# 入力ファイル名
file = "rmsd.txt"
# スペース区切りのテキストファイルからデータセットを生成
# "#"で始まる行はコメント行として無視する。
# 1列目をX軸、2列目をY軸とみなしている。
x, y = np.loadtxt(file, comments="#", unpack=True)
# 1000フレーム分のデータを100nsに置き換えるため、x軸の値に0.1をかける
x = 0.1 * x
# plot
fig = plt.figure()
fig1 = fig.add_subplot(111)
fig1.plot(x, y, label="rmsd")
# x軸の設定
fig1.set_xlim(0, 100) # 表示範囲
fig1.set_xlabel("Simulation time [ns]") # ラベル
# y軸の設定
fig1.set_ylim(0, 10) # 表示範囲
fig1.set_yticks([0, 2, 4, 6, 8, 10])
fig1.set_ylabel(
r"RMSD [$\AA$]",
)
fig1.legend()
plt.show()
# save
fig.savefig("rmsd.png", bbox_inches="tight", pad_inches=0.05)
これを実行するとrmsd.pngという画像ファイルが出力されます。
ちなみに"text.usetex": False,をTrueに、fig.savefig時に.pdf拡張子を指定するとベクター画像形式になります。後でこのファイルをIllustratorなどで開いて目盛りなど見た目を再調整したいときはこちらが便利でしょう。
ここから見た目の変更などを行いたい場合は、matplotlibのチュートリアルとかAPIリファレンスとかを読み、上で何をしているかを理解して思い通りのフォーマットを作ってあげましょう。
また、RMSFを同様にしてプロットしてみます。$x$軸は残基番号、$y$軸がRMSF[Å]になります。
#!/usr/bin/env python3.11
# %%
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update(
{
"text.usetex": False,
"font.family": "sans-serif",
"font.sans-serif": "Arial",
"font.size": 18,
"lines.color": "b",
"lines.linewidth": 1.0,
}
)
# 入力ファイル名
file = "rmsf.txt"
# スペース区切りのテキストファイルからデータセットを生成
# "#"で始まる行はコメント行として無視する。
x, y = np.loadtxt(file, comments="#", unpack=True)
# plot
fig = plt.figure()
fig1 = fig.add_subplot(111)
fig1.plot(x, y, label="rmsf")
# x軸の設定
fig1.set_xlim(0, 170) # 表示範囲
# np.arrangeで0から170まで20刻みの配列を作成
fig1.set_xticks(np.arange(0, 170, step=20))
fig1.set_xlabel("Residue Number") # ラベル
# y軸の設定
fig1.set_ylim(0, 10) # 表示範囲
fig1.set_yticks([0, 2, 4, 6, 8, 10])
fig1.set_ylabel(
r"RMSF [$\AA$]",
)
fig1.legend()
plt.show()
# save
fig.savefig("rmsf.png", bbox_inches="tight", pad_inches=0.05)
結果はこんな感じ。
N末端とC末端の残基のゆらぎがそれぞれ大きいことは、PyMOLで表示させたトラジェクトリの動きと対応していることでしょう。また、アミノ酸Asp114からLys118までの"DEDKK"の周辺でRMSFの値が大きくなっていることがわかります。これは結晶構造中でミッシング残基になっていた範囲の残基と一致しており、この領域が本質的にゆらぎが大きく、そのために結晶構造中でも残基の座標が決定できなかったことと一致しています。