2014/5/21
2014年バイオスーパーコンピューティング・ソフト講習会シリーズ「CafeMol」粗視化分子モデル
計算ソフト講習会

検崎博生
理化学研究所情報基盤センター

CafeMol (www.cafemol.org)• Features are;– Various CG models– protein/DNA/RNA– multiple basin model– accurate CG model– Simulating protein-at-work– “switching”• Under development
lipid• Developer
検崎博生、古賀信康、藤原慎司、堀直人、金田亮、李 文飛、岡崎圭一、姚 新
秋、寺川剛、伊藤真志保、高田彰二
京都大学高田研究室

CafeMol 2.1 (2013/7)
source & manual released

Overview of CafeMol• General-purpose coarse-grained (CG)biomolecular modeling and simulation software– Protein: 1 bead / 1 amino acid
– Nucleic acid:
3 beads (sugar, base, phosphate) / nucleotide– Lipid: ~3 beads / lipid• Written by FORTRAN90 with MPI and Open MP• Large-scale simulation– ～”millisecond” event by K-computer• Version 1.0 is released (only protein) (2010/12/27)– Version 2.0 (protein, DNA, RNA) (2012/5/31)– Version 2.1 (protein, DNA, RNA) (2013/7/1)

Menu

1. Molecular dynamics simulation
2. Coarse-grained models
3. Simulation methods & Implementation
4. Selected applications

Menu

1. Molecular dynamics simulation
2. Coarse-grained models
3. Simulation methods & Implementation
4. Selected applications

Biomolecular simulations• Molecular orbital(MO) method– quantum mechanics• chemical reaction• Molecular dynamics(MD) method– classical mechanics• conformational change of biomolecules• Continuum method– elastic model• muscle contraction– fｌuid dynamics• blood flow

Space & time in protein dynamics

MO

10000

Mol motor
Allostery

1000
100

Protein Foldingin Hmgb
Conf. Relax
in Hmgbα

Reactiondynamics
Heme
Relax

- 15
fs

Folding Protein G

βHairpin
Form

Helix
Form

CG MD

AA MD

10

#atoms (excluding solvent)

F1

-12
ps

-9
ns

-6
μs

time

-3
ms

0
s

3
h

log10tMolecular dynamics• Calculate molecular movementm– Classical mechanics– Numerically integrate Newton’s equation of motion– m(mass) x a(acceleration) = F(force)– interaction parameters are derived from MO• Interaction– bond length, bond angle, dihedral angle– van der Waals interactiondihedral angle– electrostatic interaction
bond angle
bond lengthFConstant temperature dynamics• Berendsen thermostat– rescale velocity to control temperature– not canonical distribution• Nose-Hoover thermostat– rescale velocity using artificial particles– canonical distribution• Langevin dynamics– random force and viscosity– fluctuation-dissipation theorem– canonical distribution

Generalized ensemble methods• Generalized ensemble methods replica1– accelerate conformational sampling replica2
replica3– only obtain equilibrium variables• Replica exchange molecular dynamics(REMD)– simulate many copies at different temperatures– exchange configuration at proper probability– extension to interaction parameters• Multicanonical method– sample “flat histogram” on energy– iterative learning of “the density of states”
canonical
ensembleE

multicanonicalensembleEMenu

1. Molecular dynamics simulation
2. Coarse-grained models
3. Simulation methods & Implementation
4. Selected applications

Why coarse-grained model ?
Reducing the number of particles to 1/100• 1 amino acid residue →1 particle (Caatom)• No water molecule(included in potential energy)
Protein folding andfunctional conformational changecoarsegrained

10^6 times speed up• Reducing the number of particles• Enlarge time step
• Low friction

AA: 100PC 1week

CG: 1PC 1minute!!

Models and energy functionsCG
1 beads / 1 amino acid

Atomica.a.

A. Go-like model
B. Atomic interaction based CG(AICG) model
C. Flexible local potential(FLP) model
D. AICG2+ model
E. Multiple basin model
F. DNA/RNA model
G. Elastic network model
H. Electrostatic and hydrophobic interactions
I. Explicit and implicit ligands

Go-like model
C. Clementi, H. Nymeyer, and J.N. Onuchic, J. Mol. Biol. (2000)

Based on the energy landscape theory
Structure based

Vprotein  Vlocal  Vgo  Vex






 K  1  cos      K  1  cos 3   

Vlocal  K b  ri ,i1  r0 i ,i 1  K    i   0 i
2i2

: bond angle
: dihedral angle(0 means native state)

i

1

3



i

0i

i



i

i

10
  r 12

r0 ij 0ij
  6 
Vgo  go  5 rij r  
i, j 

 ij  nativeVex   ex

nonnative


i, j

12

 
 
r 
 ij 

0i

Kb = 100
K = 20
K1 = K3 = 0.5go  0.18
ex  
= 4A
=1.0kcal/mol

Atomic interaction based CG (AICG) model
W. Li, P. Wolynes, S. Takada, PNAS (2011)V   k bi ( r i  r0i )CG  k ai ( i   0i )

2

i

2

i

  {  ,1[1  cos(  0i )]   i ,3 [1  cos 3( i  0i )]}iii

Atomica.a.

native   [ 5( r / r )  6 ( r / r ) ] 
ij

ij
0

ij 12

ij
0

non  native

ij 10

i  j 3  (C / r )ij 12i  j 3

Wenfei Li

1) Contact energy εij from pairwise all-atom (AA) energy
E IJ ( RIJ )   u AA ( rij )SAu AA ( r )  V ( r )  G GB
(
r
)


G
(r)
pol

iI jJ

2) Coefficients fitted by AA-derived fluctuation (23 proteins)
param

Kb

KaG

kaH

kaE

kaT

kaC

G

H

E

T

C

nloc

Av.

109.94

13.40

40.0
3

17.3
2

19.35

11.7
0

0.29

1.76

1.32

0.82

0.81

0.37

Test for fluctuation,
structural change, &folding

Protein G

0.694
0.758
Protein L

Flexible local potential(FLP) model
T. Terakawa, and S. Takada, Biophys. J. (2011)

Probability distributions to CG potential
for flexible local potentials


Cα

Cα

Cα

Cα

a

Pa (a )
Pd ( )

Loop segment library

Boltzmann inversion; prob. to pot.

ln Pa (a )Va   k B T
sin a

Vd   k BT ln Pd ( )

AICG2+ model
W. Li, T. Terakawa, W. Wang, S. Takada, PNAS (2012)

improved AICG model + FLP modelV   kb ( r  r )
I

I 2
0

I
I
  VaI ( I )   Vdih( I )
I

I

IJ
IJ 2
I
I 2
(
r

r
)
(



IJ
IJ
0
0 )
   1，3 exp( )+   1，4 exp( 
)+
2
2
2w
2 w
J =I + 2
J =I +3
nativenon  native

I  J 3I  J 3IJ
   nloc
[5( r0IJ / r IJ )12  6( r0IJ / r IJ )10 ]  (C / r IJ )12

Multiple-basin model for proteinsK. Okazaki, N. Koga, S. Takada, J.N. Onuchic, and P.G. Wolynes, PNAS (2006)

Use of 2 references V ( R | R1)
  c1  c1 
   VMB  V ( R | R 2 )   V  c 2  c2 V ( R | R1)  V ( R | R 2 )  V
VMB 2
V ( R | R1)  V ( R | R 2)  V
 
 
2
 c2   log  c1 

CG DNA model
T.A. Knotts IV, N.Rathore, D.C. Shwartz, and J.J. Pablo, J. Chem. Phys. (2007)•Three interactions sites•Phosphate•Sugar•Base•Reproduce various DNA
behavior•Salt-dependent
melting•Bubble formation•Mechanicalproperties

3SPN.1 force field
E.J. Sambrisiki, D.C. Schwartz, and J.J. de Pablo,Knotts, Biophys. J. (2009)

Vdna  Vlocal  Vstack  Vbp  Vex  Vqq  Vsolv







Vlocal  K b1  ri ,i1  r0 i ,i 1  K b 2  ri ,i1  r0 i ,i 1i K   i   0 i
i

2



i

4

  K  1  cos   
2



i

0i

i

  12    6 
0 ij 
0 ij 

Vstack  4 1  
 rij  r  i, j 

 ij  : bond angle
: dihedral angle(0 means B-type DNA)

N st

Structre-based interaction10
  r 12
r
0 ij 
0 ij 
Vbp   4bpi 5 6
r  
  rij i, j ij  Nbp 12   6 Vex  4 1   0    0     1 (if rij  dcut ), r  
i , j  rij  ij  
N ex 0 (if rij  dcut )Kb1 = 1Kb2 = 100K = 1400
K = 28
bpGC  2.532
bpAT  2.0
 =0.1839kcal/mol

3SPN.1 force field(electrostatic and solvation interaction)  rij / D
q iq jeVqq  
i , j  4 0 (T , C ) rij  (T , C )   (T ) a (C )N

Debye-Huckel theory

=78

 (T )  249.4  0.788T / K  7.20  10  4 (T / k ) 2

Debye length  0RT D  
22
 2N e I 
A q a (C )  1.000  0.2551C / M
 5.151  10  2 (C / M ) 2  6.889  10 3 (C / M ) 3
Nsolv



Vsolv    s 1  e  a ( rij  rs ) 2

s

i j

s   N AI
e N  e0 (1  [1.40418  0.268231N nt ]1 )AI  0.474876 (1  {0.148378  10.9553[ Na  ]}1 )a-1 = 5.333A
rs = 13.38A
0  0.504982

CG RNA model
N. Hori, and S. Takada, JCTC (2012)

1 nucleotide = 3 beads
Phosphate (phosphorus atom)
Sugar (center of ribose ring)
Base (pyrimidine: N1 atom)(purine: N3 atom)

CG RNA model (local)

Bond length

Local

Bond angle

Dihedral angle
5’

5’
B(Y)

B(Y)

5’

S

B(Y)
P

B(Y)

B(Y)

Base type
Y : Pyrimidine (Uracil, Cytosine)
R : Purine (Adenine, Guanine)

B(Y)
S
P
S

3’

P

S

B(R)

B(R)

S

P

B(R)
S

3’

3’

CG RNA model (nonlocal)
Nonlocal

Reference (X-ray structure’s) value

Base pairBase pairBase stack

Contact（Other）
Excluded volume

Base stackElastic network modelEK ( rij  rij0 ) 2

ij , s .t .rij0  rc

Atomic fluctuation
reproduced by ENM(Tirion96)

Electrostatic and hydrophobic interactions
Debye-Huckel form for electrostatics
Debye length  0RT D  
 2 N 2e2 I 
A q HP interactions analogous to ASACount coordination number
for each hydrophobic particle

Explicit and Implicit ligands
Explicit ligand: as a rigid molecule

Implicit ligand:
MD-MC scheme with ligand-mediated contact

V protein

k on
koff

V protein  Vimp lig

Menu

1. Molecular dynamics simulation
2. Coarse-grained models
3. Simulation methods & Implementation
4. Selected applications

Simulation method• Dynamics– Newtonian dynamics with Berendsen thermostat– Langevin dynamics– Multi-Particle Collision dynamics (MPC)• Time integration– velocity Verlet algorithm• Run mode
– Constant temperature simulation– Simulated annealing– Auto-search of Tf– Replica exchange method– Potential “switching”• Useful option– anchor, bridge, pulling, fix, box

Units in CafeMol• The length unit is A
– 1A = 10-10m• The energy unit is kcal/mol– 1kcal/mol = 6.9478pNnm = 0.04337eV– 300kBT = 0.6kcal/mol• The mass unit is cafemol-mass-unit (cafemu)– each amino acid has the mass of 10cafemu– 1cafemu = 13.7amu = 2.275x10-26kg• The unit of charge is elementary-electric charge (e)• The unit time is café time
– 1cafe time= 1.809x10-13s ~ 200fs– intrinsic dynamics is accelerated by coarse-graining
the energy landscapeInput and output files• Input files– Input file(.inp)– PDB file(.pdb)– native-info file(.ninfo)– parameter file(.para)– AICG files• Output files– data file(.data)– time-series file(.ts)– native-info file(.ninfo)– coordinate and velocity• PDB format(.pdb)• CARD format(.crd, .vdcd)– trajectory

red: essential files
green: frequently used
black: optional files• PDB format(.movie)• DCD coordinate(.dcd, .vdcd)– PSF file(.psf)– restart file(.rst)– replica information file(.rep)

Example of input file
(folding simulation of src SH3)
interaction

Input
block

input/output<<<< filenames
path = ./data<<<< energy_function
filename = sh3
LOCAL(1) L_GO
OUTPUT psf pdb dcd rst
Go modelNLOCAL(1/1) GO EXV
path_pdb = ./pdb>>>>
path_ini = ./pdb<<<< md_information
path_natinfo = ./ninfo
simulation type n_step_sim = 1>>>>
n_tstep(1) = 300000<<<< job_cntl
tstep_size = 0.2
i_run_mode = 2
Const T MD
n_step_save = 100
i_simulate_type = 1 Langevinn_step_neighbor = 100
From random
i_initial_state = 1
tempk = 300.0>>>>
n_seed = 1<<<< unit_and_state>>>>i_seq_read_style = 1
i_go_native_read_style = 1
detail/option
1 protein 1SRL.pdb

sequence/structure

Native-info

one-by-one styleall-in-one style<<<< native_info_sim1
NINFO(all/all) f1atp_all.ninfo>>>><<<< native_info_sim1
NINFO(1/1) 1 Intra-mol 1…
NINFO(3/6) 13 Inter-mol 3/6
1= f1atp_alpha_E.ninfo…
13= f1atp_alphDP_betaDP.ninfo>>>>

native-info file (alpha_E subunit)
bond 1 1 1 1 2 1 2 3.8132 1.0000 1.0000 100.0000…
angl 1 1 1 1 2 3 1 2 3 93.2170 1.0000 1.0000 20.0000…
dihd 1 1 1 1 2 3 4 1 2 3 4 67.0855 1.0000 1.0000 1.0000…
contact 1 1 1 1 5 1 5 5.9973 1.0000 1 0.3000…

pair ij

nat- distancecoefficient contact interactionUseful option1
 redefine parameters
dfcontact = 4.5
cdist_rep12 = 6.0
rneighbor_dist = 20.0
fric_const = 0.02

definition of native contact(default 6.5A)
reference distance in excluded interaction(default 4.0A)
truncation distance for neighbor list(default 24.0A)
friction constant(default 0.25)

 delete interaction
DEL_LGO(mpini-mplas)
DEL_GO(mpini,i-mplas,i/mpini,j-mplas,j) box: fix some particles or units
xbox = boxsize_x
ybox = boxsize_y
zbox = boxsize_z
boxsigma = sigma( d  3 ),Vbox  0 
 k box  d

12(0.5  d  3 ),0.5  d  
 k box ) ( d  0.5 ) (1  120.5
 0.5 
12

Useful option2
 anchor: particle i constrain to some position
ANCH i ki l0 x0 y0 z0

Vanchor  k i ( ri 0  l0 ) 2 ( ri 0  l 0 ),
 0 ( ri 0  l 0 ) bridge: bind I and j particles by a harmonic spring
Vbridge  k ij ( rij  l0 ) 2 ( rij  l 0 ),BRIDGE i j kij l0 0 ( rij  l 0 ) pulling: particle I is pulled by constant force or constant velocity
PULL_CF i fx fy fz
PULL_CV i kd vx vy vz x0 y0 z0

 fix: fix some particles or units
FIX_UNIT(unitini-unitlas)
FIX_MP(mpini-mplas)(fx, fy, fz)
kd(vx, vy, vz)

CafeMol codeCafeMol code parallization• Time integral(MPI+OpenMP)– neighboring list
• each node calculates distances between assigned pairs andmakes neighboring list

– force, energy• each node calculates force or energy including theirneighboring list
• these calculations are parallelized by OpenMP• Temperature/Hamiltonian REMD(MPI)– replica_1, replica_2, …. , replica_nPerformance of parallelization
at K Computer
Efficiency of parallelization:
99.996%

Before 2010
RICC at Riken (MPI+Open MP)
After 2011
at K-computer (SIMD)#node ~73,000 ⇒～1PFLOPS !?

MFLOPS

1.00E+09

20 nucleosomes(35918 particles)

1.00E+08

1.00E+07

MFLOPS

1.00E+06

1.00E+05

#Node

12

12

24

24

48

48

96

192

384

768

1536

3072

6144

96 192 384 768 1536 3072 6144

実行効率(%) 16.33 15.99 16.02 15.89 15.11 15.58 15.67 15.57 15.58 13.14

Menu

1. Molecular dynamics simulation
2. Coarse-grained models
3. Simulation methods & Implementation
4. Selected applications

Folding simulation of src SH3 domain
N. Koga, and S. Takada, J. Mol. Biol. (2001)

computed values experimental values

Native fluctuation by Go-like model
Src SH3 domain

Red: Go-like model(CG)Green: AMBER GB/SA(AA)
Blue: AMBER+ TIP3P(AA)
Pink: X-ray B-factor(Experiment)

Folding temperature of src SH3(Auto-search of Tf)
Bi-section method<<<< job_cntl
i_run_mode = 4
i_simulate_type = 1
i_initial_state = 1>>>><<<< searching_tf
tempk_upper = 500.0
tempk_lower = 100.0>>>>*********************************
tf_out tempk n_state d_state p_transtf_out 300.000995
51*********************************
tf_out tempk n_state d_state p_transtf_out 400.000
1 10000*********************************
tf_out tempk n_state d_state p_transtf_out 350.000
166
835
78
*********************************
tf_out tempk n_state d_state p_transtf_out 325.000
953
48
19
*********************************…
…
*********************************
tf_out tempk n_state d_state p_transtf_out 341.406
638
363
98
*********************************

Folding temperature of some proteinsProtein

Number ofamino acid

Folding
temperature(K)albumin binding
domain53

380.4

src SH3 domain

56

342.9

protein G

56

338.2

a-spectrin SH3
domain

57

360.1

Sso7d

64

332.0

protein L

78

374.2

Im9

86

382.0

cytochrom B562

106

352.2

“Switching” simulation<<<< native_info_sim1
NINFO(1/1) 1 Intra-mol 1…
switching
NINFO(3/6) 13 Inter-mol 3/6
1= f1atp_alpha_E.ninfo…
13= f1atp_alphDP_betaDP.ninfo>>>><<<< native_info_sim2
NINFO(1/1) 1…
NINFO(3/6) 13
1= f1atp_alpha_TP.ninfo…
13= f1atp_alphE_betaE.ninfo>>>>

native-info file (alpha_E subunit)
bond 1 1 1 1 2 1 2 3.8132 1.0000 1.0000 100.0000…
angl 1 1 1 1 2 3 1 2 3 93.2170 1.0000 1.0000 20.0000…
dihd 1 1 1 1 2 3 4 1 2 3 4 67.0855 1.0000 1.0000 1.0000…
contact 1 1 1 1 5 1 5 5.9973 1.0000 1 0.3000…

pair ij

nat- distancecoefficient contact interaction

Rotation mechanism of
F1-ATPase by switching Go model
N. Koga, and S. Takada, PNAS (2006)

Conformational change by
K. Okazaki, N. Koga, S. Takada, J.N. Onuchic,
Multi-basin model and P.G. Wolynes, PNAS (2006)<<<< unit_and_statei_seq_read_style = 1
i_go_native_read_style = 1
1a protein 1GGG_2.pdb
1b protein 1WDN_2.pdb>>>><<<< energy_function
NLOCAL(1a/1a) GO EXVNLOCAL(1b/1b) GO EXVMULTIGO_SYSTEM(1a) 1a/1a
MULTIGO_SYSTEM(1b) 1b/1b>>>><<<< multiple_gobdemax_mgo = 100.0
baemax_mgo = 1.0
dihemax_mgo = 0.5
ENEGAP(1)(1) 0.0 -1.8
DELTA(1ab) 28.0>>>>

Conformational change using
3-state multi-basin model
X. Yao, H. Kenzaki, S. Murakami, and S. Takada, Nature Comm. (2010)

Multidrug transporter
AcrB•Largely responsible formultidrug resistance•Asymmetric homo-trimer
structure

Suggesting mechanismof drug exportation andfunctional rotation!

Sliding movement of KIF1A
R. Kanada, T. Kuwata, H.Kenzaki, S.Takada, PLOS Comput. Biol. (2013)

1 phase: multiple-basin (T, D)
2 phase: go(D )3 phase: multiple-basin(D, phi)
4 phase: go(phi)5 phase: go(T)KIF1A:blue
tubulin:green
cargo:yellow

CG DNA simulation••••

30 bp DNA duplex
Langevin dynamics (300K)[Na+] = 69mM
cutoff length 20D(Debye length)<<<< unit_and_statei_seq_read_style = 2
i_go_native_read_style = 3
1-2 dnasequence>>>><<<< energy_functionDH
LOCAL(1-2) L_BDNA
NLOCAL(1-2/1-2) DNA ELE>>>>
3SPN.1<<<< electrostatic Intra mol 1,2
cutoff_ele = 20.0 Inter mol 1-2
ionic_strength = 0.069
diele_water = 78.0>>>><<<< in_box
xbox = 120.0
ybox = 120.0
zbox = 120.0
boxsigma = 4.0>>>>

Diffusional search of p53 on DNA
T. Terakawa, H. Kenzaki, and S. Takada, JACS, (2012)•Tumor suppressor•Transcription factor

1. Nonspecific diffusional
search on DNA
2. Specific binding

protein:AICG and FLP model
DNA:CG DNA model
protein-DNA:
excluded and electrostatic interaction

Core domain hopping
C-terminal domain sliding

Simulation of nucleosome• Electrostatic interaction+ Go potential

1KX5.pdb

go pro-dna  0.8go pro[Na+] = 200mM
go pro-dna:coefficient ofprotein-DNA Go potential<<<< energy_function
LOCAL(1-2) L_BDNALOCAL(3-10) L_GO
NLOCAL(1-2/1-2) ELE DNANLOCAL(1-2/3-10) GO EXV ELENLOCAL(3-10/3-10) GO EXV>>>><<<< electrostatic
cutoff_ele = 5.0
ionic_strength = 0.2
diele_water = 78.0>>>>
H. Kenzaki, et al unpublished data

Acknowledgement
CafeMol development has been supported by Researchand Development of the Next-Generation Integrated
Simulation of Living Matter, a part of the Development
and Use of the Next-Generation Supercomputer Projectof the Ministry of Education, Culture, Sports, Scienceand Technology.