新学術領域研究「コンピューティクスによる物質デザイン：複合相関と非平衡ダイナミクス」
平成25年度第1回研究会
2013年7月8-9日

公募研究（2013-2014)：

拡張アンサンブル密度汎関数分子動力学法によるＡＴＰ加水分解の研究（Cf. 2010-2012: 倭班の連携研究員）
岡本祐幸（OKAMOTO, Yuko）
名古屋大学・大学院理学研究科・物質理学専攻（物理系）
名古屋大学・大学院理学研究科・構造生物学研究センター
名古屋大学・大学院工学研究科・計算科学連携教育研究センター
名古屋大学・情報基盤センター

e-mail: okamoto{a}phys.nagoya-u.ac.jp
URL: http://www.tb.phys.nagoya-u.ac.jp/

お詫び：
新学術領域研究（理工系）「生命分子システムにおける動的秩序形成と高次機能発現」（領域代表：加藤晃一 [岡崎統合バイオサイエンスセンター]）
計画研究：「生体分子集団および人工分子集団の相互作用と大規模構造転換」が先日新規採択された。
新学術領域研究の重複執行は許されないので、これからは、コンピューティクスでは、班友として活動させて頂きたい。

Generalized-Ensemble Algorithm (拡張アンサンブル法)
Simulation Methods that Do Not Get Trapped in States of Local
Energy Minima[e.g., Multicanonical Algorithm, Simulated Tempering,Replica-Exchange Method, etc.]
Based on Non-Boltzmann Weight Factor
Realize random walks in potential energy space and/or any other physicalquantities (OR their conjugate parameters)
Histogram Reweighting Techniques
Can obtain thermodynamic quantities for a wide range of temperature froma single simulation run by reweighting
REVIEWS: U.H.E. Hansmann & Y.O., Curr. Opin. Struct. Biol. 9, 177 (1999);
A. Mitsutake, Y. Sugita, & Y.O., Biopolymers 60, 96 (2001);
Y.O., in Water and Biomolecules, K. Kuwajima, Y. Goto, F. Hirata, M. Kataoka,& M. Terazima (eds.) (Springer-Verlag, Berlin, 2009) pp. 61-95;
H. Okumura, S.G. Itoh, & Y.O., in Practical Aspects of Computational Chemistry II:
An Overview of the Last Two Decades and Current Trends,
J. Leszczynski and M.K. Shukla (eds.) (Springer, Dordrecht, 2012) pp. 69-101;
A. Mitsutake, Y. Mori, and Y.O, in Biomolecular Simulations: Methods and Protocols,L. Monticelli and E. Salonen (eds.) (Humana Press, New York, 2012) pp. 153-195.
岡本祐幸、計算科学講座第９巻「超多自由度系の最適化」第２章（共立, 2013) pp. 119-241.

Multicanonical Algorithm
MC: B. Berg & T. Neuhaus, Phys. Lett. B267, 249 (1991); Phys. Rev. Lett. 68, 9 (1992).
MD: U. Hansmann, Y.O. & F. Eisenmenger, Chem. Phys. Lett. 259, 321 (1996);
N. Nakajima, H. Nakamura & A. Kidera, J. Phys. Chem. B 101, 817 (1997).

Pmu (E) = n(E)Wmu (E) = const
uniform (flat) distribution inenergy-1W mu (E) = n(E)EEmin

Random Walk in Potential
Energy SpaceCanonical Ensemble
MC version:Multicanonical Ensemble
MC version:

Generalized-Ensemble Algorithms have beendeveloped in MC algorithms

MD version:

Canonical Ensemble

s
E s
 mqi   q  s mqi  f i  s mqi
i

2
s


Qs  s   mqi2  3 Nk BT   Q
s
 i


Multicanonical Ensemble
MD version:

1Wmu (E )  exp(  0Emu (E ))n(E )Emu sEmus mqi   q  s mqi  E f i  s mqii
2s2Qs  s   mqi  3 Nk BT0   Qs i

U. Hansmann, Y.O. & F. Eisenmenger, Chem. Phys. Lett. 259, 321 (1996);
N. Nakajima, H. Nakamura & A. Kidera, J. Phys. Chem. B 101, 817 (1997).

Single-Histogram Reweighting Techniques
A. Ferrenberg & R. Swendsen, Phys. Rev. Lett. 61, 2635 (1988).

 A( E ) P  E; T   A( E )n  E e A PE;
TnEe     ECEE



TECE

EHere, the density of states n(E) is obtained from the histogram of theenergy distribution Nmu(E) that was obtained from the production run ofthe multicanonical simulation:

N mu ( E )nE , where N mu ( E )  n  E  Wmu ( E ) .
Wmu ( E )

Single-Histogram Reweighting Techniques
A. Mitsutake, Y. Sugita & Y.O., J. Chem. Phys. 118, 6664 (2003).

When the physical quantity A cannot be written asa function of E, we use the following equation:

MULTICANONICAL ALGORITHM
MC：B. Berg & T. Neuhaus, Phys. Lett. B267, 249 (1991); PRL 68, 9 (1992).
MD：UHE Hansmann, Y. O. & F. Eisenmenger, Chem. Phys. Lett. 259, 321 (1996)
N. Nakajima, H. Nakamura & A. Kidera, J. Phys. Chem. B 101, 817 (1997).

Step 1: Iterations of Short Preliminary Runs toDetermine the Multicanonical Weight Factor Wmu (E)Step 2: One Long Production Run
Step 3: Analyze the Data to Obtain:* Global-Minimum Energy Configuration* Thermodynamic Quantities for Desired Temperatures(by Ferrenberg-Swendsen Single-HistogramReweighting Techniques)

WB  E; T 
Pc  E;T   Pmu  E Wmu  E Multicanonical
Canonical 50K
Canonical 1000K

30E20100-10050000

100000MC Sweeps

150000200000Enk-ave

17-Residue Helical Peptide (120000-300000 MC Sweeps)
Simulation and movie by A. Mitsutake

Canonical MC: T = 200 K

Multicanonical MC

Simulated Tempering （焼き戻し法）
A.P. Lyubartsev, et al., J. Chem. Phys. 96, 1776 (1992).
E. Marinari and G. Parisi, Europhys. Lett. 19, 451 (1992).
Temperature is a dynamical variable： Sample temperature uniformly

WST ( E;T )  exp(  E  a(T ))PST (T )   dEn ( E ) exp (   E  a (T ))  const
Random Walk in Temperature Space→ Random Walk in Energy Space
Discretize Temperature:

WST ( E;Tm )  exp(  m E  am )exp(am ) dEn(E) exp(m E)
am : Dimensionless Helmholtz free energy at temperature Tm
am is determinded by iterations of short ST runs
A. Irback & F. Potthast, J. Chem. Phys. 103, 10298 (1995).
U. Hansmann & Y.O., J. Comput. Chem. 18, 920 (1997).

Simulated Tempering （焼き戻し法）
A.P. Lyubartsev, et al., J. Chem. Phys. 96, 1776 (1992).
E. Marinari and G. Parisi, Europhys. Lett. 19, 451 (1992).

Step 1: Canonical MC/MD Simulations at Temperature Tm
for a Few Steps
Step 2: Temperature is Updated to a Neighboring Value Tm±1
a la Metropolis with Conformations Fixed

whereRepeat These 2 Steps
Canonical Distribution at Any Temperature
by Multiple Histogram Reweighting Techniques (WHAM)

Multiple-Histogram Reweighting Techniques(Weighted Histogram Analysis Method: WHAM)
A. Ferrenberg & R. Swendsen, Phys. Rev. Lett. 63, 1195 (1989).
S. Kumar, D. Bouzida, R. Swendsen, P. Kollman & J. Rosenberg, J. Comput. Chem. 13, 1011 (1992).

 A(E)nE e
 A  nE e

 EE

TEE

Given M set of histograms Nm(E), which were obtained at Tm, thefollowing WHAM equations are solved iteratively for density of states n(E)and dimensionless Helmholtz free energy f m : (nm are the total numberof samples obtained at Tm)
Mn E   N (E)mm 1M nmem 1fm   m E, where e fm   nE e  m E .
E

Multiple-Histogram Reweighting Techniques(Weighted Histogram Analysis Method: WHAM)
A. Mitsutake, Y. Sugita & Y.O., J. Chem. Phys. 118, 6664 (2003).

When the physical quantity A cannot be written as
a function of E, we first obtain the dimensionlessHelmholtz free energy fm (m = 1, …, M) by solving
the WHAM equations. We then use the following
equation:Replica-Exchange Method（レプリカ交換法）(also referred to as Parallel Tempering)
MC: K. Hukushima & K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).
MD: Y. Sugita & Y.O., Chem. Phys. Lett. 314, 141 (1999).

1. System
M Non-Interacting Replicas of the Original System at M Different
Temperatures

2. Replica-Exchange
Step 1: Independent Canonical Simulations Performed for
Each Replica
Step 2: A Pair of Replicas Corresponding to Neighboring
Temperatures are Exchanged a la Metropolis
Repeat These 2 Steps

3. Canonical Distribution at Any Temperature
by Multiple Histogram Reweighting Techniques (WHAM)

Multidimensional Replica-Exchange Method (MREM)
Y. Sugita, A. Kitao & Y.O., J. Chem. Phys. 113, 6042 (2000).

T131415169101112Temperature Exchange
5678Parameter Exchange

1234



Multidimensional ST
A. Mitsutake & Y.O., Phys. Rev. E 79, 047701 (2009);
J. Chem. Phys. 130, 214105 (2009);
A. Mitsutake, J. Chem. Phys. 131, 094105 (2009).

From Multidimensional REM toMultidimensional MUCA and ST
A. Mitsutake & Y.O., Phys. Rev. E 79, 047701 (2009);
J. Chem. Phys. 130, 214105 (2009);
A. Mitsutake, J. Chem. Phys. 131, 094105 (2009).

MMUCA: random walk in multidimensional energy
MST: random walk in multidimensional parameterMREM: random walk in multidimensional parametere.g.,
WHAM eqns.
Mn E ,V Nm ,n 1
Mm ,n nm,nem ,n 1(E ,V )

fm ,n  m (E  nV ), where e m ,n   n E ,V e
fE ,V m (E  nV ).

Examples of Multidimensional REM, MUCA, and ST
A. Mitsutake & Y.O., Phys. Rev. E 79, 047701 (2009);
J. Chem. Phys. 130, 214105 (2009);
A. Mitsutake, J. Chem. Phys. 131, 094105 (2009).

 = h = external field
V = M = magnetization1. Simulated Tempering and Magnetizing
random walk in temperature T and external field h＊Ising Model
T. Nagai & Y.O., Phys. Rev. E 86, 056705 (2012).

＊3-state Potts Model
T. Nagai, Y.O., & W. Janke,
J. Stat. Mech. (2013) P02039.

Examples of Multidimensional REM, MUCA, and ST
A. Mitsutake & Y.O., Phys. Rev. E 79, 047701 (2009);
J. Chem. Phys. 130, 214105 (2009);
A. Mitsutake, J. Chem. Phys. 131, 094105 (2009).

2. Isobaric-Isothermal Ensemble（定圧定温アンサンブル）＊ MUCA: Multibaric-Multithermal Algorithm (MUBATH)
random walk in energy E and volume V
H. Okumura & Y.O., Chem. Phys. Lett. 383, 391 (2004). (MC version)
H. Okumura & Y.O., Chem. Phys. Lett. 391, 248 (2004). (MD version)＊ REM: random walk in temperature T and pressure P
Y. Sugita & Y.O., Lect. Notes in Computational Science & Engineering,ed. by T. Schlick and H.Gun (2002) pp. 304-332.
T. Okabe, M. Kawata, Y.O. & M. Mikami, Chem. Phys. Lett. 335, 435 (2001).

＊ ST: random walk in temperature T and pressure P
Y. Mori & Y.O., J. Phys. Soc. Jpn. 79, 074003 (2010).

Examples of Multidimensional REM, MUCA, and ST
A. Mitsutake & Y.O., Phys. Rev. E 79, 047701 (2009);
J. Chem. Phys. 130, 214105 (2009);
A. Mitsutake, J. Chem. Phys. 131, 094105 (2009).

3. Umbrella Sampling

LH k q, p  K p  E0 q    Vl q, where Vl   Kl  q  dl  .
l2l 1＊ REM: Replica-Exchange Umbrella Sampling (REUS)
random walk in reaction coordinate Y. Sugita, A. Kitao & Y.O., J. Chem. Phys. 113, 6042 (2000).

＊ ST: Simulated Tempering Umbrella Sampling (STUS)
random walk in reaction coordinate Y. Mori & Y.O., Physical Review E 87, 023301 (2013).

Cf.
＊ MUCA: Metadynamics
random walk in reaction coordinate A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. USA 99, 12562 (2002).

拡張アンサンブル法によるタンパク質とリガンドの結合構造予測
Prediction of protein-ligand binding structures bygeneralized-ensemble simulations
H. Kokubo, T. Tanaka & Y.O., J. Comput. Chem. 32, 2810-2821 (2011).

Replica-Exchange Umbrella Sampling (REUS)Y. Sugita, A. Kitao & Y.O., J. Chem. Phys. 113, 6042 (2000).
束縛4束縛3束縛2束縛11234
Enhanced sampling about thedistance between the protein
pocket and ligand

By using the reweighting techniques, the ensemble average properties
without the constraint can be recovered.

Five test systems [ligand (protein)] (T = 300 K, P = 1 atm)
H. Kokubo, T. Tanaka & Y.O., J. Comput. Chem. 32, 2810-2821 (2011).
benzodiazepine(protein: MDM2)

deoxythymidine 3’,5’-bisphosphate (pdTp)
tyrosine(protein: staphylococcal nuclease)(protein: aldolase)

PDB code  PDB code
cytidylic acid (2’-CMP)(protein: ribonuclease A)

2-aminopyrimidine(protein: heat shock protein HSP90)

heat shock protein and 2aminopyrimidine
Asp-78 and ligand

Simulation and movie
by H. Kokubo

Results•

H. Kokubo, T. Tanaka & Y.O., J. Comput. Chem. 32, 2810-2821 (2011).
Starting from the configurations in which the protein and the ligand are faraway from each other in each system, our method predicted the ligand binding
structures as the global minima in free energy (or, potential of mean force) inexcellent agreement with the experimental data from PDB.

crystal predicted
potential of mean force shows
the most stable distance

We remark that for 1ROB and 1SNC, there areattempts by a popular existing docking program,
GOLD, but they failed in the predictions(classified as significant errors).

ATPの加水分解
生体内におけるエネルギー源
１．ATP＋水分子の拡張アンサンブルDFT-MDシミュレーション
２．ATP＋ATP加水分解酵素（GroEL、F1モーターなど）＋水分子の拡張アンサンブルDFT-MDシミュレーション

Y. Mori & Y.O., Physical Review E 87, 023301 (2013).

焼き戻し傘サンプル法
Simulated Tempering Umbrella Sampling (STUS)
trialUtrial

trialUreject

U'trialaccept
自由エネルギー
free energy

アンブレラポテンシャル Umbrella Potential

反応座標
reaction coordinate
遷移確率 Transition Probability

反応座標
Reaction coordinate
無次元化された自由エネルギー Dimensionless Free Energy

時間 timeマロンアルデヒドの分子内プロトン移動
Proton Transfer in Malonaldehyde
Y. Mori & Y.O., Physical Review E 87, 023301 (2013).
反応座標の定義
Definition of reaction
coordinate

マロンアルデヒド malonaldehyde： 1 分子 molecule
水 water： 71 分子 molecules
総原子数 Total No. of Atoms： 222

シミュレーション条件
Simulation Conditions
Y. Mori & Y.O., Physical Review E 87, 023301 (2013).
計算条件
プログラム Program

CP2K (ver. 2.1) *交換相関汎関数 Functional

BLYP

基底 Bases

Gaussian and plane wave

擬ポテンシャル Pseudo Potential

Goedecker-Teter-Hutter

アンサンブル Statistical Ensemble

NVT

温度制御法

Thermostat

Nosé-Hoover chains (# of chains = 3)

温度

Temperature

300 K

体積

Volume13.82×13.82×13.82 (Å3)

時間ステップ Time Step

0.5 fs

総シミュレーション時間 Total Sim. Time200 ps (50 ps × 4 runs)

Umbrella Potentials
アンブレラポテンシャル（k）

0.01 (hartree･bohr-2)

アンブレラポテンシャル（ξ0）-1.0, -0.8, -0.6, -0.4, -0.2, 0.0,0.2, 0.4, 0.6, 0.8, 1.0 (Å)

11個のパラメータを使用* http://cp2k.berlios.de/シミュレーション結果
Y. Mori & Y.O., Physical Review E 87, 023301 (2013).

Simulation and movie by Y. Mori

反応自由エネルギー Potential of Mean ForceY. Mori & Y.O., Physical Review E 87, 023301 (2013).

平均力ポテンシャル
Potential of Mean Force: PMF

平均力ポテンシャルを再重法により計算
Multistate Bennett Acceptance Ratio estimator:
MBAR** M. R. Shirts and J. D. Chodera, J. Chem. Phys. 129 (2008) 124105.

反応自由エネルギー Potential of Mean ForceY. Mori & Y.O., Physical Review E 87, 023301 (2013).

赤 red：STUS
緑 green：従来の US

まとめ
拡張アンサンブル法（特に、マルチカノニカル法、焼き戻し法、レプリカ交換法およびそれらの多次元／多変数への拡張版）は古典系ばかりでなく量子系における分子シミュレーションにおいてとても有効であることが示された。

COLLABORATORS
Ayori MITSUTAKE [Keio Univ.]（光武亜代理）
Yoshiharu MORI [Nagoya Univ. → IMS]（森義治）
Shingo ITO [Nagoya Univ.]（伊東真吾）
Stephan IRLE [Nagoya Univ.]
Yoshitake SAKAE [Nagoya Univ. & IMS]（榮慶丈）
Tetsuro NAGAI [Nagoya Univ. → Ritsumeikan Univ.]（永井哲郎）
Yuji SUGITA [RIKEN]（杉田有治）
Hironori KOKUBO [Takeda Pharm.] （小久保裕功）
Toshimasa TANAKA [Takeda Pharm.] （田中稔祐）