GENESIS Hands-on Part 3:
Generalized-ensemble simulations using GENESIS
Shingo Ito
RIKEN Center for Computational Science
2024/07/02
IUPAB2024 Hands-on Training Program CHARMM-GUI/GENESIS MD Tutorial0Schedule of GENESIS parts (6/30-7/2)
06/30 Part 1
13:30 – 15:00

GENESIS basics and GENESIS on Fugaku (Kobayashi)
LectureHands-on tutorial on Fugaku

07/01 Part 2
14:30 – 15:30

Coarse-grained simulations in GENESIS (Tan)

15:30 – 16:30

High-performance computation with GENESIS (Jung)

07/02 Part 3
13:30 – 15:00

Generalized-ensemble simulations using GENESIS (Ito)1

Contents• What are Generalized-Ensemble simulations?
• Demonstration of Temperature Replica-Exchange MD

2

Molecular dynamics (MD) simulation• Molecular dynamics (MD) simulation is useful for understanding
biological functions such as protein folding, membrane transport, ligand
binding, etc.
PTI (1977)
58 res: < 1000 atoms

Bacterial cytoplasm (2016)> 100 M atoms
40 years later...

J. A. McCammon, et. al., Nature 267, 585–590 (1977)

I. Yu, et. al., eLife 5, e19274 (2016)3

Common issues of MD simulation
Folding Funnel• The folding funnel hypothesis suggests that protein
conformation reaches the global minima (= thermalequilibrium) after infinite-time scale MD simulations.

But...
• We can perform only finite-time (ns–μs) scale MD simulations.

K. E. Dill, et al., Science. 338, 1042-1046 (2012)4

Common issues of MD simulation
Folding Funnel
Free energy• Sometimes, finite time-scale MD simulations get trapped at
one of the local minima.
• The sampling in the limited conformational space does notsatisfy the ergodicity.
BarrierLocal minima• We can not estimate physical quantities.

How can the finite time-scale MD simulation
overcome the barrier?
Global minima= Native structure

5Generalized-ensemble simulations• Generalized-ensemble simulations is one of the solutions= Conventional MD + enhanced sampling method• Generally, there are two types of enhanced sampling methods
1.

Using thermal fluctuations at high temperatures (> 300.0 K)• Simulated tempering• Temperature replica exchange

2.

Smoothing the potential/free energy surface using
an additional potential• Umbrella sampling• Metadynamics
YouTube movie by Giovanni Bussi
https://www.youtube.com/watch?v=IzEBpQ0c8TA6

Simulated Tempering MD

Folding simulationof Chignolin
Figures from H. Oshima,et al., JCTC (2019)• In the Simulated Tempering (ST) MD simulation,
the temperature changes from 𝑇! to 𝑇!±# .
• The protein, RNA/DNA etc., become more
flexible at the high temperature.
A. Irbäck, et al., J. Chem. Phys. 103, 10298 (1995)
U.H.E. Hansmann, et al., J. Comput. Chem. 18, 920 (1997)7

Simulated Tempering MD
Free energy• Simulated tempering (ST) performs aHigh temperature

random walk in temperature space.

𝑊,-./.),-0 𝐸 = exp(−𝛽𝐸)

Room temperature

𝑊1$ 𝐸; 𝑇 = exp(−𝛽𝐸 + 𝛼(𝑇))discretizationConformationMDP(E)𝑃$%&''

𝑃1$ (𝑇) = 0 𝑑𝐸 𝑛 𝐸 𝑊1$ 𝐸; 𝑇 = constant

ST-MD
P(E)

𝑃$%&''

𝑊1$ 𝐸; 𝑇! = exp(−𝛽! 𝐸 + 𝛼(𝑇! ))𝑃$%()*+𝑷𝑺𝑻
A. Irbäck, et al., J. Chem. Phys. 103, 10298 (1995)
U.H.E. Hansmann, et al., J. Comput. Chem. 18, 920 (1997)

Potential energy

Potential energy

8Simulated Tempering MD
Free energy• This is a procedure of the ST-MD simulation.
High temperature• 𝛼! is decided by the short MD simulations.

Room temperature

MD𝑃$%&''

Perform MD at 𝑇! (𝑚 = 1, … , 𝑀)

2.

Change 𝑇! to 𝑇!±# each fixed MD stepwith following probability, w.

Conformation

P(E)

1.

𝑤 𝑇! → 𝑇!±# = A

𝑃$%&''

if ∆≤ 0if ∆> 0∆= 𝛽!±# − 𝛽! 𝐸 − 𝛼!±# − 𝛼! , 𝛽! = 1H𝑘 𝑇

ST-MDP(E)1,exp −∆ ,4 !

𝑃$%()*+𝑷𝑺𝑻• It is difficult to define hyperparameter, 𝛼! ,for large/complicated system.
A. Irbäck, et al., J. Chem. Phys. 103, 10298 (1995)
U.H.E. Hansmann, et al., J. Comput. Chem. 18, 920 (1997)

Potential energy

Potential energy

9Temperature replica-exchange MD (T-REMD)
Replica 1300 K

310 K

Replica 2310 K

300 K

Replica 3

320 K

330 K

Replica 4

330 K

320 K

Time
P(E)

change → exchange• T-REMD is a parallelization model of simulated
tempering (ST) MD (= parallel tempering MD).
• T-REMD consists of M non-interacting copies (𝑖 =1, … , 𝑀) at M different temperatures, 𝑇! (𝑚 = 1, … , 𝑀).
A𝑖 = 𝑖(𝑚) ≡ 𝑓(𝑚)𝑚 = 𝑚(𝑖) ≡ 𝑓 5# (𝑖)f: permutation function• The weight factor for the state X, W(X), is given by
8𝑊678 𝑥 = exp − N 𝛽! 𝐻(𝑞 ) ! , 𝑝 ) ! )!%#

H: Hamiltonian

Potential energy

K. Hukushima, et al., J. Phys. Soc. Jpn. 65, 1604 (1996)
Y. Sugita, et al., Chem. Phys. Lett. 314, 141 (1999).
10Temperature replica-exchange MD (T-REMD)
Replica 1300 K
310 K

Replica 2310 K
300 K

RejectReplica 3320 K

Replica 4P(E)Accept320 K330 K330 KTimeAccept
𝐸 𝑥! − 𝐸 𝑥 ">0𝑥9
𝐸(𝑥 9 )𝑥)𝑥)

𝐸(𝑥 ) )𝐸(𝑥 ) )𝑥9
𝐸(𝑥 9 )and j, corresponding to neighboring temperatures,
𝑇! , and 𝑇. (𝑇! < 𝑇. ), are exchanged.
𝑋=9)… , 𝑥! , … , 𝑥.

≪0→𝑋

:

[9]!
[)]!
= (… , 𝑥! , … , 𝑥. )• According to the Metropolis criteria, satisfying thedetailed balance, the replica-exchange is “accepted”
or “rejected”.
𝑤 𝑋 → 𝑋: = A

Reject𝐸 𝑥! − 𝐸 𝑥 "• Every fixed number of MD steps, a pair of replicas, i1,exp −∆ ,if ∆≤ 0if ∆> 0∆= 𝛽. − 𝛽! (𝐸(𝑥 ) ) − 𝐸(𝑥 9 ))• Succeeding replica exchanges must have overlapping
distributions.

11Estimating unbiased physical quantities
PickupReplica 1• Physical quantities are obtained by following

300 K

310 K

Replica 2

310 K

300 K

1. Collect data at “300 K” from each replica

Replica 3

320 K

330 K

2. Estimate time averages at 300 K by following

procedure.

?"#$%&'

Replica 4

330 K

320 K

𝐴 =)!> =N 𝐴 $%&''@ (𝑡)=%#TimeP(E)• It means that we waste all the data except for thedata at T=300 K.

Wasteful!!
Potential energy

12Multistate Bennett Acceptance Ratio (MBAR)
M. R. Shirts, et al., J. Chem. Phys. 129, 124105 (2008)• MBAR can estimate dimensionless free energy differences by the following equation.
8?(𝑓V) = ln N N
9%# .%#exp −𝛽) 𝑈(𝑥9. )V∑8F%# 𝑁F exp 𝑓F − 𝛽F 𝑈(𝑥9. )• Once 𝑓V) is obtained by using MBAR equation, average of physical quantities at target
temperature, 𝐴 ABCDEA , can be derived as follows,
8?(𝑓\ABCDEA = ln N N

exp −𝑈ABCDEA (𝑥9. )V∑8F%# 𝑁F exp 𝑓F − 𝑈F (𝑥9. )9%# .%#8 ?(𝑓\ABCDEA, H = ln N N
9%# .%#𝐴 ABCDEA =A(𝑥9. )exp −𝑈ABCDEA (𝑥9. )V∑8F%# 𝑁F exp 𝑓F − 𝑈F (𝑥9. )exp −𝑓ABCDEA, H
exp −𝑓ABCDEA

exp −𝑓\ABCDEA, H~
exp −𝑓\ABCDEA

fi/k: Dimensionless free energy oftemperature label i/k
M : Total # of replicas
Nj/k: Total # of snapshot oftemperature label j/k
U : Potential energy
x : Coordinate
A : Physical quantities

13Multistate Bennett Acceptance Ratio (MBAR)• 𝑓V) is obtained by a simple self-consistent iteration as follows.
8?)𝑓V) = ln N N
F%# .%#exp −𝛽) 𝑈(𝑥F. )∑8I# 𝑁0 exp 𝑓V0 − 𝛽0 𝑈(𝑥F. )0%#unknow value

1.

Set initial (t=0) 𝑓\)=%' to 0.0 from i=1 to M

2.

Calculate each 𝑓\)=I# by using MBAR equation

3.

Check ∆𝑓\) = 𝑓\)=I# − 𝑓\)=
• If ∆𝑓\) > threshold (~1.0×105J ), go back step 2• If ∆𝑓\) ≤ threshold, go to step 4

4.

Estimate weight factors at each snapshot and temperature
14

Estimating free energy differences• Weight factor, 𝑤F. , is obtained as follows.
𝐴 ABCDEA =exp −𝑓ABCDEA, H
exp −𝑓ABCDEA

8 ?)exp −𝑓\ABCDEA, H~= N N 𝑤F. 𝐴(𝑥F. )\exp −𝑓ABCDEA
F%# .

• Since 𝑤F. is independent of A, it is available for estimating other quantities, B, C, etc.
• We can estimate the free energy differences along a reaction coordinate, 𝜉, named “potential ofmean force (PMF)”, 𝐹(𝜉), using 𝑤F. .
𝑃=-K*>= 𝑥F. = 𝑤F. 𝑃L)-M>N (𝑥F. )𝐹 𝜉 = −𝑘4 𝑇=-K*>= ln 𝑃=-K*>= 𝜉 = −𝑘4 𝑇=-K*>= ln 𝑤𝑃L)-M>N 𝜉15

Appendix

16

Potential smoothingEnergy• Potential smoothing is useful if we already know which

Biassing potential

collective variables (CVs) are important.
• CVs are quantities that describe the behavior or stateof a system.
1.

Approach 1: Decrease the potential energy barrier• Multicanonical

Biassing potential• Umbrella sampling

Original
2.

Collective variable

Approach 2: Fill the free energy surface• Metadynamics

17