(Note: These tutorials are meant to provide
illustrative examples of how to use the AMBER software suite to carry out
simulations that can be run on a simple workstation in a reasonable period of
time. They do not necessarily provide the optimal choice of parameters or
methods for the particular application area.)
Copyright McGee, Miller, and Swails 2009
AMBER ADVANCED TUTORIALS
TUTORIAL 3 - SECTION 4
MMPBSA.py
Dwight McGee, Bill Miller III, and Jason Swails
1) 出発構造の構築と平衡化したシステム作成シミュレーションの実行
このシミュレーションでモデル化するシステムは、ヒトH-Rasタンパク質と、シグナル伝達カスケードの中心であるC-Raf1（RAS-RAF）のラス結合ドメインとの複合体である。RAS-RAF複合体の部分的に平衡化し、事前に準備されたPDBは下記に置いてある。
この構造には、rasおよびRAFタンパク質と、生理的に必要なGTPヌクレオチドが含まれている。 :
出発構造を構築し、平衡化したシステムを得るシミュレーションを実行する方法は、Section 1 and Section 2を参照する。
MMPBSA.pyを使用して、結合自由エネルギーを計算するための重要なファイルは、トポロジーファイルとmdcrdファイル(ras-raf_top_mdcrd.tgz)である。
2) 並列計算でRAS-Rafの結合自由エネルギーを計算する
複合体、受容体とリガンドの相互作用エネルギーと溶媒和自由エネルギーを計算し、結合自由エネルギーの推定値を得るために、結果を平均する。チュートリアルのこの部分では結合に寄与するエントロピーの計算を実行しないので、厳密には、その結果は真の自由エネルギーを示さないが、同様のシステム間での比較に使用することはできる。Normal Mode Analysis (Nmode)の例としては、Section 5 にはAMBERのptrajモジュールを使用しQuasi-Harmonicエントロピー計算を実行するために、システムのエントロピーの寄与を計算したり、下記入力ファイルの&general namelistの最後の行をコメントアウトしている。
MM-GBSA法とMM-PBSA法を比較するため両方を使用して結合エネルギー計算しSection 3.1で得た結果と比較する。MMPBSA.py.MPIは、各スレッド（プロセス）に同数のフレームを割り当てることによって計算を並列化する。処理されたフレームの数が、スレッドの数の倍数で開始されると、最も効率的に動作する。しかし、これは必要条件ではなく、フレーム数が、スレッドの数の倍数でない場合、「残りの」フレームが均等に最初のスレッドのサブセット間で分配される。例えば、3のスレッドで50フレームを実行すると、2つのスレッドが17フレームで最後のスレッドが16で計算される。従って、第三のスレッドは、それらが進行する前に、最初の2つの計算終了を待つ必要がある。（取るように）、例えば、5スレッドは、賢明な選択であり、各スレッドが10フレームを使う。しかし、スレッドの数は、処理されるフレーム数を超えることは出来ず、MMPBSA.py.MPIはエラーメッセージを表示して終了する。MMPBSA.py.MPIの入力ファイルはMMPBSA.pyの入力ファイルと同じである。:
| mmpbsa.in |
Input file for running PB and GB &general endframe=50, verbose=1, # entropy=1, / &gb igb=2, saltcon=0.100 / &pb istrng=0.100, / |
MMPBSA.py用の入力ファイルは、AMBERのsanderモジュールのmdinファイルの設定と同様になるように設計されている。各名前リストは、ampersand (&)で開始し次に名前リストが続く。また、名前リストの終了はbackslash (/)または'&end'で表す。すべての変数の完全なリストは、ユーザーズマニュアルhereを参照する。入力ファイルは、general, pb, と gbの3つの名前リストに分割されている。一般的な名前リストは、計算の特定の部分に変数を指定するのではなく、すべての部分に変数を指定するように設計されている。このセットアップでは、受容体をRAS、リガンドはRAFと定義した。'endframe'変数はmdcrdの何フレームで終了するかをセットする。MPIのコマンドの詳細については、マニュアルまたはSection 3.4を参照する。'&gb'と '&pb' ネームリストマーカーは、スクリプトでそれらの名前リスト内で定義され指定された値でMM-GBSAとMM-PBSA計算を実行するかを決めている。
mpirun -np 4 $AMBERHOME/bin/MMPBSA.py.MPI -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp ras-raf_solvated.prmtop -cp ras-raf.prmtop -rp ras.prmtop -lp raf.prmtop -y *.mdcrd > progress.log 2>&1
または、PBSのようなキューイングシステムに下記のスクリプトを入力すると
qsub parallel.job
スクリプトparallel.jobは、bashシェルの場合以下のようになる。:
| parallel.job |
#!/bin/sh
#PBS -N rasraf_parallel
#PBS -o parallel.out
#PBS -e parallel.err
#PBS -m abe
#PBS -M email@address.com
#PBS -q brute
#PBS -l nodes=1:node:ppn=4
#PBS -l pmem=900mb
cd $PBS_O_WORKDIR
mpirun -np 4 MMPBSA.py.MPI -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp ras-raf_solvated.prmtop \
-cp ras-raf.prmtop -rp ras.prmtop -lp raf.prmtop -y bigprod.mdcrd > progress.log 2>&1
|
このスクリプトの途中経過結は、計算の進行ファイルprogress.logに印刷される。計算中のすべてのエラーがこのファイルに出力される。最後に、計算が終わると、計算は、スクリプトの各ステップの間に要した時間を示す。
| progress.log |
MMPBSA.py.MPI being run on 4 processors ptraj found! Using /scr/arwen_3/swails/i686/amber11/exe/ptraj sander found! Using /scr/arwen_3/swails/i686/amber11/exe/sander Preparing trajectories with ptraj... 50 frames were read in and processed by ptraj for use in calculation. Starting calculations Starting gb calculation... Starting pb calculation... Calculations complete. Writing output file(s)... Timing: Processing Trajectories With Ptraj: 0.126 min. Total GB Calculation Time (sander): 4.782 min. Total PB Calculation Time (sander): 28.407 min. Output File Writing Time: 0.053 min. Total Time Taken: 33.379 min. MMPBSA Finished. Thank you for using. Please send any bugs/suggestions/comments to mmpbsa.amber@gmail.com |
スクリプトのコマンドラインで '-y *.mdcrdは' '.mdcrd'で終わる作業ディレクトリ内のすべてのファイルを読み軌跡の解析にそれらを使用することを指示する。
「keep_files」はデフォルト値として１がセットされ、出力ファイルは: Parallel_output.tgz.
スクリプトはptrajを使用して、GBとPBの計算中に分析された座標である3つの非溶媒和mdcrdファイル（複合体、受容体、およびリガンド）を作成する。。*.mdoutファイルは指定されたすべてのフレームのエネルギーを含んでいる。quasi-harmonicエントロピーの計算がptrajで行われる場合、平均PDBファイルは、平均構造として作成される。MMPBSA.py.MPIによって作成されたすべてのファイルは、最終的な出力ファイルを除き、接頭辞「_MMPBSA_ 'で始まる。
| FINAL_RESULTS_MMPBSA.dat |
| Run on Sun Feb 14 19:10:43 EST 2010 |Input file: |-------------------------------------------------------------- |Input file for running PB and GB in serial |&general | endframe=50, verbose=1, | mpi_cmd='mpirun -np 3', nproc=3 |/ |&gb | igb=2, saltcon=0.100 |/ |&pb | istrng=0.100, |/ |-------------------------------------------------------------- |Solvated complex topology file: ras-raf_solvated.prmtop |Complex topology file: ras-raf.prmtop |Receptor topology file: ras.prmtop |Ligand topology file: raf.prmtop |Initial mdcrd(s): bigprod.mdcrd | |Best guess for receptor mask: ":1-166" |Best guess for ligand mask: ":167-242" |Calculations performed using 50 frames. |Poisson Boltzmann calculations performed using internal PBSA solver in sander. | |All units are reported in kcal/mole. ------------------------------------------------------------------------------- ------------------------------------------------------------------------------- GENERALIZED BORN: Complex: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -1863.7944 17.1704 2.4283 EEL -17200.7297 75.9366 10.7391 EGB -3249.6511 65.2075 9.2217 ESURF 91.3565 1.3938 0.1971 G gas -19064.5240 77.8536 11.0102 G solv -3158.2946 65.2224 9.2238 TOTAL -22222.8186 51.0216 7.2155 Receptor: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -1268.1888 14.2342 2.0130 EEL -11557.0773 71.7127 10.1417 EGB -2532.0669 57.7003 8.1600 ESURF 64.2843 1.1143 0.1576 G gas -12825.2661 73.1118 10.3396 G solv -2467.7826 57.7110 8.1616 TOTAL -15293.0487 35.3527 4.9996 Ligand: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -529.3090 9.4198 1.3322 EEL -4684.4720 36.1449 5.1117 EGB -1688.9631 26.5353 3.7527 ESURF 37.0493 0.6185 0.0875 G gas -5213.7811 37.3522 5.2824 G solv -1651.9138 26.5425 3.7537 TOTAL -6865.6949 25.8878 3.6611 Differences (Complex - Receptor - Ligand): Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -66.2966 4.2751 0.6046 EEL -959.1803 34.9190 4.9383 EGB 971.3789 33.0497 4.6739 ESURF -9.9770 0.3759 0.0532 DELTA G gas -1025.4769 35.1797 4.9752 DELTA G solv 961.4018 33.0518 4.6742 DELTA G binding = -64.0750 +/- 6.3729 0.9013 ------------------------------------------------------------------------------- ------------------------------------------------------------------------------- POISSON BOLTZMANN: Complex: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -1863.7944 17.1704 2.4283 EEL -17200.7297 75.9366 10.7391 EPB -3207.7160 66.4023 9.3907 ECAVITY 67.8762 0.7818 0.1106 G gas -19064.5240 6061.1875 857.1813 G solv -3139.8399 66.4069 9.3914 TOTAL -7686.8660 52.5400 7.4303 Receptor: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -1268.1888 14.2342 2.0130 EEL -11557.0773 71.7127 10.1417 EPB -2483.7242 56.4551 7.9840 ECAVITY 47.1495 0.4737 0.0670 G gas -12825.2661 5345.3320 755.9441 G solv -2436.5747 56.4571 7.9842 TOTAL -5250.2060 38.5188 5.4474 Ligand: Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -529.3090 9.4198 1.3322 EEL -4684.4720 36.1449 5.1117 EPB -1670.4169 27.6694 3.9131 ECAVITY 28.0328 0.4133 0.0584 G gas -5213.7811 1395.1865 197.3092 G solv -1642.3841 27.6725 3.9135 TOTAL -2350.3020 25.1197 3.5525 Differences (Complex - Receptor - Ligand): Energy Component Average Std. Dev. Std. Err. of Mean ------------------------------------------------------------------------------- VDWAALS -66.2966 4.2751 0.6046 EEL -959.1803 34.9190 4.9383 EPB 946.4251 34.5128 4.8808 ECAVITY -7.3062 0.3004 0.0425 DELTA G gas -1025.4769 1237.6138 175.0250 DELTA G solv 939.1189 34.5141 4.8810 DELTA G binding = -86.3579 +/- 8.3264 1.1775 ------------------------------------------------------------------------------- ------------------------------------------------------------------------------- WARNINGS: igb=2 should be used with mbondi2 pbradii set. Yours are modified Bondi radii (mbondi) |
統計ファイルの先頭には、日付/時刻、入力ファイルのコピー、リガンドと受容体のマスクのための最適の推測、スクリプトで使われたファイル、分析されたフレームの数、使用されたPBソルバ（もしあれば）を含んでいる。統計ファイルの残りの部分は、すべての平均エネルギー、標準偏差、およびPBに続いたGBの平均の標準誤差を含んでる。各セクションの後、結合のΔGは、エラー値と一緒に与えられる。ファイル内の別の用語の意味は次の通りである:
VDWAALS = van der Waals contribution from MM.
EEL = electrostatic energy as calculated by the MM force field.
EPB/EGB = the electrostatic contribution to the solvation free energy calculated by PB or GB respectively.
ECAVITY = nonpolar contribution to the solvation free energy calculated by an empirical model.
DELTA G binding = final estimated binding free energy calculated from the terms above. (kCal/mol)
受容体とリガンドの接合ポテンシャル項の値が、単一のtrajectoryアプローチを使用した複合体接合ポテンシャルを正確にキャンセルする必要があるので、全ガス相エネルギーが報告されていないことに注意する。エネルギーが許容公差内でキャンセルしない場合は、エラーメッセージが発生する。
一般的に複合体形成に有利な静電エネルギーと不利な溶媒和自由エネルギーを予測する。これは、結合粒子の脱溶媒和することと, それらを結合インタフェースに整列させるために使用しなければならないエネルギーを象徴している。
負の全結合自由エネルギー-86.36キロカロリー/モルから、純水中ではこのタンパク質-タンパク質複合体は形成されるが、複合体形成に不利なエントロピーの寄与を推定しなかったので、結果は実際の結合自由エネルギーに等しくないことを心に留めておく。GBのアプローチは、わずかに低い結合エネルギーを与えるが、それでもこれは、複合体を形成する傾向にあることを示唆していることに注意する。
CLICK HERE TO GO TO SECTION 3.5
CLICK HERE TO GO TO MMPBSA.py HOME
CLICK HERE TO GO BACK TO THE INTRODUCTION
(Note: These tutorials are meant to provide
illustrative examples of how to use the AMBER software suite to carry out
simulations that can be run on a simple workstation in a reasonable period of
time. They do not necessarily provide the optimal choice of parameters or
methods for the particular application area.)
Copyright McGee, Miller, and Swails 2009