計算生命科学の基礎10（奥村 久士）

生体分子動力学シミュレーションの 基礎と応用 Hisashi OKUMURA Exploratory Research Center on Life and Living Systems Institute for Molecular Science The Graduate University for Advanced Studies Japan 奥村久士 自然科学研究機構 生命創成探究センター 分子科学研究所 計算科学研究センター 総合研究大学院大学 構造分子科学専攻

Books to read 参考書 Computer Simulation of Liquids M. P. Allen, D. J. Tildesley Oxford Science Publications (2017/8) Understanding Molecular Simulation (Third Edition) From Algorithms to Applications Daan Frenkel and Berend Smit Academic Press (2023/7) 分子シミュレーション ―古典系から量子系手法まで 上田 顕 裳華房 (2003/10) コンピュータ・シミュレーションの基礎（第2版）: 分子のミクロな性質を解明するために 岡崎 進, 吉井 範行 化学同人 (2011/7)

分子動力学シミュレーションの講義ビデオ 「分子動力学シミュレーション」で動画検索 https://www.youtube.com/watch?v=6B3BE7-iIPk テキストも公開

• https://www.youtube.com/watch?v=6B3BE7-iIPk

内容 １．生体分子動力学シミュレーションの基礎 ２．温度・圧力の制御 ３．タンパク質の分子動力学シミュレーションの例

内容 １．生体分子動力学シミュレーションの基礎 ２．温度・圧力の制御 ３．タンパク質の分子動力学シミュレーションの例

二体問題 力学 ２つの質点の運動を解く問題。 独立な２つの一体問題として解くことができる。 ニュートンの運動方程式 F12 Fi = mi ri F21＝ - F12 m1 F=F 2 dr d r , r= 2 上付き点・は時間微分 r = dt dt 太字はベクトル m2

二体問題 力学 ２つの質点の運動を解く問題。 独立な２つの一体問題として解くことができる。 ニュートンの運動方程式 F12 = m1r1 m1r1 + m2 r2 = F12 + F21 ( m1 + m2 ) r重心 = 0 r重心 m1r1 + m2 r2 = m1 + m2 重心は慣性運動 Fi = mi ri F21 = m2 r2 F12 F21  1 1  r12 = r1 − r2 = − = +  F12 m1 m2  m1 m2   r12 = F12 換算質量 = 1  1 1    +  m1 m2 

三体問題 力学 ３つの質点の運動を解く問題。 特殊な条件の場合しか解析的に解くことはできない。 ニュートンの運動方程式 F12 m1 Fi = mi ri F21 m2 F23 F13 F32 F31 m3

分子動力学シミュレーション 計算機を使って原子、分子の動きを 時々刻々数値的に解く。 Fi = mi ri ニュートンの運動方程式 v2 m1 v1 v1 m m1 1 m3 mv33 v3 m2 m3 v2 ri = vi Fi vi = mi m2 m2

分子動力学シミュレーション 計算機を使って原子、分子の動きを 時々刻々数値的に解く。 初期条件  座標 r1 , r2 , , rN  速度 v1 , v2 , , v N 運動方程式 ri = vi Fi vi = mi Fi = mi ri

分子動力学 シミュレーションの例 • N = 108 • V = 5.0×5.0×5.0 • T0*= 0.1, 1.0 （LJ単位系） • Lennard-Jones model アルゴン Ar 固体 液体 低温 T0 = 12 K 高温 T0 = 120 K

周期境界条件 一方向から粒 子が飛び出た ら反対側から 箱の中に戻す。 rc

タンパク質 アミノ酸が多数つながったひも状の分子 H H N H O Cα H H H C O R1 H アミノ酸 N C O R2 アミノ酸 H H H Cα N O Cα R1 C O H H N Cα R2 ペプチド結合 O C O H H

20種類のアミノ酸 有坂文雄著 「バイオサイエンスのための蛋白質科学入門」 裳華房 2004 より

タンパク質の構造（形）の例 アルファヘリックス ベータヘアピン アルファヘリックス + ベータヘアピン

Example of a biomolecule: alanine dipeptide 生体分子の例：アラニンジペプチド

生体分子のポテンシャルエネルギー E = Ebond +Eangle +Edihedral+Eelec+ELJ Ebond =  kr ( r − r0 ) 2 r bond Eangle =  k ( −  ) 2  0 angle Ediheadral vn =  1 + cos ( n −  )  diheadral 2 AMBER, CHARMM, OPLS, etc. 

生体分子のポテンシャルエネルギー E = Ebond +Eangle +Edihedral+Eelec+ELJ Eelec qi q j 1 = i , j 4 0 rij  Aij Bij  ELJ =   12 − 6    r r i , j  ij ij  AMBER, CHARMM, OPLS, etc. r

ポテンシャルエネルギーのカットオフ MD シミュレーション ではポテンシャルの カットオフを行う。   12   6  u ( r ) = 4   −     r    r  rc ≦ L/2 いくつかのカットオフ ・不連続 ・連続、微分不可能 ・１階微分可能 ・２階微分可能 ・カットオフなし T. Sakaguchi and H. Okumura: J. Phys. Soc. Jpn. 82 (2013) 034001.

MDシミュレーションの手順 1．初期条件 座標と運動量の初期条件を設定する。 2．平衡化 温度（時には圧力も）を制御しながら、平衡状態に至らせる。 3．サンプリング 長いMDシミュレーションを実行し、多くの微視的状態を サンプルする。統計サンサンブル平均を計算する。

シンプレクティック分子動力学法 ハミルトニアン pi2 H = + V (q) i 2 mi 正準方程式 H qi = , pi H pi = − qi 物理量 Z (q, p) の時間変化 dZ   H  H    =   −   Z = DH Z dt  i  pi qi qi pi   形式解 Z ( t + t ) = e DH t Z ( t ) 時間発展演算子

シンプレクティック分子動力学法 時間発展演算子の分割 ハミルトニアンが H = K ( p ) + V ( q )と分割される場合  H  H    K  V   DH =   − −  =   = DK + DV qi pi  i  pi qi qi pi  i  pi qi 近似 : e DH t =e DV t 2 e DV t = 1 + DV t + DV pi = − e DK t e DV t 2 1 2 2 DV t + 2 V  pi = Fi qi pi + O ( t 3 ) 1 2 2 DK t + 2 K  p DK qi = qi = i pi qi mi e DK t = 1 + DK t + V  D pi = − Fi = 0 qi pi K  pi D qi = =0 pi qi mi V  DV qi = − qi = 0 qi pi K  DK pi = pi = 0 pi qi 2 V 2 K pi2 K = i 2mi

シンプレクティック分子動力学法 近似 : e DH t = e DV t 2 速度ベルレ法 t e DK t e DV t 2 + O ( t 3 ) t e operation : p = p + Fi 2 n +1/2 p e DK t operation : qin +1 = qin + i t mi DV e DV 2 t 2 operation : n +1/2 i n i n pin +1 = pin +1/2 + Fi n +1 t 2 影のハミルトニアン   2 H H H 1  2 H H H H = H −   −2  24 i j  pi p j qi q j qi q j pi p j Hに近い保存量 H が存在する.  2 t +  

シンプレクティック分子動力学法 近似 : e DH t = e DK t 2 e DV t e DK t 2 + O ( t 3 ) 位置ベルレ法 e DK t 2 operation : qin +1/2 e DV operation : e DK t 2 n p t = qin + i mi 2 pin +1 = pin + Fi n +1/2 t n +1 i operation : q n +1/2 i =q pin +1 t + mi 2 影のハミルトニアン   2 H H H 1  2 H H H H = H −   −2  24 i j  qi q j pi p j pi p j qi q j Hに近い保存量 H が存在する.  2 t +  

例：水 シンプレクティック 能勢-Poincaré 熱浴 非シンプレクティック 能勢-Hoover 熱浴 能勢-Poincaré 熱浴ではハミルトニアンはよく保存する。 影のハミルトニアン：Hに近い保存量   2 H H H 1  2 H H H  2 H = H −   −2 t +  24 i j  pi p j qi q j qi q j pi p j 

シンプレクティック条件 • 正準変換 q Q  =   →  =   ,  =  ( )  p P • ハミルトンの運動方程式 =J H H ,  =J    0 1 J =  -1 0   • ζの時間発展  H  H T H  ( ) =  = M = MJ = MJ = MJM   • シンプレクティック条件 ζ=ζ(η) が正準変換 ⇔ MJMT=J    M ij =  i  j

誤差の評価 MDシミュレーションプログラムのデバッグにはエネル ギー保存則を使うのが便利。エネルギーの誤差ΔH の時間ステップ幅Δt 依存性が積分法から予想される 次数になっているかチェックする。 ΔH ∼ (Δt)k もし次数k が予想通りなら力とエネルギーがコンシステ ントに計算されている。 ローカルエラー 1ステップでの誤差 δH = |Hn+1 − Hn| Hn · · · nステップ目でのハミルトニアンの値 ベルレ法では δH = O((Δt)3)

グローバルエラー 時刻 t = 0 と後の固定した時刻 t = nΔt とでの誤差 ΔH = |H (nΔt) − H (0)| nΔt を一定に保ちながら、Δt を細かくする。(Δt に 反比例するようにn を増やす) 例えば、速度ベルレ法では ΔH ∝ nδH = nO ((Δt)3) = O ( (Δt)−1 ) O ((Δt) 3) = O ((Δt) 2) 通常グローバルエラーの次数はローカルエラーよりも 1次低い。

t = nΔt を一定に保ちながら、 Δt を細かくしてMDシミュレーシ ョンを行い、ハミルトニアンのグローバルエラー ΔH を調べる。 ΔH と Δt の対数を取り、傾きから誤差の次数を計算する。 ΔH = Δt 2 logΔH = log (Δt 2) = 2 log (Δt)

内容 １．生体分子動力学シミュレーションの基礎 ２．温度・圧力の制御 ３．タンパク質の分子動力学シミュレーションの例

能勢修一教授 Prof. Shūichi Nosé (1951-2005) 能勢-Hoover 熱浴 Nosé-Hoover thermostat pi ri = mi pi = Fi −  pi  1 pi2  =   − 3NkT0  Q  i mi 

能勢の熱浴 H NVT p'i2 PS 2 = + E (r ) + + gkT0 ln s 2 2Q i 2ms 運動量と時間をスケール 仮想時間t’での運動方程式 ・・・正準方程式 p'i  dri  dt' = ms 2 変数変換  p'i  dp'i = F p = i i  dt' s  dt' P ds dt =  = s s  dt' Q   p'i2  dPs = 1  − gkT0  2  dt' s    i ms   S. Nosé: Mol. Phys. 52 (1984) 255 S. Nosé: J. Chem. Phys. 81 (1984) 511 p'i dt' pi = , dt = s s 実時間tでの運動方程式 ・・・正準方程式ではない pi   ri = m  p =F −s p i i  i s   s = s Ps  Q  2 p P = i  s  m − gkT0 i  能勢-Hoover熱浴   ri = pi m   p = F − p i i  i   pi2 1 − gkT0   =   Q i m  

瞬間温度を用いると gkT (t ) =  i pi2 m 最後の式は以下のように書き換えられる。 gk  = T (t ) − T0  Q フィードバックメカニズム T (t )  T0 →  増加 → p 減少 → T (t ) 減少 T (t )  T0 →  減少 → p 増加 → T (t ) 増加 熱浴の質量 Q 大 →  のダイナミクスは遅くなる Q 小 →  のダイナミクスは速くなる

Lennard-Jones 流体 への応用 • 周期境界条件 • N = 500, ρ = 0.8 • rc = 4.0 • T0 = 1.0 • T (t = 0) = 0.5 Lennard-Jones単位系 (ϵ = σ = 1).   12   6  u ( r ) = 4   −     r    r 

能勢-Hoover 熱浴における時間発展 物理量 A (r, p, ) の時間発展 A A A A ( r, p, ) =  ri  +  pi  +  ri pi  i i 時間発展演算子（Liouville 演算子）    D =  ri  +  pi  +  ri pi  i i 運動方程式から   pi   1 pi2 D =   +  ( Fi −  pi )  +   − gkT0   pi Q  i m i mi ri i   A ( r, p, ) = DA 形式解 A ( t + t ) = e Dt A ( t )

時間発展演算子の分割 D = D1 + D2 + D3   pi  1  pi2 D1 =   +  − gkT0   Q i m i mi ri    D2 =  F  pi i  D3 = −   pi  pi i 鈴木-Trotter 分割 Approximation : e Dt = e 例えば、 e D1t 1 2 2 = 1 + D1t + D1 t + 2 D3 t 2 e D2 t 2 e D1t e D2 e t 2 D1t e D3 t 2 + O ( t 3 ) pi ri = ri + t mi 2   1 p D1t i e  =  +   − gkT0  t Q i m 

D2 についても同様に e D2 t pi = pi + Fi t D3 の演算についてはΔtの高次項がゼロにはならないが 公式 ex = 1 + x + 1 2 x + 2 , を用いると、以下の形になる e D3 t 1 2 2   pi = 1 + D3t + D3 t +  pi 2!   1 2  = pi 1 + ( −t ) + ( −t ) + 2!  = pi e −t   

能勢-Hoover 熱浴の時間発展 t   pi  pi exp  −  2   t pi  pi + Fi 2 pi ri  ri + t mi  1 pi2    +   − gkT0  t Q i m  t pi  pi + Fi 2 t   pi  pi exp  −  2   “←”: プログラムにおける代入 G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M. L. Klein, Mol. Phys., 87, 1117 (1996).

圧力の制御 Andersenの方法 pi V ri = + ri mi 3V V pi = Fi − pi 3V PV V= M  1  pi2 PV = +  ri  Fi  − P0  3V  i mi i  瞬間圧力 P ( t ) 圧力浴 体積の変化

Andersen の方法 H NPH pi2 H. C. Andersen: J. Chem. Phys. 72 (1980) 2384 PV 2 = + PV 2 + E (V r ) + 0 2 M 3 i 2 mV i 1 3 座標と運動量のスケール 1 3 ri = V ri , pi = V − 1 3 pi 正準方程式 pi  d ri = 2  dt 3 mV i  変数変換 dp 1 1 i 3 = V Fi  ri = V 3 ri  dt 1 −  dV PV pi = V 3 pi =   dt M  dPV  pi2 1  = 2 +  Fi  ri  − P     0  dt 3 V i mV 3 i  i    pi V r = + ri  i mi 3V   V  pi = Fi − 3V pi   V = PV  M   1    pi2  PV =    +  Fi  ri  − P0   i  3V  i mi    瞬間圧力 P ( t ) 

フィードバック機構 P (t )  P0 → V  0 → P (t ) が増える P (t )  P0 → V  0 → P (t ) が減る ピストンの質量 M 大 → V のダイナミクスは遅くなる。 M 小 → V のダイナミクスは速くなる。 ハミルトニアン H は保存されるので H = 一定 平衡状態では体積の運動エネルギー PV2/2M は H より はるかに小さい。 PV 2/2M << H ゆえにエンタルピー H はほぼ一定: H ≡ H0 + P0V ≈ 一定 アンダーセンの方法は近似的に NPH 一定（定圧定エン タルピー）のアンサンブルを生成する

Andersen の方法における時間発展 物理量 A ( r , p,V , PV ) の時間発展 A ( r , p,V , PV ) =  ri  i A A A A +  pi  +V  + PV  , ri pi V PV i 時間発展演算子 （Liouville 演算子）     D =  ri  +  pi  +V  + PV  , ri pi V PV i i 運動方程式より pi   PV   1 D= + V  Fi  +  + 2  pi M V  3V i mV 3 ri i i 1 3 A ( r , p,V , PV ) = DA     pi2 , 2 +  Fi  ri  − P    0  3 i  PV  i mV i  形式解 A ( t + t ) = e Dt A ( t )

時間発展演算子の分割 D = D1 + D2 + D3  pi2  D1 =  + 2  5  r 3 i mV i 3mV 3 PV i i i pi PV  D2 =  M V   1 D3 = V  Fi  + pi  3V i 1 3   i Fi  ri − P0  P V 鈴木-Trotter 分割 近似 : e Dt = e D3 t 2 e D2 例えば 1 e D1t = 1 + D1t + D12 t 2 + 2 t 2 e D1t e D2 t 2 e D3 e t 2 D1t + O ( t 3 ) ri = ri + pi mV i 2 3 e D1t PV = PV +  i t pi2 3mV i 5 3 t

Andersen の方法における時間発展 t pi  pi + V Fi 2  t  1 P − r  F PV  PV +   0 i i V 3 i   2 1 3 PV t V V + M 2 pi ri  ri + 2 t 3 mV i PV  PV +  i t pi  pi + V Fi 2  t  1 Fi  ri − P0  PV  PV +    3V i  2 1 3 “←”: プログラムにおける代入 2 i p 3mV i PV t V V + M 2 5 3 t

温度・圧力の制御 定温定圧分布(NPT) 能勢・Andersenの方法 pi V ri = + ri mi 3V s V  pi = Fi −  +  pi  s 3V  P s=s s Q pi2 Ps =  − 3 NkT0 i mi 瞬間温度 3 NkT (t ) 熱浴 エネルギーの出入り 圧力浴 体積の変化 PV V =s M  1    pi2 PV = s    +  ri  Fi  − P0  i  3V  i mi   瞬間圧力 P (t )

内容 １．生体分子動力学シミュレーションの基礎 ２．温度・圧力の制御 ３．タンパク質の分子動力学シミュレーションの例

アミロイドーシス タンパク質が間違った形にフォールディングし、 アミロイド線維を形成することにより引き起こされる病気． ・ 例：アルツハイマー病 認知症の１つ．脳の萎縮．老人斑 老人斑で見られるアミロイドb β1 ペプチドのアミロイド線維 Lührs et al., Proc. Natl. Acad. Sci. USA 102, 17342 (2005) β2 Abアミロイド線維の離合集散 のメカニズムを調べる ↓ 分子動力学シミュレーション

神経細胞毒性を持つものはどれか？ 単量体 オリゴマー アミロイド線維 アミロイド仮説 オリゴマー仮説 ＝ 最近はオリゴマーの 方がより毒性が高い と報告されている。 神経細胞毒性

超音波によるアミロイドβ線維破壊の分子シミュレーション サインカーブ状の圧力 P (t) P (t) = P + P sin(2 t/t) 0 period t = 1 ns 300 MPa P0 P time t -100 MPa • Nosé-Hoover熱浴の質量 P0 = 100 MPa Q = 10 (kcal/mol) ps2 P = 200 MPa t = 1 ns (⇔ 1 GHz) • Andersen圧力浴の質量 M = 10-7 (kcal/mol) ps2 /Å6 20個の初期速度

動画 圧力が正の時は何も起 きない。 気泡は通常、b2の疎水 性残基周辺で生成。 気泡中でもアミロイドは 形状保つ。 気泡が崩壊するとき、 水分子がb1の親水性 残基に衝突し、アミロイ ドが破壊される。 H. Okumura and S. G. Itoh: J. Am. Chem. Soc. 136 (2014) 10549-10552

アミロイド線維の長さ依存性 何回目の負圧で気泡ができたか？ 初期条件20個 Ab×12量体 Ab×6量体 Ab×3量体 気泡生成1回 だけおきた。 アミロイド線維が短いほど気泡できにくい。壊れにくい。 ⇔超音波でアミロイド線維の断片が均質化 阪大後藤研の実験 E Chatani, Y-H Leea, H Yagia, Y Yoshimura, H Naikib, and Y Goto: Proc. Natl. Acad. Sci. USA 106 (2009) 11119.

まとめ １．生体分子動力学シミュレーションの基礎 分子動力学シミュレーションの概要 生体分子のモデル、力場 時間発展手法 ２．温度・圧力の制御 能勢の方法、アンダーセンの方法 ３．タンパク質の分子動力学シミュレーションの例 タンパク質凝集体のシミュレーション