粗視化生体分子シミュレータ
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 1.0
source & manual released• Features are;– Various CG protein models
multiple basin model
accurate CG model– Simulating protein-at-work“switching”• Under development
DNA/RNA lipid
DNA/RNA,

Overview of CafeMol• General-purposep pcoarse-grainedg((CG)) biomolecularmodeling and simulation software– Protein: 1 bead / 1 amino acid
– Nucleic acid: 3 beads (sugar, nucleotide,
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)

Menu
ModelsSimulation methods
Implementation
Selected applications
I
In-progressmodelsd l & methods
th d

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

Atomicaa
a.a.

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

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
V  K  r  r   K      K  1  cos     K  1  cos 3   
2

local

b

2



0 i ,i 1

i ,i 1ii

0i

bond angle
dihedral angle(0 means native state)

i

1

3



0i

i



i

i

i

  r 12  r 10 
0 ij
0 ij
Vgo  go  5   6   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 
ex 
= 4A
=1.0kcal/mol

Atomic interaction based CG (AICG) modelV   kbi (r i  r0i )CG  k ai ( i   0i )

2

i

2

i

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

Atomica.a.

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

ij 12

ij
0

i  j 3

ij 10

ij
0

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

Wenfei Li

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

2) Coefficients
C ffi i t fitted
fitt d b
by 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
0

0.29

1.76

1.32

0.82

0.81

0.37

GNM
MSM

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. .7575∼

0. 0.665 5∼

00. .55
55

0.0
<
0
0. .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 2
V ( R | R1)  V ( R | R 2)  V
 
 
2
 c2   log  c1 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

HP interactions analogousgto 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

kon
koff

V protein  Vimp lig

Menu
ModelsSimulation methods
Implementation
Selected applications
I
In-progressmodelsd l & methods
th dSimulation method• Dynamics– Newtonian dynamics with Berendsen thermostat– Langevin dynamics• Time integration– velocity Verlet algorithm• Run mode
– Constant temperature simulation– Simulated annealing– Auto-search of Tf– Replica exchange method– Potential “switching”

Menu
ModelsSimulation methods
Implementation
Selected applications
I
In-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 Rik
Riken

spee
ed up

100

10

1
speed upope(s)
mpi(s)

8

16

32

64

128

7.85

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 hybrid parallelization(MPI and Open MP)
20 nucleosomes (35918 particles)
RICC at Rikentime(s)

45
40
35
30
25

ope(s)
mpi(s)

20

main_loop(s)
15
10
5
0

OMP(force)
MPI(force)
MPI(replica)
speed up

1

8
1
64
1

2

8
2
64
1.75

4

8
4
64
3.19

8

8
8
64
5.26

15

8
15
64
5.64

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/SA
Blue; 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) 2 3
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 = ./pdbpppath_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
Bi
section 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.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
tf_outt ttempkk n_statet t dd_statet t p_transttf_out 325.000
953
48
19
*********************************…
…
*********************************
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
domain53

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

78

374.2

Im9

86

382.0

cytochrom B562

106

352 2
352.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
f
l h
f
…
13= f1atp_alphDP_betaDP.ninfo>>>><<<< native_info_sim2
native info sim2
NINFO(1/1) 1…
NINFO(3/6) 13
1= f1atp_alpha_TP.ninfo
f
l h
f
…
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 of
F1-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)23
NLOCAL(1b/1b)( b/ b)23
MULTIGO_SYSTEM(1a) 1a/1a
MULTIGO_SYSTEM(1b) 1b/1b
i use atom protein = 0
i_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 unpublished data

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
carge:yellow

Wenfei Li

Frustration, specificity & nonlinearity
Frustrationin large-amplitude motion of allostericproteins
Li Wolynes,
Li,
W l& Takada,
T k d PNAS 2011

Off-lattice Go model(merge of Go model & ENM)
Clementi Nymeyer,
Clementi,
Nymeyer & Onuchic 2000
V (R | R0 )  Kr (bi  bi0 )2  K (θi θi0 ) 
2

bonds

angles

１
3
{
K
[1
－
cos(
φ
－
φ
)]

K
φi－φi0 ))]}
 φii0
φ[1－cos(3(dih d ldihedral  r 12  r 10  nonnat  D 12
0ij
0ij



    ε2  

5
6

ε

1 
 r   i j 3  r   rij 
i j 3
 ij   ij 


nat contact

Subscript 0 means
the value at native

FragilebondPermanent
Kr = 100ｅ Kθ1 = 20ｅ Kφ3= 1.0ｅ Kφ= 0.5ｅbond
ε１=ε２=0.17ｅ D=0.45（Å）

Atomic interactions in allosteric proteins:
Ene g decomposition
EnergyAnalysisy of residue-pairwisepenergygy from atomic force fieldE ( RIJ )  uAA,ij (r)IJiI jJ

Gohlke et al., JMB, 2003, 330, 891

CG; I

CG; J

EIJ: coarse grained contact

uenergyAA,ij : All atom energy between atom
pair (i, j)
Using
Usg LCPOC O (pair(pa wise)se)SA
U AA ( r )  V ( r )  G GB
(
r
)


G
(
r
)
pol
AA energy include vacuum part andsolvation part, by AMBER99SB force field.

Case ett al.
C
l AMBER10,
AMBER10 2008
Weiser et al, JCC, 1999, 20,217

Contact energies in single-domain proteins
E
Energyddecomposition
iti
Exponential law!!

Contact energies in allosteric proteins:
Ene g decomposition
EnergyAdenylate kinase

1ake
close

close
l

open

4ake
open
•Diverse•Diff bet. classes

Atomic interactions in allosteric proteins:
Ene g decomposition
EnergyAdenylate kinase
Red: oppositely charged
Pink: charged+polar
Green:charged+apolar
Cyan:polar+polar
Green:polar-apolar or polar+apolar
Black: identically charged

D54, K200, D158, E167,
E162 K50 and D61
E162,

Atomic interactions in allosteric proteins:
Ene g decomposition
Energy41 allosteric p
proteins in pdb
pComp of class-average for each protein

sharedspecific

41 pairspAtomic interaction based CG (AICG) model
Deriving parameters
V ( R | R0 )  kbi (ri  r0i )2

Multiscale strategies:

Chu and Voth, 2006
Li et al., 2010,
See also
Trylska et al
Gohlke et al., 2003,i

 kai ( i  0i )2i {i [1  cos( i  0i )]i

 i [1  cos3( i  0i )] / 2}
native  ij [5(r0ij / rij )12  6(r0ij / rij )10 ]i  j 3nonnative

Energy decomposition →relative weight of contact energies
 Match rmsf bet. AA and CG→
scale local & tertiary interactions

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

For 23 training proteins--- Flow chart
CG parameterskb

CG parameters updating
C

No
Converged?

ka

X

 X

 ij   nlocwij

w

CG MD

IJ

AA-FF
EnergydecompositionuAA  E IJ W IJ

yes
stop

CG: coarse grained
AA: all atom
FF fo
FF:
force
ce field

Fb

FaX

FX

Fnloc

Fb

FaX

FX

Fnloc

AA MD

AA-FF

Fluctuation
matchinghi

Chu and Voth, Biophys. J., 2006,90,1572

F: mean square fluctuations

AICG model for mean fluctuations
Near native fluctuation (RMSFs)

Comp against AA MD

Stat for 30 proteins

AICG model for mean fluctuations
30 proteins
Fujitsuka et al, 2006, Proteins
Yang et al, 2007, Structure

Table 3. Average correlation coefficients (C.C.) and standard errors (S.E.) between the rmsfs
derived by AA model and by different CG models based on the proteins of testing set.

ENM

AICG model

Models

GNM

ANM

hete

homo

hete nloc
hete-nloc

C.C.

0.694

0.648

0.758

0.722

0.738

S.E.

0.018

0.031

0.021

0.027

0.031

hete: heterogeneous model.
homo: homogeneous model.
hete-nloc: only the nonlocal interactions are heterogeneous.

AICG model for predicting structural change
adenylatedlkinase
ki


o c
o c
d  (R  R ) / | R  R |i 
overlap (i )  v  d
V ( R | R open )PCA for Go model,
NM for ANM

V ( R | R close )AICG model for predicting structural change
C
Correlationl i w. structure changeh41 struct pairs of allosteric proteins
Gerstein NAR1998
Gerstein,

ENM

New model

ANM

hete

homo

hete-nloc

MO
M.O.

0 480 (0.037)
0.480(0 037)0 540 (0.037)
0.540(0 037)

0 511 (0.038)
0.511(0 038)0 518 (0.040)
0.518(0 040)C.O.

0.600 (0.040)

0.657 (0.039)

0.628 (0.040)

0.626 (0.042)

M.O.

0.421 (0.032)

0.517 (0.036)

0.475 (0.036)

0.512 (0.038)

C.O.

0.556 (0.037)

0.638 (0.038)

0.608 (0.038)

0.631 (0.039)

Modelsopenclose l

V ( R | R open )closeopen

V (R | R

close)

Both the interaction heterogeneity and anharmonicity are important for predicting theconformational change.

Large-amplitude fluctuation
adenylatedlkinaseki

V ( R | R close ) Fluctuation in close statecloseLarge-amplitude fluctuation in AKE
harmonic

homo

Hete-nloc

LID domaindi has
h higherhi h mobility
bilit thanth theth NMP domain.
diInteraction heterogeneity enhances the collective motions.

Large-amplitude fluctuation in AKE

hete-nloc

Very rare event in a long trajectoryMenu
ModelsSimulation methods
Implementation
Selected applications
I
In-progressmodelsd l & methods
th d

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       K  1  cos   

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

i

i

2



i

i

0i



i

0i

i

  12    6 
0 ij 
0 ij 

Vstack  4 1  
 r  i , j  rij
ij



 


4

bond angle
dihedral angle(0 means B-type DNA)

Nst

Go type 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 
bpAT 
 
 =0.1839kcal/mol

3S
3SPN.1force
o ce fielded(electrostatic and solvation interaction)  rij / D
qiqj

e
Vqq  
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 ) 2
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
i 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 

DNA duplex• 30 bp oligomer of DNA• Langevin dynamics (300K)• [Na+] = 69mM<<<< unit_and_state
unit and statei_seq_read_style = 2
i_go_native_read_style = 3
1-2 dnasequence>>>>DH<<<< energy_function
NLOCAL(1-2/1-2) 7 11
i use atom protein = 0
i_use_atom_protein
3SPN.1
i_use_atom_dna = 0>>>><<<< electrostatic Intra mol 1,2
I t moll 1
Inter
1-2
2
cutoff_ele
ff l = 20
20.00ionic_strength = 0.069
diele_water = 78.0>>>><<<< in_box
xbox = 120.0
ybox = 120.0
zbox
b = 120.0
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  go ppro
[Na+] = 50mM<<<< energy_function
NLOCAL(1-2/1-2) 7 11
NLOCAL(1-2/3-10)
NLOCAL(12/3 10) 2 3 7
NLOCAL(3-10/3-10) 2 3
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CG lipid model
H. Noguchi, and M. Takasu, Phys. Rev. E (2001)• i th lipid molecule:• 1 hydrophilicy ppparticle (j(j=1,, blue))• 2 hydrophobic particles (j=2,3, gray)Self-assemble
assemble bilayer structure• Self

CG lipid model
Vlip  Vlocal  Vrep  Vhydrobond angle

Vlocal  K b  r( i , j ),( i , j1)  r0 (i , j ),(i , j 1)  K  cos  i
i , j 1, 2

2

 r ( i , j ), ( i ' , j ' ) 
 1Vrep    exp  20
i  i ', j
Nst

Kb = 500
iK = 500
 = 0.6594kcal/mol
0 6594k l/ l
= 7.5A
2* = 10c2 = 4.75
4 75(  i , j   j * 1)
3* = 14(  j * 1   i , j   j *) c3 = 6.75


 0.5 
Vhydro    0.25(  i , j   j *) 2  c j
i , j  2,3 (  j *  i, j )c j
1Density dependenti , j   rattractioni  i ', j ' 2 , 3
exp 20 ( i , j ),( i ', j ')  1.9   1
s

 

Formation of vesicle
S. Fujiwara, et al unpublished data<<< job_cntl
i_run_mode = 2
i simulate type = 1
i_simulate_typei_initial_state = 5>>>><<<< unit_and_statei_seq_read_style = 3
i_go_native_read_style = 3
1 lipidsequence>>>><<<< initial_lipid
nmp_transverse_lipid = 20
nmp_longitudinal_lipidg= 20
nlayer_lipid = 1
grid_size_lipid = 1.075
z_coord_lipid(1) = 1.0>>>><<<< energy_function
NLOCAL(1/1) 17 19
i_use_atom_protein = 0
i_use_atom_dna = 0>>>>

Acknowledgement
CafeMolC
f M lddevelopment
l
th
has b
been supported
t db
by R
Researchhand 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.