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

Size: px

Start display at page:

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

ゆきさやぶき
4 years ago
Views:

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

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

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

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

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

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

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

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

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

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

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

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

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

14 球状の配置 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 の球の外側に出て行こうとすると系の中心に向けて束縛力をかける系の表面に位置する水分子は中心付近の水とは異なる環境に置かれる 4

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

16 圧力の計算 F df = E TS, = de TdS de = PdV SdT F P = = kt V T N NkT = + ri V 3V i= N NkT = + V 3V f = PdV ln Z V N ij i= j= i+ i r f + TdS ij T = SdT kt Z Z V ビリアルの定理 T 周期境界条件ではこちらを使う相互作用のない系 ( 理想気体 ) では PV = NkT = nrt 6

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

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

19 水溶液中のシミュレーション () 5. Tools Amber Write Prmtop を選択し Folderが C: Users iu Desktop になっていることを確認し File nameに 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.inpcrd ala4- wat.pdbを生成したala4-watフォルダに移動 9

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

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

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

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

24 課題平衡化 (eq eq) とプロダクション (prod) における温度 (TEMP) と圧力 (PRESSURE) 体積 (VOLUME) の時間変化をプロットせよ時間刻み Δt は fs 平衡化 (eq eq) とプロダクション (prod) における水素結合長の変化をプロットせよこれらのプロットから何が言えるか考察せよ 4

25 溶媒和自由エネルギーの近似 () 以下のような熱力学過程を考える q q q 3 q5 q q q 3 q5 q 4 q 4 q q q 3 q5 q 4 電荷 0 の溶質を溶媒に溶かす G np 電荷を移動する G pol 溶媒和自由エネルギー G solv = G np + G pol 5

26 溶媒和自由エネルギーの近似 () 非極性項 ( G np ) は炭化水素の溶媒和自由エネルギーの実験データから溶媒接触表面積 (solvent-accessible surface area, SASA) に比例すると近似できる G = σ A + b 極性項は溶媒を連続誘電体とみなして電磁気学の理論を用いて求める G np pol = A: SASA σ, b: 経験的パラメータ [ φ( r) φ ( r) ] ρ( r) dr vac 静電ポテンシャル溶質の電荷分布 6

27 誘電体コンデンサーに比誘電率 ε の誘電体を挿入すると誘電体の表面に電荷が現れ極板間の電場を打ち消す静電ポテンシャルは /ε となる + +δq δq +δq δq +δq δq 水溶液中では水分子が配向して誘電体として働き静電相互作用を弱める 7

28 連続誘電体モデル分子表面にプローブ球 ( 水の場合半径.4 Å) を転がした時球の中心が作る軌跡溶媒接触表面 (solvent-accessible surface, SAS) SAS からプローブ球の半径分内側の点がつくる表面分子表面 (molecular surface, MS) MS の内側を低誘電率 (ε = ~4) 外側を溶媒の誘電率 ( 水の場合 ε = 80) の誘電体とみなす Leach, A.R. Molecular modeling principles and applications nd Ed. Pearson, England 8

29 Poisson-Boltzmann 方程式連続誘電体モデルにおいて静電ポテンシャルを与える塩がない場合 Poisson 方程式 ε r φ r = 4π ρ r [ ( ) ( )] ( ) 静電ポテンシャル溶質の電荷分布塩が存在する場合塩の電荷分布は Boltzmann 分布に従う [ ε ( r) φ( r) ] = 4π [ ρ( r) + ρ ( r) ] ion 塩の電荷分布 9

30 Generalized Born モデル Poisson-Boltzmann 方程式の問題点力の計算ができない計算コストが高い Generalized Born モデルの特徴イオンの溶媒和自由エネルギーの式を拡張 G pol = q a ε G GB ij 計算コストが低い力を解析的に求めることが可能 f pol = = r ε + a a i j exp N q q f i i, j= GB ( r 4a a ) ij j i j 30

31 非極性項のモデル横軸に溶媒接触表面積縦軸にモル溶解度の対数をプロット * モル溶解度 sと自由エネルギー µ µ = RT ln s µ w w o : 水溶液中での標準化学ポテンシャル µ o 現在では比例定数 σ に 5 cal mol - Å - が使われる ** この項に極性項を合わせて GB/SA(PB/SA) モデルと呼ばれる : 炭化水素の標準化学ポテンシャルアルカンシクロアルカンアルキルベンゼン *Hermann J. Phys. Chem. 76, 754 (97). **Sitkoff et al. J. Phys. Chem. 98, 978 (994). 3

32 参考 : ソフトウェア DelPhi Poisson-Boltzmann 方程式を解き静電ポテンシャルを計算する無料 Software:DelPhi AMBER Generalized Born モデルを使用した分子動力学シミュレーションが可能有料 CHARMM Generalized Born モデルを使用した分子動力学シミュレーションが可能有料 3

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

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

35 初期構造からのずれ (RMSD). ChimeraのMD Movieでmin.dcd min.dcd eq.dcd eq.dcd prod.dcdを開く. Select Atom Specifier を選択し Atom Specifier to と入力し OK Cα 原子を選択 3. MD Movieのメニューの Analysis Plot RMSD を選択 4. Ignore hydrogensを false にし Plot RMSD [Å] Step 35

36 課題初期構造からの Cα 原子のずれ (RMSD) の時間変化をプロットせよ平衡化 (eq eq) とプロダクション (prod) における温度 (TEMP) と圧力 (PRESSURE) 体積 (VOLUME) の時間変化をプロットせよ時間刻み Δt は fs CRN 分子内の水素結合をつ以上についてこの水素結合距離の時間変化をプロットせよ ( どの残基のどの原子間か明示すること ) これらのプロットから何が言えるか考察せよ 36

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

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

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

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

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

42 NAMD における設定 () SHAKEを使う場合は以下の設定を行う rigidbonds all usegrouppressure yes 4

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

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

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

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

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

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

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

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

51 Particle Mesh Ewald 法 () 点電荷を以下のつの電荷分布に分ける = + 点電荷 ρ σ ガウス分布に従う電荷分布 ( ) i r = q i exp i πσ σ 3 残りの電荷分布 r r 5

52 Particle Mesh Ewald 法 () ガウス分布に従う電荷分布はなめらか高速 Fourier 変換を用いてPoisson 方程式と解き静電ポテンシャルを求める σ ϕ r = 4π ρ r ( ) ( ) ~ 4π ϕ ~ k Fourier 変換 σ ( k) = ρ ( k) 発散を防ぐため全電荷は 0 にする必要がある 5

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

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

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

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

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

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

58 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 δ 位にプロトン化 ε 位にプロトン化 δ, ε 位にプロトン化 His 側鎖の pk a は中性付近であるためつの窒素原子とも水素原子が結合した状態も十分にとりうる His 周りの水素結合ネットワークからプロトン化状態がわかる 58

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

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

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

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

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

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