2015年11月30日 バイオスーパーコンピューティング神戸2015

神戸大学大学院システム情報学研究科
計算科学専攻 計算生物学講座
田中 成典







リガンドドッキングやタンパク質フォールディン
リガンドド
キングやタンパク質フルデ ン
グにおける長時間の構造緩和ダイナミクス平衡状態（ボルツマン分布）にどのように接近していくか？その数学的構造を議論
量子－古典対応に基づき、ランジュバン・ダイ
ナミクスを拡散量子モンテカル 法で解析
ナミクスを拡散量子モンテカルロ法で解析
ポテンシャル場に応じて、系のンタルピ 、ポテンシャル場に応じて、系のエンタルピー、エントロピー、自由エネルギーがどのように時間変化するか？

Structure--Based Drug Design (SBDD)
Structureタンパク質とリガンド(薬剤候補)の立体構造に基づく薬剤設計
近年の技術的進歩
 コンピュータの計算能力の急速な進化 タンパク質の立体構造の解明
 X線立体構造解析やNMRといった測定技術の進歩
 インフラの整備（SPring-8、京）タンパク質－リガンド系の分子シミュレーションによる定量的な設計が可能となった

計算機
計算機シミュレーションを用いたSBDDのメリットを用 たリ ト タンパク質とリガンドの相互作用（自由）エネルギーを計算できる
 立体的な結合構造が利用あるいは予測できるため、詳細な分子間の立体的な結合構造が利用あるいは予測できるため 詳細な分子間の相互作用を解析できる

ドッキング・シミュレーション（スコア関数）
ドッキングシミュレ ション（ス ア関数）


分子動力学シミュレーション（古典力場）


タンパク質のリガンド結合サイトにある薬剤分子の結合様式の分布や熱力学的性質を調べる。

量子化学的シミュレーション（電子状態）
量子化学的シミュレション（電子状態）

離れているリガンド分子がタンパク質のどのサイトにど ようなポどのようなポーズで結合していくか？
で結合して くかタンパク質と結合したリガンド分子の相互作用や化学反応の詳細を電子レ ルで解析する。
学反応の詳細を電子レベルで解析する。

さらに、よりマクロな、細胞・臓器レベルのシミュレ ションもある。
レーションもあるArtificial Bee Colony (ABC) アルゴリズム

AutoDockRas
Rasタンパク質と薬剤探索
タンパク質と薬剤探索

従来の抗がん剤

分子標的薬

成人の死因第一位であるがんの治療
 従来の抗がん剤は嘔吐や脱毛など大きな副作用を持つ

分子標的薬 の期待
分子標的薬への期待 がん細胞特有の分子を標的とするため副作用を抑えることが出来る
最も有望な分子標的：Rasタンパク質
がん 原因となるタ パク質 がんの原因となるタンパク質 最も普遍的ながん遺伝子産物と呼ばれている
最も普遍的な遺伝子産物 呼ばれる

Rasタンパク質を分子標的とした薬剤開発に注目が集まっている（神戸大学医学研究科・片岡研究室との共同研究）

Rasタンパク質とは？とは？
Rasファミリー低分子量Gタンパク質の一種



細胞増殖を制御する分子 イ チ 役割を持細胞増殖を制御する分子スイッチの役割を持つ



Rasタンパク質はGTPと結合している
質
結合



（GTPの構造式）

GTP結合型Rasには不活性型のState1と活性型のState2が存在する

Rasは不活性な状態と活性な状態を行き来する？

？

Rasタンパク質の特徴
Rasタンパク質の特徴
State1(不活性)突然変異によりState2が維持されState2(活性) ると細胞がガン化

薬剤が結合するポケットが存在する

薬剤が結合するポケットが存在しない

この状態遷移の機序の解明は分子シミュレーションの重要な課題



Molecular Mechanics (MM)

力場（force field）を利用 → エネルギー極小状態の探索

Molecular Dynamics (MD)

生体分子系を構成する全ての原子に対して古典的なNewtonの運動方程式を解くの世界→ μs ～ ms の世界へ
 F i = m ia i

リガンドドッキング・シミュレーション→ 有効に結合する阻害剤の探索「与えられたポテンシャル場での最適化問題」と位置づけることができる。
「与えられたポテンシャル場での最適化問題」と位置づけることができる
平衡分布（ボルツマン分布）にどのように近づいていくか？

Transfer operator for probability density
pinto
Decomposition
eigenstates

Operator that propagates the probability density ρAnalogy
gy to variational pprinciple
p in qquantum chemistryy （LCAO-MO）
Expansion by basis functions
Eigenvalue problem(cf., tICA, Relaxation mode analysis)

U(r)) : potential
U(t ti l fieldfi ld
D: diffusion coefficient

Fokker-Planck-Smoluchowski方程式

平衡状態：ボルツマン分布（P. Faccioli et al., Phys. Rev. Lett. 97 (2006) 108101.）(Villin headpiece)

Classical-quantum correspondence
Imaginary-time
S h diSchrodinger
equationti
Diff iDiffusion+B
BranchinghiTime-independent Schrodinger eq.
Ground state with E0 = 0

Effective action: 分布関数 P(r,t) のpropagator

この非平衡過程をシミュレートできないか？
→ 拡散量子モンテカルロ（DMC）法

：基底状態t→∞ζ: 一様乱数様乱数[0,1]

平均ポテンシャルエネルギー

エントロピ
エントロピー

λ：未定乗数変分原理 → ボルツマン分布 P0

F0 は最小自由エネルギは最小自由エネルギーを与えるを与える（等号は（等号は） より）

結局、
拡散モンテカルロ法により、波動関数の時間発展をシミュレートすれば その長時間極限として基底状態がシミュレートすれば、その長時間極限として基底状態が得られ、量子－古典マッピングにより、それが自由エネルギー最小の分布（ボルツマン分布）を与える。
ギ 最小 分布（ボ ツ分布）を与える。

[0, 1]の一様乱数はメルセンヌ・ツイスター法で生成
正規乱数分布χはボックス・ミュラー法で生成
kBT = 1.0, D = 1.0, Δt = 10-4, M = 20,000 walkers( )(0) = 0x(m)0.00 ((m = 11,2,2 …, M)
に対する波動関数の時間変化

kBT = 1.0, D = 1.0, Δt = 10-4,
M = 20,000
20 000 walkers
lkx(m)(0) = 22.00 (m = 11,2,2 …
…, M)
に対する波動関数の時間変化

kBT = 0.1,0 1 D = 1.0,
1 0 Δt = 10-44,
M = 20,000 walkers
x(m)(0) = 2.0 (m = 1,2, …, M)
に対する波動関数の時間変化

kBT = 1.0, D = 1.0, Δt = 10-4,
M = 20,000
20 000 walkers
lk[-5, 5] に一様に広がった分布からスタートした波動関数の時間変化
タ トした波動関数 時間変化

エントロピー減少kBT = 1.0, D = 1.0, Δt = 10-4,
M = 20,000,walkersx(m)(0) = -1.0 (m = 1,2, …, M)
に対する波動関数の時間変化

D = 1.0, Δt = 10-4, M = 20,000 walkers
初期配置 x(m)(0) = -1.0
1 0 ((m = 11,2,2 …, M)
で温度を変えてみる。（kBT = 0.1 – 1.0）

Walkersの1/3以上がW
lk
1/3以上が初めて x > 0 に移る
ステップ数を I とする。

アレニウス型で良くフィット
Δ = 0.57（U(x)のポテンシャル障壁は 1.0）

kBT = 1.0, D = 1.0, Δt = 10-3,
M = 20,000 walkers

U( ) Veff(x),U(x),( ) Φ0(x)の形状( )の形状

x(m)(0) = -1.0 (m = 1,2, …, M)
に対する波動関数の時間変化（3411ステップで全体の1/3のwalkersが x > 0 に移動）

エントロピーは最初増加、しばらくントロピ は最初増加、しばらく安定した後、分布が x > 0 に移行するのに伴い、減少する。
分布の x > 0 への移動はエンタルピーの減少が駆動→ ordered state

遷移過程で、ポテンシャルエネルギー（エンタルピー）は大きくゆらぎながら 平均は大きくゆらぎながら、平均としては減少ポテンシャル U(x,y) と基底状態 Φ0(x,y)=C exp[-U(x,y)/2kBT] の形状

U(x,0), U(0,y), Φ0(x,0), Φ0(0,y)の形状(x,y)=((xy)=(-11,0)に浅いminimum0)に浅いminimum(x,y)=(1,0)に深いminimum

基底状態 Φ0(x,y) のy積分、x積分を行って、x軸、y軸に射影したもの→x
x=1,1, yy=0に分布が集中0に分布が集中

kBT = 1.0, D = 1.0, Δt = 10-3, M = 10,000 walkers

700ステップ

800ステップ

kBT = 1.0, D = 1.0, Δt = 10-3, M = 10,000 walkers

900ステップでの射影 Φ(x), Φ(y)1100ステップ

1000ステップ

1200ステップ

1300ステップ

2000ステップ

1400ステップ

500ステップ、10万ステップでの射影
Φ(x), Φ(y) と基底状態の比較

等温古典確率系緩和ダイナミクス ⇔ 量子系基底状態探索
 Langevin eq., Fokker-Planck eq. ⇔Schrodinger eq., Feynman path integral⇔ Diffusion Monte Carlo U, S, F の時間的変動の追跡
 Entropy dominance vsvs. Enthalpy

急峻なポテンシャルの谷間に落ちて行くとき、エントロ
ピ の減少が見られる（ピーの減少が見られる（explicitly
li itl illillustrated）。
t t d）

多自由度系への拡張も容易 有限温度のoptimizationの一般的枠組