計算生命科学の基礎

分⼦動⼒学計算における拡張アンサンブルシミュレーション
神⼾⼤学⼤学院⼯学研究科 応用化学専攻
中津井雅彦

コンテンツ• 先週の補足– KBDD: GPCR

• 統計アンサンブル– 運動方程式の一般化– NVTアンサンブル(カノニカルアンサンブル)– NPTアンサンブル(T-Pアンサンブル)• 拡張アンサンブル– 必要性– マルチカノニカル法– レプリカ交換法– 実例

MP-CAFEE法のGPCRへの適用
GPCRのMP-CAFEE計算手順構築
complex
平衡化計算

膜＋GPCR＋リガンド＋溶媒のプロトコル構築→力場問題→リガンドの部分電荷

complex
非平衡計算

パラメータファイルの作成（&検討）
安定性の解析→RMSD→膜厚→Area per Lipid

ligand
平衡化計算

ligand
非平衡計算

実行スクリプトの作成• パラメータの検討

dGの計算

力場：gromos43a1
リガンド：PRODRG＋GAMESS＋RESP
3

三井情報 佐藤氏のスライドより

運動方程式の一般化
ニュートンの第二法則( )=213=i54ラグランジアン=−ポテンシャル運動エネルギー

ラグランジュの運動方程式

dd

一般化座標

−=0=( ,,…,)

ハミルトニアン
ハミルトニアン( ) : 系の全エネルギー,=− ( , ): 一般化座標: 一般化運動量
ハミルトニアンの正準方程式=,=−d−ddd=

全エネルギーの保存

d=d−

全エネルギーの時間変化が0なので、エネルギー保存則が成⽴する「コンピュータ・シミュレーションの基礎」(第2版), 岡崎進, 吉井範行, 化学同人 2011=0

代表的なアンサンブル
アンサンブル
個々の座標や運動量は異なるが、熱⼒学的には同じ状態にある座標・運動量の集団
分⼦動⼒学計算では、運動方程式を解いて得られるタンパク質構造のスナップショットの集合(トラジェクトリ)• ミクロカノニカルアンサンブル (NEVアンサンブル)– 粒⼦数(N), エネルギー(E), 体積(V)一定

通常の分⼦動⼒学計算はNEVアンサンブル
ニュートンの運動方程式, ハミルトニアンの正準方程式はエネルギー保存則が成り⽴つため• カノニカルアンサンブル (NVTアンサンブル)– 粒⼦数(N), 体積(V), 温度(T)一定• T-Pアンサンブル (NPTアンサンブル)•粒⼦数(N), 圧⼒(P), 温度(T)一定

温度一定のMD (NVTアンサンブル)
カノニカルアンサンブル
運動エネルギーと温度の間の関係)

3" = $ − 1 &' (22速度スケーリング法
温度を制御するためには、速度を制御すればいい" = *"+*=

3$ − 1 &' (,∑"「タンパク質計算科学 基礎と創薬への応用」, 神谷成敏, 肥後順一, 福西快文, 中村春木, 共立出版, 2009

温度・圧⼒一定のMD (NPTアンサンブル)
能勢・アンダーソンの方法(温度制御)(圧⼒制御)

圧⼒を制御する方法(ピストン)

外界からの圧⼒./ を考える
ポテンシャルエネルギー01 = ./ 034

全ての座標を、0 でスケーリング
座標

速度= 02= 02「コンピュータ・シミュレーションの基礎」(第2版), 岡崎進, 吉井範行, 化学同人 2011

温度・圧⼒一定のMD (NPTアンサンブル)

温度を制御する方法
熱浴

時間の進み方が異なる(t -> t’)と考える=*+運動量の関係=+*0 2「コンピュータ・シミュレーションの基礎」(第2版), 岡崎進, 吉井範行, 化学同人 2011

温度・圧⼒一定のMD (NPTアンサンブル)
仮想系のハミルトニアン+ =+() 2

0 2*++02 ).6++ 3$ ' (, ln * + .1 : + ./ 027

現実系の運動方程式; =−0; =102−20 *+30 **; =*2* 3$&'( − (/ +*7 2**. − ./ + 0<*

0, =, *, 6 に関する項を消去
0, =に関する項を消去: ミクロカノニカルアンサンブル (NEV): カノニカルアンサンブル (NVT)「タンパク質計算科学 基礎と創薬への応用」, 神谷成敏, 肥後順一, 福西快文, 中村春木, 共立出版, 2009

評価データセット
共結晶構造が存在し、比較可能なKiや値が登録されているも共結晶構造が存在し、比較可能な やKd値が登録されているものを検証データとのを検証データとして使用
ドッキングスコア開発用データーセット• CHK1 : 6 リガンド
• CDK2 : 12 リガンド
• ERK2 : 6 リガンド
• uPa :7 リガンド

1セットあたり 約130,000時間(約15年) / コア
合計 5,330,000時間 ( 約608年) / コア「京」の30%を使用して約約 6日
文献調査により選択したGPCR• GPCR : 13 リガンド

CHK1： MP-CAFEE法による予測と実験結果比較
109

実験値 (pKi)

強弱876

実験では、新規な化合物を化学合成し生物活性評価するのに2ヶ月程度かかる
実験では、新規な化合物を化学合成し生物活性評価するのに ヶ月程度かかるのにヶ月程度かかるのに対しのに対し「京」では、3日程度で人間の労働力無しで、結果の取得が可能。
「京」では、 日程度で人間の労働力無しで、結果の取得が可能。
5-14.5-12.5-10.5-8.5

京が予測した結合(kcal/mol)
京が予測した結合自由エネルギー
結合自由エネルギー
数多くの化合物を計算するためにには、もっと高速な計算機が必要
数多くの化合物を計算するため
強 ⇒ ポスト「京」への強い期待！ 弱MP-CAFEE法による予測と実験結果比較
CDK2

ERK

-3

8.5-487.5-6

実験値（実験値（pKi））

実験値
実験値（
（ ΔG kcal/mol））-5-7-876.5-96-10-11-22-17-12
計算値（）
計算値（ΔG kcal/mol）-75.5-18.5-18-17.5-17-16.5-16計算値（ΔG計算値（kcal/mol））-15.5-15-14.5

平衡化したドッキングポーズとは？
十分な平衡状態に達している必要がある(λ1＝＝0)Equilibrium step:
5 sampling velocities × 50 ns MDs

状態遷移が一瞬の変化である必要がある

λ1＝＝0

λ2＝＝0.1λ3＝＝0.2λ32＝＝1.0

Unequilibrium step:
6 sampling velocities × 32 λ points × 2 ns MDs* H. Fujitani, et. al., Physical Review E, 79, 021914, (2009)

平衡化したドッキングポーズとは？
E85
ポケット残基のポケット残基のCαで重ね合わせた時のリガンドの残基の で重ね合わせた時のリガンドのRMSD
で重ね合わせた時のリガンドのC87

RMSD

平衡化○R95E85平衡化×平衡化×

RMSD

C87

■x-ray■トラジェクトリー40ns

R95Time(ps)同じ初期構造でも初速度によりドッキングポーズが変化する

非平衡計算開始構造の違いによる
MP-CAFEE法の結合自由エネルギー値への影響
E85
C87

R95平衡化×
平衡化×

E85
C87

平衡化状態におけるドッキングポーズを取得することが重要課題

R95平衡化○

現在の課題①～平衡化した構造とは？～プロテアーゼ（uPa)阻害剤の結合自由エネルギー計算結果
各化合物に対してΔGを幾つか算出してプロット各化合物に対してΔGの最小値をプロット

-15
5.5

66.5

77.5

8-17dG(計算値計算値)
計算値-197-21-23

6818-25

計算値 ΔG(kcal/mol)

5Abbott pKi-16-18-205678r = -0.40-22-24-26-27-29-28
Abbott pKi

開始構造（平衡化の選択時間）
開始構造（平衡化MDの選択時間）
平衡化によってΔGによって Gのバラつきが大きい実験値(実験値(Ki)を再現する傾向がみられるKi)を再現する傾向がみられる

平衡化MDのエネルギープロファイル

どのようにして、平衡化プロファイルから最安定構造を選択するか？
理化学研究所 計算科学研究機構 荒木氏

現在の課題②～結合ポーズの正確な予測～
E85
C87

R95 間違った結合ポーズ

E85
C87

あらかじめ正しいドッキングポーズを取得することが重要課題

R95
正しい結合ポーズ
理化学研究所 計算科学研究機構 荒木氏
三井情報(株) 佐藤氏

CHK1のドッキング計算（Redocking）
共結晶構造のRedockingでも結晶構造を当てられない場合もある
でも結晶構造を当てられない場合もある
共結晶構造のCSARID#Pose

GSCORE TOP
DOCKING POSE

RMSD TOP
DOCKING POSE

RMSD

RANK

RMSD

5370.29

10.29

19901.32

110.4740970.81

20.3170540.31

10.3172

67.90

63.998987

9.26

40.51

＊ドッキングソフト：GLIDEを使用を使用＊ドッキングソフト：
三井情報(株) 佐藤氏

Docking計算法で予測された結合ポーズとX線構造との比較（CSAR: CHK1_89の例）の例）
反転■x-ray■docking pose
三井情報(株) 佐藤氏

直行

並行平衡化MD計算によるドッキングポーズ変化
rank1

rank69rank72

rank56■x-ray■docking pose
■md50ns

50nsMDでは結晶構造と同じ結合ポーズ
では結晶構造と同じ結合ポーズはでは結晶構造と同じ結合ポーズは現れない⇒ レプリカ交換MDなどサンプリング手法を検討する必要あり？などサンプリング手法を検討する必要あり？
レプリカ交換三井情報(株) 佐藤氏

分⼦動⼒学シミュレーションの課題(1)
シミュレーション時間
１ヶ月程度で到達可能なシミュレーション時間（タンパク質、水系）～100ns: 通常の計算機サーバ＜6～10μs:
GPU搭載型計算機（20日≒10 秒）～10μs: 京（並列計算の合計値）～1ms: Anton（データ通信がボトルネック）
Antonを用いた計算例
を用いた計算例
ユビキチンのfolding,
ユビキチンのunfolding
Non-native
native生命化学反応の時定数（ミリ秒～秒）

リガンドのタンパク質結合過程のシミュレーション
Non-native

native

ユビキチンのfolding,unfolding
ユビキチンのPiana et al., PNAS 2013, 110, 5915–5920,
Fig. 1 A

リガンドのタンパク質結合過程のシミュレー
リガンドのタンパク質結合過程のシミュレーション
Src kinase(チロシンキナーゼ)とdasatinibの結合
Shan et al., J. Am. Chem. Soc. 2011, 133,9181–9183,
Supplementary Material (1_si_001 (AVI))

Piana et al., PNAS 2013, 110, 5915–5920

理化学研究所 計算科学研究機構 荒木氏のスライドより

Shan et al., J. Am. Chem. Soc. 2011, 133, 9181–9183

リガンドのタンパク質結合過程のシミュレーション
リガンドのタンパク質結合過程のシミュレーション
Src kinase(チロシンキナーゼ)とdasatinibの結合
Shan et al., J. Am. Chem. Soc. 2011, 133, 9181–9183,
Supplementary Material (1_si_001 (AVI))

Src kinase(チロシンキナーゼ)とdasatinibの結合

分⼦動⼒学シミュレーションの課題(2)分⼦数タンパク質のシミュレーションの場合
システム構成：タンパク質：１個 （水：104-105個）

アボガドロ数（1023個)

＜＜

…
…
…
:

:

:

:

:

:

…
～10nm

…

１分子シミュレーションから算出された物理量で熱力学量を推定しなければならない…

京：88000ノードを使用して、1mMのタンパク質溶液1mLの挙動を1秒間シミュレートしようとすると、、、⇒ 1s / (10ns/day) ×1017個/105個 = 1020年かかる!!
年かかる
理化学研究所 計算科学研究機構 荒木氏のスライドより

準安定構造からの長時間MDで最安定構造にたどり着けるか？
で最安定構造にたどり着けるか？
準安定構造からの長時間
⻑時間MDによる化合物結合ポーズの追跡(CDK2_CS18)

結晶構造 (equil_2)

docking rank10 (equil_3)

水素結合数55000

Protein-Ligand energy (kJ/mol)

0100100

エネルギー値が異なる
200
20000.80.2
1μs

0.40.6

エネルギー値が異なる。
0.81μs00.20.40.6

x_ray
rank10

結合ポーズは1μs後も同じにならない⇒正しい結合ポーズ（最安定構造）を予測するには温度レプリカ等の拡張サンプリングが必要
理化学研究所 計算科学研究機構 荒木氏のスライドより

トラップからの脱却
通常のMDシミュレーションでは、タンパク質の構造がエネルギー極⼩の領域にトラップされることがある( )エネルギーが上がれば、トラップから抜け出すことができる原⼦の速度を上げる、温度を上げるカノニカルアンサンブルでのエネルギー分布
温度一定のMD.(>)

高エネルギーでの探索
トラップから抜けられるが、求めるエネルギー(温度)での構造は得られない

この分布で探索したい>
このエネルギー分布内で座標が移動
なかなか高エネルギーにはならない

この分布での構造が欲しい拡張サンプリング
マルチカノニカル法.(>)

カノニカル(低温)

カノニカル(高温)

マルチカノニカル>

一様のエネルギー分布空間において、構造探索を⾏う

レプリカ交換法.(>)>
複数のエネルギー分布でMDを実⾏し、それらの間で温度または構造の交換を⾏いつつ、構造探索を⾏うマルチカノニカル法.(>)

一様分布になるように、人⼯的なエネルギー分布.?@ (>)を定義.?@ > = A > <?@ > =定数<?@ > = A > B

A(>): 状態密度>

系の状態密度に反⽐例するように、マルチカノニカル重み因⼦<?@ (>)を定義

任意の温度(, におけるマルチカノニカルポテンシャルエネルギー>?@ >; (, = &' (, ln A(>) = (, G(>)
マルチカノニカル法の運動方程式>?@ >; (,EFC =>ニュートンの運動方程式
C = EC※事前に、短いMDでマルチカノニカル重み因⼦<?@ (>)を決定しておく必要がある「生体系のコンピュータ・シミュレーション」, 岡崎進, 岡本祐幸 [編], 化学同人, 2002

レプリカ交換法
複数の温度で、MD計算を並列して実⾏する温度を交換する

交換成⽴(((2(H

同じ系のレプリカ
温度が異なるシミュレーション

交換不成⽴

シミュレーション時間

レプリカ交換法の手順
1. 互いに相互作用しない系のコピー(レプリカ)を複数用意する
2. 各レプリカで、異なる温度において独⽴にカノニカルシミュレーションを実⾏する
3. 温度値が隣接した二つのレプリカを、詳細つり合いの条件に基づいて交換する∆= S? − S>

T−>

1
遷移確率I J → J + = Lexp(−∆)∆≤ 0∆> 0

4. 2. 3.を繰り返す。
事前に重み因⼦を決定する必要がない。
独⽴したMDシミュレーションを複数実⾏するため、並列コンピュータに適する。
「生体系のコンピュータ・シミュレーション」, 岡崎進, 岡本祐幸 [編], 化学同人, 2002

実例

• KBDDプロジェクト– CSAR CHK1_89 re-docking #1からの、X線結晶構造の探索

Docking計算法で予測された結合ポーズとX線構造との比較（CSAR: CHK1_89の例）の例）
反転■x-ray■docking pose
三井情報(株) 佐藤氏

直行並行

様々な温度でのエネルギー分布〜通常MD(NPT)でのtotal energyを使用〜

lnP(E)

298.0 299.4 300.8 302.3 303.7 305.1 306.6 308.0 (K)系のエネルギー(kJ/mol)lnP(E)

エネルギー分布が重なるようにレプリカの温度を設定する必要がある
298K（1.4K刻み）系のエネルギー(kJ/mol)
理化学研究所 計算科学研究機構 荒木氏のスライドより

398K

レプリカ交換法による、タンパク質-化合物相互作用エネルギーの経時変化
相互作用エネルギー(kJ/mol)

タンパク質-化合物相互作用エネルギーの経時変化
通常MD（25℃）

1μs
相互作用エネルギー(kJ/mol)

シミュレーション時間（ps）

T-REMD（ 25℃のレプリカ）
広いエネルギー空間が探索できる

13ns
シミュレーション時間（ps）

理化学研究所 計算科学研究機構 荒木氏のスライドより

準安定構造からの温度レプリカによる構造サンプリング
T-REMD 91レプリカ(25℃～175℃)

for CDK2_CS18

主成分分析結果(25℃, replica#0)
4

ドッキング構造からの通常MD（
ドッキング構造からの通常 （1μs)

3X線構造からの通常線構造からの通常MD（線構造からの通常 （1μs)

2
1PC10-4-3-2-10

12345-1

X線構造からの線構造からのremd(25℃℃, 0-5ns)
線構造からの-2-3

ドッキング構造からのremd(25℃℃, 0-13ns)
ドッキング構造からの-4-5PC2
通常MDよりは広い構造空間を探索できており、X線構造のリガンド配座近傍の探索は行われているものの、たどり着いてはいない（現在温度幅を拡大して検証中）
理化学研究所 計算科学研究機構 荒木氏のスライドより

バイオグリッドHPCIプロジェクト「新薬開発を加速する「京」インシリコ創薬基盤の構築」
プロジェクト「新薬開発を加速する「京」インシリコ創薬基盤の構築」
バイオグリッド
KBDD (K supercomputer-based drug discovery project by Biogrid pharma consortium)「京」産業利用枠: 新薬開発を加速する「京」インシリコ創薬基盤の構築

申請主体: NPO法人バイオグリッドセンター関西,研究代表者 京都大学医学研究科 奥野恭史
製薬企業(２３社): アステラス製薬(株)、アスビオファーマ(株)、エーザイ(株)、小野薬品工業(株)、科研製薬(株)、(株)カネカ、キッセイ薬品工業(株)、杏林製薬(株)、協和発酵キリン(株)、参天製薬(株)、塩野義製薬(株)、千寿製薬(株)、大正製薬(株)、大日本住友製薬(株)、田辺三菱製薬(株)、帝人ファーマ(株)、東レ(株)、日産化学工業(株)、日本新薬(株)、日本たばこ産業(株)、マルホ(株)、MeijiSeikaファルマ(株)、持田製薬(株)IT企業(2社): (株)京都コンステラ・テクノロジーズ, 三井情報(株)大学等: 京都大学大学院医学研究科, (独)産業技術総合研究所, 理研HPCI企画調整グループ

謝辞
KBDDコンソーシアムメンバー
NPO法人 バイオグリッドセンター
奥野恭史 教授, 荒木望嗣 様 (理化学研究所計算科学研究機構), J. B. Brown 先生 (京都⼤学医学研究科),
白石慧 様, 佐藤美和 様 (三井情報(株)), 広川貴次 研究チーム⻑ ((独)産業技術総合研究所),
⾦井千⾥ 主任研究員((株)京都コンステラテクノロジーズ)
製薬企業 (２３社):
アステラス製薬(株)、アスビオファーマ(株)、エーザイ(株)、⼩野薬品⼯業(株)、科研製薬(株)、(株)カネカ、キッセイ薬品⼯業(株)、杏林製薬(株)、協和発酵キリン(株)、参天製薬(株)、塩野義製薬(株)、千寿製薬(株)、⼤正製薬(株)、⼤日本住友製薬(株)、⽥辺三菱製薬(株)、帝人ファーマ(株)、東レ(株)、日産化学⼯業(株)、日本新薬(株)、日本たばこ産業(株)、マルホ(株)、MeijiSeikaファルマ(株)、持⽥製薬(株)IT企業 (2社):(株)京都コンステラ・テクノロジーズ, 三井情報(株)(独) 産業技術総合研究所
理研・計算科学研究機構
高度情報科学技術研究機構

http://sciencechannel.jst.go.jp/M120001/detail/M120001027.html

JST サイエンスチャンネルにて放送中