2015/9/1
2015年バイオスーパーコンピューティング・ソフト講習会シリーズ「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• 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. Biomolecular simulation and CafeMol
2. Model and applications of CafeMol
3. Recent development of CafeMol

Menu

1. Biomolecular simulation and CafeMol
2. Model and applications of CafeMol
3. Recent development of CafeMol

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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!!

Molecular dynamics• Calculate molecular movementmF– Classical mechanics– Numerically integrate Newton’s equation of motion– m(mass) x a(acceleration) = F(force)– interaction parameters are derived from MO (AA MD)• Interaction– bond length, bond angle, dihedral angle– van der Waals interactiondihedral angle– electrostatic interaction

bond angle
bond length

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

Menu

1. Biomolecular simulation and CafeMol
2. Model and applications of CafeMol
3. Recent development of CafeMol

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

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

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

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

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 

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 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 DNA simulation••••

30 bp DNA duplex
Langevin dynamics (300K)[Na+] = 69mM
cutoff length 20kD(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, Plos Comput. Biol.
to be publishedNucleosome wrapping simulation
For demonstration
Na+=200mM

Dinucleosome wrapping simulation
For demonstration
Na+=100mM

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 kD  
 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

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)

Menu

1. Biomolecular simulation and CafeMol
2. Model and applications of CafeMol
3. Recent development of CafeMol

Recent development of CafeMol• Partial charge– RESPAC by Terakawa• CG GBSA– CG Accessible Surface Area (ASA) by Li– CG Generalized Born (GB) by Chang• Multicanonical simulation by Ito• New DNA models– 3SPN.2 by Terakawa and Tan– 3SPN.2C by Tan• Refined hydrophobic interaction by Tanaka andKanada• Implicit hydrodynamics by Hori
• 2-beads model by Li

Determination of partial charge
in Cα CG model
Glu

Ala

Conventional

partial charge

-1

-0.3

0

-0.7

Terakawa
JCTC 2013

protein specific partial charge

RESP-like procedure from Atom to CG (RESPAC)
Minimize χ w.r.t. q’sDetermination of partial charge
in Cα CG model
Lac repressor

Proteininterior

Arg22

AA, PBCG, DH+RESPAC

Proteinexterior

CG, DH+integer

50

CG Accessible Surface Area （ASA）
POPS: Fraternali & Cavallo, NAR 2002,
Based on Wodak & Jenin, PNAS(1980)

Li

SASA   Aii

pi pij bij  rij  Ai ( r )  Si   1 



Si
j 1

j i 
NN

bij  0；when rij  Ri  R j  2 Rsolv
bij    Ri  Rsolv   Ri  R j  2 Rsolv  rij 1  ( R  R )r  ; when r  R  R  2 R
j

i

1
ij

Si  4 ( Ri  Rsolv ) 2

ij

i

j

solv

CG Accessible Surface Area （ASA）
Total SASA by POPS and Naccess
POPS-AA vs. Naccess

POPS-CG vs. Naccess
10000

R=0.99slope=1.10
8000

6000
45 proteins
4000
4000

6000

8000

SASA_NACS_total

10000

SASA_POPS-CG_total

SASA_POPS-AA_total

10000

R=0.99slope=1.04

8000

6000
45 proteins
4000
4000

6000

8000

SASA_NACS_total

10000

Debye-Huckel model (ionic solution)

Ionic strength

CG Generalized Born Model
Still’s formulae for GB

Born radius Ri is defined by

Based on Dominy-Brooks formulation of GB

Ri   166 / G pol , i

1
1
1


dV
4

Ri RvdW , i in , r  RvdW , i ri

55

CG Generalized Born Modelvia multiscale fitting of Born self energy
We fit self-energy of point charge
based on Dominy-Brooks formulation of GB~ G pol , iParameter optimization-1/ Born radius (CG, GB)

Test against AA PB

AA , PB

fitting-1/ Born radius (AA, PB)

Solvation free energy (CG, GB)

G pol , i

CG ,GBLe Chang

Solvation free energy (AA, PB)

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.