分子シミュレーションの展開
スーパーコンピューティングに向けて
次世代生命体統合シミュレーション研究開発
分子スケールチーム

理化学研究所 次世代計算科学研究開発プログラム
横浜市立大学大学院生命ナノシステム科学研究科
木寺 詔紀

タンパク質は動く
1．このように動く
1－1．データベースから見えるタンパク質の応答的運動
1－2．応答的タンパク質機能像
2．そのように動くタンパク質機能の分子シミュレーション
2－1．もっと速さを！
2－2．次世代生命体統合シミュレーション 分子スケール(1) 多剤排出トランスポーターを例に分子動力学計算
粗視化モデル計算(2) 量子化学計算(3) マルチスケールシミュレーション

1．このように動く
1－1．データベースから見えるタンパク質の応答的運動
タンパク質は一般的に「動く」と言えるものなのか？
「動く」とは意味のあることなのか？
「動く」、どのように？
Protein Data Bank

リガンド非結合状態

リガンド結合状態

因果関係
リガンド結合 → 構造変化

Single Molecules

Classification

T. Amemiya et al. J.Mol.Biol. 408, 568 (2011)

この程度動くCαRMSD ＞ 1Å立体構造変化あり

立体構造変化なし

528311

動かないことの意味？
酵素機能の分布EC6 EC1
EC5EC4EC2EC3EC3
Domain Motion

EC2Local Motion

EC2EC3No Significant
MotionR. Koike et al. J.Mol.Biol. 379, 397 (2008)

このように動く

Domain Motion

Local Motion

Classification

Buried Ligand
CαRMSD ＜ 1Å

104
動かなければリガンドは入れない！

Buried Ligand Classの特徴
結合状態

非結合状態

結晶水
284

10-Formyltertrahydrofolate dehydrogenase

側鎖177
Trp operon repressor

リガンド結合との因果関係

原因がリガンド結合である
221

原因がリガンド結合でない
181

other than buried ligand

Classification

結晶場の影響

結晶場によってひずみを受けている

網羅的シミュレーション

PDB ID1

PDB ID2

RMSDs from initialstructures at 10 ns [Å]Pair-wise RMSD [Å]Protein 1

Protein 2

InitialMinimum(time [ns])1n0u_A

1n0v_D

38.6

28.1

34.7

18.6 (5.2/0.4)

2qyu_A

2qza_A

30.6

11.4

6.1

9.0 (0.1/0.3)

1eut_A

1w8n_A

23.1

14.8

2.7

2.0 (0.1/0.1)

2zgv_A3c3b_A

3.6

2.3

7.0

1.8 (1.1/6.0)

1e4f_T

1e4g_T

1.6

3.3

2.7

1.7 (5.2/5.4)

2v8i_A

2v8j_A

1.9

1.1

2.1

1.0 (7.4/3.3)

2ps0_A

2ps3_A

3.22.62.0

1.6 (7.9/3.7)

3b8s_B

3b9d_A

5.4

3.0

2.4

1.3 (2.6/1.3)

2ddb_A2epf_D

3.3

4.9

1.8

1.2 (0.7/0.9)

2o1c_A

2o1c_D

2.0

2.1

1.5

0.8 (2.6/1.1)

1rif_A

1rif_B

4.7

2.9

1.1

1.3 (0.7/0.2)

SinglelinkerMultiple
linker寺田透（東京大学、理研）

Correlation with the number of linkers
coupled domain
motions
independent
domain motions

1．このように動く
1－2．応答的タンパク質機能像 予測可能性

FoldingFunction

機能シミュレーション＝ 応答的非平衡シミュレーション

線形応答理論

Molecular Dynamics Simulation

fj

Elastic Network Model果因Linear Response Theory of Protein Structural Changes

r   Σf

Maltose-binding periplasmic protein

ドメイン運動：線形応答理論の予測
閉じる 開く閉運動

45

閉運動 : ●開運動 : ▲再配置運動: ■

14

再配置運動

開運動

ローカル運動：線形応答理論の予測

121

38
閉運動

閉運動 : ●開運動 : ▲開運動ローカル運動の多くは非線形的Order-disorder transition

Helix-coil transition

リガンド結合に伴うローカル運動の約60%に見られた

2．タンパク質機能の分子シミュレーション
2－1．もっと速さを！

Molecular dynamics simulations for biomolecules

Bovine pancreatic trypsin inhibitor(BPTI) 454 non-hydrogen atoms
In vacuum~8.8 ps

IBM 360/91：~1MFLOPSAnalyses of Trajectory Data
RMSD/RMSF
Autocorrelationfy
Crosscorrelation

Protein’s vibrational dynamics

BPTI

> 1ms
~10,000 atoms including
explicit water

108 time longer simulation

D. E. Shaw et al. Science 330, 341 (2010)

Bacteriophage lambda protein W (gpW)
58 residues

Lignocellulosic Biomass
A system was constructed oflignocellulosic biomass containing 52
lignin molecules each with 61 monomers,
the same cellulose fibril as describedabove and 1,037,585 TIP3P watermolecules, totaling 3,316,463 (or 3.3million) atoms.

7,300 times larger system(3 x 107 times more interactions)
JaguarPF, a Cray XT5 massively parallel processing
computer with 30,000 Opteron 2.3 GHz cores
200TFLOPS
2x108 times larger power

R Schulz, et al. J. Chem. Theory Comput 5, 2798 (2009)33 years: 1977

2010

Quantitative Measure
Computer Power2 x108Target Systems
Size (Interaction)
Time Scale

3 x107
1 x108

Qualitative Measure
Physics: Vibrational Dynamics

Biology: Biomolecular Functions
Real Simulations

Next-Generation Supercomputer 「京」10 PFLOPS
6.4 x 105 cores
1PB memory

Next Generation Integrated Simulationof Living Matter2 タンパク質機能の分子シミュレーション
2－2 次世代生命体統合シミュレーション 分子スケール
量子化学 QM

林重彦（京大）

ProteinDF

中村春木（阪大）

佐藤文俊
（東大）

高精度
反応

全原子分子
動力学 MMMarble

池口満徳
（横市大）

「動き」の中にタンパク質機能を見いだす分子シミュレーション

マルチスケール
Platypus

杉田有治
（理研）
木寺詔紀
（理研）

高効率
運動粗視化モデル
CG CafeMol
高田彰二
（京大）

(1) 多剤排出トランスポーターを例にMultidrug Efflux Transporter AcrB
Multidrug resistance:
A serious problem in the
chemotherapy and in the
antibiotic treatment

H+

AcrB is a principal multidrug
exporter.

Functionally rotating mechanism

Proton Transfer(QM)
Cooperative Structural Transitions(MM-MM/CG)
Drug Transport(MM/CG-CG)
S Murakami S, et al. Nature 443, 173 (2006)

Proton transfer → Structure changea.
Protonation states & Gate opening
All atom molecular dynamics simulation
T. Yamane, M. Ikeguchi (YCU)
Structure change → Drug transportb.
Ligand binding sites
3D-RISM calculation
T. Imai, A. Kidera (Riken), F. Hirata (IMS)c.

Cooperative motions
Triple Go model simulation
X. –Q. Yao, S. Takada (Kyoto U)

a.

All atom molecular dynamics simulation

T. Yamane, M. Ikeguchi (YCU)

Extrusion state
D407C - K940Nz
D408C - K940Nz
T978C - K940Nz-

D407 D408

-

T978
D407
D408

D407H D408H-

D407 D408H

Time (ns)

K940

b.

3D-RISM Calculation

3D-RISM/KH equationsh (r )    c (r ')  wvv' ( r '- r )   hvv' ( r '- r )  dr 'exp  d (r )  - 1 for d (r )  0h (r )  d (r )for d (r )  0
where d (r )  -  u (r )  h (r ) - c (r )Theoretical prediction

Thermolysin

Isopopanol+ water

T. Imai, et al. J. Am. Chem. Soc. 131, 12430 (2009)

Crystal structure

Isopopanol
AcetoneAcetonitrile
PhenolAppllication to AcrB

Fragment Molecules Used as Ligands

T. Imai, et al., J.Phys.Chem B in press

“Multi-functional” Sites
AccessBinding

MFS-3
MFS-1

MFS-2

Extrusion

c.

Triple Go-model Simulation

Ca CG model of Porter domain
Multiple-basin model

∆12

∆23

Okazaki, et al. PNAS 103, 11844 (2006)∆13

Trimer interface
Three types of multiple LJ type

∆13

∆V2Binding∆V3 Access
Extrusion

LigandMinocycline

Repulsive potentialXin-Qiu Yao, et al. Nature Commun. 1, 117 (2010)

?

BE?
?

EA?

AB?

?

(2) 量子化学計算
Have you seen MO of protein?
HOMO

Fe-d orbitalmodel
CO myoglobin

localized

LUMO

cytochrome cdelocalized

All-electron calculations (not divided) using ProteinDF
佐藤文俊（東大)

(3) マルチスケールシミュレーション
Coupled Simulation

Speed

Boundary conditions

Average quantities

Precision

Homogeneous Systems

Coupled Simulation in Biomolecules
Parallel Algorithm

NSerial AlgorithmMN >GS Ayton, WG Noid, GA Voth, Current Opinion in Structural Biology 17, 192 (2007)

MQMMMCoupling between Enzymatic Reaction
and Protein Motions
1. On-the-fly QM/MM MD simulation - S. Hayashi (Kyoto U)
H QM/MM  rQM , rMM 

2. QM/MM umbrella sampling - H. Nakamura (Osaka U)
H QM/MM  rQM , rMM 

01rQM   rQM
 1 -   rQM

3. QM/MM free energy - S. Hayashi (Kyoto U)
F  rQM   -1

ln  drMM exp -  H QM/MM  rQM , rMM On-the-flyMD Simulation

Bovine rhodopsin
Hayashi S, Tajkhorshid E, Schulten K, Biophys J. 96, 403 (2009)umbrella sampling
Cis-trans isomerization of asolvated Proline dipeptide
01rQM   rQM
 1 -   rQM

Hydrogen bondHydration

Y. Yonezawa et al., J. Am. Chem. Soc., 131, 4535 (2009)free energy
F  d, rQM   -1

ln  drMM exp  -  H QM/MM  d, rQM , rMM 

Complete separation of QM/MM optimization and MD sampling
Conventional

PresentMMQM

 MM  rMM QM  d, rQM iteration

MM  d , rQM , rMM  

S. Hayashi (Kyoto U)

MMQM

QM  d, rQM  MM  rMM 

reweighting0exp -   H QM/MM  d , rQM , rMM  - H QM/MM  d 0 , rQM, rMM 00
drexp
Hd,r,r
Hd,r, rMM  
MMQM/MMQM
MMQM/MMQM

0MM  d 0 , rQM, rMM 

Parallel algorithm
3D-RISM

Serial algorithm
Coupled Hamiltonian
Path Search – Path Ensemble
Sequential Monte Carlo Simulation

MS-MM/CG Simulation

MMCG drives MMMM is pulled by CGCGH  x; z   V  x   U  z   k 2   x  - z MM CGxzcoupling

2Removing the influence from the couplingH  x; z   V  x   U  z   k 2   x  - z MM CGxz
Hamiltonian
Exchange(multicopy)k→0

H  x; z   V  x MMMultiscale Enhanced
Sampling (MSES)

2couplingk→∞U z  0

Path Search(multicopy)- V x- F  z   ln  dxe     x  - z potential of mean force of MM
potential on z

On-the-fly String Method
L. Maragliano, E. Vanden-Eijnden,Chem. Phys. Lett. 446, 182 (2007).

2-2. Multiscale Enhanced Sampling (MSES)
K. Moritsugu (Riken)

k1

k3

k2

Hamiltonian
Exchange

k= 0

K Moritsugu, T. Terada, A. Kidera, J. Chem.Phys. 133, 224105 (2010)

Enhanced Sampling

Metropolis Criterion
Temperature REmn   m - n  V  xm ; m  - V  xn ; n min 1,exp   mn  MSESmn    km - kn  W  xm , z m  - W  xn , z n W (x, z) 1.
2.
3.

21
  x  - z 2

The constraint W can be set independently of V(x)The number of degrees of freedom in W can be small
The potential is harmonic
Excellent scalability has been attained in MSES.

Intrinsically Disordered Protein: Sortase A

Enhance the disordered regions

N Suree, et al. J. Biol. Chem. 284, 24465 (2009)Dynamic and Disordered Loops
20 replicas

Disordered Loop only

12 replicas

calcium+peptide

2-3. Path Search: On-the-fly String MethodH  x; z   V  x   U  z   k 2   x  - z MM CGstartk→∞U z  02coupling

Path Search(multicopy)

path

potential of mean force on z
end

No guiding potential → Entropy driven
Require some constraints → Path search
Path Search Methods
Minimum energy paths
Intrinsic reaction coordinate (K Fukui, 1981)
Self-penalty walk (R. Elber, 1987)
Conjugate-peak refinement (M. Karplus, 1992)
Nudged-elastic band (H. Jonsson, 1994)
String method (E. Vanden-Ejinden, 2002)
Finite tempearture paths
Max flux method (J.E. Straub, 1997)
Finite tempearture string method (E. Vanden-Ejinden, 2005)
On-the-fly string method (E. Vanden-Ejinden, 2007)

Y. Matsunaga(Riken)

Alanine-dipeptide in GB waterz  f ,y 
8 replicas
FinalInitial[kcal/mol]

Minimum Free Energy Path

Adenylate kinase in water(65,000 atom system)z: 20 principal components
65 replicas

タンパク質は動く
1．このように動く
1－1．データベースから見えるタンパク質の応答的運動
1－2．応答的タンパク質機能像
2．そのように動くタンパク質機能の分子シミュレーション
2－1．もっと速さを！
2－2．次世代生命体統合シミュレーション 分子スケール(1) 多剤排出トランスポーターを例に分子動力学計算
粗視化モデル計算(2) 量子化学計算(3) マルチスケールシミュレーション