名古屋大学 大学院理学研究科 物理学教室 最終講義
再収録 November 18, 2022

計算機シミュレーションによる研究を振り返って（その４）岡本 祐幸（おかもと ゆうこう）
名古屋大学
情報基盤センター 大規模計算支援環境研究部門
国際本部 グローバル・エンゲージメントセンター
e-mail: okamoto{a}cc.nagoya-u.ac.jp
URL: https://yuko-okamoto.github.io/homepage/index.shtml

研究内容の話に戻ります。
研究内容については、今年（2022年）３月１４日（月）に日本生物物理学会のサブグループ「生体分子シミュレーション・モデリング」の第２回研究会で、「拡張アンサンブル法による生体分子シミュレーションによる研究を振り返って」というタイトルで講演させて頂き、その動画が既に以下のYoutubeサイト
招待講演 岡本祐幸先生（2022/3/14）
https://www.youtube.com/watch?v=qmPxJG_gT10にアップロードされています。
今回、もし時間切れとなり、研究内容を十分お話しできませんでしたら、このYoutube動画をご覧頂けると幸いです。

Met-enkephalin in aqueous solution (RISM) (Level 3)
In gas phase, the peptide is compact due tointra-chain hydorgen bonds (b-hairpin structures).
In aqueous solution, however, the peptide is extended:
In agreement with NMR experiments.

8 low-energy conformations from
Simulated Annealing (SA) runs
are superposed.

M. Kinoshita, Y. O. & F. Hirata, J. Am. Chem. Soc. 120, 1855 (1998).

Met-enkephalin in aqueous solution (RISM) (Level 3)•

Multicanonical Simulation

Simulation and movie by A. Mitsutake

In aqueous solution, the peptide is extended:
In agreement with NMR experiments.
A. Mitsutake, M. Kinoshita, Y. O. & F. Hirata,
Chem. Phys. Lett. 329, 295 (2000).
See also: A. Mitsutake, M. Kinoshita, Y. O. & F. Hirata,
J. Phys. Chem. B 108, 19002 (2004).

Average End-to-End Distance of Met-enkephain
in Gas Phase and in Aquesou Solution (RISM): ECEPP/2

A. Mitsutake, M. Kinoshita, Y. O. & F. Hirata,
Chem. Phys. Lett. 329, 295 (2000).
See also: A. Mitsutake, M. Kinoshita, Y. O. & F. Hirata,
J. Phys. Chem. B 108, 19002 (2004).

Hpoly-1

Helix Propensities of Homopolymers: (Ala)10, (Val)10 & (Gly)10
Average Helicity (Level 2)

Level 0: Y.O. & U. Hansmann, J. Phys. Chem. 99, 11276 (1995).
Level 2: A. Mitsutake & Y.O., Chem. Phys. Lett. 309, 95 (1999).
A. Mitsutake & Y.O., J. Chem. Phys. 112, 10638 (2000).

Hpoly-2

Average Energy of (Ala)10
Level 0 (Gas Phase)

Level 2

A. Mitsutake & Y.O., Chem. Phys. Lett. 309, 95 (1999).

Hpoly-3

Specific Heat of (Ala)10

Helix-Coil Transition Temperatures:
T = 420 K (Level 0)
T = 340 K (Level 2)

A. Mitsutake & Y.O., Chem. Phys. Lett. 309, 95 (1999).

Importance of Solvent Effects• Level 1 Distance-Dependent Dielectric Function e(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

We compared the lowest-energy conformations obtained from simulations
with the native conformations (more accurate solvation theory gives
better agreement).

C-Peptide (N = 13): X-Ray

C-Peptide of Ribonuclease A
Amino-Acid Sequence:
Lys-Glu-Thr-Ala-Ala-Ala-Lys-Phe-Glu-ArgGln-His-Met

C-Peptide (N = 13): X-Ray

Level 0 (Gas Phase)

RMSD = 1.9 A

H. Kawai, Y.O., M. Fukugita, T. Nakazawa & T. Kikuchi, Chem. Lett. 1991, 213.
Y.O., M. Fukugita, T. Nakazawa & H. Kawai, Protein Eng. 4, 639 (1991).

C-Peptide (N = 13): X-Ray

Level 0 (Gas Phase)

RMSD = 1.9 A

Level 1RMSD = 1.4 A

U. Hansmann & Y.O., J. Phys. Chem. B 102, 653 (1998).
U. Hansmann & Y.O., J. Phys. Chem. B 103, 1595 (1999).

C-Peptide (N = 13): X-Ray

Level 0 (Gas Phase)

RMSD = 1.9 A

Level 1RMSD = 1.4 A

Level 2

RMSD = 0.8 A
M. Masuya & Y.O., unpublished.

BPTI(16-36) (N = 21): X-Ray

Fragment of Bovine Pancreatic Trypsin Inhibitor
BPTI(16-36)
Amino-Acid Sequence:
Ala-Arg-Ile-Ile-Arg-Tyr-Phe-Tyr-Asn-AlaLys-Ala-Gly-Leu-Cys-Gln-Thr-Phe-Val-TyrGly

BPTI(16-36) (N = 21): X-Ray

Level 0 (Gas Phase)

T. Nakazawa, H. Kawai, Y.O. & M. Fukugita, Protein Eng. 5, 495 (1992).

BPTI(16-36) (N = 21): X-Ray

Level 0 (Gas Phase)

Level 1

T. Nakazawa & Y.O., J. Peptide Res. 54, 230 (1999).

BPTI(16-36) (N = 21): X-Ray

Level 1

Level 0 (Gas Phase)

Level 2

Y.O., M. Masuya, M. Nabeshima & T. Nakazawa, Chem. Phys. Lett. 299, 17 (1999).

C-Peptide of Ribonuclease A (N=13) in Water (Level 3)
Amino-Acid Sequence:
ACE-Ala-Glu-Thr-Ala-Ala-Ala-Lys-Phe-Leu-Arg-Ala-His-Ala-NME
AMBER 99 Force Field (peptide)
TIP3P (water)

13 residue peptide
1387 Water Molecules(R = 22 Å)
Total 4365 atoms

NMR Studies (Osterhout 1988)
50-60% helix in solution
First three residues are
distorted
Formation of a salt bridge
between Glu2 and Arg10

Initial Configuration

C-Peptide of Ribonuclease A (N=13) in Water
Y. Sugita & Y. O., Biophys. J. 88, 3180 (2005).
1. REM (32 replicas)
100 ps X 32

2. MUCAREM (4 replicas)
1 ns X 4# of tunneling events = 0(-12850, -8250)# of tunneling events = 5(-12850, -8250)

C-Peptide of Ribonuclease A (N=13) in Water
Y. Sugita & Y. O., Biophys. J. 88, 3180 (2005).
3. REMUCA (1 replica)
15 ns# of tunneling events = 46(-12850, -8250)

Time series ofpotential energy

C-Peptide of Ribonuclease A (N=13) in Water

Canonical (T = 275 K)

Simulation and movie by Y. Sugita

REMUCA

C-Peptide of Ribonuclease A (N=13) in Water
REMUCA Simulation with AMBER (parm99)

Simulation and movie by Y. Sugita

C-Peptide of Ribonuclease A (N=13) in Water
Y. Sugita & Y. O., Biophys. J. 88, 3180 (2005).

RMSD from “native structure”
Time Series of RMSD

Rg: radius of gyration

Time Series of Rg

ヘルムホルツ自由エネルギー

F (T ) = - kBT log Z =< E >T - TS (T )1
1 - b E ( x)Peq ( x) = WB ( x; T ) = e
ZZ
where Z =1
dx WB ( x; T ) and b =k BTFree Energy Landscape (C-peptide)
T = 300 K

At T = 300 K, 3 characteristic structural groups exist, among
which the group of native-like structure with Glu2-Arg10 salt
bridge is the global minimum in free energy.
Y. Sugita & Y. O., Biophys. J. 88, 3180 (2005).

G-Peptide (N=16) of Protein G in Water
Amino-Acid Sequence:
Gly-Glu-Trp-Thr-Tyr-Asp-Asp-Ala-Thr-Lys-Thr-Phe-Thr-Val-Thr-Glu

16 residue peptide
2328 TIP3P water molecules(R = 26 Å)
Total: 7231 atoms• G-peptide is a C-terminal fragment(41-56 residues) of protein G B1
domain.
• This b-hairpin structure is stable inaqueous solution (F. Blanco, G. Rivas& L. Serrano, Nature Struct. Biol. 1,584 (1994); N. Kobayashi, S. Honda, H.
Yoshii, H. Uedaira & E. Munekata,
FEBS Lett. 366, 99 (1995)).

“Native Structure” of G-Peptide

G-Peptide of Protein G (N=16) in Water
MUCAREM Simulation with CHARMM (ver. 22)

Simulation and movie by T. Yoda

Snap-shot Structures of G-peptide
observed in the Simulation

T. Yoda, Y. Sugita & Y.O., Proteins 66, 846 (2007).

Comparisons of force field parameters for
all-atom protein models with explicit solvent
T. Yoda, Y. Sugita & Y.O., Chem. Phys. Lett. 386, 460-467 (2004);
T. Yoda, Y. Sugita & Y.O., Chem. Phys. 307, 269-283 (2004).

Force-field parameters for all-atom models
Commonly used conformational potential energy

Econf = å K r (r - req ) 2bonds+ å Kq (q - q eq ) 2
anglesVn[1 + cos(nf - g )]dihedrals 2

Bond-stretching
Bond-bending+ å

Dihedral angle

é Aij Bij
qi q j ù+ å ê 12 - 6 + 332
úe Rij úûi< j ê
ë Rij Rij

Non-bonding interactions(Lennard-Jones and electrostatic)

These energy terms include some force-field parameters (blue color)

The existing force fields have different force-field parameters

Typical example of folding simulations using different force fields
C-peptide(13 residues)
Lys-Glu-Thr-Ala-Ala-Ala-Lys-Phe-Glu-Arg-Gln-His-Met

AMBER ff94

AMBER ﬀ96

Simulation and movie by Y. Sakae
Method：Simulated annealing, Simulation time : 1.0 nsec, Temperature : 700～200 K, Solvent model : GB/SA

α-helix

Force-field dependency of secondary-structure properties

C-peptide(13 residues)

Helicity of C-peptide
T. Yoda, Y. Sugita & Y.O.,
Chem. Phys. Lett. 386, 460 (2004).

AMBER parm94
AMBER parm96
AMBER parm99
CHARMM22
OPLS-AA/L
GROMOS96

7. 分子科学研究所 (1995-2005)
2001年8月、米国カリフォルニア州San Diegoの隣のLa JollaのScripps研究所Charlie Brooksを杉田有治さんと訪問しました。
I gave a seminar talk.
In one of the slides I had the unpublishedfigure on the right. Namely, π-helix
formation was found in C-peptide by
MUCAREM simulations with CHARMM22.

T. Yoda, Y. Sugita & Y.O., Chem. Phys. Lett. 386, 460 (2004);
T. Yoda, Y. Sugita & Y.O., Chem. Phys. 307, 269 (2004).

I think this influenced the development of CMAP of CHARMM force field:
M. Feig, A.D. MacKerell Jr. & C.L. Brooks III, J. Phys. Chem. B 107, 2831 (2003).
A.D. MacKerell Jr., M. Feig & C.L. Brooks III, J. Am. Chem. Soc. 126, 698 (2004).

Force field refinement for all-atom protein models
Y. Sakae & Y.O., Chem. Phys. Lett. 382, 626-636 (2003);
Y. Sakae & Y.O., J. Phys. Soc. Jpn. 75, 054802 (9 pages) (2006);
Y. Sakae & Y.O., Mol. Sim. 36, 138-158 (2010);
Y. Sakae & Y.O., J. Chem. Phys. 138, 064103 (2013), etc.
Review:
Y. Sakae & Y.O., in Computational Methods to Study the Structure and Dynamics ofBiomolecules and Biomolecular Processes – from Bioinformatics to Molecular
Quantum Mechanics, A. Liwo (ed.) (Springer-Verlag, Berlin Heidelberg, 2019) pp.
203-256.

Conformational Energy

Econf = EBL + EBA + Etorsion + Enonbond
Backbone-torsion energy termVnEtorsion = åå [1 + cos(nF - g n )]F n 2
= E (j ,y ) + Erest
backbone dihedral
angles φ and ψ

Φ
n
γn
Vn

: all dihedral angles: number of waves: phase: Fourier coefficient

rest of the torsion terms

Side chain

316

Backbone-torsion energy surfaces of some force fields

180

180

90

90

0

-180
-90

y

-90
0
90

f

-90

y

-90
0
90

f

AMBER parm94

90

0

-180

-180
180

0
f

90

-180
180

CHARMM22

y

y

-90
-180

90

f

AMBER parm96

180

AMBER parm99

180

180

90

0
-90

-90

-180
180

90

-90

0

-180
0

180

-180

180

0

-180
-90

-90
0

f

90

-180
180

OPLS-AA

y

90
0

-180
-90

-90
0

f

90

-180
180

OPLS-AA/Ly1. Proposal of new functional form
Y. Sakae and Y.O., J. Phys. Soc. Jpn. 75, 054802 (9 pages) (2006).

conventional energy term

Vm
Vn
E (j ,y ) = å [1 + cos(mj - g m )] + å [1 + cos(ny - g n )]m 2
n 2

New backbone-torsion-energy term
¥

¥

m =1

n =1

e (j ,y ) = a + å (bm cos mj + cm sin mj ) + å (d n cos ny + en sin ny )¥¥+ åå ( f mn cos mj cos ny + g mn cos mj sin ny
m =1 n =1+ hmn sin nj cos ny + imn sin mj sin ny )a,bm,cm,dn,en,fmn,gmn,hmn,imn: Fourier coefficient

New torsion energy term

Etorsion = e (j ,y ) + Erest318

We consider the truncated Fourier series:e (j ,y ) = a + b1 cos j + c1 sin j + b2 cos 2j + c2 sin 2j+ d1 cosy + e1 siny + d 2 cos 2y + e2 sin 2y+ f11 cos j cosy + g11 cos j siny
+ h11 sin j cosy + i11 sin j sinyThis function has 13 Fourier-coefficient parameters.

319

Example of the application of new backbone-torsion energy term
Ramachandran plot
Carl Branden and John Tooze著（勝部幸輝ほか訳）「タンパク質の構造入門第2版」
ニュートンプレス（2000）

180
90

-90

-90

0
f

90

-180
180

α-helix region

90

90

0

-180

180

180

y

0

-180
-90

-90

0
f

90

-180
180

0

-180
y

-90

-90

0
f

90

-180
180

β-structure region

Energy surface of new energy term can represent Ramachandran space directlyyApplication to AMBER parm94 and AMBER parm96
AMBER parm94

180
90
0

-180
-90

-90
0
90

f

y

-180
180

~

AMBER parm96~

α-helix~~

β-structure

180
90
0

-180
-90

-90
0

f

90

y

-180
180

~

~

α-helix~~

321
β-structureC-peptide

Results of folding simulations
AMBER parm94

AMBER parm96

Original

Original

α-helix region

α-helix region

β-structure region

β-structure region

Simulated annealing simulation

Simulation time : 1ns (1,000,000 MD steps × 1.0fs × 60 times)
Temperature : 2,000K to 250K (Berendsen’s method) 322
Solvent : GB/SA modelG-peptide

Results of folding simulations
AMBER parm94

AMBER parm96

Original

Original

α-helix region

α-helix region

β-structure region

β-structure region

Simulated annealing simulation

Simulation time : 1ns (1,000,000 MD steps × 1.0fs × 60 times)
Temperature : 2,000K to 250K (Berendsen’s method) 323
Solvent : GB/SA modelOptimization method of force-field parameters
{ m}
!
¶Etot
f im = - !
¶xim

Nm

Number of atoms in molecule mf imForce acting on atom i

{m}
Etot

Total potential energy for molecule m

Atom if im

Molecule m

If force-field parameters are of ideal values, all native structuresare stable without any force acting on each atom in molecules.
(we expect F=0)
In reality, F≠0, and because F ≥0, we can optimize the force-fieldparameters by minimizing F with respect to these parameters.
Y. Sakae and Y.O., Chem. Phys. Lett. 382, 626-636 (2003).

Structures of selected 100 proteins

Resolution: >= 2.0Å, Sequence similarity of aminoacid: >= 30%, Number of residue: < 200

Results: Optimized force-field parameters and its energy surface
Optimized ff
180
90
0

-180
-90
0
90

@

@

-90

e (j ,y ) = 0.835 cos j - 0.327 sin j- 0.088 cos 2j + 0.100 sin 2j+ 0.287 cosy - 0.160 siny
+ 0.019 cos 2y - 0.054 sin 2y- 0.427 cos j cosy + 0.247 cos j siny
+ 0.114 sin j cosy + 0.603 sin j siny

-180
180

AMBER parm96

AMBER parm94

180

180
90

90
0

-180
-90

-90
0

f

90

-180
180

y

0

-180
-90

-90
0

f

90

-180
180

y

Test simulations

C-peptide(13 residues)

G-peptide(16 residues)

blocked by COCH3- and –NH2
Replica-Exchange Molecular Dynamics (REMD) simulations
Simulation time : 5.0 ns (32 replica)
Temperature : 700K to 250K (Nosé-Hoover method)
Solvent : GB/SA model
Program : Modified TINKER program package
327Comparison of helicity and strandness
C-peptide

300 K

Comparison of helicity and strandness
G-peptide

300 K

Potential of mean force

DG (f ,y ) = -k BT ln PB (f ,y )at 300 K

Optimized ff

C-peptide
AMBER ff94

AMBER ff96

Potential of mean force

DG (f ,y ) = -k BT ln PB (f ,y )at 300 K

Optimized ff

G-peptide
AMBER ff94

AMBER ff96

膜タンパク質の立体構造予測
Prediction of Membrane Protein Structures byReplica-Exchange Monte Carlo Simulations
H. Kokubo & Y.O., Chem. Phys. Lett. 383, 397-402 (2004).
H. Kokubo & Y.O., J. Chem. Phys. 120, 10837-10847 (2004).
H. Kokubo & Y.O., J. Phys. Soc. Jpn. 73, 2571-2585 (2004).
H. Kokubo & Y.O., Chem. Phys. Lett. 392, 168-175 (2004).
H. Kokubo & Y.O., Biophys. J. 96, 765-776 (2009).
R. Urano, H. Kokubo & Y.O., J. Phys. Soc. Jpn. 84, 084802 (2015).
R. Urano & Y.O., J. Chem. Phys. 143, 235101 (2015).

8/ 201
31 1
/2
01
2

10

20

09

08

20

20

07

20

06

20

05

20

04

03

20

194,259 entries as ofAugust 17, 2022, but only
about 2 % of them aremembrane protein structures .

All PDBAll TM

20

Motivation:
Protein Data Bank (PDB) :

90000
80000
70000
60000
50000
40000
30000
20000
10000
0

It is estimated that 20-30 % of all genes in most genomes encode
membrane proteins. However, only a small number of detailed
structures have been obtained for membrane proteins because oftechnical difficulties such as high quality crystal growth.
No. of Membrane Proteins■ Eubacteria◆ Archaebacteria▲ Eucaryotes

No. of ORF

S. Mitaku, Biophysics 42, (2002) 104-109 (in Japanese).
A. Krogh et al., J. Mol. Biol. 305 (2001) 567-580.

We can know Membrane
Spanning Regions
by Prediction Tool on the Web:SOSUI,HMMTOP,TMHMM..etc.

Predict transmembrane
helix configurations bymolecular simulations

Sequence
MNGTEGPNFYVP、、、、

No. N terminaltransmembrane region

C terminal

type

length1

10

IWLALGTALMGLGTLYFLVKGMG

32

SECONDARY

23

2

45

TTLVPAIAFTMYLSMLLGYGLTM

67

SECONDARY

23

3

77

IYWARYADWLFTTPLLLLDLALL99

SECONDARY

23

4

106

TILALVGADGIMIGTGLVGALTK128PRIMARY

23

5

136

WWAISTAAMLYILYVLFFGFTSK158

PRIMARY

23

6

176

VTVVLWSAYPVVWLIGSEGAGIV198

PRIMARY

23

7

203

ETLLFMVLDVSAKVGFGLILLR224

PRIMARY

22

Tertiary StructureTwo-stage model
Individual helices of a protein are postulated to bestable separately as domains in a lipid bilayer.
And then side-to-side helix association is driven,
resulting in a functional protein.
J.-L. Popot and D.M. Engelman, Annu. Rev. Biochem. 69, 881 (2000).

lIndividual helix stability as domains is a consequence of the hydrophobic effect andmain-chain hydrogen bonding.
lSpecific folding energy is provided mainly by packing of the preformed helices witheach other, by loop structures, and by interactions with prosthetic groups.
lAdditionally, ion pairs and hydrogen bonds between helices are sometimes found, andgeneral contributions are made by interactions with the lipid environment.

Simulation Conditions
Harmonic constraint are set only
when helices are apart at somedistance from membrane boundary
planes.

We keep backbone structures rigid
in order to perform the simulation
with helix structures.

We neglect loop region aminoacids, lipids, and water.

Harmonic constraint are set only
when helices are apart at somedistance.

We construct ideal a-helices
from amino-acid sequence.

CHARMM force-field (ver. 19)
Replica-exchange Monte Carlo simulation(rigid translation, rigid rotation of each helix andtorsion angle rotation of side chains)DIMERIC TRANSMEMBRANE DOMAIN OF HUMAN
GLYCOPHORIN A
Simulation and movie by H. Kokubo

Principal Component Analysis:
case for ε=1.0①
cluster1②
cluster2③
cluster3

200 K

342 K

585 K

1888 K

Structure of Each Cluster:case for ε=1.0

Viewed from different angles

Native(NMR)

Cluster1

Cluster2

Cluster3

H. Kokubo & Y.O., Chem. Phys. Lett. 392, 168-175 (2004).
H. Kokubo & Y.O., Biophys. J. 96, 765-776 (2009).

Bacteriorhodopsin （Case of 7 Helices）
Native Structure

Initial structures of the simulation

From the top
From the side

Replica1

Replica2
Initial structures

Replica14

32 random initial configurations obtainedby simulations at the highest temperature
The replica-exchange MC simulation of 32 replicas

MC time seriesT

Potential Energy

RMSD

MC time series

Replica６Replica１４Replica１９Replica２５

Snapshots (Replica 16)

REMC(from above)

Replica 14

⑤
④
③⑥
②⑦
①

Native structure

Simulation and movie by H. Kokubo

REMC(from above)

Replica 14

⑤⑥⑦
④③②①

Native structure
Simulation and movie by H. Kokubo

H. Kokubo & Y.O., Biophys. J. 96, 765-776 (2009).

Local Minimum Structure

RMSD 4.8 Å(all Cα)

Native Structure

Replica 14

H. Kokubo & Y.O., Biophys. J. 96, 765-776 (2009).

Local Minimum Structure

RMSD 4.8 Å(all Cα)

Native Structure

Replica 14

Principal
Component
Analysis at(a) T = 500 K(b) T = 976 K(c) T = 5000 K

Typical
Structuresof the 7
Clusters
Top view

Side view

H. Kokubo & Y.O.,
Biophys. J. 96, 765-776 (2009).

Bacteriorhodopsin （Case of 7 Helices）
R. Urano & Y.O., J. Chem. Phys. 143, 235101 (2015).

Introduced bending of TM helices⑥⑪

Native Structure⑩⑮

Simulation and movie by R. UranoNative Structure (left) andGlobal-Minimum Free-Energy
Structure (right) from REM Simulations

R. Urano & Y.O., J. Chem. Phys. 143, 235101 (2015).

Bending of TM helices reproduced.
R. Urano & Y.O., J. Chem. Phys. 143, 235101 (2015).

RMSD (Å)

0.5 0.8

1.0

1.1

0.9

1.8

0.7

http://spider.art.coocan.jp/biology2/cell2012.html

T. Nagai, R. Ueoka & Y.O., J. Phys. Soc. Jpn. 81, 024002 (2012).

Phase behavior of DPPC bilayer
studied by REMD with a coarse-grained
model
MARTINI 2.0
a coarse-grainedmodel.
four atoms are treatedas one site, so thatlarger system size
and/or time scale
can be studied.

lipid

phosphategrycerolwater
cholinepalmitoyl(C1-tail)palmitoyl (C2-tail)Simulation of a 32-lipid system• system– 32 DPPC and 500 water particles ( including 50 anti-freeze particles)• temperature– Thermostat: Nose-Hoover method in the production run– Temperature distribution: 127 points between283 K and 390 K• pressure– Reference pressure: １atom– Barostat: Parrinello-Raman method in the production run• Frequency of replica exchange– Every 100 step

• Time step– 20 fs

Results
Animation
A trajectory of one of
the replicas.
Water is suppressed
for clarity.

Simulation and movie by T. Nagai

Results
T. Nagai, R. Ueoka & Y.O., J. Phys. Soc. Jpn. 81, 024002 (2012).
Change of Potential of Mean Force (PMF)

Simulation and movie by T. Nagaib f (x ,h )= - ln ò dqd (x - x (q))d (h - h (q))e - b E ( q )Results
T. Nagai, R. Ueoka & Y.O., J. Phys. Soc. Jpn. 81, 024002 (2012).

Potential of Mean ForceThree states are found.
untilted gelsoltilted gel
295.76 K

T. Nagai, R. Ueoka & Y.O., J. Phys. Soc. Jpn. 81, 024002 (2012).

Results
Case for a bigger system: 128 lipids

sol

gel
(un-tilt)

gel(tilt)

Structural Characteristics of Monomeric Aβ42 on
Fibril in the Early Stage of Secondary Nucleation
Mechanism
K. Noda, Y. Tachi, & Y.O.,ACS Chemical Neuroscience 11, 2989-2998 (2020).

Path for amyloid β fibril formation
A

C

Primary nucleationB

362/ 25

Elongation

Secondary nucleation[1] M. Törnquist et al., Chem. Commun. 54, 8667 (2018).
[2] S. I. A. Cohen et al., PNAS 110, 9758 (2013).

3. 方法REST2 (Replica Exchange with Solute Tempering 2)363/ 25

REST2:[1] P. Liu et al., PNAS 102, 13749 (2005).
[2] T. Terakawa, T. Kameda, and S. Takada, J. Comput. Chem. 32, 1228 (2011).
[3] L. Wang et al., J. Phys. Chem. B 115, 9431 (2011).
[4] S. Jo and W. Jiang, Comp. Phys. Commun. 197, 304 (2015).

Simulation 1｜Aβ42 monomer in solution

364/ 25[1] W. L. Jorgensen et al., J. Chem. Phys. 79, 926 (1983).

365/ 25

Simulation 2｜Aβ42 fibril and monomer in solution[1] W. L. Jorgensen et al., J. Chem. Phys. 79, 926 (1983).
[2] Y. Xiao et al., Nat. Struct. Mol. Biol. 22, 499 (2015).

366/ 25

Test of simulations
J coupling constant
カップリング定数（J-coupling）

RMSD from experimentThis study

0.66 Hz

Sugourakis et al. [2]0.73 HzSong et al. [3]0.82 Hz

J-coupling is in good agreement[1] J. Roche, et al., Biochemistry 55, 762 (2016).
[2] N. G. Sugourakis, et al., J. Mol. Bio. 405, 570 (2011).
[3] W. Song, et al., Sci. Rep. 5, 11024 (2015).

Radius of gyration and number of contacts

367/ 25#$%
Radius of𝑟!"Gyration
In solution ＜ on fibril

Intracontacts
残基A ＞残基BIn solution
on fibril

Unfolded on fibril

368/ 25

5. 結果と考察

線維上でのヘリックス構造

酸性
塩基性

C-terminus
Lys16

Glu22

N-terminus–+

Helix formation

Monomer on fibril

Helix formation：in solution ＜ on fibril

369/ 25

Development of New Generalized-Ensemble Algorithms
REVIEWS: A. Mitsutake, Y. Sugita, & Y.O., Biopolymers 60, 96 (2001);
Y.O., J. Mol. Graphics Modell. 22, 425 (2004);
Y. Sugita, A. Mitsutake, & Y.O., in Lecture Notes in Physics,
W. Janke (ed.) (Springer-Verlag, 2008) pp. 369-407.
Combination of Replica-Exchange Method (REM) andMulticanonical Algorithm (MUCA)
1. Replica-Exchange Multicanonical Algorithm (REMUCA)
Y. Sugita & Y.O., Chem. Phys. Lett. 329, 261 (2000).
A. Mitsutake, Y. Sugita, & Y.O., J. Chem. Phys. 118, 6664; 6676 (2003).

• Multicanonical weight factor is determined from a short REM
simulation by multiple histogram techniques
2. Multicanonical Replica-Exchange Method (MUCAREM)
Y. Sugita & Y.O., Chem. Phys. Lett. 329, 261 (2000).
A. Mitsutake, Y. Sugita, & Y.O., J. Chem. Phys. 118, 6664; 6676 (2003).

• Multicanonical simulations are performed for each replica and apair of replicas are exchanged

Development of New Generalized-Ensemble Algorithms
REVIEWS: A. Mitsutake, Y. Sugita, & Y.O., Biopolymers 60, 96 (2001);
Y.O., J. Mol. Graphics Modell. 22, 425 (2004);
Y. Sugita, A. Mitsutake, & Y.O., in Lecture Notes in Physics,
W. Janke (ed.) (Springer-Verlag, 2008) pp. 369-407.
Combination of Replica-Exchange Method (REM) and
Simulated Tempering (ST)
3. Replica-Exchange Simulated Tempering (REST)
A. Mitsutake & Y.O., Chem. Phys. Lett. 332, 131 (2000).

• Simulated tempering weight factor is determined from ashort REM simulation by multiple histogram techniques
4. Simulated Tempering Replica-Exchange Method (STREM)
A. Mitsutake & Y.O., J. Chem. Phys. 121, 2491 (2004).

• Simulated tempering simulations are performed for each replica
and a pair of replicas are exchanged
See also:
M.K. Fenwick & F.A. Escobedo, J. Chem. Phys. 119, 11998 (2003).

KEY ELEMENT: Multiple-Histogram
Reweighting Techniques (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).

Given M set of histograms Nm(E), the following WHAM equations are
solved iteratively for density of states n(E) and dimensionless Helmholtz
free energy f m :Mn (E ) =å N (E )

m =1
M

m

å nme

m =1

fm - bm E, where e- fm= å n ( E )eE

- bm E.

Estimating Density of States of Complex
Systems by REWL-MUCAREM Simulations* T. Hayashi and Y.O.,“An Efficient Simulation Protocol for Determining theDensity of States: Combination of Replica-Exchange
Wang-Landau Method and Multicanonical ReplicaExchange Method,”Phys. Rev. E 100, 043304 (2019).
* T. Hayashi, C. Muguruma, and Y.O.,“Calculation of the Residual Entropy of Ice Ih by Monte
Carlo Simulation with the Combination of the ReplicaExchange Wang-Landau Algorithm and MulticanonicalReplica-Exchange Method,”
J. Chem. Phys. 154, 044503 (11 pages) (2021).

MUCAREM and REWL
* MUCAREM (Multicanonical Replica-Exchange Method)
Y. Sugita and Y.O., Chem. Phys. Lett. 329, 261 (2000).
A. Mitsutake, Y. Sugita, and Y.O.,
J. Chem. Phys. 118, 6664 (2003).
A. Mitsutake, Y. Sugita, and Y.O.,
J. Chem. Phys. 118, 6676 (2003).
In MUCAREM, we use the multiple-histogram
reweighting techniques extensively.
* REWL (Replica-Exchange Wang-Landau)
T. Vogel, Y. W. Li, T. W¨ust, and D. P. Landau,
Phys. Rev. Lett. 110, 210603 (2013).
T. Vogel, Y. W. Li, T. W¨ust, and D. P. Landau,
Phys. Rev. E 90, 023302 (2014).

REWL-MUCAREM protocol:
T. Hayashi and Y.O., Phys. Rev. E 100, 043304 (2019).

We first perform REWL, then MUCAREM.
We tested the performances of MUCAREM only, REWL only,
and REWL-MUCAREM with the 2-dimensional Ising Model.

REWL-MUCAREM protocol:
T. Hayashi and Y.O., Phys. Rev. E 100, 043304 (2019).

Exact Solutions:P.D. Beale, Phys. Rev. Lett. 76, 78 (1996);
A.E. Ferdinand and M.E. Fisher, Phys. Rev. 185, 832 (1969).

Calculation of residual entropy of ice
by multicanonical simulations
B.A. Berg, C. Muguruma & Y.O., Phys. Rev. B 75, 092202 (2007);
B.A. Berg & W. Yang, J. Chem. Phys. 127, 224502 (2007);
C. Muguruma, Y.O. & B.A. Berg, Phys. Rev. E 78, 041113 (2008);
B.A. Berg, C. Muguruma & Y.O., Mol. Sim. 38, 856 (2012);
T. Hayashi, C. Muguruma & Y.O., J. Chem. Phys. 154, 044503 (2021).

S0experiment = 0.82(5)cal/deg/mole

W.F. Giauque and M. Ashley, Phys. Rev. 43 81 (1933).

S0experiment = 0.815(26)cal/deg/mole

S0 = k ln(W )where

S0Pauling = 0.80574...

O. Haida, T. Matsuo, H. Suga, and S. Seki,
J. Chem. Thermodynamics 6, 815 (1977).

W = (W0 ) NPauling0W

cal/deg/mole

3=
2

L. Pauling, J. Am. Chem. Soc. 57 2680 (1935).
Onsager pointed out this is just a lower bound.

L. Onsager and M. Dupuis, J. Rend. Sc. Inst. Fis. “Enrico Fermi” 10, 294 (1960).
Onsager’s student Nagle improved by series expansion:

S

Nagle
0

= 0.81480(20)cal/deg/mole

W0Nagle = 1.50685(15)J.F. Nagle, J. Mat. Phys. 7, 1484 (1966).

S0MUCA = 0.81550(21)cal/deg/mole

W0MUCA = 1.50738(16)B.A. Berg, C. Muguruma & Y.O., Phys. Rev. B 75, 092202 (2007).

S

MUCA
0

= 0.81527(17)cal/deg/mole

W0MUCA = 1.50721(13)C. Muguruma, Y.O. & B.A. Berg, Phys. Rev. E 78, 041113

Ice Rules
1. There is one hydrogen atom on each bond(then called hydrogen bond)
2. There are two hydrogen atoms near each oxygen atom(these constitute a water molecule)

Lattice structure of one layer of ice I

B.A. Berg, C. Muguruma & Y.O., Phys. Rev. B 75, 092202 (2007).
We proposed the following two Hamiltonians in which
only the groundstate satisfies both of the Ice Rules.
6-state H2O molecule model

E = -å h(b, s , s )
1
b

b

where

2
b

[Ice Rule (2) is always enforced]

ì 1, for H-bondh(b, s , s ) = í
î0, otherwise
1b

E = -å f (s, bs1 , bs2 , bs3 , bs4 )s2
b2-state H-bond model[Ice Rule (1) is always enforced]

ì2, for 2 H nuclei close to s
ï
1
2
3
4
f ( s, bs , bs , bs , bs ) = í1, for 1 or 3 H nuclei close to s
ï 0, for 0 or 4 H nuclei close to s
î

where

Single-Histogram Reweighting Techniques
A. Ferrenberg & R. Swendsen, Phys. Rev. Lett. 61, 2635 (1988).
The number of states n(E) is obtained from the histogram of the energy
distribution Nmu(E) that was obtained from the production run of themulticanonical simulation:

N mu ( E )n(E) =, where N mu ( E ) = n ( E )Wmu ( E ) .
Wmu ( E )Then the residual entropy is obtained by

S0 = k ln(n( E0 )), where E0 is the energy of the groundstate
where the normalization for n(E) is set byì6 N for the 6-state modelåE n( E ) = í22 N for the 2-state modelîat b = 0

Cf. B.A. Berg and T. Celik, Phys. Rev. Lett. 69, 2292 (1992).

B.A. Berg, C. Muguruma & Y.O., Phys. Rev. B 75, 092202 (2007).

S0 = k ln(W )where

S0MUCA = 0.81550(21)W = (W0 ) N
cal/deg/mole

W0MUCA = 1.50738(16)SUMMARY up to 2012
B.A. Berg, C. Muguruma & Y.O., Mol. Sim. 38, 856 (2012).

We have improved the accuracy of residual entropyof ordinary ice to beN
S0 = k ln(W )(1935)
Pauling
0

W

3
=
2

where

W = (W0 )(1966)(2008)

W0Nagle = 1.50685(15)W0MUCA = 1.50721(13)S0Pauling = 0.80574... S0Nagle = 0.81480(20)S0MUCA = 0.81527(17)But the accuracy of experiment is just one digit.
S0experiment = 0.815(26)O. Haida, T. Matsuo, H. Suga, and S. Seki,
J. Chem. Thermodynamics 6, 815 (1977).

We would like the calorimetry experimentalists toimprove the values.

New Calculations
In 2012, we used up to N = 2880 water molecules.
For more water molecules, Wang-Landau did not work
any more.
B.A. Berg, C. Muguruma & Y.O., Mol. Sim. 38, 856 (2012).

Our new method, REWL-MUCAREM (combination of
Replica-Exchange Wang-Landau method
T. Vogel, Y.W. Li, T. Wust & D.P. Landau, Phys. Rev. Lett. 110, 21 (2013).

and Multicanonical Replica-Exchange Method)
Y. Sugita & Y.O., Chem. Phys. Lett. 329, 261 (2000).

worked for N = 4704.
T. Hayashi, C. Muguruma & Y.O., J. Chem. Phys. 154, 044503 (2021).

New Calculations
T. Hayashi, C. Muguruma & Y.O., J. Chem. Phys. 154, 044503 (2021).

New Calculations
T. Hayashi, C. Muguruma & Y.O., J. Chem. Phys. 154, 044503 (2021).

Latest Value

Comparisons with Other Calculations
T. Hayashi, C. Muguruma & Y.O., J. Chem. Phys. 154, 044503 (2021).

Problem with Random Numbers:
T. Hayashi, C. Muguruma & Y.O., J. Chem. Phys. 154, 044503 (2021).

Marsaglia
G. Marsaglia, A. Zaman, and W.W. Tsang, Stat. Probab. Lett. 9,35 (1990).

vs.
Mersenne Twister
M. Matsumoto and T. Nishimura, ACM Trans. Model. Comput.
Simul. 8, 3 (1998).

Mersenne Twister is a better random number generator
than Marsaglia.

Problem with Random Numbers:
T. Hayashi, C. Muguruma & Y.O., J. Chem. Phys. 154, 044503 (2021).

In the range (0 〜 5)×10-7, Mersenne Twister gives auniform distribution of random numbers, while
Marsaglia gives only 9 values (green horizontal lines in
the left figure).

球状タンパク質の立体構造予測
T. Yoda, Y. Sugita & Y.O., Biophys. J. 99, 1637 (2010).
T. Yoda, Y. Sugita & Y.O., Proteins 82, 933-943 (2014).

Challenge the prediction of the 3-dimensional
structure of a small protein by MUCAREM.
ヴィリン ヘッドピース サブドメイン(36 残基; 596 原子)
半径 30 Åの水の球 (3513 個の水分子);
全原子数は 11,135 個

1

20

10

30

MLSDEDFKAVFGMTRSAFANLPLWKQQNLKKEKGLF
helix1helix2

helix3

Computational Details
T. Yoda, Y. Sugita & Y.O., Biophys. J. 99, 1637 (2010).

(Force Field = CHARMM22/CMAP for protein& TIP3P for water)(1) REMD with 96 replicas in implicit solvent (GB/SA);
initial conformation was fully extended(2) Unfolded protein w/o any secondary structures
was soaked in a sphere of radius 30Å (with 3513
TIP3P water molecules)(3) REMD with 128 replicas (T = 250 K ～ 700 K)(4) Determine multicanonical weight factors by WHAM
(iterate several times to refine weight)(5) Two production runs of MUCAREM with 8 replicas(MUCAREM1: 1.127 µs in total covering T = 269 K ～ 699 K
MUCAREM2: 1.157 µs in total covering T = 289 K ～ 699 K)

Villin headpiece subdomain(36 amino acids; 596 atoms)
in sphere of water of radius 30 Å
(3513 water molecules);altogether 11,135 atoms

MUCAREM simulation

Potential Energy Time Series

RMSD Time Series

MUCAREM2 (Replica 5)

動画提供：依田隆夫准教授（長浜バイオ大学）

Snapshots from MUCAREM2

T. Yoda, Y. Sugita & Y.O., Biophys. J. 99, 1637 (2010).

Native-Like Structures Obtained from
MUCAREM2 (Left) and MUCAREM1 (Right)

Main-Chain RMSD = 1.1 Å(residues 2 to 35) [Replica 5]

Main-Chain RMSD = 3.3 Å(residues 2 to 35) [Replica 8]
RMSD = 1.2 Å (residues 9 to 32)Native-Like Structures Obtained from MUCAREM
T. Yoda, Y. Sugita & Y.O., Biophys. J. 99, 1637 (2010).

Main-Chain RMSD = 1.1 Å (residues 2 to 35) [Replica 5]
灰色：自然の構造（PDB ID: 1YRF）、緑色：シミュレーションの結果

a-Helix Formation
Nn ≤ 7
Nn : no. of non-local side-chain
native contacts(=18 for X-ray structure)
Solid: 300 K,
Dashed: 350 K,
Short Dashed: 400 K,
Dotted: 500 K,
Dotted Dash: 600 K

Nn = 8, 9

Nn ≥ 10

T. Yoda, Y. Sugita & Y.O., Biophys. J. 99, 1637 (2010).

Salt Effects on Folding of a Small GlobularProteinT. Yoda, Y. Sugita & Y.O., Proteins 82, 933-943 (2014).
Challenging the prediction of the 3-dimensional structureof a small protein by MUCAREM.
Villin headpiece subdomain(36 amino acids; 596 atoms)
sphere of salted water with radius
30 Å(3494 water molecules,
11 K+, 13 Cl- ≈ 0.2 M KCl);Primary Sequence of HP-36
1

20

10

30

MLSDEDFKAVFGMTRSAFANLPLWKQQNLKKEKGLF
helix1helix2

helix3

Native-Like Structure (Global Minimum in Free
Energy) Obtained from MUCAREM Simulation(Left)
T. Yoda, Y. Sugita & Y.O., Proteins 82, 933-943 (2014).

Main-Chain RMSD = 1.25 ÅExperimental Structure(PDB ID: 1YRF)

Free Energy Landscape
T. Yoda, Y. Sugita & Y.O., Proteins 82, 933-943 (2014).

in Pure Water

in 0.2 M Salted Water遺伝的アルゴリズムによるタンパク質の立体構造予測

Prediction of Protein Structures by Genetic
Algorithms

T. Hiroyasu, M. Miki, M. Ogura & Y. O., J. IPS Jpn. 43, 70-79 (2002).
Y. Sakae, T. Hiroyasu, M. Miki & Y. O., Pacific Symposium on Biocomputing
16, 217-228 (2011).
Y. Sakae, T. Hiroyasu, M. Miki & Y. O., J. Comput. Chem. 32, 1353-1360(2011).
Y. Sakae, T. Hiroyasu, M. Miki, K. Ishii & Y.O., J. Phys.: Conf. Ser. 487,
012003 (2014).
Y. Sakae, T. Hiroyasu, M. Miki, K. Ishii & Y.O., Molecular Simulation 41,
1045-1049 (2015).
Y. Sakae, J. Straub, & Y.O., J. Comput. Chem. 40, 475-481 (2018).

Parallel Simulated Annealing MC/MD with Genetic Crossover(PSAMC/GAc2 and PSAMD/GAc2)

T. Hiroyasu, M. Miki, M. Ogura & Y. O., J. IPS Jpn. 43, 70-79 (2002).
Y. Sakae, T. Hiroyasu, M. Miki & Y. O., J. Comput. Chem. 32, 1353-1360 (2011).

High T

Low T

Parallel Simulated Annealing Simulations

Step 1. Parallel simulated annealing simulations for a certain MC/MD steps
Step 2. Genetic crossover
Repeat these two steps.
A

B

C

D

GA crossover
GA crossover
GA crossover
GA crossover

E

FGenetic Crossover (2-Point Crossover)Exchange a Randomly Chosen Pair of Dihedral Angle Sets
Structure A

Dihedral Angle

Structure B

All dihedral angles in n randomly selectedconsecutive amino-acid residues areexchanged. n = 2 – 10, motivated by thefragment assembly method (D. Baker).

Instead of simulated annealing, we
now repeat the following two steps:
Step 1. [Local Update] Parallel MC/MD
simulations at temperature T
Step 2. [Global Update]
Genetic crossoverHence, the detailed balance conditionis assured and canonical ensemble attemperature T is obtained.
Cf. it is similar to hybrid MC.
By introducing M temperature values,
we can easily generalize the method to
REMC/REMD.

Y. Sakae, T. Hiroyasu, M. Miki, K. Ishii & Y.O.,
Molecular Simulation 41, 1045-1049 (2015).
Y. Sakae, J. Straub, & Y.O.,
J. Comput. Chem. 40, 475-481 (2018).

Example 1
Trp-Cage (PDB ID: 1L2Y)
20 residues

method: MD (Langevin dynamics)
temperature: 282 K
solvent: GB/SA
force field: AMBER ff03
no. of individuals: 16
simulation time per individual: 100ps×100（10ns）

Simulation（Individual No.2）

PDB Structure（NMR）

RMSD : 0.809 Å

Y. Sakae, T. Hiroyasu, M. Miki, K. Ishii & Y.O., Molecular Simulation 41, 1045-1049 (2015).

Example 2
Villin Headpiece (PDB ID: 1YRF)
36 residues

method: MD (Langevin dynamics)
temperature: 300 K
solvent: GB/SA
force field: AMBER ff03
no. of individuals: 32
simulation time per individual: 200ps×100（39.4ns）

Simulation（Individual No. 11）

PDB Structure (X-ray)

RMSD : 2.234 Å

Y. Sakae, T. Hiroyasu, M. Miki, K. Ishii & Y.O., Molecular Simulation 41, 1045-1049 (2015Example 3
Protein A (PDB ID: 1BDD)
46 residues (10-55 out of 60)

method: MD (Langevin dynamics)
temperature: 300 K
solvent: GB/SA
force field: AMBER ff03
no. of individuals: 32
simulation time per individual: 1.0ns×90（90ns）

PDB Structure（NMR）
Simulation（Individual No. 5）
RMSD : 1.707 Å
Y. Sakae, T. Hiroyasu, M. Miki, K. Ishii & Y.O., J. Phys.: Conf. Ser. 487, 012003 (2014).

機械学習による蛋白質の立体構造予測（AlphaFold2）
with Ryuga Sugano

最近、AlphaFold,AlphaFold2というディープラーニングに基づく手法が、CASP(Critical Assessment of protein Structure Prediction)というタンパク質立体構造予測オンテストで圧倒的な正答率を挙げている。

CASP : タンパク質構造予測コンペ

https://doi.org/10.1126/science.370.6521.1144

https://www.itmedia.co.jp/news/articles/2012/
01/news053.html

評価指標
評価指標としてRMSD(平均二乗偏差)を用いる
min
平行移動、回転N : 比較したい全原子数
δi : i番目の原子の予測結果と対応する正解構造の距離(Å)
RMSD=0で完全一致、0に近いほど良い精度

結果

青：PDBの立体構造
緑：AlphaFold2の予測結果

結果(7LTW…二量体)

RMSD…1.485

青：PDBの立体構造
緑：AlphaFold2の予測結果

AlphaFold2
ブラウザ版とインストール版(名大不老)
ブラウザ版

インストール版(不老)7E1C
355残基.

1時間55分
RMSD:1.406

54分
RMSD:1.591

7FH8
676残基

2時間32分
RMSD:0.914

1時間34分
RMSD:0.854

7WHP_A
1058残基out of memory

1時間50分
RMSD:1.658

結果(7WHP_A,名大の不老を用いて)
7WHP_chainA
1058残基RMSD…1.658

計算時間1時間50分

青：AlphaFold2予測結果
緑：PDB構造応用例：新型コロナウイルスのオミクロン株のReceptor Binding Domain (RBD)の構造)

223残基のうち15残基に変異RMSD:0.822

ピンク：オリジナル：新型コロナウイルスのRBDのPDB立体構造
緑：オミクロン株のRBDのAlphaFold2の予測結果

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

Hamiltonian Replica-Exchange Method (HREM)
H. Fukunishi, O. Watanabe & S. Takada, J. Chem. Phys. 116, 9058 (2002).

T

13

14

15

16

9

10

11

12

5

6

7

8

1

2

3

4

Temperature Exchange
Parameter ExchangelCf. Two-dimensional MUCA:
J. Lee, M.A. Novotny, & P.A. Rikvold, Phys. Rev. E 52, 356 (1995).
J. Higo, N. Nakajima, H. Shirai, A. Kidera & H. Nakamura, J. Comput. Chem.
18, 2086 (1997).

Multidimensional Replica-Exchange Method(MREM)（also referred to as Hamiltonian REM）
Y. Sugita, A. Kitao & Y.O., J. Chem. Phys. 113, 6042 (2000).

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

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

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

Van der Waals Replica-Exchange Method(vWREM)
S.G. Itoh, H. Okumura & Y.O., J. Chem. Phys. 132, 134105 (2010).
HamiltonianH l (q, p) = K (p) + E l (q )Potential Energy

E l (q ) = E P (q P ) + E PS (q P , qS ) + ES (qS )where P stands for protein and S solvent

EP (qP ) = Vl (q P ) + Erest (q P )and

N P -1Vl (qP ) = å
k =1

éæ ls ö12 æ ls ö6 ù
kl
kl
ê
ú
4
e
ç
÷
ç
÷
å
kl ç
ç r ÷ úêè rkl ÷ø
l =1
è kl ø ûëNP

Namely, the van der Waals radii of the atoms only in the P-P interactions are
scaled by l.

Van der Waals Replica-Exchange Method(vWREM)
S.G. Itoh, H. Okumura & Y.O., J. Chem. Phys. 132, 134105 (2010).

Alanine dipeptidefyDecrease atom size
Original atom size

Move easily!

Time series of dihedral
angles f

Van der Waals
REMD>

Conventional
REMD

f

y

>

Canonical
MD

Effectiveness of samplings

Multi-Overlap Algorithm
MC version: B. Berg, H. Noguchi & Y.O., Phys. Rev. E 68, 036126 (2003).
MD version: S.G. Itoh & Y.O., Chem. Phys. Lett. 400, 308 (2004);
S.G. Itoh & Y.O., J. Chem. Phys. 124, 104103 (2006).

Multi-overlap algorithm realizes a random walk in theconfigurational space between two appropriate reference
configurations.

・dihedral-angle distance
A dihedral-angle distance is defined by
1 n
0
d=
d
(
v
,
v
å
a
i
i ).
p n i =1

where vi is the dihedral angle i, vi0 is the dihedral angle i0of the reference configuration, and the distance d a (vi , vi )is given byd a (vi , vi0 ) = min(| vi - vi0 |, 2p - | vi - vi0 |)0 £ d a (vi , vi0 ) £ p .

Therefore,0 £ d £ 1.

If d = 0, all dihedral angles are coincident with those of thereference configuration.
dihedral-angle distance is an indicator of how similar theconfiguration is to the reference configuration.

In the canonical ensemble at a constant temperature, the probability
distribution with respect to two dihedral-angle distances d1 and d 2
is given by
Pc (d1 , d 2 ) = ò dEn(d1 , d 2 ; E ) exp( - b E )In the multi-overlap ensemble at a constant temperature, utilizing
a function f (d1 , d 2 ) , the probability distribution with respect totwo dihedral-angle distances is constant:Pmuov (d1 , d 2 ) = ò dEn(d1 , d 2 ; E ) exp(- b E + f (d1 , d 2 )) = const
Pc (d1 , d 2 )0Pmuov (d1 , d 2 )1 d1 , d 2

0

1 d1 , d 2

The multi-overlap weight factor is defined by
Wmuov = exp(- b Emuov ) = exp[- b ( E - k BTf (d1 , d 2 )] .

The multi-overlap MD simulation is performed by solving
the following modified Newton equation:
dpi
¶E
¶E
¶f (d1 , d 2 )= - muov = + k BT
.
dt
¶qi
¶qi
¶qi

The function f (d1 , d 2 ) is determined by the following process.
At first, we set f (d1 , d 2 ) = 0 . Then, if d1 = c1 , d 2 = c2 after ith
MD step, the function f (d1 , d 2 ) is updated byfi (d1 = c1 , d 2 = c2 ) = f i -1 (d1 = c1 , d 2 = c2 ) - a .

Here a is an appropriately chosen positive constant. This
procedure is iterated until the probability distribution Pmuov isflat.
(S. G. Itoh and Y. Okamoto, Chem. Phys. Lett. 400 (2004) 308)

The density of states is defined byn( E ) = ò dd1dd 2 N (d1 , d 2 ; E ) exp( b E - f (d1 , d 2 )),where N (d1 , d 2 ; E ) is the histogram of the probability
distribution.
The expectation value of a physical quantity A at any
temperature T ¢ is give by

A T¢ =- b ¢E
dEA
(
E
)
n
(
E
)
e
ò

ò dEn( E )e- b ¢Esimulation：canonical MD simulation
multi-overlap MD simulation
system：Met-enkephalin (Tyr-Gly-Gly-Phe-Met) in vacuum.
the N-terminus and the C-terminus were blocked with theacetyl group and the N-methyl group, respectively.
temperature：300K (Gaussian constraint method)force field：CHARMM param 22 (all-atom model)
time step：0.5fs (leap-frog method)
total simulation time：24ns (+1ns equilibration)

We performed a multi-overlap MD simulation for Met-enkephalinwith two reference configurations.

Reference configuration 1

Reference configuration 2

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

d2

d1

0.8

0.3

0.3

0.2

0.2

0.1

0.1

0
0

5000

4

1 10

4

1.5 10time (ps)

4

2 10

4

2.5 10

0

0

5000

4

1 10

time (ps)The time series of the dihedral-angle distances in aconventional canonical MD simulation.

4

1.5 1042 10

4

2.5 10

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

d2

d1

0.8

0.3

0.3

0.2

0.2

0.1

0.1

0

0

5000

4

1 10

4

1.5 10time (ps)

4

2 10

4

2.5 10

0

0

5000

4

1 10

4

1.5 10time (ps)The time series of the dihedral-angle distances in themulti-overlap MD simulation.

4

2 10

4

2.5 10Canonical MD

Multicanonical MD

r1 < 0.5 Þ
r2 < 0.5 Þ
Simulation and movie byS.G. Itoh
Multioverlap MD

Multicanonical MD

Multioverlap MD

400

400

350

350

350

300

300

300

250

250

250
200
150

E(kca/mol)

400

E(kcal/mol)

E(kcal/mol)Canonical MD

200
150

200
150

100

100

100

50

50

50

0
0

5000

4

1 10

4

1.5 10time (ps)

4

2 10

4

2.5 10

0

0

5000

4

1 10

4

1.5 1042 10

4

2.5 100
0time (ps)Time series of potential energy

5000

4

1 10

4

1.5 10time (ps)

4

2 10

4

2.5 10Average potential energy and specific heat asfunctions of temperature:
Solid curve: MUOV, Dashed curve: MUCAMUOV results are reliable only in the range～250 K < T < ～350 K

A1ref. 1

A2

ref. 2BC
The free-energy landscape obtained
from the multi-overlap MD
simulation at T=300K with respect
to the RMSD axes r1 , r2 .

C

Multicanonicaｌ-Multioverlap Algorithm
S.G. Itoh & Y.O., Phys. Rev. E 76, 026705 (2007).

In canonical ensemble at a constant temperature, the probability
distribution with variables E , d is given by
Pc ( E , d ) µ n( E; d )Wc ( E; T0 ) º n( E; d ) exp( - b 0 E ) ,where n( E; d ) is the density of states
We employ Emco ( E , d ) so as to obtain a constant probability
distribution with E , d
Pmco ( E , d ) = n( E; d )Wmco ( E; d )º n( E; d ) exp(- b 0 Emco ( E , d ))= const・再重法
状態密度は次のように表されるn( E ) = å N ( E; d ) exp( b 0 Emco ( E , d ))d

ここで N ( E ; d ) はシミュレーションにより得られたヒストグラム.
この状態密度を用いて,物理量 A の温度T での統計平均を
次のように得ることができる
-b E
A
(
E
)
n
(
E
)
e
å

AT = E

-b E
n
(
E
)
e
å
E

・シミュレーション条件
シミュレーション：マルチカノニカル‐マルチオーバーラップ
MDシミュレーション
マルチカノニカルMDシミュレーション
マルチオーバーラップMDシミュレーション
系：真空中のMet-enkephalin (Tyr-Gly-Gly-Phe-Met)
N末はアセチル基,C末はN-メチルアミド基でそれぞれ修飾.
系の温度：300K (ガウス束縛法)
力場：CHARMM param 22 (all-atom model)
時間刻み：0.2fs (leap-frog法)
全シミュレーション時間：4ns (+1nsの平衡化)

Reference
conformation

350

300

300

300

250

250

250

200
150

200
150

0

500 1000 1500 2000 2500 3000 3500 4000

time step (ps)Multicanonical MD

50

200
150
100

100

100
50

E (kcal/mol)

350

E (kcal/mol)

E (kcal/mol)

350

0500 1000 1500 2000 2500 3000 3500 4000

time step (ps)Multioverlap MD

50

0

500 1000 1500 2000 2500 3000 3500 4000

time step (ps)Multicanonicalmultioverlap MD

Time series of the potential energy for each simulations
MUCA-MUOV now allows accurate reweighting in awide temperature range.

0.8

0.7

0.7

0.7

0.6

0.6

0.6

0.5

0.5

0.5

0.4

0.4

0.4

d

0.8

d

d

0.8

0.3

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0

0
0

500 1000 1500 2000 2500 3000 3500 4000

time step (ps)Multicanonical MD

0

500 1000 1500 2000 2500 3000 3500 4000

time step (ps)Multioverlap MD00

500 1000 1500 2000 2500 3000 3500 4000

time step (ps)Multicanonicalmultioverlap MD

Time series of the dihedral-angle distance for each simulations

Average potential energy and specific heat as functionsof temperature:
Blue and black curves: MUCA & MUCA-MUOV
Red curve: MUOVMUCA and MUCA-MUOV results are reliable in a wide
temperature range.

S.G. Itoh & Y.O., Phys. Rev. E 76, 026705 (2007).

Multicanonical MD

Multioverlap MD

Multicanonicalmultioverlap MD

Histogram with E and d axes for each simulations

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 parameter
MREM: random walk in multidimensional parameter
e.g.,
WHAM eqns.
M()

n E ,V =åNm ,n =1
Mm ,nå nm,nem ,n =1(E ,V )

fm ,n - bm (E + lnV ), where e- fm ,n()= å n E ,V e
E ,V- bm (E + lnV ).

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).

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).

l = h = external field
V = M = magnetization
1. 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.

T. Nagai & Y.O., Physics Procedia 34, 100 (2012).
T. Nagai & Y.O., Phys. Rev. E 86, 056705 (2012).

Application of Two-dimensional
ST method to Two-dimensional
Ising model and Its CrossoverSystem and method
Hamiltonian: H = E − hM,

where E = − ∑ <i , j> σ iσ j , M = ∑σ ii

Two-dimensional
Ising model (4X4)− β (E−hM )+a(T ,h)Weight: w({σ i },T ,h) = e“Simulated Tempering and Magnetizing” (STM).
（Both temperature T and external field h are updated）ℎ

critical point
~2.27

ß

ℎ

ℎ=0

𝑇
First-order

à𝑇

Second-orderReweighting methods
The WHAM or MBAR method enable one to obtain thedensity function and free energy of systems.
- This gives free energy necessary for the ST methods.
- Knowing the density function let one calculate thermal
average at any point.

WHAM

S. Kumar, et al.: J. Comput. Chem. 13, 1011 (1992).
A. Mitsutake and Y.O.: J. Chem. Phys. 130, 214105 (2009).

Reweighting methods
MBAR

M. Shirts and J. Chodera: J. Chem. Phys. 129, 124105 (2008).

Simulation and movie by T. Nagai

Results and Discussion

𝑇, ℎ• In this simulation, thetemperature andexternal field realized
the random walks.
Therefore thesimulations are
performed properly.

緑：σ=+1赤：σ=-1

Temperature: 20 points in [1:5]
External field: 21 points in [-1.5, 1.5]
System: 20X20

T. Nagai & Y.O., Phys. Rev. E 86, 056705 (2012).

Results and Discussion• To investigate the effect of multidimensional generalization,– We performed “Simulated Magnetizing (SM),”
where 𝑇 is fixed and ℎ is allowed to change.

T. Nagai & Y.O., Phys. Rev. E 86, 056705 (2012).

Simulation and movie by T. Nagai

Results and Discussion
T is fixed at 1.97 (below Tc)

T is fixed at 3.88 (above Tc)

Below the Tc, H
changes in only one
side. Do not go over 0.

Simulation and movie by T. Nagai

Results and Discussion
This animation illustrates how the temperature and external
field moved during the STM simulation.
ℎ

ß

ℎ=0

à𝑇
At every frame,
5000 MC sweeps.
Change by 1000
MC sweeps.

ST-Freq: 50
TvsH=0.5
L=20

Results and Discussion
T. Nagai & Y.O., Phys. Rev. E 86, 056705 (2012).

M as a function of Lt when h=0.

Lt

We applied the finite-size scaling form.
Here, M ~ t1/8
Ref) D.P. Landau: Phys. Rev. B 13, 2997 (1976).

Results and Discussion
T. Nagai & Y.O., Phys. Rev. E 86, 056705 (2012).
The Binder cumulant for h = 0 .
The curves were obtained with the reweighting method.

K. Binder, Z. f. Physik B 43, 119 (1981).

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. Grandcanonical Ensemble

λ = -μ = - chemical potential

V = N = no. of particles＊ REM: random walk intemperature T and chemical potential μ
D. Matsubara & Y.O., J. Chem. Phys. 152, 194108 (2020).
Lennard-Jones Fluid System:→
Reduced unit: r* = r/σ , E* = E/ε , μ* = μ/ε , T* = kBT/ε , V* = L*3, L* =L/σ ,= <N>/V* , V* = L*3

ρ* = ρσ3Grandcanonical MC(GCMC)
Translation:

Insert:Delete:Metropolis Criterion

Grandcanonical REM (GCREM)
M Replicas with

D. Matsubara & Y.O., J. Chem. Phys. 152, 194108 (2020).

Tm : T1 , T2 , ! , TM00µ m : µ1 , µ2 , ! , µ M
1

1

for T-exchangefor μ-exchange
Repeat these two steps.

RESULTS

Grandcanonical REM (GCREM)
D. Matsubara & Y.O., J. Chem. Phys. 152, 194108 (2020).
Random walk in T-μ space

Acceptance ratios of T-exchange (a,b) and for μ-exchange (c,d)

Grandcanonical REM (GCREM)RESULTSD. Matsubara & Y.O., J. Chem. Phys. 152, 194108 (2020).

Average number of particles <N>

Particle number fluctuations CμGrandcanonical REM (GCREM)RESULTSD. Matsubara & Y.O., J. Chem. Phys. 152, 194108 (2020).

Binder Cumulant

where the number density is used for m:m = ρ* =<N>
V*

These results give: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. Isobaric-Isothermal Ensemble（定圧定温アンサンブル）＊ MUCA: Multibaric-Multithermal Algorithm (MUBATH)
random walk in potential 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., in Lect. Notes in Computational Science & Engineering,ed. by T. Schlick and H.Gun (2002) pp. 304-332; cond-mat/0102296.
T. Okabe, M. Kawata, Y.O. & M. Mikami, Chem. Phys. Lett. 335, 435 (2001).
Also, see D. Paschek & A. Garcia, Phys. Rev. Lett. 93, 238105 (2004).

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

Liquid-Gas Phase Transitions Studied byMultibaric-Multithermal Simulations
MUBATH

Conventional NPT

H. Okumura & Y.O., J. Phys. Soc. Jpn. 73, 3304 (2004).

Liquid-Gas Phase Transitions Studied byMultibaric-Multithermal Simulations
MUBATH

Conventional NPT

H. Okumura & Y.O., J. Phys. Soc. Jpn. 73, 3304 (2004).

Liquid-Gas Phase Transitions Studied byMultibaric-Multithermal Simulations

H. Okumura & Y.O., J. Phys. Soc. Jpn. 73, 3304 (2004).

Liquid-Gas Phase Transitions Studied byMultibaric-Multithermal Simulations
H. Okumura & Y.O.,
J. Phys. Soc. Jpn. 73, 3304 (2004).

MUBATH

NPT (reweighted)

Liquid-Gas Phase Transitions Studied byMultibaric-Multithermal Simulations

H. Okumura & Y.O., J. Phys. Soc. Jpn. 73, 3304 (2004).

ST Simulation in Isobaric-Isothermal Ensemble
Y. Mori & Y.O., J. Phys. Soc. Jpn. 79, 074003 (2010);
Y. Mori & Y.O., J. Comput. Chem. 38, 1167-1173 (2017);
G. La Penna, Y. Mori, R. Kitahara, K. Akasaka, & Y.O., J. Chem. Phys. 145,
085104 (2016);
Y. Mori et al., in preparation.

temperature and pressure become dynamical variables.
random walk in temperature and pressure
CF. Multibaric-Multithermal Algorithm (MUBATH)
random walk in potential 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)
CF. Y. Mori & H. Okumura, J. Phys. Chem. Lett. 4, 2079–2083 (2013).
Cf. previous works:
D. Paschek & A. Garcia, Phys. Rev. Lett. 93, 238105 (2004);
T. Imai, S. Ohyama, A. Kovalenko, & F. Hirata, Protein Sci. 16, 1927 (2007);
T. Imai & Y. Sugita, J. Phys. Chem. B 114, 2281 (2010).

Pressure-Induced Unfolding of Ubiquitin
NMR Experiments
30 bar – 3000 bar

Simulated System• 76 amino acids• 6232 water molecules• 19985 atoms

PDB: 1V80
R. Kitahara, S. Yokoyama & K. Akasaka, J. Mol. Biol. 347, 277 (2005).

Time series of pressure P, potential energy E, andvolume V
for ubiquitin

Y. Mori & Y.O., J. Comput. Chem. 38, 1167 (2017).

Large structural fluctuations• Fluctuations of distance d between pairs of
Ca atoms.

f [Å]Amino Acid Residues
Large fluctuations observed in agreement with experiments.

Structural changes under high pressurer
Simulation and movie by Y. Mori

distribution
[Å-1]

r
r [Å]

Ubiquitin and water molecules

at low pressure
Simulation and movie by Y. Moriat high pressure

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).
4. Umbrella SamplingL

Hk (q, p) = K( p) + E0 (q) + å l Vl (q), where Vl (x) = Kl (x (q) - dl ) .
l

2

l =1
＊ REM: Replica-Exchange Umbrella Sampling(REUS)
random walk in reaction coordinate x

Y. Sugita, A. Kitao & Y.O., J. Chem. Phys. 113, 6042 (2000).

＊ ST: Simulated Tempering Umbrella Sampling (STUS)
random walk in reaction coordinate x

Y. Mori & Y.O., Phys. Rev. E 87, 023301 (2013).

＊ MUCA: Metadynamics (Wang-Landau in reaction coordinate)
random walk in reaction coordinate x

A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. USA 99, 12562 (2002).
see also:
T. Huber, A.E. Torda, W.F. van Gunsteren, J. Comp. Aid. Mol. Des. 8, 695 (1994);
H. Grubmuller, Phys. Rev. E 52, 2893 (1995).

Prediction of Protein-Ligand Binding Structures byReplica-Exchange Umbrella Sampling
H. Kokubo, T. Tanaka & Y.O., J. Comput. Chem. 32, 2810-2821 (2011);
H. Kokubo, T. Tanaka & Y.O., J. Chem. Theor. Comput. 9, 4660-4671 (2013);
H. Kokubo, T. Tanaka & Y.O., J. Comput. Chem. 34, 2601-2614 (2013);
Y.O., H. Kokubo & T. Tanaka, J. Chem. Theor. Comput. 10, 3563-3569 (2014);
S. Tsukamoto, Y. Sakae, Y. Itoh, T. Suzuki & Y.O., J Chem. Phys. 148, 125102 (2018);
H. Mishima, Y. Itoh, T. Kurohara, T. Suzuki, N. Asada, K. Kusakabe, and Y.O., J. Comput. Chem.
44, 1604-1609 (2023).

Replica-Exchange Umbrella Sampling (REUS)
Y. Sugita, A. Kitao & Y.O., J. Chem. Phys. 113, 6042 (2000).
L

Hk (q, p) = K( p) + E0 (q) + å l Vl (q), where Vl (x) = Kl (x (q) - dl ) .
.

l =1

M

- b( E+Vl (x ))
N
(E,
x
)e
å m

Potential of Mean ForcePb ,l (E, x ) = m =1M

1

Wb (x ) = - ln Pb ,l = 0 (x ).

b

2

l

ån e
m

WHAM

(

)

f m - bm E+ Vl m (x ),m=1

e - f m = å å Pb m ,l m (E,x ).

E
xSee also:
E.M. Boczko & C.L. Brooks, J. Phys. Chem. 97, 4509 (1993).
B. Roux, Comp. Phys. Commun. 91, 275 (1995).

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 and2-aminopyrimidine
Asp-78 and ligand

Simulation and movie
by H. Kokubo

Results: REUS with 24 replicas•

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

crystalpredicted

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

potential of mean force shows
the most stable distancePrediction of Protein-Ligand Binding Structures by
2-dim H-REMD: Replica-Exchange Umbrella
Sampling (REUS) and Replica-Exchange Solute
Tempering (REST): REUS/REST
H. Kokubo, T. Tanaka & Y.O., J. Comput. Chem. 34, 2601-2614 (2013) .
Cf. REUS/gREST: S. Re, H. Oshima, K. Kasahara, M. Kamiya & Y.Sugita,
PNAS 116, 18404-18409 (2019).

Potential Energy:

bn
bn
[i ]
Emn (q ) = U ll (q ) +U ls (q[i ] ) + U ss (q[i ] ) + Vm (q[i ] )
b0
b0
[i ]

l: ligands: protein/water

Total no. of replicas: M × N
REST parameters (n=1, 2, …, N)
REUS parameters, i.e., umbrella potentials(m=1, 2, …, M)

REST: P. Liu, B. Kim, R.A. Friesner, & B.J. Berne, PNAS 102, 13749 (2005).
T. Terakawa, T. Kameda, and S. Takada, J. Comput. Chem. 32, 1228 (2011).
REUS: Y. Sugita, A. Kitao & Y.O., J. Chem. Phys. 113, 6042 (2000).

2-dim H-REMD: REUS (M=24) + REST (N=8)(192 replicas)
H. Kokubo, T. Tanaka & Y.O., J. Comput. Chem. 34, 2601-2614 (2013) .

REUS (1T4E)

No. ofRandom Walk
Cycles

REUS/REST (1T4E)

Prediction of Ligand Binding Affinity by 2-dim HREMD: REUS/REST simulation
Y.O., H. Kokubo & T. Tanaka, J. Chem. Theory Comput. 10, 3563-3569 (2014).

Studies for Ligand Binding Selectivity ofHistone Deacetylase Inhibitor by Molecular
Dynamics Simulation
S. Tsukamoto, Y. Sakae, Y. Itoh, T. Suzuki, Y.O.,
J. Chem. Phys. 148, 125102 (2018).

Epigenetic regulations of gene
expression
Transcription isaccelerated

Gene activation

Me

Me

Condensed
chromatin
histone(HAT) acetylation

Deacetylation (HDAC)(KDM) H3K4 methylation

H3K9 methylation (G9a)Ac

acetyl group

Ac

Relaxedchromatin

DNA

Ac

Me

Me
Me

methyl group

Transcription isrepressed

A
c

Gene silencing

HDAC isozymes
There are 18 HDAC isozymes. They
are divided into 4 classes by
sequence homology and functions.
The HDACs which belong to thesame class have very similar
sequences and structures. When weuse HDAC inhibitors for therapy, we
should consider selectivity of theinhibitors. Unnecessary inhibition
may cause critical side effect.

• Class 1: HDAC1, HDAC2, HDAC3, HDAC8• Class 2: HDAC4, HDAC5, HDAC7, HDAC9,HDAC6, HDAC10• Class 3: SIRT1, SIRT2, SIRT3, SIRT4,
SIRT5, SIRT6, SIRT7• Class 4: HDAC11

Yellow: Crystal structureof HDAC2 (pdbID: 3max).
Blue: Crystal structure of
HDAC3 (pdbID: 4a69)

Purpose of this study
T247
T. Suzuki, Y. Kasuya, Y. Itoh, Y. Ota, P.
Zhan, K. Asamitsu, H. Nakagawa, T.
Okamoto, N. Miyata PLoS ONE, 8,
e68669 (2013)

1. To reproduce the selectivity of T247 byusing molecular dynamics simulations
2. To elucidate the reason why T247 has
the selectivity for HDAC3
3. To propose new selective inhibitors

REUS MD Simulation (One of the Replicas)
Simulation and movie by S. Tsukamoto

Potential of mean force
We succeeded in reproducing the selectivity.
HDAC2

HDAC3

The most
stable statePotential of mean force
S. Tsukamoto, Y. Sakae, Y. Itoh, T. Suzuki, Y.O., J. Chem. Phys. 148,
125102 (2018).

HDAC2

Docking state

HDAC3

Docking stateHydrophobic structures at binding site
S. Tsukamoto, Y. Sakae, Y. Itoh, T. Suzuki, Y.O., J. Chem. Phys. 148, 125102 (2018).

HDAC3

HDAC2

Blue: T247 Red: hydrophobic structures

Weak hydrophobic
interaction among T247
and hydrophobic residues
in HDAC2 causes weakRadii of gyration of the hydrophobic structures
S. Tsukamoto, Y. Sakae, Y. Itoh, T. Suzuki, Y.O., J. Chem. Phys. 148, 125102 (2018).

205

HDAC2 G E - Y F P
HDAC3 G N Y F F PStructure of wild type and mutant of HDAC2

Phe:206

Gly:150
Phe:200

His:17
9

Tyr:205

Gly:143

Phe:151

Leu:272

HDAC2-WT
Hydrophobic structure was
openedPhe:144

Phe:199
His:172
Leu:26
6

HDAC2-Y198F

Hydrophobic structure wasclosed焼き戻し傘サンプル法
Simulated Tempering Umbrella Sampling (STUS)
Y. Mori & Y.O., Phys. Rev. E 87, 023301 (2013).
trial

trial

trial

U

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

U

trial

U'

accept
自由エネルギー
free energy

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

反応座標
Reaction coordinate
無次元化された自由エネルギー Dimensionless Free
Energy時間 timeマロンアルデヒドの分子内プロトン移動
Proton Transfer in Malonaldehyde
Y. Mori & Y.O., Phys. Rev. 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., Phys. Rev. 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/シミュレーション結果 Results
Y. Mori & Y.O., Phys. Rev. E 87, 023301 (2013).

Simulation and movie by Y. Mori

Potential of Mean Force （平均力ポテンシャル）
Y. Mori & Y.O., Phys. Rev. E 87, 023301 (2013).

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

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

Potential of Mean Force （平均力ポテンシャル）
Y. Mori & Y.O., Phys. Rev. E 87, 023301 (2013).

red：STUS
green：conventional US

S. Sadhukhan, D. Munoz, C. Adamo, & G.E. Scuseria, Chem. Phys. Lett. 306(2008) 83.

Formation of one phthalocyanine from fourphthalonitriles and one iron atom with REUS and
Density Functional Tight Binding (DFTB) method
S. Ito, S. Irle, and Y.O., Comput. Phys. Commun. 204, 1-10 (2016);
S. Ito, Y. Wang, Y.O. & S. Irle, J. Chem. Phys. 149, 072332 (2018).

Initial coordinates• A simulation was already tried with only DFTB, but thephthalocyanine was not formed.
• Therefore, we carried out the simulation in this mechanism with
REUS and DFTB by DFTB+/REM.

Movies of trajectories
Simulation and movie by S. Ito

Product: Structure M

Product: Structure P(Phthalocyanine)Conjecture on scheme for forming phthalocyanine••

S. Ito, S. Irle, and Y.O., Comput. Phys. Commun. 204, 1-10 (2016);
S. Ito, Y. Wang, Y.O. & S. Irle, J. Chem. Phys. 149, 072332 (2018).

Implementation of REUS to GAMESS•

S. Ito, Fedorov, Y.O. & S. Irle, Comput. Phys. Commun. 228, 152 (2018).

DFTB-REMD Simulations of Deca-Alanine
with/without D3 Dispersion Corrections•

S. Ito, Fedorov, Y.O. & S. Irle, Comput. Phys. Commun. 228, 152 (2018).

Red: DFTB;

Blue: DFTB with D3 dispersion corrections

SUMMARY
We have shown thatgeneralized-ensemble algorithms
are particularly suitable forbiomolecular simulations.

COLLABORATORS (Former Group Members)
Ayori MITSUTAKE [IMS à Keio Univ. à Meiji Univ.]
Ulrich H.E. HANSMANN [IMS à Michigan Tech. Univ. à Univ. of Oklahoma]
Masato MASUYA [IMS à Kagoshima Univ.]
Yuji SUGITA [IMS à Univ. Tokyo à RIKEN]
Takao YODA [IMS à Nagahama Inst. Bio-Science]
Hiroshi NOGUCHI [IMS à Juelich à Univ. of Tokyo]
Hisashi OKUMURA [IMS à Nagoya Univ. à Rutgers Univ. à IMS]
Satoru G. ITOH [IMS à Nagoya Univ. à NIH à IMS]
Yoshitake SAKAE [IMS à Hiroshima Univ. à Nagoya Univ. à
Temple Univ. à RIST]
Hironori KOKUBO [IMS à Univ. of Houston à Takeda Pharm. à
Chugai Pharm.]
Yoshiharu MORI [Nagoya Univ. à IMS à Kitasato Univ. à Kobe Univ.]
Tetsuro NAGAI [Nagoya Univ. à Ritsumeikan Univ. à Univ. of Tokyo à
Fukuoka Univ.]
Ryo URANO [Nagoya Univ. à Boston Univ. à Nagoya Univ. à Okayama Univ.]
Shuichiro TSUKAMOTO [Nagoya Univ.]
Shingo ITO [Nagoya Univ. à Boston Univ. à RIKEN]
Yuhei TACHI [Nagoya Univ. à Deloitte Analytics]
Daiki MATSUBARA [Nagoya Univ. à RIKEN]
Takuya HAYASHI [Nagoya Univ. à Kureha]
Kohei NODA [Nagoya Univ. à JSR]
Ryuga SUGANO [Nagoya Univ. à NHK]COLLABORATORS
Hikaru KAWAI [Cornell Univ. à Univ. of Tokyo à KEK à Kyoto Univ. à
National Taiwan Univ.]
Takeshi KIKUCHI [Cornell Univ. à Univ. of Tokyo à Ciba-Geigy à
Ritsumeikan Univ.]
Masataka FUKUGITA [Kyoto Univ. à Univ. of Tokyo]
Takashi NAKAZAWA [Nara Women’s Univ.]
Fumio HIRATA [IMS à Ritsumeikan Univ. à Toyota RIKEN]
Masahiro KINOSHITA [Kyoto Univ.]
Jose ONUCHIC [Univ. of San Diego à Rice Univ.]
Bernd A. BERG [Florida State U.]
Chizuru MUGURUMA [Chukyo Univ.]
Tomoyuki HIROYASU [Doshisha Univ.]
Mitsunori MIKI [Doshisha Univ.]
Katsuya ISHII [Nagoya Univ.]
Takahiro KATAGIRI [Nagoya Univ.]
Satoshi OHSHIMA [Nagoya Univ.]
Ichiro TAKAHASHI [Nagoya Univ.]
Toshimasa TANAKA [Takeda Pharm.]
Ryuichi UEOKA [Sojo Univ.]
Wolfhard JANKE [Univ. of Leipzig]
John STRAUB [Boston Univ.]
Takayoshi SUZUKI [Kyoto Prefect. Univ. of Medicineà Osaka Univ.]
Yukihiro ITOH [Kyoto Prefect. Univ. of Medicineà Osaka Univ.]
Stephan IRLE [Nagoya Univ. à Oak Ridge National Lab.]
Dmitri G. FEDOROV [AIST]