分子シミュレーション効率化のためのアルゴリズム

Size: px

Start display at page:

Download "分子シミュレーション効率化のためのアルゴリズム"

しゅんすけいせき
5 years ago
Views:

1 平成 6 年 6 月日分子モデリングと分子シミュレーション分子動力学法実習東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム寺田透

2 本日の講義内容ペプチドの分子動力学シミュレーション水溶液環境のモデルペプチドの分子動力学シミュレーション課題タンパク質の分子動力学シミュレーション課題シミュレーションの高速化シミュレーション実行上の注意点

3 ペプチドの生成 (). UCSF Chimera.8. を起動. Tools Structure Editig Build Structure を選択 3. Start Structure peptide sequece を選択し Peptide Sequece に A を 4 個入力し Apply 3

4 ペプチドの生成 () 4. Add Peptide Sequece ウインドウで主鎖の二面角を指定し OK ( ここでは α helix 構造を作るのでデフォルトのままで良い ) 5. メイン画面に生成されたペプチドの構造が現れるので Actios Atoms/Bods show Actios Ribbo hide 6. File Save PDB で ala4.pdb としてデスクトップに保存 4

5 力場パラメータの割り当て. Tools Structure Editig AddH を選択し OK 水素原子を付加. Tools Structure Editig Add Charge を選択し Stadard residues の力場に AMBER ff99sb を指定し OK 3. Tools Amber Write Prmtop を選択し Folder が C: Users iu Desktop になっていることを確認し File ame に ala4 Select force field type に AMBER ff99sb を指定し Save 5

6 分子動力学シミュレーションの手順. 講義のページから amd.exe と tcl85t.dll をダウンロードしデスクトップに保存. 講義のページから ala4.zip をダウンロードしデスクトップに保存 3. ala4.zip をダブルクリックして解凍し生成されたフォルダ (ala4) に先に作成した ala4.prmtop と ala4.ipcrd を移動 4. ru.bat をダブルクリックして実行 6

7 参考 : ソフトウェア NAMD 本講義で使用無料 AMBER CHARMM 力場に対応 Gromacs 無料 AMBER CHARMM GROMOS 力場に対応他に AMBER CHARMM など 7

8 シミュレーション結果の表示 (). Chimera を起動する ( 起動済みの場合は File Close Sessio を選択 ). Tools MD/Esemble Aalysis MD Movie を選択 3. Trajectory format に NAMD (prmtop/dcd) Prmtop に ala4.prmtop を DCD に mi.dcd eq.dcd prod.dcd を指定し OK 8

9 シミュレーション結果の表示 () 4. MD Movie ウインドウにある再生ボタンをクリックし最初はリボンモデルのまま運動の様子を観察せよ 5. Actios Atoms/Bods show Actios Ribbo hide でスティックモデルに変更せよ 6. 原子の色分けの不具合を File Ope で color.com を開いて修正し運動の様子を観察せよ 9

10 初期構造からのずれ (RMSD). Select Atom Specifier を選択し Atom Specifier to Select と入力し OK Cα 原子を選択. MD Movie のメニューの Aalysis Plot RMSD を選択 3. Igore hydroges を false にし Plot 0

11 参考 : 水素結合距離の測定. MD Movie ウインドウの Aalysis Plot Distaces を選択. 原子間距離を測りたい原子のペアの一方を Ctrl キーを押しながら左クリックで選択しもう一方を Ctrl キーと Shift キーを押しながら左クリックで選択 3. MD Plots ウインドウの Plot ボタンをクリック a helix では i 番目のカルボニル酸素と i+4 番目のアミド窒素が水素結合を形成する Brade & Tooze タンパク質の構造入門第版より引用

12 水素結合距離 [Å] RMSD [Å] シミュレーションの結果初期構造最終構造ステップシミュレーションによって α helix 構造が壊れていることに注意ステップ

13 水溶液環境のモデル () 今回のシミュレーションは真空中で行われており水分子による溶媒効果は考慮されていない生体分子のシミュレーションにおいては水溶液環境を適切なモデルを用いて再現する必要がある 3

14 水溶液環境のモデル () 現在以下の方法がよく用いられている水分子を陽に配置球状に配置直方体状に配置周期境界条件溶媒和自由エネルギーを近似的に求める非極性項溶媒接触表面積に比例極性項連続誘電体モデル Poisso-Boltzma 方程式 Geeralized Bor モデル 4

15 球状の配置 r cap r E cap x x y y z z k 0 0 r cap r rcap r rcap 0 r 0 水分子の蒸発を防ぐため分子が半径 r cap の球の外側に出て行こうとすると系の中心に向けて束縛力をかける系の表面に位置する水分子は中心付近の水とは異なる環境に置かれる 5

16 周期境界条件中央のセルと同じものが無限に繰り返すセルから出て行った分子はそのセルの反対側から入るどの分子も同じ環境系が隣接セルからの影響を感じないように系のサイズを十分に大きくする必要がある 6

17 圧力の計算 F df E TS, de TdS F P V NkBT V NkBT V T 3V 3V de PdV SdT k B l Z T V N i N r i f N ij i ji PdV i r f TdS ij T SdT kbt Z Z V ビリアルの定理 T 周期境界条件ではこちらを使う相互作用のない系 ( 理想気体 ) では PV = Nk B T = RT 7

18 圧力の制御周期境界条件におけるセルの大きさを変化させることで圧力を制御する r i 圧力減圧力増 αr i 分子の重心位置も同様にスケールされる分子内の原子の相対位置は変化しない瞬間的に P < 0 となることがある P NkBT V 3V N N r ij i ji f ij r r ij f i r ij f i f j r ij f f 0 i j f j r f f 0 ij i j 8

19 水溶液中のシミュレーション (). Chimera を起動し ala4.pdb を開く. Stick 表示に変更し水素原子を付加する 3. Tools Structure Editig Solvate を選択し Solvate method に Box Solvet Model に TIP3PBOX Box size に 6 を入力し OK 4. Tools Structure Editig Add Charge で Stadard residues に AMBER ff99sb を指定し OK 9

20 水溶液中のシミュレーション () 5. Tools Amber Write Prmtop を選択し Folderが C: Users iu Desktop になっていることを確認し File ameに ala4-wat Select force field typeに AMBER ff99sb を指定し Save 6. File Save PDB を選択しデスクトップに ala4-wat.pdbとして保存 7. 講義のページから ala4-wat.zipをダウンロードしデスクトップに保存し解凍 8. ala4-wat.prmtop,ala4-wat.ipcrd ala4- wat.pdbを生成したala4-watフォルダに移動 0

21 水溶液中のシミュレーション (3) 9. ala4-wat フォルダを開き restrait.pl をダブルクリック ala4-wat_rest.pdb が生成 0.mi.i を以下の通り修正セルのサイズの情報を転記 ala4-wat.ipcrd をワードパッドで開き末尾を表示セルサイズの整数値に近く 3 5の積で表せる数

22 水溶液中のシミュレーション (4). ru.bat をダブルクリックしシミュレーションを以下の順に実行 ( 約 4 分 ) エネルギー最小化 ( 水分子のみ )(mi) エネルギー最小化 ( 全体 )(mi) 3 平衡化 (0 300 K)(eq, 0 ps) 4 平衡化 ( 定圧 )(eq, 0 ps) 5 プロダクション (prod, 0 ps)

23 Temperature [K] Pressure [bar] Volume [Å 3 ] 結果の解析. コマンドプロンプトを起動し以下を実行 > cd Desktop ala4-wat > eergy.pl eq.log eq.log prod.log. eergy.csv が生成されるので Excel で開く Time [ps] Time [ps] Time [ps] 3

24 平衡化における体積の変化水を配置する際少数の水分子を小さな系で平衡化したモデルタンパク質の周囲にあてはめているがタンパク質の原子と衝突する水分子は機械的に取り除いているため配置した水分子とタンパク質の間に隙間ができる定温定圧シミュレーションを行い水分子の配置を最適化すると隙間が埋まり体積が減少する水のモデルタンパク質の周りに隙間ができる定温定圧シミュレーション体積が減少 4

25 課題平衡化 (eq eq) とプロダクション (prod) における温度 (TEMP) と圧力 (PRESSURE) 体積 (VOUME) の時間変化をプロットせよ時間刻み Δt は fs エネルギー最小化 (mi mi) 平衡化 (eq eq) プロダクション (prod) における Cα RMSD の変化をプロットせよこれらのプロットから何が言えるか考察せよ 5

26 計算時間 () 対象 : 球状に配置した水分子 (TIP3Pモデル) Amber のSaderモジュール使用計算にはItel Xeo Processor 8コア使用時間刻みDtは0.5 fs psの計算にかかる時間 ( 単位は秒 ) を計測原子数 T total [s] 比率 T b [s] T b /T total * * s あたり 4.9 日かかる 6

27 分子シミュレーションの効率化時間刻み Dt を長くする SHAKE 法多重時間積分法非共有結合相互作用の計算の近似カットオフ法多重極子展開法 Particle mesh Ewald (PME) 法本講義では赤枠の 3 つの方法について解説 7

28 SHAKE 法時間刻みは最も速い運動の周期の 0 分のから 0 分の最も速い運動は X H 伸縮運動周期は約 0 fs Δt = 0.5~ fs 次に速い運動は X X 伸縮運動周期は約 0 fs SHAKE 法により X H 結合長を固定長い時間刻み (Δt = fs) の使用が可能 8

29 SHAKE の適用例 Methaol の分子動力学シミュレーション ( 温度制御なし ) における全エネルギーの誤差 ( 初期値との差 ) の推移全エネルギーの誤差 [kcal mol - ] SHAKEなし (Dt = fs) SHAKEあり (Dt = fs) 時間 [ps] 9

30 SHAKE の適用例平均絶対誤差 [kcal mol - ] SHAKE なし SHAKE あり時間刻み Δ t [fs] SHAKE を用いると時間刻み fs でも SHAKE なしの 0.5 fs に匹敵する精度が得られる 30

31 NAMD における設定 () SHAKEを使う場合は以下の設定を行う rigidbods all usegrouppressure yes 3

32 非共有結合相互作用の扱い非共有結合相互作用は原子のペアについて計算する必要がある N 原子系では N(N )/ のペア非共有結合相互作用は距離が離れるほど弱くなる (va der Waals 引力は r 6 に比例静電相互作用は r に比例 ) 離れている原子同士は相互作用しないとみなすカットオフ法 3

33 カットオフ法 r c 原子 i から半径 r c の範囲内にある原子との非共有結合相互作用の計算を行うこの範囲にある原子の平均個数を M とすると非共有結合相互作用の計算量は N(N )/ から NM に減少する 33

34 ペアリストの作成 r l r c カットオフ半径 r c 以内にある原子ペアのリストを作成する必要があるこの計算量は N(N )/ カットオフ半径 r c より外側の半径 r l の範囲でリストを作っておき原子の最大移動度が r l r c を超えた時にリストを更新するようにすると計算量を削減できる 34

35 周期境界条件の場合 () 周期境界条件では基本セルのコピーが無限に続くので全ての原子ペアについて相互作用を近似せずに直接計算することは不可能

36 カットオフ法の適用カットオフ半径によっては基本セルの周辺のイメージセルも考慮する必要がある ( 左の例では 6N + N(N )/ ペアの計算が必要 ) y x 36

37 Miimum image covetio カットオフ半径 r c を最も短い基本セルの辺の長さの分の以下にすれば考慮すべきペア数は N(N )/ でよい miimum image covetio y x 37

38 カットオフ法の問題点 Va der Waals 相互作用は遠距離では r 6 の項が支配的 va der Waals 相互作用はカットオフ法で十分な精度で計算可能静電相互作用は r に依存カットオフ法では精度良く評価することが困難原子がカットオフ半径の範囲から出入りする際にエネルギーが変動するため全エネルギーは保存しない r 6 r Dt 後相互作用なし相互作用あり 38

39 カットオフしない計算法 () r 中央の基本セル内の原子同士だけでなく基本セル内の原子と周囲のイメージセル内の原子との間の相互作用も計算する位置 rにおける静電ポテンシャル : r j r r q j j 39

40 カットオフしない計算法 () 電荷は周期的に分布している電荷分布をフーリエ級数で表すことができる 40 x x i x i x dx x i x x i x exp ~ exp ~ exp ~ exp ~ 0

41 カットオフしない計算法 (3) 3 次元の場合 4 z y x z y x z z y y x x d i i z i y i x i r r r r r, exp ~ , exp ~ exp exp exp ~ cell

42 カットオフしない計算法 (4) Poisso 方程式を解く 4 r r r r r r r r r r cell cell exp ~ ~ exp exp 4 ~ 4 i d i d i ポテンシャルエネルギー関数発散を防ぐために全電荷 ( ) は 0 でなければならない ~ 0

43 Particle Mesh Ewald 法 () x y z の範囲はマイナス無限大から無限大 r を離散化することでこの範囲を限定できる exp ~ exp ~ exp ~,, x x y y z z x x y y z z K K K K k K K K z z z y y y x x x i i i K k K k K k k K k K K k K k k k 高速フーリエ変換で計算できる

44 Particle Mesh Ewald 法 () 電荷分布は点電荷からなる一般には格子点に乗らないガウス関数でぼかす 44

45 Particle Mesh Ewald 法 (3) 本来の電荷分布からの差を求める σ = 点電荷ガウス関数でぼかした電荷分布差の電荷分布静電ポテンシャルをガウス関数でぼかした電荷分布がつくる静電ポテンシャルと差の電荷分布がつくる静電ポテンシャルの和で表す 45

46 Particle Mesh Ewald 法 (4) r σ 小 σ 大差の電荷分布では点電荷のまわりにこれを打ち消す反対の符号の電荷が分布静電ポテンシャルはr より速く0に減衰カットオフ法でも精度よく計算できる 46

47 Elapsed time [s] 計算時間 () Computatioal time for ps 000 近似なし 00 PME 0 PME+SHAKE Number of atoms 水分子の系で計算時間を計測近似なしでは原子数 N の乗に比例 PME を使用することでほぼ NlogN に比例 SHAKE を併用することで時間刻みを 4 倍 ( fs) にでき計算速度は 3. 倍程度高速化した 47

48 NAMD における設定 () PME 法を使う場合は以下の設定を行う cutoff 0.0 switchig off cellbasisvector cellbasisvector cellbasisvector PME yes PMEGridSizeX 45 PMEGridSizeY 40 PMEGridSizeZ 45 extededsystem XSC_file_ame とはいずれかを記載 48

シミュレーション実行上の注意点 () 立体構造の取得 PDB(http://www.

49 シミュレーション実行上の注意点 () 立体構造の取得 PDB( からダウンロード通常生物学的に機能しうる単位である biological uit 構造に対してシミュレーションを行う例 : Ribouclease T (PDB ID: I0X) 非対称単位の構造 Biological Uit の構造 49

50 シミュレーション実行上の注意点 () 欠失残基はモデリングなどで補う N 末端 C 末端が欠失している場合は欠失残基の前後の残基を acetyl 基 N-methyl 基でブロックしても良い水素原子付加 SS 結合の有無 His のプロトン化状態に注意 50

51 His のプロトン化状態 O H N CH C O H N CH C O H N CH C CH CH CH HN N HN N NH NH d 位にプロトン化 e 位にプロトン化 d, e 位にプロトン化 His 側鎖の pk a は中性付近であるためつの窒素原子とも水素原子が結合した状態も十分にとりうる His 周りの水素結合ネットワークからプロトン化状態がわかる 5

52 シミュレーション実行上の注意点 (3) リガンドの力場パラメータは分子動力学ソフトウェアに含まれていないので自分で作成するか Amber Parameter Database * 等から取得する PME を利用する場合は電荷を中性にするためにカウンターイオンを配置平衡化は十分に時間をかけて行う少なくとも s 程度初期構造からあまりずれないように束縛し平衡化の過程で束縛力を徐々に弱めるのが良い 5 *

53 タンパク質の MD シミュレーション (). Chimera で PDB ID CRN の構造を開く. Stick 表示に変更する 3. 水素原子を付加する 4. 水分子を直方体状に配置する 5. 電荷を付加する ( 標準残基の力場パラメータに AMBER ff99sb を指定 ) 6. パラメータファイルを保存する ( ファイル名は CRN 力場パラメータは AMBER ff99sb) 53

54 タンパク質の MD シミュレーション () 7. PDB ファイルを保存 ( ファイル名 :CRN.pdb) 8. 講義のページから CRN.zip をダウンロードしデスクトップに解凍 9. 生成されたフォルダを開き先程保存した CRN.prmtop CRN.ipcrd CRN.pdb を移動 0.restrait.pl を実行 CRN_rest.pdb が生成.mi.i のセルのサイズを修正.ru.bat をダブルクリックし実行 ( 約分 ) 54

55 課題初期構造からの Cα 原子のずれ (RMSD) の時間変化をプロットせよ平衡化 (eq eq) とプロダクション (prod) における温度 (TEMP) とポテンシャルエネルギー (POTENTIA) の時間変化をプロットせよ時間刻み Δt は fs これらのプロットから何が言えるか考察せよ 55

56 課題の提出課題課題の結果と考察をつの PowerPoit ファイルにまとめて寺田宛に送ること ( エクセルファイルはサイズが大きいので送らないこと ) その際件名は分子モデリング課題とし本文に氏名と学生証番号を明記すること 56

分子シミュレーション効率化のためのアルゴリズム

分子シミュレーション効率化のためのアルゴリズム平成 5 年 5 月 7 日分子モデリングと分子シミュレーション分子動力学法実習東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究プログラム寺田透本日の講義内容ペプチドの分子動力学シミュレーション水溶液環境のモデルペプチドの分子動力学シミュレーション課題タンパク質の分子動力学シミュレーション課題シミュレーションの高速化シミュレーション実行上の注意点ペプチドの生成

分子シミュレーション効率化の ためのアルゴリズム

分子シミュレーション効率化のためのアルゴリズム