Microsoft PowerPoint - molsim-okamoto-07I

Size: px

Start display at page:

Download "Microsoft PowerPoint - molsim-okamoto-07I"

ひさともたかひ
7 years ago
Views:

1 第 11 回分子シミュレーション夏の学校 August 1-3, 2007 拡張アンサンブル法と生体分子シミュレーション Yuko OKAMOTO( 岡本祐幸 ) 名古屋大学大学院理学研究科物理学教室 okamoto{a}phys.nagoya-u.ac.jp URL:

2 Contents (I) Protein Structures and Genomes (II) Computer Simulations in Canonical Ensemble and in Generalized Ensemble (III) Generalized-Ensemble Simulations of Protein Folding

3 (I) Protein Structures and Genomes

4 CELL Genetic information is stored in the chromosomes in the nucleus. Chromosomes are composed of double helices of DNA. From: J.E. Darnell, Jr., The Processing of RNA, Copyright 1983 by Scientific American, Inc.

5 uman Genome 3 billion base pairs in total April, 2003 Declaration of End of uman Genome Reading About 30,000 genes in total [one gene corresponds to one protein] 22nd Chromosome (1999) 34 million base pairs; 545 genes 21st Chromosome (2000) about 34 million base pairs; 225 genes

7 Genetic Code Amino Acid Codon Alanine アラニン (Ala) GCU, GCC, GCA, GCG Glutamate グルタミン酸 (Glu) GAA, GAG Glutamine グルタミン (Gln) CAA, CAG Aspartate アスパラギン酸 (Asp) GAU, GAC Asparagine アスパラギン (Asn) AAU, AAC Leucine ロイシン (Leu) UUA, UUG, CUU, CUC, CUA, CUG Glycine グリシン (Gly) GGU, GGC, GGA, GGG Lysine リジン (Lys) AAA, AAG Serine セリン (Ser) UCU, UCC, UCA, UCG, AGU, AGC Valine バリン (Val) GUU, GUC, GUA, GUG Arginineアルギニン (Arg) CGU, CGC, CGA, CGG, AGA, AGG Threonine トレオニン (Thr) ACU, ACC, ACA, ACG Proline プロリン (Pro) CCU, CCC, CCA, CCG Isoleucine イソロイシン (Ile) AUU, AUC, AUA Methionine メチオニン (Met) AUG(Start codon) Phenylalanine フェニルアラニン (Phe) UUU, UUC Tyrosine チロシン (Tyr) UAU, UAC Cysteine システイン (Cys) UGU, UGC Tryptophan トリプトファン (Trp) UGG istidine ヒスチジン (is) CAU, CAC UAA, UAG, UGA (Stop codon)

8 One of the Most Important Problems in the Post Genome Era From the DNA sequence that has been determined by the human genome problem, predict the three-dimensional structure of the corresponding protein and its biological function.

9 AMINO ACID N C R C O O R = SIDE CAIN

10 EXAMPLES OF AMINO ACID SIDE CAINS C χ 1 χ 3 O C C C C C χ 2 C C χ 1 a) b) GLY ALA TYR c)

11 EXAMPLES OF AMINO ACID SIDE CAINS C S χ 4 χ 4 O C O χ 2 S C C χ 2 C χ 3 C χ 2 C χ 3 χ 1 χ 1 χ 1 d) e) f) CYS MET GLU

12 N C R 1 C O O + N C R 2 C O O 2 O N C R 1 O C N C R 2 C O O PEPTIDE BOND

13 PROTEIN R 1 O R 2 O R N N 2 C C N C C N C φ 1 ψ 1 ω 1 φ 2 ψ 2 ω 2 φ Ν ψ Ν C O O Protein

14 MYOGLOBIN (N=153) Mb-1

15 MYOGLOBIN (N=153) α-elix Mb-2

16 IMMUNOGLOBULIN (N=216) β-seet Igg

17 TPI TRIOSE POSPATE ISOMERASE (N=248) α-elix and β-seet TPI

18 Three-Dimensional Structure of Protein (Shape) Closely Related to Its Biological Function

19 3D Structure of Proteins Their Biological Functions ( 酵素反応酸素の運搬貯蔵等 ) E.g., Myoglobin O 2 or CO function

20 Purpose Prediction of the Three-Dimensional Structures of Proteins from the First Principles (Try to predict the native protein structures by a simulation from a random initial conformation, using only the amino-acid sequence information)

21 Interest of Basic Research Understanding Protein Folding Mechanism (Protein Folding Problem) Understanding the Mechanism of ow Proteins Carry Out Their Biological Functions Applications New Drug Design De Novo Design of Artificial Proteins with Specific Functions Study of Diseases that are Caused by Misfolding of Proteins(e.g., Mad Cow Disease, Altzheimer s Disease)

22 Af-1 Anfinsen s Experiment (1960 s)

23 Af-2 Anfinsen s Experiment (1960 s) Add Denaturant

24 Af-3 Anfinsen s Experiment (1960 s) Add Denaturant Remove Denaturant

25 ANFINSEN S DOGMA Three-Dimensional Structures of Proteins are Determined by Their Amino-Acid Sequence Information (The State of the Lowest Free Energy) Thus, once the correct energy function of the protein system is given, one should be able to predict the native three-dimensional structure of protein by computer simulations A-dogma

26 Our Strategy Solvation Theory Electronic Structure Theory Simulation Algorithm Energy Function of Protein System Protein Tertiary Structure Prediction from the First Principles

27 Strategy One Use accurate energy function Conformational energy of protein itself (bond stretching + bond bending + torsion + Lennard-Jones + electrostatics, etc.) E total = K R ( R R 0 ) 2 + K θ ( θ θ 0 ) 2 + bonds angles V n [ ] 1 + cos ( nφ γ ) + 2 dihedrals i< j A ij 12 B ij 6 R ij R ij + q iq j εr ij + Interaction energy with solvent Strg-1

28 Solvent Effects Level 1 Distance-Dependent Dielectric Function ε(r) (Crude but computationally cheap) Level 2 Term Proportional to Solvent Accessible Surface Area Level 3 Explicit Solvent Molecules Included (Accurate but time-consuming) *Statistical Mechanical Theory of Liquids, etc. Reference Interaction Site Model (RISM) Scaled-Particle Theory + Poisson-Boltzmann *Direct Inclusion of Water Molecules

29 Strategy Two Use powerful simulation algorithms that do not get trapped in states of energy local minima Strg-3

30 質量 m の質点の運動ある時刻力 f 座標 q (t = 0) t 秒後座標 q (t = t) ニュートンの運動方程式 E mq = = f q

31 力が重力の場合ニュートンの運動方程式ある時刻 2 dq m 2 重力 f dt E = = q f 加速度 at () 2 dq = = g= dt GM R 2 2 t 秒後速度座標 dq vt () = () t = gt dt 1 2 qt () = gt 2

32 力が電気力の場合ニュートンの運動方程式ある時刻電気力 f m 2 dq 2 dt E = = q f 加速度 Q+ t 秒後電荷 Q at dq dt r = q q 2 + () = = kqq r

33 計算機シミュレーションニュートン方程式を数値的に解く dq q( t+ Δt) q( t) q( t+δt) q( t) vt () = () t = lim = dt Δ t 0 Δt Δt 2 dq vt+δt vt qt+ Δt qt+δ t+ qt ( ) ( ) ( 2 ) 2 ( ) ( ) at () = = lim = 2 Δ 0 2 dt t Δt ( Δt) ニュートンの運動方程式の数値解 2 dq 2 dt E = = q f () t qt ( + 2 Δ t) = 2 qt ( +Δt) qt ( ) + ( Δt) 2 m m f

34 (II) Computer Simulations in Canonical Ensemble and in Generalized Ensemble

35 Microcanonical Ensemble ミクロカノニカルアンサンブル Isolated System: E tot = const 孤立系 : E tot = 一定

36 Canonical Ensemble カノニカルアンサンブル System in eat Bath (Exchange Energy w/ eat Bath) 熱浴中の系 ( 熱浴とエネルギーをやりとり ): T = const = 一定

37 n(e) Canonical Ensemble at Temperature T P B (E) = n(e)w B (E) Density of States E W B (E) = exp(- β E ) E E Boltzmann Factor Canonical Probability Distribution cano

38 MONTE CARLO State x x = q1, q2, q3,, qn; p1, p2, p3,, pn { } Generate states one after another. ( 1) ( 2) ( 3) ( v) ( v+ 1) x x x x = x x = x Suppose at the ν-th step the state was x j and the candidate for the (ν+1)-th step is x k. Suppose also that the transition probability w satisfies the following: Markov Chain ( ( v) ( v 1) ; ( 1), ( 2),, ( v 1) = ) ( ) j + = k = j k w x x x x x x x w x x j k

39 By definition, the transition probability w satisfies ( v+ 1 ) ( ) ( v ) ( ) ( ) P x = P x w x x k j j k j where P (ν) (x) is the probability distribution of state x at the ν-th step. The equilibrium probability distribution P eq (x) should satisfy P x = P x w x x ( ) ( ) ( ) eq k eq j j k j A necessary condition for this to be satisfied is to impose a detailed balance condition: ( ) ( ) = ( ) ( ) P x w x x P x w x x eq j j k eq k k j

40 One solution to this detailed balance condition is (Metropolis method): wx ( x) j 1, if P ( x ) P ( x ) eq k eq j P ( x ) = k eq k < eq k P ( x ) eq j P ( x ) min 1, eq k = P ( x ), if P ( x ) P ( x ) eq j eq j

41 In Canonical Ensemble, where the equilibrium probability distribution is the Boltzmann weight factor, we have 1 1 Peq ( x) = WB ( x; T) = e Z Z eq eq ( βδe ) βe( x) where Z = dxw ( x; T) and β = P P ( x ) k ( x ) j 1 kt βδe = e where ΔE Ex ( ) Ex ( ) 1, if ΔE 0 wx ( j xk) = βδe e, if ΔE > 0 = min 1,e B B k j

42 MOLECULAR DYNAMICS Newton equations of motion: Microcanonical Ensemble mq i E = = q i Nose s method: Canonical Ensemble at temperature T E s s mq = = i mqi fi mqi qi s s 2 s Qs = s mq i 3NkBT + Q i s f i 2

43 Example of MC simulations: Lennard-Jones fluid E LJ σ σ = 4ε i< j r ij r ij r = q q where x = q, q,, q { } 1 2 N 12 6 ij i j where qi qi +Δqi Δ q = ξρ i ρ= ( ρ,ρ,ρ ) x random numbers 0.5 < ρ,ρ,ρ < 0.5 x y y z z Trial Move: q old q new

44 Met-Enkephalin: Global Minimum Structure in Gas Phase Tyr-Gly-Gly-Phe-Met (N = 5) E = -12 kcal/mol

45 Canonical 1000K E MC Sweeps

46 Canonical 600K E MC Sweeps

47 Canonical 300K E MC Sweeps

48 Canonical 50K E MC Sweeps

49 徐冷法 (Simulated Annealing) S. Kirkpatrick, C. Gelatt, Jr. & M. Vecchi, Science 220, 671 (1983). Reproduce a Crystal-Making Process on a Computer SA-2

50 Application of Simulated Annealing to Systems of Biopolymers. Kawai, T. Kikuchi & Y.O., Protein Eng. 3, 85 (1989). See also: S. Wilson, W. Cui, J. Moskovitz & K. Schmidt, Tetrahedron Lett. 29, 4373 (1988). C. Wilson & S. Doniach, Proteins 6, 193 (1989). A. Brunger, J. Mol. Biol. 203, 803 (1988). M. Nilges, G. Clore & A. Gronenborn, FEBS Lett. 229, 317 (1988). For a review see: Y.O., Recent Res. Devel. In Pure & Applied Chem. 2, 1 (1998).

51 Simulated Annealing( 徐冷法 ) P B (E) = n(e)w B (E) P B (E) igh T E E P B (E) Low T E E min SA-2

52 Simulated Annealing: T = 1000 K 50 K Canonical 1000K Canonical 50K E MC Sweeps

53 C-peptide-1 C-Peptide (N = 13): X-Ray C-Peptide of Ribonuclease A Amino-Acid Sequence: Lys-Glu-Thr-Ala-Ala-Ala-Lys-Phe-Glu-Arg- Gln-is-Met

54 M. Masuya & Y.O., unpublished. Cp-first

64 M. Masuya & Y.O., unpublished. Cp-last

65 C-peptide-4 C-Peptide (N = 13): X-Ray Level 0 (Gas Phase) RMSD = 1.9 A Level 1 Level 2 RMSD = 1.4 A RMSD = 0.8 A

66 分子動力学シミュレーション (Molecular Dynamics) ニュートン方程式を数値的に解く dq q( t+ Δt) q( t) q( t+δt) q( t) vt () = () t = lim = dt Δ t 0 Δt Δt 2 dq vt+δt vt qt+ Δt qt+δ t+ qt ( ) ( ) ( 2 ) 2 ( ) ( ) at () = = lim = 2 Δ 0 2 dt t Δt ( Δt) ニュートンの運動方程式の数値解 m 2 dq 2 dt E = = q f 2 f() t qt ( + 2 Δ t) = 2 qt ( +Δt) qt ( ) + ( Δt) m ( 速度 ) ベルレ法などいろいろな数値解 ( 積分法 ) がある

67 2 f() t = mω q() t qt ( + Δ t) = qt ( ) +Δtvt ( ) 調和振動子 f() t vt t vt tat vt t vt tω qt m η qt () 1 t () t ( t t) M () t () t 2 vt () η η Δ = +Δ = = ω Δt 1 η 2 ( +Δ ) = () +Δ () = () +Δ = () Δ () det M 1 ( ) 2 2 = + ω Δt だんだん発散 qt ( + Δ t) = qt ( ) +Δtvt ( + Δt) vt ω 2 ( +Δ t) = vt ( ) Δt qt ( ) det M = 1 M ω ( Δt) Δt = 2 ω Δt 1 安定な周期運動 J 0 1 = -1 0 シンプレクティック条件 MJM T =J 位相空間体積保存 det M=1

68 剛体分子の定温シンプレクティック分子動力学シミュレーション法の開発. Okumura, S.G. Itoh & Y.O., J. Chem. Phys. 126, (17 pages) (2007).

69 カノニカル MD 剛体分子 MD の発展非シンプレクテック位相空間体積非保存非シンプレクテック位相空間体積保存シンプレクテック位相空間体積保存カノニカル能勢の熱浴 (1984) 予測子修正子法能勢 oover 熱浴 Martyna ら (1996) 能勢 Poincaré 熱浴 Bond ら ( 陰的 1999) 能勢 ( 陽的 2001) 剛体分子 4 元数予測子修正子法 4 元数松林中原 (1999) 4 元数 Miller ら (2002) 本研究能勢 Poincaré 熱浴 Miller らの剛体 MD 組み合わせ初めて

70 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD 初めてハミルトニアン 2 2 p' i 1 T T PS = s + ( ) ( ) (, ) ln 2 π' 2 i Sqi DS i qi π' i + E rq+ + gkt s i 2ms i i 8s 2Q 0 0 運動方程式 pi r i = m 1 q = S ( q ) ω 2 Ps s = s Q (4) i i i s p i = Fi pi s s Iωi i = Ni ωi ( Iiωi) Iiωi s P p = + ω Ι ω gkt 2 i T s i ι i i m i 0

71 能勢 Poincaré 熱浴 +シンプレクティック剛体 MD 初めてハミルトニアンの分割 2 p' 1 ( T ) 2 i = s + ln 2 π' 2 i P qi + gkt s i 2ms i 8I1s 2 i ( T ) 2 1 ( T ) 2 sp + π' i P2qi + π' i P3qi + se(, r q) + 8Is 8Is 2Q i 時間発展演算子の分割 Δt Δt Δt Δt exp exp exp exp exp exp Δt Δt Δt Δt exp D exp D exp exp 2 D 3 D 4 + O Δt [ DΔ t] = D D D D [ D Δt] S 3 ( )

72 3. 水分子系への応用水分子 : 80 分子 TIP3Pモデル温度 :300 K, 密度 :0.997 g/cm 3 周期境界条件静電力の計算 : Ewald 法計算時間 : 1.5ns 時間ステップ : Δt = 2fs, 3fs, 4fs, 5fs a) 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD b) 能勢 Poincaré 熱浴 + 非シンプレクティック剛体 MD c) 能勢 oover 熱浴 + シンプレクティック剛体 MD d) 能勢 oover 熱浴 + 非シンプレクティック剛体 MD

73 熱浴結果ンプレクティック非シンプレクティック剛体非シンプレクティックシΔt = 2 fs 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD ハミルトニアンよく保存

74 熱浴結果ンプレクティック非シンプレクティック剛体非シンプレクティックシΔt = 3 fs 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD ハミルトニアンよく保存

75 熱浴結果ンプレクティック非シンプレクティック剛体非シンプレクティックシΔt = 4 fs 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD ハミルトニアンよく保存

76 熱浴結果ンプレクティック非シンプレクティック剛体非シンプレクティックシΔt = 5 fs 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD ハミルトニアンやや増加

77 剛体系のシンプレクティック法のまとめカノニカル (NVT) アンサンブルにおける剛体分子の陽的シンプレクティック MD 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD 初めてシンプレクティック MD の利点ハミルトニアンに近い保存量が存在ハミルトニアンよく保存 :Δt=4 fs まで O.K. 従来の非シンプレクティックな方法典型的な時間ステップ Δt=0.5~2 fs +Andersen の定圧法定温定圧アンサンブルでも可. Okumura, S.G. Itoh, and Y.O., J. Chem. Phys. 126, (2007)

Perl + α. : DNA, mrna,,

2009 Perl + α. : DNA, mrna,, DNA .. DNA A C G T DNA 2 A-T, C-G DNA NH 2 NH 2 O - O O N P O - O CH 2 O N N O - O P O CH 2 O N O - O O P O NH 2 O - O - N CH 2 O N O OH OH OH DNA or RNA (U) (A) (G) (C)