A coarse-grained simulator:
CafeMol
検崎博生1,古賀信康1, 藤原慎司1, 堀直人1, 金田亮1,
1 2, 岡崎圭一
12
李 文飛1,2
岡崎圭 1, 姚 新秋1, 高田彰二1,2
1京都大学理学研究科生物物理学教室
2 JST-CREST
JST CREST

CafeMol (www.cafemol.org)
Kenzaki• CafeMol 2.1 (2013/7)
source & manual released• Features are;– Various CG models
protein/DNA/RNA
multiple basin model
accurate CG model– Simulating protein-at-work“switching”• Under development
lipidOverview of CafeMol• General-purposep pcoarse-grainedg((CG)) biomolecularmodeling and simulation software– Protein: 1 bead / 1 amino acid
– Nucleic acid:
3 beads ((sugar,g base, phosphate)) / nucleotide– Lipid: ~3 beads / lipid• Written byy FORTRAN90 with MPI and Openp MP
• Large-scale simulation– ～”millisecond” event byy K-computerp• Version 1.0 is released (only protein) (2010/12/27)– Version 2.0 (protein, DNA, RNA) (2012/5/31)

Menu
ModelsSimulation methods
Implementation
Selected applications
I
In-progressmodelsd l & methods
th dMenu
ModelsSimulation methods
Implementation
Selected applications
I
In-progressmodelsd l & methodsth dModels and energy
gy functions
CG1 beads / 1 amino acid

Atomicaaa.a.

A. Off-lattice Go model
B Atomic interaction based CG modelB.
C. Multiple basin model
D DNA/RNA modelD.
E. Elastic network model
F Electrostatic and hydrophobic interactions
F.
G. Explicit and implicit ligands

Off-lattice Go 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 i2i2

θ: bond angle
φ: dihedral angle(0 means native state)

i

13

φ0iii

φii

⎡ ⎛ r ⎞12 ⎛ r ⎞10 ⎤0 ij0 ij
Vgo = εgo ∑ ⎢5⎜ ⎟ − 6⎜ ⎟ ⎥⎜ rij ⎟⎜r ⎟ ⎥i, j ⎢⎝⎠⎝ ij ⎠ ⎦⎣nativeVex = ε ex

nonnative

∑
i, j

12

⎛σ ⎞
⎜ ⎟⎜r ⎟⎝ ij ⎠0iKb = 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,W.
Li P.
P Wolynes,
W lS.
S Takada,
T k d PNAS (2011)
V = ∑ kbi (r i − r0i )CG+ ∑ k ai (θ i − θ 0i )2i2i+ ∑ {ε φi ,1[1 − cos(φ i − φ0i )] + ε φi ,3 [1 − cos 3(φ i − φ0i )]}i

Atomica.a.

native+ ∑ ε [5(r / r ) − 6(r / r ) ] +ij

ij0i > j −3ij 12ij
0non − nativeij 10∑ ε (C / r )ij 12i > j −3

Wenfei Li1) Contact energy εij from pairwise all
all-atomatom (AA) energy
E IJ ( RIJ ) = ∑∑ u AA ( rij )SAu AA ( r ) = V ( r ) + ΔG GB
pol ( r ) + ΔG ( r )i∈I j∈J2) Coefficients
C ffi i t fitted
fitt d bby AA-derived
AA d i d fluctuation
fl t ti (23 proteins)
t i )
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
00.29

1.76

1.32

0.820.81

0.37

GNMMSM

0.4

Protein G

0.694
0.758

0.2

Protein L
C.C.

00

1.

∼

∼

0. 0.995 5∼

0. 0.885 5∼

00. .75
75∼

0. 0.665 5∼

00. .55
55

0.0<00. .4
45 5

Proobability

Test for fluctuation,
structural changechange, &folding

0.6

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 =2V ( R | R1) − V ( R | R 2) − ΔV− (+Δ )2⎛ c2 ⎞
χ = log⎜ ⎟⎝ c1 ⎠

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

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) + K ∑ (r − r )V = K ∑ (r − r2localb1i ,i +1i(+ Kθ ∑ θ i − θ 0 ii

40 i ,i +1b20 i ,i +1i ,i +1i) + K ∑ (1 − cos(φ − φ ))2φi0ii⎡⎛ σ ⎞12 ⎛ σ ⎞ 6 ⎤0 ij ⎟0 ij ⎟⎥−⎜Vstack = 4ε 1 ∑ ⎢⎜⎜⎟⎜ r ⎟ ⎥i , j ⎢ rij⎝⎠⎝ ij ⎠ ⎦⎣θ: bond angleφ: dihedral angle(0 means B-type DNA)

Nst

Structre-based interaction⎡ ⎛ r ⎞12 ⎛ r ⎞10 ⎤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 < dcutt ),)
⎜
⎟
⎟
⎜
⎥
⎢
i , j ⎝ rij ⎠⎝ rij ⎠ ⎦⎣Nex

= 0 (if rij > dcut )Kb1 = 1εKb2 = 100εKθ = 1400εKφ = 28εεbpGC = 2.532εεbpAT = 2.0ε20ε =0.1839kcal/mol

3S3SPN.1forceo ce fielded(electrostatic and solvation interaction)⎛⎞ − rij / κDq iq j⎜⎟eVqq = ∑⎜⎟i , j ⎝ 4πε 0ε (T , C ) rij ⎠ε (T , C ) = ε (T )a (C )N

Debye-Huckel theory

ε (T ) = 249.4 − 0.788T / K + 7.20 ×10 − 4 (T / k ) 2a (C ) = 1.000 − 0.2551C / M+ 5.151×10 − 2 (C / M ) 2 − 6.889 ×10 −3 (C / M ) 3Nsolv[Vsolv = ∑ ε s 1 − ei< j] −ε− a ( rij − rs ) 2sεs = ε N AI
eN = e0 (1 − [1.40418 − 0.268231N nt ]−1 )AI = 0.474876(1 + {0.148378 + 10.9553[ Na + ]}−1 )α-1 = 5.333A
rs = 13.38A
ε0 = 0.504982εRNA model

N. Hori and S. Takada (2012)

1 nucleotide = 3 beads

Phosphate
Phh t ((phosphorushhatom)t )Sugar (center of ribose ring)
Base (pyrimidine: N1 atom)(purine: N3 atom)

RNA model (local)()
Bond lengthLocal

Bond angle

Dihedral angle5’

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

3’PS

B(R)

B(R)

SS
P

B(R)
S

3’3’

RNA model ((nonlocal))
Nonlocal

Reference (X-ray structure’s) value

Base pairBase pairBase stack
Contact（Other）
E l d d volume
Excludedl

Base stackElastic network model∑0 2
E=K (rij − rij )ij , s.t .rij0 < rc

Atomic fluctuation
reproduced by ENM(Tirion96)

Electrostatic and hydrophobic interactionsDebye-Huckel form for electrostatics

HP interactions analogousgto ASACount coordination numberfor each hydrophobic particle

Explicit and Implicit ligands
Explicit ligand; as a rigid molecule

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

V protein

konkoff

V protein + Vimp −ligMenuModelsSimulation methods
Implementation
Selected applications
I
In-progressmodelsd l & methods
th dSimulation method• DDynamics
i– 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Menu
ModelsSimulation methods
Implementation
Selected applications
IIn-progressmodelsd l & methods
th dCafeMol code

• Parallelization– neighboring list, force, energy→hybrid(MPI+Open
hybrid(MPI+Open MP)– replica exchange→MPI(temperature/Hamiltonian REMD)

Performance of MPI parallelization
1300 base pairs DNA (7798 particles)
BG/L at RikRiken

speeed up

100

10

18speed up 7.85
ope(s)
mpi(s)

163264128
15.38 29.64 54.76 94.65

328.62 164.64
5.44
5.90

82.65
5.84

42.18
5.71

22.01
5.70

High parallelization efficiency

Performance of parallelization at K-Computer
BeforeBf2010RICC at Riken (MPI+Open MP)
After 2011
at K-computer (SIMD)

Efficiency of parallelization
parallelization:
99.996%#node#d ~73,000
73 000 ⇒～1PFLOPS !?

MFLOPS

20 nucleosomes(35918 particles)

1.00E+09

1.00E+08

1 00E+07
1.00E+07

MFLOPS

1 00E+06
1.00E+06

1.00E+05#Node12244896192

384768153630726144

12244896 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
ModelsSimulation methods
Implementation
Selected applications
I
In-progressmodelsd l & methods
th dNative fluctuation by off-lattice Go model
Src SH3 domain

Red; Go modelGreen; AMBER GB/SABlue; AMBER+ TIP3P
Pink; X-ray B-factor

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

computed φ−values experimental φ−values

Example of input file(f ldi simulation(foldingil ti off src SH3)

Input
block

Input/output<<<< filenames<<<< energy_function
Go model NLOCAL(1/1) GO EXV
path = ./data
filename = sh3
i use atom protein = 0
i_use_atom_protein
Interaction
OUTPUT pdb dcd
i_use_atom_dna = 0
path_pdb = ./pdb>>>><<<< md_information
path_ini = ./pdb
pppath_natinfo = ./ninfo
n_step_sim = 1>>>>
n_tstep(1) = 300000<<<< job_cntl
tstep_size = 0.2
Simulation typei run mode = 2
i_run_mode
n step save = 100
n_step_save
Const T MD
i_simulate_type = 1
n_step_neighbor = 100
Langevin
i_initial_state = 1
i_com_zeroing = 0
From random>>>>i_no_trans_rot = 0<<<< unit_and_state
tempk = 300.0
D t il/ ti
Detail/option
i_seq_read_style = 1
n_seed = 1
i_go_native_read_style = 1>>>>
1 protein
1SRL.pdb
1SRL pdb>>>>Sequence/structureFoldingo d g temperature
te pe atu e ofo src
s c SH3
S 3(Auto-search of Tf)
Bi-section
Bisection method<<<< job_cntl
job cntl
i_run_mode = 4
i_simulate_type = 1
i_initial_state = 1>>>><<<< searching_tf
tempk_upper = 500.0
tempk lower = 100.0
tempk_lower>>>>*********************************
tf_out tempk n_state d_state p_transtf_out 300.000
99551*********************************
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
16683578*********************************
tftf_outt ttempkk n_statet t dd_statet t p_transttf_out 325.000953
4819*********************************…
…
*********************************
tf_out tempk n_state d_state p_transtf out 341.406
tf_out341 406
638
363
98*********************************

F ldi temperature
Foldingtt
off some proteinst iProtein

Number ofamino acid

Folding
temperature(K)albumin binding
domain53380.4

src SH3 domain

56

342.9

protein G

56

338.2

α-spectrin SH3
domain

57

360 1
360.1

Sso7d

64

332.0

protein L

78374.2

Im986382.0

cytochrom B562

106352 2352.2

“Switching”
g simulation<<<< native_info_sim1
native info sim1
Intra-mol 1 NINFO(1/1) 1…
Inter-mol 3/6 NINFO(3/6) 13
1= f1atp_alpha_E.ninfo
fl hf…
13= f1atp_alphDP_betaDP.ninfo>>>><<<< native_info_sim2
native info sim2
NINFO(1/1) 1…
NINFO(3/6) 13
1= f1atp_alpha_TP.ninfo
fl hf…
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
5 9973 1.0000
1 0000 1 0.3000
0 3000…

pair ij

nat- distanceRotation mechanism ofF1-ATPase by switching Go model
N KogaN.
Koga, and S
S. Takada
Takada, PNAS (2006)

Conformational change by MBP
K. Okazaki, N. Koga, S. Takada, J.N. Onuchic, and P.G. Wolynes, PNAS (2006)<<<< unit_and_statei_seq_read_style = 1
i_go_native_read_style = 1
1aprotein
1GGG_2.pdb
1bprotein
1WDN_2.pdb
1WDN 2 pdb>>>><<<< energy_function
NLOCAL(1a/1a) GO EXVNLOCAL(1b/1b)( b/ b) GO EXVMULTIGO_SYSTEM(1a) 1a/1a
MULTIGO_SYSTEM(1b) 1b/1b
i use atom protein = 0i_use_atom_protein
i_use_atom_dna = 0>>>><<<< multiple_gobdbdemax_mgo= 100.0
100 0baemax_mgo = 1.0
dihemax_mgo = 0.5
ENEGAP(1)(1) 0.0 -1.8
1.8
DELTA(1ab) 28.0>>>>Sliding movement of KIF1A
R. Kanada, et al PLOS Comp Biol (in press)

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

DNA duplex• 30 bp oligomer of DNA• Langevin dynamics (300K)• [Na+] = 69mM<<<< unit_and_statei_seq_read_style = 2
i go native read style = 3
i_go_native_read_style
1-2 dnasequence>>>><<<< energy_functiongyDH
LOCAL(1-2) L_BDNA
NLOCAL(1-2/1-2) DNA ELE
i_use_atom_protein = 0
3SPN.1
i use atom dna = 0
i_use_atom_dna>>>><<<< electrostatic Intra mol 1,2Inter mol 1-2cutoff_ele = 20.0
ionic_strength = 0.069
diele_water = 78.0>>>><<<< in_box
in boxxbox = 120.0
ybox = 120.0
zbox = 120.0
boxsigma = 4.0>>>>Simulation of nucleosome•

Electrostatic interaction +Go potential

1KX5 db
1KX5.pdb
εgo ppro-dna:coefficient ofprotein-DNA Go potential

εgo ppro-dna = 0.5ε0 5εgo ppro
[Na+] = 50mM<<<< energy_function
NLOCAL(1-2/1-2) ELE DNANLOCAL(1-2/3-10)
NLOCAL(12/3 10) GO EXV ELENLOCAL(3-10/3-10) GO EXV
i_use_atom_protein = 0
i_use_atom_dna = 0>>>><<<< electrostatic
cutoff_ele = 20.0
ionic strength = 0.05
ionic_strength
diele_water = 78.0>>>>H. Kenzaki, et al unpublished data

Menu
ModelsSimulation methods
Implementation
Selected applications
IIn-progressmodelsd l & methods
th dAcknowledgement
CafeMolCf M lddevelopmentlthhas bbeen supported
t db
by R
Research
hand Development of the Next-Generation Integrated
Simulation of Living Matter
Matter, a part of the Development
and Use of the Next-Generation Supercomputer Projectof the Ministry of Education, Culture, Sports, Scienceand Technology.