1. No category

生体機能における化学反応の役割
エネルギー変換、機能制御の主役
物理化学特論
膜間プロトン
ポテンシャル
光エネルギー
化学的エネルギー
• 発色団の光異性化
• プロトン移動
林 重彦
力学エネルギー
（Fo及びγの回転）
ADP+Pi
膜間プロトン
ポテンシャル
化学的エネルギー
（ATPの合成）
ATP
京都大学大学院理学研究科
化学反応
Catalysis: Stabilization of
Reaction Transition State (TS)
Lock-and-Key mechanism for substrate in TS
機能発現過程
電子状態変化
タンパク質の構造変化
局所的
大域的
反応基質とタンパク質
環境の相互作用
TS Stabilization in Enzymes
Ka of various host-guest complexes
Enzyme-TS
Ka of TS substrate
kcat
log K a   log
Kmkuncat
Exceptionally large
Reactant
TS
Product
Pauling Nature (1948)
“enzymes are molecules that are complementary
in structure to the activated complexes of the
reactions that they catalyze, ..., [rather than]
entering into reactions”
Zhang, Houk, ACR (2005)
Enzymesubstrate
Zhang, Houk, ACR (2005)
Catalysis: Stabilization of
Reaction Transition State (TS)
Molecular Motor Function Controlled by
Enzymatic Activity
F1Fo ATPsynthase
ADP
Lock-and-Key mechanism for substrate in TS
local
Pi
~10 kcal/mol
ATP
Reactant
TS
Product
Preorganization:
Binding pocket structure of enzyme
is prepared to stabilize substrate in
the TS conformation
reversible!
Noji et al.
motor
Kamerlin & Warshel, Proteins (2009)
Protein dynamics?
Noji, Takeuchi, et al.
global
織り込まれる様々な現象
‐F1 分子モーターを例に‐
分子モーター
化学反応
の因果律
織り込まれる様々な現象
‐F1 分子モーターを例に‐
分子モーター
化学反応
の因果律
サブユニットの
構造変化
サブユニットの
構造変化
複合体の
構造変化
複合体の
構造変化
含まれる自由度の数
mechanical
ATP generator
合成酵素の
因果律
含まれる自由度の数
量子
力学
古典
力学
合成酵素の
因果律
分子シミュレーションの手法
分子シミュレーションの手法
電子
電子
MD
小分子
静的（ps）
タンパク質一分子
~ s
複合体
> s
シミュレーション効率（系の大きさ・計算時間）
分子シミュレーションの手法
ギャップ
QM
化学反応
MD
粗視化
小分子
静的（ps）
タンパク質一分子
~ s
Methodologies
電子
分子
MD
モデル
分子の記述の精度
• Empirical force fields (MM)
• No explicit description of electronic states
QM/MM
モデル化
粗視化
小分子
静的（ps）
タンパク質一分子
~ s
複合体
> s
シミュレーション効率（系の大きさ・計算時間）
複合体
> s
シミュレーション効率（系の大きさ・計算時間）
1. Molecular dynamics (MD)
simulations
古典力学
• Large protein systems
QM
モデル化
モデル
モデル
粗視化
分子
分子の記述の精度
分子
分子の記述の精度
QM
2. Molecular orbital or density
functional theories (QM)
• Explicit description of electronic states
• Huge computational cost
3. QM/MM simulations
• Reactive moieties are described by ab
initio QM methods.
• Protein environment is treated by MM.
Methodologies
Methodologies
1. Molecular dynamics (MD)
simulations
1. Molecular dynamics (MD)
simulations
古典力学
• Large protein systems
• Large protein systems
• Empirical force fields (MM)
• No explicit description of electronic states
• Empirical force fields (MM)
• No explicit description of electronic states
2. Molecular orbital or density
functional theories (QM)
• Explicit description of electronic states
• Huge computational cost
量子力学
2. Molecular orbital or density
functional theories (QM)
transition state
• Explicit description of electronic states
• Huge computational cost
量子力学
3. QM/MM simulations
3. QM/MM simulations
• Reactive moieties are described by ab
initio QM methods.
• Protein environment is treated by MM.
• Reactive moieties are described by ab
initio QM methods.
• Protein environment is treated by MM.
2013年ノーベル賞研究
SH and Ohmine (2000) (QM/MMopt)
SH, Tajkhorshid, Schulten (2003)
(excited state QM/MM-MD)
講義の構成
1. 分子間相互作用のおさらい
2. 量子化学（QM）の簡単な説明
3. 分子力場（MM）法の簡単な説明
4. QM/MM の説明
Martin Karplus
Michael Levitt
Arieh Warshel
“for the development of multiscale models for
complex chemical systems”
5. 酵素反応経路解析
6. タンパク質の柔らかさと酵素活性
7. QM/MM 法を用いたタンパク質機能設計
分子の相互作用
タンパク質構造と機能
部分電荷の間の相互作用
DNA塩基対の相補性
静電マップ
真空中のクーロン相互作用
V
ペプチドの部分電荷
q1q2
4 0 r
 ヘリックス
媒質中のクーロン相互作用
G
C(1997)
Journal of Chemical
Education, 74, 771
反応性
分子認識
分子集合体の物性
物性
H2O
q1
q7
電荷集団を代表する点 R
（重心など）を定義し、そこか
らのポテンシャルを考える
ri = ri + R
q2
q5
q6
q3
q4
(X)
X = (x, y, z)
点 X に電荷 Q があるとき
の相互作用エネルギー V
V  Q( X )
qi
( X )  
i 1 4 0 X  ri
qi, ri : 原子の電荷及び座標
 媒質の誘電率
r =  / 0 : 相対誘電率
例 : 水分子の r = 78
r の異なるクー
ロン相互作用
酵素活性
ポテンシャルの多極子展開
q8
q1q2
4 r
塩は熱では溶けにくい
水にはよく溶ける
CO2
複数の部分電荷が作るクー
ロンポテンシャル (X)
V
qi
i 1 4 0 X  R  ri
( X )  

μ  X  R
qtot

4 0 X  R 4 0 X  R 3

（|X – R| の逆数の高次の項）
距離が遠いと小さくなる
qtot   i 1 qi
全電荷
μ   i 1 qi ri 双極子モーメント
Journal of Chemical Education, 74, 771 (1997)
電気双極子
分子や基が遠く離れている
とき、部分電荷の集団として
現れる相互作用を考える
有極性分子 : 永久（部分電荷から
生じる）電気双極子モーメントを持つ
無極性分子 : 永久電気双極
子モーメントを持たない
電気双極子
電気双極子モーメント = ql
q
q
l : 部分電荷間の距離
単位 : デバイ （D）
1 D = 3.33564 x 1030 Cm
多原子分子の双極子モーメント
双極子モーメントはベクトル量
電荷分布と座標から双極子
モーメントを計算
 = (x, y, z)
分子の分極
外部電場 
双極子–誘起双極子相互作用
*
誘起双極子は永久双極子の向
きの変化に追随する
双極子モーメントの大きさ 
 = || = (x2 + y2 + z2)1/2
q1
q8
R
q7
q2
q5
q6
q3
q4
誘起双極子モーメント *
外部電場  により分子の電子状
態が歪み、無極性分子であっても
電気双極子が誘起される
分子内のある一点
R を基準にとる
* = 
ri = ri – R
多原子分子の双極子モーメント
μ   i 1 qi ri
分散相互作用
強い角度
依存性
X、Y : 窒素、酸素、フッ素
1. 分子内の電子の位置の揺らぎ
により瞬間的に双極子が生じる
静電モデル（最も単純な記述）
2. 相手の双極子を分極させること
により二つの双極子が相関する
1. 電子分布の揺らぎやすさ
2. 誘起分極の大きさ
相互作用エネルギー
（クーロン相互作用）
X – H ・・・ Y
分散相互作用 （ロンドン力）
相互作用の大きさは分極率
に依存する
水素結合
電気陰性度の高い原子の間に
水素原子が介在し結びつける
希ガス等のような無極性分子
同士でも弱い相互作用が働く
瞬間双極子の間の相互作用
 : 分極率
  x*    xx  xy  xz   E x 
 * 
 
 y     yx  yy  yz   E y 
  *  
 
 z   zx  zy  zz   Ez 
一般に異方性がある
ロンドンの式
II
2   
V    1 62  12
3
r
I1  I2
I1、I2 : イオン化ポテンシャル
電子求引性の X 原子に結合して
生じたプロトンの部分正電荷と、
原子 Y にある孤立電子対の部分
負電荷の間のクーロン相互作用
ルイス酸塩基複合体
X – H + :Y → X – H ・・・ Y
分子オービタル
 = c1x + c2H + c3Y
分子オービタル
反結合性
非結合性
結合性
Lennard-Jones 相互作用
量子化学（QM）計算
原子間の分散引力と近距離の斥力を併せて記述
分子の電子状態を計算する
VLJ (R ) 
Natom

i j
 
4 ij  ij
 Rij


  ij  



 
R

 ij  
12
分散力
パウリの排他原理により、近
距離で強い斥力がはたらく
波動関数の重なりなので本当は指数
関数が良いが、簡便な計算のために、
r-12 の斥力とすることが多い
量子化学計算の方法
原理・精度・計算コストが異なる
TD-DFT（励起状態）
電子相関○
分散力×
多配置×
電荷移動×
現在の趨勢
密度汎関数法（DFT）
電子相関○
分散力×
多配置×
CASSCF法
電子相関×
分散力×
多配置○
電子相関 : エネルギー
分散力 : 弱い相互作用
多配置 : 金属、励起状態
原子核
N
MP2 摂動法
電子相関○
分散力○
多配置×
CAS-PT法
電子相関○
分散力○
多配置○
他にも色々な手法あり
（SAC-CI とか）
基底状態には DFT 法（分散力が必要な場合は MP2）
励起状態には TDDFT 法（うまくいかないこと多し）か CAS/CASPT
Ze
N
Nelec
1 2 Nelec QM Z A Nelec 1 QM Z A ZB
ˆ
 
H    i   
i 2
i
A riA
i  j rij
A B RAB
電子の運動
エネルギー
近距離斥力
Hartree-Fock法
電子相関×
分散力×
多配置×
e電子
電子状態のハミルトニアン
6
電子
e-
原子核から 電子間反発 核間反発
の引力
様々な手法がある
• Hartree-Fock（HF） （変分法、平均場近似）
• CI、MCSCF （変分法、電子相関）
• MP2、MRMP、CASPT2 （摂動法、電子相関）
• CC、SAC/SAC-CI （クラスター展開法、電子相関）
• DFT （変分法、密度汎関数にて電子相関）
エネルギーの期待値の表式
基底関数  (r) の一電子・二電子積分で記述される
E   Hˆ 
エネルギーの期待値
NQM
Z A ZB
A B rAB
 Eel (D, d)  
NBF NBF
Eel (D, d)   h D 


一項目は電子のエネルギー
二項目は核間反発
1 NBF NBF NBF NBF
 g  d 
2    
NQM
 1
Z 
h   dr1* (r1 )   12   A   (r1 )
A r1A 
 2
1
g    dr1dr2* (r1 ) (r1 )  * (r2 ) (r2 )
 r12 
一電子積分
二電子積分
DFT 法は二電子項が異なる
密度行列
変分法と Roothaan 方程式
波動関数により与えられる電子密度の行列
エネルギーを変分パラメータで変分し波動関数を決定
エネルギーは電子密度行列要素を含む
NBF NBF
1 NBF NBF NBF NBF
Eel (D, d)   h D   g  d 
2    
 
d 二粒子密度行列
D 一粒子密度行列
独立粒子モデル
HF 法（閉殻）の場合
（スレーター行列式）
  Aˆ  (r ) (r )
(r
)
(r )
1
1
1
NBF
2
a (r )    (r )Ca

Nel /2
Nel 1
Nel /2
Nel
分子軌道
Nel / 2
1
d   D D  D D
2
a
Ca : LCAO 係数 これを計算で決定する
D  2  CaC*a
NBF NBF
Eel (D, d)   h D 

HF 法の場合
密度行列の微分
dD  CP法などで解
dEel (D, d) NBF NBF  dh
D  h
  

かれる
dRA
dRA 
   dRA
dd  
積分の微分
1 NBF NBF NBF NBF  dg 
d   g 




解析的に
dRA 
2      dRA
計算可能
重なり積分の微分
HF 法の場合
Nel / 2
dS 
dEelHF (D, d) NBF NBF  dh
D   2 aCaC*a
  

dRA
dRA 
a
   dRA

1
d 

2     dRA
C に関する非線形方程式なの
で反復的に解く（SCF）
F(C)C  SCε
NBF NBF
1


F  h   D  g   g   Fock 行列
2


 
S   dr1* (r1 ) (r1 )
重なり積分
a : 軌道エネルギー
CI や MCSCF 法では、CI 展開係
数が（も）、変分パラメータとなる
分子の形
核に働く力を計算することにより構造最適化などを行う
dg 
分子軌道 a で汎関数変分
Roothaan 方程式 : Ca に対する行列方程式（固有値問題）
核座標の微分
NBF NBF NBF NBF

1 NBF NBF NBF NBF
 g  d 
2    
Roothaan 方程式の微
分からの条件より密度
行列の微分がなくなる
分子のポテンシャルエネルギー e-
ZBe
電子
E  Eel   el (r1 , r2 ,  , rN ;R1,  R NZ )
NZ
Z Z e2
 A B
核間反発
電子のエ
A  B 4 0 RAB
R
ネルギー
 V (R1 ,  R NZ )
e電子
ZAe
ある与えられた核の位置の電子の波動
関数及びエネルギーを決定することによ
り、核座標の関数として与えられる
分子の形
ポテンシャルエネルギーが極小となる構造
エネルギー勾配法 g(R ,  R )  V (R1,  R N )
1
N
R
により決定される
Z
Z
平衡結合長
エネルギー成分
分子間相互作用 : エネルギー分割法
分子間相互作用を分割して特徴をとらえる
Kitaura and Morokuma
二つの分子 A と B の間の
分子間相互作用エネルギー
分子A
VA B  E AB  E A0  EB0
複合体（超分子）
のエネルギー
分子B
HF 波動関数に対して考える
孤立分子の HF エネルギー E A0   0A Hˆ A  0A
孤立分子の
HF 波動関数
0
0
E   Hˆ B  0B E0  E A  EB
複合体分子の HF エネルギー E AB   AB Hˆ AB  AB
Hˆ AB  Hˆ A  Hˆ B  VˆAB
0
B
0
B
各寄与の水－水間
距離に対する依存
性の摂動論による
見積り
分子A
分子B
A
B
反対称化演算子
交換を許す
CT
交換反発エネルギー EEX  E3  E1 PL
（4）   Aˆ   *  *  相手分子の空
4
A B
BA
軌道を混ぜる
PL
ES
EX
電荷移動エネルギー ECT  E4  E3
VAB  EES  EPL  EEX  ECT  EMIX
分散力 : 量子的な誘起双極子相互作用（HF法では記述できない）
タンパク質中の水分子の相互作用
エネルギー分割の例
水２量体
0
0
（1） 1   A  B E1  1 Hˆ AB 1
静電エネルギー EES  E1  E0
分子B
（2）  2   A  B 分子 B（A）の影響下での A（B）の HF 波動関数
分極エネルギー EPL  E2  E1
（3）   Aˆ   0  0  A と B 間の電子の
3
孤立分子の
エネルギー
分子A
直積
bacteriorhodopsin
FT-IR (Kandori et al.)
Ab initio QM/MM normal mode analysis
K-BR
Kandori, D2O
O-D and N-D
?
水‐水間距離
水素結合が非常に強い
Wat402 とその周りとの相互作用の分割 (kcal/mol, DZV(O+))
Escf
EES
EEX
EPL
ECT
EMIX
-27.4
-64.7
62.8
-14.1
-18.2
7.9
分子力学（MM）
分子動力学（ＭＤ）シミュレーション
• 簡便な分子力場（MM）ポテンシャ
ルエネルギー関数を用いて、系の
エネルギーと力を計算する。
• Newton 方程式に従って、原子
の時間発展を計算する。
• アンサンブル平均などを計算し統
計量などを得る。
水
構成する原子の座標 R
= (R1, R2, …, RNatom)
の簡単な関数
分子力場 （MM force field）
V MM (R ) 

1 r
kij Rij  Rij0
2
bond ( ij )
Nbond


Nangle

2
1 
kijk  ijk  ijk0
angle ( ijk ) 2



2
分子内相互作用
分子の結合による
相互作用
計算は
n
N
N
V
QM に比 
  ijkl 1  cos  nijkl   ijkln 
べると圧 dihedral ( ijkl ) n( dihedral ) 2 
倒的に軽 N
 A
B 
分子間相互作用
   ijLJ  ij12  ij6 
い
非結合性の相互作用
dihedral
d
atom
タンパク質の
折りたたみ
AMBER 
force field
MMの分子内相互作用ポテンシャル
結合距離 Rij、結合角 ijk、及び二面角 ijkl の簡単な関数
MM
Vintra
(R ) 


1 r
kij Rij  Rij0
2
bond ( ij )
Nbond

Nangle


2
1 
kijk  ijk   ijk0
angle ( ijk ) 2


Ndihedral
Nd


dihedral ( ijkl ) n ( dihedral )
n
Vijkl
i

k
ijk
Rij
2
j
ijkl
1  cos nijkl  
2 
n
ijkl
Rij 
 qi q j 

  Rij 
Natom
  ijcl 
i j
解析的に微分可能
• 安定構造の最適化
• 分子動力学計算
MMの分子間相互作用ポテンシャル
 Aij
Bij  レナードジョーンズ
 12  6  （LJ）相互作用
Rij 
i j
 Rij
分散力と交換斥力
Natom
 qi q j 
cl
   ij 
 クーロン相互作用
i j
  Rij 
 : 近接原子間の相互作用を除く
MM
(R ) 
Vinter
Natom

LJ
ij
ij
l

 Rij
i j
生体膜中の
膜タンパク質
（ロドプシン）

結合距離及び結合角のポテンシャルは、極小値からの二次関数
平衡構造の周りの微小の変位しか考慮していない
二面角のポテンシャルは、Rij や ijk と相関しない少数のフーリエ級数
複数の（準）安定構造の間の相対的安定性（ペプチド高分子）
通常のシミュレーション
では電荷は変化せず
分極相互作用の取り込み
• 原子点誘起双極子モーメント
• 結合点誘起双極子モーメント
• 電荷平衡法 （fq-SPC）
• Charge Response Kernel
電荷移動相互作用 まず無理
ほぼ全ての原子の対の
相互作用を計算
MM 計算の律速
高速化
計算量 ~ Natom2
• カットオフ（精度がよろしくない）
• エワルド法、PME 法
• 多重極展開法
分子力場パラメータ
化学的に異なる結合は
異なる性質を持つ
分子力場パラメータ（続き）
Atom Type : 化学的に同種の原子を定義
V MM (R ) 
とはいえ、例えばタンパ
ク質の場合、アミノ酸の
数（20 種類）は限られて
いるので、それらの結合
を網羅したパラメータを
決定するのは可能


1 r
k ij Rij  Rij0
bond ( ij ) 2
Nbond


Nangle


Ndihedral
Nd


V MM (R )
 Fj (R )  
R j
ポテンシャルの原子
座標に関する微分は
解析的に求まる
数値積分 velocity Verlet 法
簡便さと丸め誤差
r (t   t )  r (t )   tv (t )  1/ 2 t 2a(t )
の精度によって
v (t  1/ 2 t )  v (t )  1/ 2 ta(t )
様々な手法がある
v (t   t )  v (t  1/ 2 t )  1/ 2 ta(t  1/ 2 t )
t + t
n
Vijkl


n 
1  cos nijkl   ijkl

2 
計算時間
安定な時間発展計算のための時間刻み  t は、
1~2 fs （物理的要請） 分子振動の周期の数十分の一
Newton 方程式
t - t t
2
Natom
q q 
   ijcl  i j 
i j
  Rij 
系の原子の時間発展を計算
dt 2

Natom
 A
B 
   ijLJ  ij12  ij6 
Rij 
i j
 Rij
Newton 方程式の時間積分
mj
2
1 
k ijk  ijk   ijk0
angle ( ijk ) 2
dihedral ( ijkl ) n ( dihedral )
AMBER force field
JACS 117, 5179-5197
(1995)
d 2R j (t )

t - t t
t + t
t - t t
t + t
t - t t
t + t
100 days with 1 fs / step = 8.64 Ms
Trajectory time
10 ps
QM-MD
GPU
Anton
864 s = 14.4 m
100 ps
86 s
1 ns
8.6 s
10 ns
MM-MD
Computational time/ step
0.86 s
100 ns
12 step / s
1 s
120 step / s
1 ms
120k step / s
MM <<< QM
Anton:
A Game Changing Technology
Molecular dynamics simulations (MM-MD)
Voltage gated K+ channel
David E. Shaw
系の分子集団の統計的振る舞いをどのような条件で
記述するか
温度の定義
運動エネルギー K(R) と
等分配の法則より
• Computer science at Columbia
Univ.
• Hedge fund “D. E. Shaw & Co.”
• “D. E. Shaw Research”:
development of a special
purpose computer for MD
simulation, “Anton”
K (R ) 
1 N
3N
mi R i 2 
kBT

2 i
2
エネルギー一定アンサンブル
• Newton 力学なのでエネルギーは一定
• 相転移などエントロピーが大きく変化するときに不具合
温度一定アンサンブル
p  F   p
全く異なるアーキテクチャ
極低レイテンシ: strong scaling
系の温度とアンサンブル
分子機能が見えてしまった
統計平均
ターゲット温度となる
ようにコントロール
圧力一定アンサンブル
Berendsen の方法
Nose-Hoover の方法
ターゲット圧力となるように系の大きさを変える
MD シミュレーションでの統計平均
カノニカルアンサンブル
分布をどのように計算するか？
温度一定アンサンブル、あるいは注目している系の周りの熱浴と
みなせる系が十分大きいとき、注目している系の振る舞いはカノ
ニカルアンサンブルで与えられる（と期待される）。
原理的には : 非常に多くの独立の系をシミュレーションしてその
間で平均をとる、、、が普通の系ですでに計算上ムリ。
カノニカル分布
系が平衡状態のときには、時間平均をする
位置 x と運動量 p の分布 P(x,p)
  N 1

P ( x,p)  Z 1  exp     
pi 2  V ( x )  
  i 1 2mi
 
  N 1

Z   dx  dp exp     
pi 2  V ( x )  


  i 1 2mi
 


1
kbT

分配関数
Nt
f  Nt 1  f  x (t i ) 
i 1
時間相関関数
（Time Correlation Function）
N
統計平均




f   dx  dpf ( x,p)P ( x,p)
任意の関数 f(x,p) の平均
f (t )f (0)  N 1  f  x(t   i )  f  x ( i ) 
i 1
タンパク質の揺らぎ
液体の構造・動径分布関数
ある原子から別の原子まで
の距離の分布
Ni2+
in water 井内ら
根平均自乗揺らぎ
V. Spiwok. J.Mol.Model. 2007. 13. 485
Root Mean Square Fluctuation
揺らぎやすさをみる
RMSFi 
ri 2
ri  ri  ri
分散共分散行列
水の動径関数
揺らぎの相関と大きさをみる
r
n(r )    4 r 2g (r )dr 
0
配位数
i 原子と j 原子の運動が相
関していなければ Cij = 0
Ni2+ に H2O が 6 配位
拡散係数と速度相関関数
三次元の拡散
r (t )  r (0)
2
D: 拡散係数
 6Dt
二乗変位が時間に比例
Cij  x i x j
CS2 の自己拡散
対角化すると主成分が求まる
Methodologies
1. Molecular dynamics (MD)
simulations
古典力学
• Large protein systems
• Empirical force fields (MM)
• No explicit description of electronic states
速度相関関数による
拡散係数の表式
2. Molecular orbital or density
functional theories (QM)
1 
D   v( )  v(0) d
3 0
3. QM/MM simulations
大雑把な見積りには便利だが、
積分の収束性に若干問題あり
青：AHA、赤：PPA
• Explicit description of electronic states
• Huge computational cost
量子力学
• Reactive moieties are described by ab
initio QM methods.
• Protein environment is treated by MM.
SH and Ohmine (2000) (QM/MMopt)
SH, Tajkhorshid, Schulten (2003)
(excited state QM/MM-MD)
QM/MM ハミルトニアン
1
Hˆ QM / MM     i 2 
i 2
Nelec
i
V
ZA
r
i
A

Nelec
iA

ij
NQM NMM
Nelec NMM

Nelec NQM
q
Z q
a r a  A a rA a
ia
Aa
MM
QM MM
V
MM
MM
QM–MM 間の静電的相
互作用を除く相互作用
• LJ相互作用
MM領域内
• 分子内相互作用
の相互作用
気相中のQM 分
Z Z
1
  A B 子の電子状態の
rij AB rAB
ハミルトニアン
QM–MM 間の静電的相互作用
• 第一項は電子と MM 原子の有
効点電荷との相互作用で、一電
子演算子（核引力項に類似）
• 第二項は原子核と MM 原子の
電荷との相互作用で、古典的
クーロン相互作用関数で表現
NQM
QM/MM エネルギー
QM / MM
Etotal
  Hˆ QM / MM 
 EelQM / MM (D, d)
NQM
QM-MM 間の相互作用
1. 静電相互作用
- QM 分子の電子と MM の電荷が直接相互作用する
- MM の静電場の中での QM 電子状態が決定される
電子と無関係の項
NQM NMM
Z q
Z A ZB
MM
MM
   A a  VQM
MM  VMM
r
r
A B
A
a
AB
Aa

NBF NBF
QM / MM
EelQM / MM (D, d)   h
D 


NMM
QM / MM
h
 h    dr1* (r1 )
a
一般に、狭義の意味での QM/MM 法では、（林の解釈では）MM 原
子の有効点電荷と QM 分子の電子の間の静電相互作用を取り込
み、環境の静電場の影響を QM 分子の電子状態に反映させる
QM 分子の電子状
態のエネルギー
気相中のエネルギー
との違いは一電子積
分項に現れる
1 NBF NBF NBF NBF
 g  d 
2    
qa
 (r1 )
r1a
MM 電荷との
相互作用の
一電子積分
一電子積分なので計算
コストが余りかからない
QM 領域の取り方
一般的な QM/MM 系
QM region
2. 分極相互作用
- QM 分子の電子は MM の静電場で分極する
- MM 分子は通常は点電荷のまま分極しない
（分極力場を用いることも原理的には可能だが計算が大変）
3. 交換反発と分散相互作用
- MM 力場のレナードジョーンズ相互作用にて計算
4. 電荷移動相互作用
- まったく記述できていない
大雑把に言って、分極と電荷移動相互作用が効きそうな分子間相
互作用がある場合は、すべて QM 分子に含めるのが望ましい
分散相互作用は QM の記述でも難しい、、、
反応が起こっている部分の周りの一層は取りたい
静電相互作用の重要性
タンパク質内反応では、反応中心が強く分極してい
ることが多い
非常に強く分極した反応中心を安定化
-
+
-
反応中心クラス
ターモデル
+
+
-
電荷移動のartifact
分極安定化の相互作用は意外
と長距離
… +
-
+
-
+
-
…
タンパク質環境がないと
プロトン移動してしまう
分極した反応中心を安定化して
いる分極した分子をさらに安定
化している分極した分子が、、、
QM 分子の外側は MM
で囲んだ方が良い
Link Atom Approach
原子を置いて安定な QM 分子にする
主に水素原子（H）
パラメトライズされた ECP を置く手法もある
問題点
• link atom と MM 電荷の相互作用
• MM 電荷と QM 分子との距離が近い
境界の QM 原子から 1–2 及び 1–3 相
互作用をする MM 電荷との相互作用を
除いてしまう
• そもそも QM-MM 境界にまたがる 1–2
や 1–3 等の相互作用を定義できない
（後述）
QM-MM境界に結合がある場合
一般的コメント
• 良い方法は余りないので、反
応中心からなるべく離れたとこ
ろに境界を持ってくる
• 反応中心のすぐそばに境界を
取らなければならないときは、
精度の悪さに左右されない量で
議論をする
しなければならないこと
境界で結合を切ることにより生じる価電子の穴を埋める
• 原子を足す（link atom）
• localized orbital などを用いる（GHOとか）
結合境界での分子内相互作用
二つの方法がある
• すべて MM 分子内相互作用関数で記述する（link
MM k
atom は何も関係させない）
最も単純な方法
QM
LA
• link atom との結合を、境界の結合と関係づける
（距離は比例、方向は両者の結合を合わせる）
より丁寧な方法（ONIOM で採用）
若干の注意点
• link atom との結合と境界との結合が同じでない
• link atom との結合の構造が QM の性質に大きく関与する
場合は、境界の取り方がそもそもあまり良くない場合もある
l
構造最適化
構造最適化のマイクロイタレーション
核座標の微分が求まるので構造最適化が可能
• QM 領域の構造最適化のイタレーションと MM 領域の構造最
適化のイタレーションを交互に行う
• MM 領域の構造最適化では、QM 分子との相互作用を古典
的なクーロン相互作用で行う。その際、QM 分子の部分電荷は
ESP 法などで決める
Zhang et al. (2000)
問題点
MM 原子座標の微分に一電子積分の微分が現れてくる
N
QM / MM

dEtotal
d  QM NMM Z Aqa
MM
MM

 VQM
 
MM  VMM 
dRa
dRa  A a rAa



 d qa 
    dr1* (r1 ) 
  (r1 ) D
  

 dRa r1a 

NBF NBF
MM 領域は一般的に原子数が多く、QM 領域より非常に
多くの最適化ステップが必要なため、積分計算が中に入
ると計算時間が非常に大きくなる
QM 構造最適化
（MM構造は固定）
静電相互作用は一
電子積分で計算
QM分子の電子密度に
基づき部分電荷を決定
• An operator form of ESP derived charges
• Electronic embedding in the classical Coulombic
representation
Ngrid
 1t a 1W(d)  Nel (d) 
1
V  (d)
q(d)  a 1W(d)  a 11
 Wa  
t 1
1a 1
 ra


Fock operator
Periodic boundary condition
NQM NMM
q
f QM / MM (d,R, X )  f 0 (d,R )      ap p qˆa (R )
rap
a
p
• Easy to implement Ewald summation
technique for PBC in QM/MM
 Long-range interaction
 Use of sophisticated MD programs
静電相互作用は点
電荷間のクーロン
相互作用で計算
問題点 : QM と MM の構造最適化に用いる
ハミルトニアンが違う
RESP Charge Operator
Ten-no et al. for RISM-SCF, SH and Ohmine for QM/MM
MM 構造最適化
（QM構造と部分電荷は固定）
分子の運動と反応
分子の運動
配置空間のポテンシャルエネルギー
曲面上の運動
ポテンシャルエネルギー曲面（ 3 N 次元）
E(R) 
ˆ (r;R )  (r;R)
 ele (r;R) H
ele
ele
 ele (r;R )  ele (r;R )
ele
ele
 ele (r;R) 電子波動関数（BO 近似）
分子運動（Newton 方程式）
   E(R)
Mi R
i
Ri
配置空間のトラジェクトリ
反応：始状態と終状態をつなぐ分子運動
反応速度の定義
反応ダイナミクス
ポテンシャル曲面の地形が反応ダイナミクスを決める
反応速度の厳密な定義（ミクロカノニカル集合）
F + H2 → FH + H
反応分割面を通るトラジェクトリの流束
F (E )  h3N  dqdp E  H(q, p) F (q, p) (q, p)
ミクロカノニカル 分割面での流束
f (q) p
F (q, p)    f (q)

q m
 (q, p)
反応性・エネルギー移動
自由エネルギー曲面
カノニカル集合での反応系を考える
反応に関わる座標（R）とそれ以外の座標（X）に分割する
F (R)  kBT ln  dX exp    E(R, X)
反応座標 s(q) を導入する
F  kBT ln  dR exp    F (R )
F (sr )  kBT ln  dq  s(q)  sr  exp    E(q)
反応座標方向の分布
(sr )  e   F ( s )
r
反応は自由エネルギー basin
（準安定状態） の間の遷移
終状態
始状態と終状態をつなぐトラジェ
クトリなら 1、そうでなければ 0
反応速度
F (E )
k(E ) 
nR (E )
n (E )  F (E )
R
反応分割面
（3N-1 次元）
f (q)  0
反応方向決定
始状態
状態密度
配置空間（ 3N 次元）
遷移状態理論（TST）
仮定
1. 反応始状態と終状態を分ける遷移状態が定義できる
2. 遷移状態は始状態と熱平衡にある
3. 遷移状態に反応分割面があり、それを超えて終状態
側に入ったトラジェクトリは再交差しない（必ずそのま
ま反応する）
帰結
1. 熱平衡の仮定により反応分割面での流束を求める際
の確率密度がエネルギー論で求まる。
2. 再交差禁止の仮定により流束が遷移状態において
のみの情報で求まる（全配置空間は見なくて良い）。
反応速度が生成物と遷移状態の情報だけで求まる
遷移状態
遷移状態理論の反応速度
遷移状態 : ポテンシャル曲面上の鞍点
ミクロカノニカル反応流束
p
FTST (E )  h3N  dqdp H(q, p)  E   qr  qr*  r  (pr )
m
鞍点：
基準振動の振動数の一つが虚数の
不安定停留点
遷移状態
（TS）
• 振動数が虚数の方向が
反応方向（ある種の反応
E *
座標となる）
• それに直交する方向が分
終状態
割面で、鞍点が分割面で
始状態

は最もボルツマン分布で
q
qr
の確率密度が高い（分割 反応座標に
面の代表点となる）
反応座標
直交する座標
（3N-1 次元）
反応速度
始状態
遷移状態

調和近似・高温（古典）近似

kTST (T )  R e  E
2
R
*
終状態

kTST (T )  ZR 1  dEnR (E )kTST (E )e   E
0

kBT  F *
e
h
F *
*
Z TS e  E
e

ZR
活性化自由エネルギ
  F *
反応を律速する段階（RDS）：反応流束のボトルネック
*
kRA 
k AB 
kBC 
kCP 

 A 

 B 

 C 

P
R 




kRA 
E *
k AB 
kBC 
各ステップでの正反応の流束
kCP 
FIJ   [I ]kIJ 
TS
FA*
kBT  kBT
F
e
 BA
k
[B]
h


 e kBT
K 
FB*
k
[A]

kBT kBT
e
h
Heaviside
分割面
qr に直交する面 階段関数
多段階反応過程
R : attempt frequency
平衡定数
FTST (E )
nR (E )
カノニカル反応速度
TST 反応速度のあれこれ
一次元の反応速度
kTST (E ) 
A
A
R
B
B
どこが RDS?
C
P
多段階反応過程（続き）
多段階反応過程（続き２）
中間状態の自由エネルギーが始状態より大きい場合
中間状態の自由エネルギーが始状態より小さい場合
FIJ   [I ]kIJ   [R]e  FIR
*
kBT  FIJ* 
kT
 [R] B e  FIJR
e
h
h
始状態からみた遷移状態の自由エネルギーが最も高
いところが律速段階となる
A
C
R
B
平衡が移動する
A
C
P
TST の仮定の一つ：始状態と遷移状態が平衡にある
遷移状態より前の最も低いエネルギーの準安定状態
から見た遷移状態の自由エネルギーが最も高いところ
が律速段階となる
逐次反応
同位体効果
B
R
P
kRI
kIP
R 
 I 
P
中間状態の自由エネルギーが始状態より低い場合
d[R]
d[I ]
d[P ]
 kRI [R]
 kRI [R]  kIP [I ]
 kIP [I ]
dt
dt
dt
kRI
中間状態は立ち上
e kRI t  e kIP t
[I]  [R]0
がってから減衰
kIP  kRI

[P ]  [R]0

kIP 1  e

 kRI t
  k 1  e 
RI
kIP  kRI
 kIP t
速度定数が大きく
異なると遅い方し
か見えなくなる
（律速段階）
反応座標が同位体原子の移
動を伴うとき、ゼロ点エネル
ギーが変わるため活性化エ
ネルギーが異なる
TS
R
反応座標が同位体原子
の移動を伴わなければで
なければあまり効かない
HD 同位体効果で
最大 10 倍ぐらい
F *
TS
F *
Product
Product
Reactant
qr
H、D

q
H、D
Reactant
qr
NTP Hydrolysis in Enzymes
プロトンのトンネル効果（量子効果）
活性化エネルギーより低いエ
ネルギー領域で、トンネル効
果により反応が起こる
ATP + H2O
反応振動数

*
~10 kcal/mol
TS
109-fold (F1-ATPase)
ADP + Pi
Diversity of proposed mechanisms
反応速度が増加
F *
反応振動数が大きいほどト
ンネル効果は大きい
D より H の方が反応速度
がより増加する
反応座標がプロトン移動でなけ
れば、同位体効果はない
Product
Reactant
qr
プロトン移動の反応障
壁が高いときに顕著に
現れる
Complex Reaction Nature
• Two elementary steps are involved:
1. Proton transfer (PT) for a water activation
2. P-O bond dissociation (POD)
OH- formation (PT)
-1
-1
High catalytic activity
of enzymes
-2
Ras + GAP (signaling, GTPase)
• Warshel et al. (1994) : a penta-covalent TS after PT
• Nemukhin et al. (2007) : a typical SN1 scheme
Myosin (motor, ATPase)
• Cui et al. (2004) : PT to Pi through a serine
• Nemukhin et al. (2007) : PT to a carboxylate
F1-ATPase (motor, ATPase)
• Hayashi et al. (2003) : PT to Pi through a water molecule
• Beke-Somfai et al. (2011) : Stepwise reactions
F1-ATPase: A Chemical-Mechanical
Energy Converter Hayashi et al., JACS (2012)
Noji et al.
(water activation)
ADP
Pi
• Highly polar electronic nature (-4 on TP)
• Difficulty in experiments due to complex
turnover cycles of functional processes
10 kcal/mol
ATP
Single molecule experiments
can detect clearly the
catalytic step out of its
complex turnover cycle.
Overall Reaction Profile
Hybrid QM/MM Simulation
Reactant
Determination of
the reaction path
in the protein
PT
17.0
(18.9)
POD
14.7
(15.0)
14.5
Equilibration and
modeling by
MD simulation
(Ito, Ikeguchi)
MD : MARBLE
(Ikeguchi et al.)
QM/MM : Hayashi and
Ohmine (JPCB 2000)
Product
SH, Shaikh, Umemura
QM/MM system
Exp. 13.5
PT
13.5
(17.0)
HBR
12.8
13.1
RDS
QM region
HBR
1.4
QM system
• 86 atoms
• B3LYP/6-31(+)g**
• # of basis function : 839
5.6
kcal/mol
HBR
4.6
-1.6
H3O+
0.9
0.0
2.3
(w/o ZPE)
• RDS is PT.
• POD occurs before PT.
Chain Activation Mechanism
Reaction Scheme
Efficient activation of the lytic water molecule (PT)
POD
PT
general
base
H3O+
RDS
activation
(dissociative PO3-)
water
relay
• Formation of a reactive (unstable) PO3- moiety
• Enhancement of PT by interaction with PO3-
Catalysis of POD
Changes of HB distances
Stabilization of
-phosphate
• Arg finger
• P-loop
Verification by Single Molecule
Experiments Ueno, Iino, Noji
Theoretically designed single molecule rate measurements
R373K
ATP (D2O, H2O)
ATPS (D2O, H2O)
H-D isotope effect
Analogue ATP
Mutant (KIE)
15
30
1,100 (1.0)
910 (0.88)
In other NTPases?
Arg finger
P-loop
Rate Limiting Step
POD or PT ?
Rate determining step
Prediction
Experiment
2.5
3.0
Sequential Order of POD and PT
H-D isotope effect
ATPS: an analogue ATP
ATP (D2O, H2O)
Proton
transfer
（PT）
H3O+
PT
ATPS (D2O, H2O)
Prediction
Experiment
The rate determining step is PT.
2.5
3.0
Order
（first）
P-O
dissociation
（POD）
PT
POD
Ratio to native
POD
PT
Prediction
Experiment
1/15
1/30
Correlation between POD and PT
Weak correlation
Strong correlation
Chain activation
mechanism
POD
PT
POD
Meta-Pi
A Mutant at Arginine Finger
R373K
PT
• Stable meta-phosphate
• High activation barrier for PT
• Unstable meta-phosphate
• POD and PT are close
both energetically and
geometrically.
R373K
native
KIE disappears!
POD
PT
Close proximity
between POD and PT
Conformational change of protein
Verification by Single Molecule
Experiments Ueno, Iino, Noji
Theoretically designed single molecule rate measurements
ATP (D2O, H2O)
Ratio to native (KIE)
Prediction
1/1,100 (1.0)
Experiment
1/910 (0.88)
ATPS (D2O, H2O)
R373K
織り込まれる様々な現象
‐F1 分子モーターを例に‐
分子モーター
化学反応
の因果律
サブユニットの
構造変化
H-D isotope effect
Prediction
Experiment
2.5
3.0
Analogue ATP
Mutant (KIE)
15
30
1,100 (1.0)
910 (0.88)
The agreements support our reaction scheme.
量子
力学
古典
力学
複合体の
構造変化
含まれる自由度の数
合成酵素の
因果律
織り込まれる様々な現象
‐F1 分子モーターを例に‐
分子モーター
化学反応
の因果律
量子
力学
Atomistic View of Slow Protein
Dynamics
Long time MD simulation
(Anton)
Shaw et al. Science (2010)
BPTI
古典
力学
サブユニットの
構造変化
揺らぎ
複合体の
構造変化
含まれる自由度の数
合成酵素の
因果律
Slow Dynamics of Protein is
Important for Enzymatic Processes
adenylate
kinase
~ms
w/o ligand
conformational
substate
(NMR, SM Spec.)
Conformationally
excited state
w/ ligand
BPTI
!
Intermittent transitions among five conformational substates
Slow ( > nano seconds) and non-linear
Slow Dynamics of Protein is
Important for Enzymatic Function
Dyhydrofolate Reductase (DHFR)
NMR
relaxation
-dissipation
Dynamic
energy
landscape
Wright and co-workers (2006)
Slow dynamics of protein
space
~ s - ms
< ps (TS relaxation)
global
local
Huge temporal and spatial mismatches
Protein relaxation influences reaction energetics.
1
  rlx
kad
<
time
Chemical reaction step
A Possible Scenario: Reorganization
Protein conformational change
Is Slow Dynamics of Protein
Important for Chemical Step?
1
  rlx
kdia
Chemical Step
If the reaction is slower than protein relaxation,
protein conf. change reduces the reaction barrier.
QM/MM-MD Simulations
Bacteriorhodopsin
Computational time of MD
Rhodopsin
100 days with 1 fs / step = 8.64 Ms
Hayashi et al.
BJ (2003)
Trajectory time
10 ps
Hayashi et al.
BJ (2009)
Prediction of time-resolved spectroscopic signals
QM/MM 100 ps
-MD
1 ns
10 ns
MM-MD
GPU
Pump-probe
Nature (2001)
Near IR pump-probe
Nature (2010)
Limited to short time trajectory calculations
Anton
Computational time/ step
864 s = 14.4 m
86 s
8.6 s
0.86 s
100 ns
12 step / s
1 s
120 step / s
1 ms
120k step / s
Ras-GAP QM
15 min / step
Computational cost
QM>>>MM
A long-time direct QM/MM-MD is difficult.
Protein conformational change
Difficulty in Determination of TS
QM/MMで自由エネルギー計算
熱揺らぎの環境の下での化学反応経路解析
MM-MD
（~ s）
MD や MC 計算により、熱揺らぎを受ける系の構造の
配置をサンプリングし、自由エネルギーを計算する
問題点 : QM/MM の高いコストにより、サンプリング数が取れない
QM/MM
（< ns）
Chemical Step
A short-time QM/MM-MD cannot reach the TS.
A MM-MD cannot describe a chemical reaction.
Basic Idea to Overcome the Difficulty
Complete separation between
• QM/MM calculation of a chemical reaction and
• MD simulation for protein conformational sampling
QM/MM free energy (FE) geometry optimization
Geometry optimization of the QM molecule on a free
energy surface defined with the MM conformational
thermal distribution obtained by MD simulations
解決法
1. static な経路を QM/MM で計算して、その経路に沿って MD で
自由エネルギーの補正をする（ただし、経路が必ずしも自由エネ
ルギーの極小経路であるとは保障されない）
2. 半経験的手法を QM に用いる（AM1、DFTB）
3. 自由エネルギー構造最適化法
QM/MM-FE Geometry Optimization
Nagaoka and co-workers
Free energy functional
F []  kBT ln   dRdX exp    E QM / MM {(R, X)}, R, X  
 : electronic wf, R : QM coordinates, X : MM coordinates
Give up sampling of R
F [, R]  kBT ln  dX exp    E QM / MM {(R, X)}, R, X  
Free energy surface of QM coordinates
Gradient on FE surface: mean gradient
F
E QM / MM

Ri
Ri
Problem: Electronic WF depends on X.
~ 1  second (> 1,000 times longer than direct QM/MM-MD)
Its computational cost is not very different from a
direct QM/MM-MD.
X
Mean Field (MF) QM/MM Method
Yamamoto
F [, R]  kBT ln  dX exp    E QM / MM {(R, X)}, R, X  
Kosugi and SH, JCTC (2012)
Kosugi and SH, JACS (2012)
Problem of QM/MM-FE
Conventional
Present method
MM
MF approximation
QM
MF-QM/MM
(R, X)   MF (R )
QM(d;R )
F [ MF , R]  kBT ln  dX exp    E QM / MM { MF (R )}, R, X  
Variational
solution
QM/MM RWFE-SCF Method
ˆt R  V Cl  d, R 
f QM / MM  d, R   f 0  d, R   q
MF of MM region is computed by MD.
MM(X)
QM(d;R )
MM(X)
Frequent iteration of
QM/MM opt. and MD
sampling is necessary.
Convergence problem!
MM (d, R; X) 
QM
MM
Yang et al.
reweighting
MM distribution is reweighted.
Much quicker convergence of
MM thermal distribution




exp    E QM  MM (d, R; X)  E QM  MM (d0 , R 0 ; X) 

   (d , R ; X)
0
0
MM
exp    E QM  MM (d, R; X)  E QM  MM (d0 , R 0; X) 

 0
Well-defined variational method
QM/MM Long-Range Interaction
Procedure
MD QM/MM MD QM/MM
i
R opt
 R 0i 1
…
MD QM/MM MD QM/MM
Distribution
of energy
difference
i
qopt
 q0i 1
Update
3 – 15 ns
~20,000 confs.




exp    E QM  MM (d, R; X)  E QM  MM (d0 , R 0 ; X) 

   (d , R ; X)
0
0
MM
exp    E QM  MM (d, R; X)  E QM  MM (d0 , R 0; X) 

 0
Conventional: sphere cluster
for continuum electron density
• Present method: RESP operator
treatment (Hayashi and Ohmine, 2000)
makes its implementation possible.
The method allows one to use
directly MD samples obtained
by sophisticated MD program
packages using PME method.
MD
MM (d, R; X) 
Ewald summation technique
• Conventional: hard to implement
Complete separation of
QM/MM and MD sampling
Present treatment:
periodic boundary condition
Test System: Amylase
33 ns
33 ns
Spectral Tuning of C1C2
Kamiya, Kato, Ishitani, Nureki,
SH, CPL (2012)
21 ns
静電相互作用マップ
•
•
正電荷を青（赤）い場所に置く
と、青方（赤方）シフト
負電荷ならばその逆
4 ns
240 ns
Glu162
Kosugi and Hayashi, JACS (2012)
Kosugi and Hayashi, JCTC (2011)
Summary
• Characteristic slow dynamics of proteins can
contribute to enzymatic catalysis.
• QM/MM simulations now almost catch up with
MD simulations with MM force fields in terms
of time scale the simulations can handle.
• The QM/MM method developed can easily
exploit advance in MD simulations with MM
force fields.
• The QM/MM method can be a powerful tool to
rationally develop protein tools with new
functionalities.
ChR2-ET
Gunaydin et al., Nat. Neuro. (2010)