第 11 回分子シミュレーション夏の学校 August 1-3, 2007 拡張アンサンブル法と 生体分子シミュレーション Yuko OKAMOTO( 岡本祐幸 ) 名古屋大学 大学院理学研究科 物理学教室 e-mail: okamoto{a}phys.nagoya-u.ac.jp URL: http://jegog.phys.nagoya-u.ac.jp/tb/
Contents (I) Protein Structures and Genomes (II) Computer Simulations in Canonical Ensemble and in Generalized Ensemble (III) Generalized-Ensemble Simulations of Protein Folding
(I) Protein Structures and Genomes
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.
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
CELL Proteins are Synthesized at Ribosome of the Cell From: J.E. Darnell, Jr., The Processing of RNA, Copyright 1983 by Scientific American, Inc.
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)
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.
AMINO ACID N C R C O O R = SIDE CAIN
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)
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
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
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
MYOGLOBIN (N=153) Mb-1
MYOGLOBIN (N=153) α-elix Mb-2
IMMUNOGLOBULIN (N=216) β-seet Igg
TPI TRIOSE POSPATE ISOMERASE (N=248) α-elix and β-seet TPI
Three-Dimensional Structure of Protein (Shape) Closely Related to Its Biological Function
3D Structure of Proteins Their Biological Functions ( 酵素反応 酸素の運搬 貯蔵等 ) E.g., Myoglobin O 2 or CO function
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)
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)
Af-1 Anfinsen s Experiment (1960 s)
Af-2 Anfinsen s Experiment (1960 s) Add Denaturant
Af-3 Anfinsen s Experiment (1960 s) Add Denaturant Remove Denaturant
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
Our Strategy Solvation Theory Electronic Structure Theory Simulation Algorithm Energy Function of Protein System Protein Tertiary Structure Prediction from the First Principles
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
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
Strategy Two Use powerful simulation algorithms that do not get trapped in states of energy local minima Strg-3
質量 m の質点の運動 ある時刻 力 f 座標 q (t = 0) t 秒後 座標 q (t = t) ニュートンの運動方程式 E mq = = f q
力が重力の場合 ニュートンの運動方程式 ある時刻 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
力が電気力の場合 ニュートンの運動方程式 ある時刻 電気力 f m 2 dq 2 dt E = = q f 加速度 Q+ t 秒後 電荷 Q at dq dt r = q q 2 + () = = 2 2 + kqq r
計算機シミュレーション ニュートン方程式を数値的に解く 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
(II) Computer Simulations in Canonical Ensemble and in Generalized Ensemble
Microcanonical Ensemble ミクロカノニカルアンサンブル Isolated System: E tot = const 孤立系 : E tot = 一定
Canonical Ensemble カノニカルアンサンブル System in eat Bath (Exchange Energy w/ eat Bath) 熱浴中の系 ( 熱浴とエネルギーをやりとり ): T = const = 一定
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
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
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
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
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
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
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
Met-Enkephalin: Global Minimum Structure in Gas Phase Tyr-Gly-Gly-Phe-Met (N = 5) E = -12 kcal/mol
Canonical 1000K 30 20 E 10 0-10 0 50000 100000 150000 200000 MC Sweeps
Canonical 600K 30 20 E 10 0-10 0 50000 100000 150000 200000 MC Sweeps
Canonical 300K 30 20 E 10 0-10 0 50000 100000 150000 200000 MC Sweeps
Canonical 50K 30 20 E 10 0-10 0 50000 100000 150000 200000 MC Sweeps
徐冷法 (Simulated Annealing) S. Kirkpatrick, C. Gelatt, Jr. & M. Vecchi, Science 220, 671 (1983). Reproduce a Crystal-Making Process on a Computer SA-2
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).
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
Simulated Annealing: T = 1000 K 50 K Canonical 1000K Canonical 50K 30 20 E 10 0-10 0 50000 100000 150000 200000 MC Sweeps
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
M. Masuya & Y.O., unpublished. Cp-first
M. Masuya & Y.O., unpublished. Cp-last
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
分子動力学シミュレーション (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 ( 速度 ) ベルレ法などいろいろな数値解 ( 積分法 ) がある
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 2 2 1 ω ( Δt) Δt = 2 ω Δt 1 安定な周期運動 J 0 1 = -1 0 シンプレクティック条件 MJM T =J 位相空間体積保存 det M=1
剛体分子の定温シンプレクティック分子動力学シミュレーション法の開発. Okumura, S.G. Itoh & Y.O., J. Chem. Phys. 126, 084103 (17 pages) (2007).
カノニカル MD 剛体分子 MD の発展 非シンプレクテック位相空間体積非保存非シンプレクテック位相空間体積保存シンプレクテック位相空間体積保存 カノニカル 能勢の熱浴 (1984) 予測子 修正子法能勢 oover 熱浴 Martyna ら (1996) 能勢 Poincaré 熱浴 Bond ら ( 陰的 1999) 能勢 ( 陽的 2001) 剛体分子 4 元数予測子 修正子法 4 元数松林 中原 (1999) 4 元数 Miller ら (2002) 本研究 能勢 Poincaré 熱浴 Miller らの剛体 MD 組み合わせ初めて
能勢 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
能勢 Poincaré 熱浴 +シンプレクティック剛体 MD 初めてハミルトニアンの分割 2 p' 1 ( T ) 2 i = s + ln 2 π' 2 i P qi + gkt s i 2ms i 8I1s 2 i 3 1 0 0 1 ( T ) 2 1 ( T ) 2 sp + π' i P2qi + π' i P3qi + se(, r q) + 8Is 8Is 2Q i 1 2 3 5 時間発展演算子の分割 Δt Δt Δt Δt exp exp exp exp exp exp 2 2 2 2 Δt Δt Δt Δt exp D exp D exp exp 2 D 3 D 4 + O Δt 5 2 2 2 2 [ DΔ t] = D D D D [ D Δt] 4 5 4 3 2 1 2 S 3 ( )
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
熱浴結果 ンプレクティック非シンプレクティック剛体非シンプレクティックシΔt = 2 fs 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD ハミルトニアンよく保存
熱浴結果 ンプレクティック非シンプレクティック剛体非シンプレクティックシΔt = 3 fs 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD ハミルトニアンよく保存
熱浴結果 ンプレクティック非シンプレクティック剛体非シンプレクティックシΔt = 4 fs 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD ハミルトニアンよく保存
熱浴結果 ンプレクティック非シンプレクティック剛体非シンプレクティックシΔt = 5 fs 能勢 Poincaré 熱浴 + シンプレクティック剛体 MD ハミルトニアンやや増加
剛体系のシンプレクティック法のまとめ カノニカル (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, 084103 (2007)