2016年度 ポスト「京」重点課題1

GENESIS講習会（講義）
理化学研究所 生命システム研究センター
分子機能シミュレーション研究チーム
李秀栄、尾嶋拓、新津藍、杉田有治（資料提供：小林千草）

2016年9月8日 FOCUS（神戸）今日の内容n 13:30 – 13:50 講義• 生体分子の分子動力学計算• GENESISの特徴

n 休憩（10分）
n 14:00 – 16:00 実習•••

GENESISのコンパイル
実行例(1) 分子動力学法計算と解析
実行例(2) レプリカ交換分子動力学計算と解析

今日の内容（改定）n 13:30 – 13:45 講義• 生体分子の分子動力学計算• GENESISの特徴

n 13:45 – 14:45 実習（１）••

GENESISのコンパイル
分子動力学法計算と解析

n 休憩（15分）
n 15:00 – 16:00 実習（２）•レプリカ交換分子動力学計算と解析

分子動力学法(Molecular Dynamics; MD)
n 粒子間の相互作用力を計算し、ニュートンの運動方程式を解く
d 2ri (t )dU (r N )mi
= Fi ( t ) = −
2
dt
dri

相互作用エネルギー粒子の座標

n 粒子の運動から、構造の安定性、構造ゆらぎ、更にアンサンブルを計算し、自由エネルギー面（結合エネルギーなど）を解析できる
n しかし、MD計算はとても時間がかかる…
• タンパク質のダイナミクスをミリ秒計算するには、ニュートンの運動方
程式を1012回計算する必要がある• 例えば1stepに1msかかるなら、約31.7年必要！

タンパク質のシミュレーション
fs
-15

ps
-12

ns
-9

µs
-6

分子の振動
タンパク質
タンパク質
＋水

二次構造変化

超巨大
生体分子
＋溶媒

size

構造緩和
HP362
(~104)

Folding of
Globular proteins8~2010

HIV-1 capsid7(~6.5x107)

粗視化モデル

log10t

BPTI 4, Ubiquitin 5(~2x104)

2000

Ribosome6(~2x106)

全原子モデル()は原子数time

タンパク質の反応

BPTI 1(~6x102)
AQP13(~105)

s
0

Folding

1983

膜タンパク質＋水＋膜分子ms
-3

Drug export
from AcrB9, 10

1. Brooks et.al, PNAS,
80,pp6571-6575 (1983)
2. Duan et al., Science,282,
pp740-3744 (1998)
3. de Groot et al., Science,294,
pp2353-2357 (2001)
4. Shaw et al., Science,330,
pp341-346 (2010)
5. Piana et al., PNAS,110,
pp5915–5920 (2013)

6. Whitford et al., PLOS Comput.
Biol., 9, e1003003 (2013)
7. Zhao et al., Nature, 497,
pp643-646 (2013)
8. Clementi et al.. JMB, 298, pp
937-953 (2000)
9. Yao et al., Nature Comm. 1:117(2010)
10. Yao et al., JACS, 135,
pp7474-7485 (2013)

could be.

important for the question of interest should be coars
Second, no general procedure to select the type andof atoms per CG bead exists, which leads to a largeof CG models being proposed [11]. For drug desig
a high resolution of the ligand and active site is
be important, atomistic/CG hybrid or multiscale mparticularly interesting [12].
In classical MD simulations, the potential energy
commonly denoted as the ‘force-field’ describes the inbetween the atoms (or generally particles). It consists obonded terms between covalently bound atoms (i.e. bon
torsions) and nonbonded terms (i.e. van der Waals inand electrostatic interactions). The parameters of theare typically fitted to reproduce data from higher-level ca
コンピュータスクリーニングの導入で
and/orexperiments. For the majority of biomolecule
J. Chem.
Phys.a 137,
230901 (2012)proteins, 合成効率はどう変わる？
DNA, lipids andsugars,relatively
small nbuilding blocks can be used, which have to be paramete
once. Traditional biomolecular force-fields such as AM
OPLS [14], CHARMM [15] and GROMOS [16] have been
ously improved and tested over the years. The diversit
organic molecules by contrast prevents
the use of buildi
コンピュータース
and ligands have to be parameterizedindividually. Para
クリーニング
the bonded terms and the van der Waals interactions aretaken from a generalized force-field,
whereas the parti
医薬品化学者
charges are derived from a QM calculation. The force-fi
commonly used for this purpose are the general AMBER(GAFF) [17], the OPLS all-atom force-field (OPLS-AA)
CHARMM general force-field (CGenFF) [19,20] and theautomatic topology builder [21].

重点課題1：創薬計算の高精度化・高速化

Computational tools

The most established computational tools available for classicalMD simulations are the software suites GROMACS [1], NAMD [2],CHARMM [3] and AMBER [4] (Table 1). MD simulations are anumerical solution of Newton’s equations of motion (Eq. (1)):@Vðr Þ
n 基質結合とタンパク質構造変化の両方を取り込む(1)@rn スパコンを活用した高並列スクリーニング
wherem is the mass of atom i, a the acceleration of atom i and f isN

f i ¼ mi a i ¼ "
i

i

i

i

the force acting on atom i given by the partial, spatial derivative of
the potential energy function V which is dependent on the posi-

創薬計算に主に用いられてい230901-2
D. L. Mobley and P. V. Klimovich
るMDソフトと力場
TABLE 1binding of another ligand, thus yielding minimal errors. While

Main software
force-fields
in MD.it seems reasonable to us
thisandcancellation
is notused
guaranteed,that for small ligand modifications,energyDeveloping group
MD suite relative freeForce
fieldcal-culations will indeed be more efficient. Also, accurate predicGROMACS [1]GROMOS [16]tion of relative binding free energies would be well suited to
CHARMM [3]CHARMM[15]
Harvard (USA) applications in lead optimization,
as noted above.
For thesereasons,
we
focus
here
on
relative
calculations.
AMBER [4]AMBER [13]San Francisco (USA)
Alchemical free energy calculations work by introducing
GROMOS [22]GROMOS [16]Zurich (CH)a series of intermediate unphysical states spanning between
Purdue/Yale University(USA)– example, for relativeOPLS[14]
the desiredend states. For
calculations,binding
free
energies
arecomparedbychanging
one
ligand
NAMD [2]–Urbana-Champaign (USA)
into another or turning off interactions of one ligand in a re[23] of another– ligand in a
Barcelona (ES) ceptor while turning onACEMD
interactions
DESMOND states here– are, roughly
New York (USA)receptor (Fig. 3). The intermediate
speaking, associated with fractional presence of each ligand,as discussed in more detail elsewhere.4, 18–21
www.drugdiscoverytoda（MorDer, J. et al.
Drug Discovery
Today, 2015,
Alchemical
free energy calculations
are built20:686-702,
on ground25
結合エネルギーの変化
and Zwanzig,
among
Kirkwork laid by Kirkwood
FIG. 1. The probability of synthesizing a compound with a particular bindMobley, D. & Klimovich,P. V. J.24Chem.
Phys.
2012,others.
137:230901.）
ing free energy change. Filled regions indicate those compounds with at least
wood described a coupling parameter approach, where a pa-

合成確率

Groningen (NL)Generalized Ensemble SimulaCon Systems (GENESIS)n 生体系において、効率的で精度の良い自由エネルギー計算
手法の開発が目的
n 高並列計算 - 京コンピュータなどでの超並列計算が可能（ただし、普通のPCクラスタでも動きます）
n 巨大な生体系も可能
n 全原子モデルのみでなく、粗視化モデル等異なる分子モデ
ルへも応用できるアルゴリズムを採用
n レプリカ交換法（拡張アンサンブル法の一つ）による
自由エネルギー計算も可能

n 開発チーム
計算科学研究機構(AICS), 粒子系生物物理研究チーム
プロジェクトリーダー： 杉田有治
主開発者：
鄭載運、森貴治、小林千草、松永康佑、安藤格士（理科大）、神谷基司、依田隆夫（長浜バイオ大）、Michael
Feig（ミシガン州立大）
公開バージョン
：Ver. 1.1.0

n ホームページ
hVp://www.aics.riken.jp/labs/cbrt/
ソースコードとマニュアは上記ページ内のリンクよりダウンロード可n 参考文献J. Jung, T. Mori, C. Kobayashi, Y. Matsunaga, T. Yoda, M. Feig, and Y.
Sugita, WIREs Comput. Mol. Sci., 5, 310-323 (2015).

生体分子MD計算のボトルネック•

分子動力学法は膨大な回数のニュートン方程式を解く•

最も時間がかかる部分は粒子間の相互作用計算の部分である
生体分子の相互作用計算は「力場(force ﬁeld)」と呼ばれる経験的な関数で記述される
force ﬁeldはおおまかに２つの部分に分けられる– 結合性相互作用： 原子間の共有結合によるもの– 非結合性相互作用：

長距離相互作用（電荷、van der Waals力)によるもの

粒子数(N) に対してO(N2)で時間がかかるため、最も時間のかかる部分

結合性相互作用• Bond• Angle• Dihedral

Bond

Angle
Dihedral

非結合性相互作用• Coulomb• van der Waals→ 非結合相互作用の高速化が必要

非結合相互作用の高速化
2)からO(N1)に時間短縮
非結合相互作用の高速化
カットオフ長を導入することでO(N1.

1. 非結合相互作用はカットオフ長を導入することでO(N2)からO(N1)に変更できる
O(N2)

O(N1)

全てのペアを計算：
O(N2)全てのペアを計算

カットオフ長以内のペアを計算：O(N1)
カットオフ長以内のペアを計算カットオフ長より長い部分のクーロン相互作用はフーリエ変換(FFT)により逆空間で計算
2. 2. カットオフ長より長い部分のクーロンはフーリエ変換により逆空間で計算
する

erfc

U
|

|

4

実空間（近距離）

3.

2

exp
G

G ⁄4
G

4

逆空間（長距離）cos G ∙

4

Self energy（定数）

Self energy
実空間
逆空間
(定数)
(近距離)
(長距離)
1,2に加えて、非結合相互作用の高速化、高並列化は近距離項、長距離項の双3. 1,2に加えて、非結合相互作用の高速化、高並列化は近距離項、長距離項の双方で更方で更に行われるに行われる
7

GENESISの高速化・高並列化
n 最も時間のかかる非結合相互作用計算を高速化するため
GENESISでは下のアルゴリズムを新規に開発、組み込む– Inverse lookup table approach (近距離相互作用計算の高速化）(Jung et al., J. Comput. Chem., 34:2412–2420 (2013))– Midpoint cell methods (近・長距離相互作用計算の高並列化）(Jung et al., J. Comput. Chem., 35:1064–1072 (2014))– ParallelizaCon of FFT (長距離相互作用計算の高並列化）(Jung et al., Comp. Phys. Comm., 200: 57-65 (2016))

n 上の３つに加えて、トラジェクトリなどの書き出し・読み込みを高速化するため並列I/Oスキームを独自に持つ

Inverse lookup table法
Jung et al., J. Comput. Chem., 34:2412–2420, (2013)n 近距離相互作用は距離の関数であるため、カットオフ長までの距離の代表点での相互作用の値を計算、メモリ(table)に記憶しておく。
メモリに記録されている相互作用値r
GENESISの方法

従来の方法
0.10

0.05

0.05

0.00

0.00

-0.05

-0.05

12

-0.10

12

1/r -1/r

6

1/r -1/r

6

0.10

-0.15

-0.10
-0.15

-0.20

-0.20

-0.25

-0.25

-0.30

1

2

3

4

5

6
2

r

7

8

9

10

-0.30
0.1

0.2

0.3

0.4

0.5

0.6
2

1/r

0.7

0.8

0.9

1.0

n 従来の方法ではr2の線形関数、３次関数で内挿していた物を、GENESISでは1/r2の線形関数で内挿し、高速で精度の良い計算を可能とするMidpoint cell 法

Jung et al., J. Comput. Chem., 35:1064–1072, (2014)

n ドメイン分割法（ほとんどのMDプログラムで採用）
•
•

•
•

系全体をカットオフ長より長い一辺をもつセルで区切る
エネルギー計算は隣接のboxのみを考えればよい(黄色のセル内の原子の相互作用は黄と緑のセル内の原子のみをカウントすればよい）
通信回数が減少されるため、並列度が良くなる
ドメイン分割法で異なるセル間の相互作用計算をどのCPUコアに割り振るかが並列度に重要な問題にな
る

•

•

7

56

63

Midpoint法1

n Midpoint cell法•

0

従来のmidpoint法では、２つの原子の中間地点のセルを受け持つコアが計算する
それぞれの原子が存在する二つのセル（８, 10とする）の中間のセル（7,11のどちらか）を受け持つコアが計算を担当する
通信回数が減少されるため、並列度が良くなる
1 KJ. Bowers et al, JCP 124, 184109 (2006)

Midpoint cell法

(a)ベンチマークデータ
Simulation time (ns/day)

128

DHFR

64
32

(b)
NAMD

32

SPDYN

16

CHARMMn ベンチマーク詳細
16

NAMD

ApoA1

SPDYN
CHARMM

8
4

ATDYN
参照：hVp://www.aics.riken.jp/labs/cbrt/benchmark/2

J. Jung, T. Mori, C. Kobayashi, Y. Matsunaga, T. 1Yoda, M. Feig, and Y. Sugita, WIREs
PC cluster
PC(2015).
cluster
Comput.
Mol. Sci., 5, 310-323
4

n例

16

WIREsComputational
ComputationalMolecular
MolecularScience
Science
WIREs

8

88
44

(c)
(c)
88
44

Simulation time (ns/day)

1616

STMVSPDYNSPDYN

4CHARMMCHARMM

2
1
128
128

STMV
STMV

0.125

1616

(e)
3232

16

6464

16
32
GENESIS
GENESIS

NAMD
NAMD

NAMD

(d)
32

32

SPDYN
SPDYN

16

SPDYN
NAMD

CHARMMCHARMMCHARMM

SPDYN

ATDYN
ATDYN
ATDYN

CHARMMPC cluster
32
3232

64
6464

PCclustercluster
PC
128
256
512
128 256
256 512
512
128

128

GENESIS
256

512

32

SPDYN

STMV

CHARMM

8
16

SPDYN

4

8

ATDYN

2
1

4
0.5

NAMD

PC cluster

16

32

2

64

128

256

512

(d)

(d)
(d)
8

1

PC cluster

64

64
NAMD

ApoA1

32
STMV
3232 STMV
NAMD
STMV
4 STMV
16
SPDYN
SPDYN
1616
SPDYN
SPDYN
2
8
88
CHARMM
1
4
NAMD
NAMD
CHARMMCHARMM
44
0.5
2
22
0.25
1
11
PC
cluster
0.125
0.5
computer
KKcomputer
PCclustercluster
PC
0.5
0.5
128
512
16
32
64
128
256
512
128
512
2048
8192 32768
32768
128
512
2048
8192
128 256
256 512
512
128

16

0.5
0.5

0.125
0.125

512
512

0.5

(b)

NAMD
NAMD

22

0.25
0.25

141
0.5
0.5 16
1616

(c)

0.25

11

256
256

512

(b)
(b)

282

PCclustercluster
PC

6464
0.5

256

16
44

ATDYN
ATDYN

3232

128

128 DHFR
3232 ApoA1
ApoA1
64
1616
32
88

NAMD
NAMD

Simulation time (ns/day)

3232

Simulation time (ns/day)

6464

y)

Simulation time (ns/day)
Simulation time (ns/day)

64

(a)

128 DHFR
128
DHFR

1616

Simulation time (ns/day)
Simulation time (ns/day)

32

WIREs Computational Molecular Science

(c)

(a)
(a)

ATDYN

8

128

Crowding (11.7 million)

256

SPDYN

0.5

128

512

512

K computer

NAMD

2048

8192

32768

京コンピュータでの計算
K computer
8
大規模資源での高並列効率を実現
2048
8192
32768

(f)

Crowding (103.7 million)

GENESISの内容物（1）n MDプログラム‒ ATDYN (ATomic decomposiDon DYNamics simulator):• atomic decomposiDonを使用• 粗視化モデル（Cα-GO, all atom GO)が計算可能• わかり易いコードでユーザによる開発可能‒ SPDYN (SPaDal decomposiDon DYNamics simulator):

今回の実習で• domain decomposiDonを使用
はこちらを使用• 超並列・高速なルーチン (Midpoint cell/3次元分割FFT/並列 I/O）

特徴

ATDYN

SPDYN

システムの分割法

原子分割

ドメイン分割New lookup table○○

レプリカ交換法○○

粗視化モデル○×

3次元分割FFT

×

○

並列I/O

×

○

r-RESPA (1.1.0)

×○GENESISの内容物（2）Isotropic, Semi-iso, An-iso,XY-ﬁxed)n ATDYN/SPDYN共通
•
•

•

最適化
² Steepest Decent法
Integrator² Leapfrog² Velocity Verlet
アンサンブル² NVE² NVT– Langevin– Bussi– Berendsen² NPT– Langevin Piston– Bussi– (Isotropy of SimulaDon box:•

拘束計算(constraint)² SHAKE (Leapfrog)² RATTLE (Velocity Verlet)² SETTLE• FFT (PME)² FFTE• Restraint funcDons² PosiDon² Bond² Angle² Dihedral angle² RMSD• Target/Steered MD• String method（赤字は1.1.0の新機能）

GENESISの内容物（3）n その他解析ツール••

トラジェクトリ変換（trjcnv）² crd_convert² pcrd_convert（並列I/O用計算用）² remd_convert（REMD計算用）
トラジェクトリ解析（trjana）² mbar_analysis² pathcv_analysis² pmf_analysis² gval_analysis² rmsd_analysis（RMSD、RMSF）² trj_analysis（距離、角度など）² wham_analysis•

リスタートファイル変換（rstcnv）² rst_convert² rst_upgrade•

PCA解析（pcaana）² avecrd_analysis² eigmat_analysis² ﬂccrd_analysis² prjcrd_analysis• 反応経路探索関連（rpath_generator）• 解析共通（libana）