2017年11月22日

遠隔インタラクティブ講義「計算生命科学の基礎（IV）」

QM/MM法を用いたタンパク質の機能解析
広島市立大学大学院情報科学研究科
鷹野 優
ytakano@hiroshima-cu.ac.jp

概要
I.

序論

II. QM/MMの理論
III. 金属タンパク質中の補欠分子の電子構造
III-1.

酸素運搬タンパク質ヘムエリスリン

III-2.

光合成反応中心のスペシャルペア

IV. 生体高分子生化学的機能解析のための分子計算技術の開発
IV-1.

半経験的分子軌道法によるQM部分の高速化

IV-2.

超並列QM/MM-MD連成計算プログラムPlatypus

V. まとめI. 序論タンパク質は生命活動を担うタンパク質 ＝ 生命活動をおこなう分子

免疫

物質輸送

構造維持

触媒

生きるのに必要な様々なはたらきをする
不具合が起きると病気に！
Molecule of the Month © David S. Goodsell and RCSB PDB licensed under CC-BY-3.0 licenseタンパク質と病気
骨形成不全症通常骨折にならないようなごく軽微の外的要因でも骨折する。
ヒトの組織の構造タンパク質、I型コラーゲンの変異が原因

骨形成不全症

正常

Gly

Gly→Alaタンパク質の一部が変化

Ala

構造を維持できなくなる→骨形成不全症

タンパク質と薬
抗インフルエンザ薬：タミフル

Structure Based Drug Design
Molecule of the Month © David S. Goodsell and RCSB PDB licensed under CC-BY-3.0 license
https://sites.google.com/a/comp-prot-sci.org/tanpakushitsu-no-rittai-kouzou-nyuumon/chapter-6

Our goal
構造機能相関
分子構造

分子のかたち電荷分布
軌道相互作用機能（例） Croタンパク質
DNAと強く結合し、RNAポリメラーゼによる転写を邪魔します
DNA

赤：負電荷
青：正電荷

Croタンパク質

タンパク質がDNAを認識する

Our goal
分子構造ー電子構造ー機能相関
分子構造

分子のかたち

電子構造

機能
電荷分布
軌道相互作用分子構造・電子構造は機能発現をどのように制御しているか？

タンパク質の電子構造
1) タンパク質は巨大で、ヘテロで、ダイナミックな系
2) タンパク質の機能の発揮は局所部位の電子状態に依存するタンパク質のつくる電場

光合成反応中心タンパク質のつくる環境は電子構造に影響

タンパク質の電子構造
1) タンパク質は巨大で、ヘテロで、ダイナミックな系
2) タンパク質の機能の発揮は局所部位の電子状態に依存する

はたらく時には形が変わる。
タンパク質のうごき(ダイナミクス)は電子構造に影響

タンパク質の電子構造
1) タンパク質は巨大で、ヘテロで、ダイナミックな系
2) タンパク質の機能の発揮は局所部位の電子状態に依存する
計算サイズ
活性中心そのもの

小活性中心をとりかこむタンパク質
溶質分子をとりかこむ溶媒分子
アンサンブルダイナミクス

大資源（計算機・人間の能力）には限りがある「どのようにしてモデル化するか?」

タンパク質の電子構造モデル化の功罪
功：計算量の減少
系の単純化
罪：情報の欠損
資源の分配をどうするか
なにを理解したいのか？
賢いモデルの選択は本質の理解に必須タンパク質の電子構造
1) タンパク質は巨大で、ヘテロで、ダイナミックな系
2) タンパク質の機能の発揮は局所部位の電子状態に依存する

QM/MM法の適用Electrosta+c
eﬀect
QM
Stericeﬀect

MMPhotoactive Yellow Protein

II. QM/MMの理論

QM/MM法・多階層連結計算法の一つ
・電子構造が重要な部分（化学反応）は量子力学(QM)を
それ以外は分子力場(MM)を用いる
・1976年にWarshelとLevittが提案
・2013年ノーベル賞（Warshel, Levitt, Karplus)

h0p://www.nobelprize.org/nobel_prizes/chemistry/laureates/2013/popular-chemistryprize2013.pdf

QM/MM法
・全系(System)をQM領域(Inner)とMM領域(Outer)に分割
・QMとMMの境界に共有結合が含まれる場合は特別な取り扱いが必要

Inner
(QM)
Outer
(MM)

SystemSystem = Inner + Outer + (Link)
共有結合

QM/MM法に使う計算手法
・基本的に制限はない
QM: 密度汎関数法(DFT)
Hartree-Fock法(HF)
post-HF法(MP2, CASSCFなど）
半経験的分子軌道法（AM1, PM6法など）
MM: AMBER
CHARMM
GROMOS
OPLS-AA

QM/MM法の種類
全系のエネルギーの表現法に関して大きく二種類ある減法表現（ONIOM法）
sub
EQM/MM(System ) =EMM (System ) + EQM ( Inner+Link ) − EMM ( Inner+Link )

+

=

‒

QM-MM
相互作用

加法表現（QM/MM法）
add
EQM/MM(System ) =EMM (Outer ) + EQM ( Inner+Link ) + EQM/MM ( Inner,Outer )

=

+

+

QM-MM
相互作用

タンパク場の効果
mechanical embedding
QMの電子密度はMMと相互作用しない
Steric eﬀectのみ考慮

electrostatic embedding
QMのハミルトニアンにMMとの静電相互作用をとりこむ
Electrostatic eﬀectとSteric eﬀectを考慮
MMの電荷がQM領域に近いときには過分極を起こす→・電荷スケーリング法
QM境界付近のMM電荷の大きさを調整
・電荷シフト法
QM境界付近のMM電荷をより外の原子にシフトする
・ガウス型電荷分布

QM/MM境界問題
QM/MM境界に化学結合が存在する場合
QM空間の切断->不対電子対の発生->QM（分子軌道）の歪み
正しいエネルギー及び力が得られない

1. どのように切るか？
2. どこを切るか？

QM/MM境界問題
1. どのように切るか？
(1)Link atom method
QM領域とMM領域の境界にある
化学結合の末端を水素原子などで置き換える

Hydrogen Atom etc…
MM
不対電子対

(2)局在化軌道法
QM領域とMM領域の境界にある
化学結合を局在化軌道(SLBOs)を用いて表現する。
Lin and Truhlar, Theor. Chem. Acc. 117, 185–199 (2007)

QM

QM/MM境界問題
2. どこを切るか？
(1) 経験的、化学的洞察
活性中心から遠方の領域を選択
sp3炭素で切る(2) 線形応答関数解析
QM/MMモデル化に起因する
QM領域の電子構造の誤差の最小化δρ (r)δρ (r) = ∫δ v(r!)dr!
δ v(r!)

線形応答関数解析
従来：QM/MMモデルによる誤差 δ v(r!) の最小化
QM/MMモデル由来のQM領域の誤差δρ (r)の最小化[ガイドライン](i) 共有結合では電子トランスファーが大きい平衡核間距離ではなく、結合不安定が起こる(結合が切れかけた)状態で最大となる。
(ii) トランスファーブロックがない状況ではサイト数が多くなる
につれフリーデル振動型の振動・減衰となる。
(iii) アラニンジペプチド, グルタチオン, アラニントリペプチド系
ではsp3接合点が伝播をブロックする。水素結合経由の伝播は約0.001(ex. 0.1 Hartreeの δ v(r!) に対し0.0001の誤差伝播)。
αヘリックス系では明白に誤差伝播が構造化される(右図)。

従来のsp3接合点での切断が有効
QM/MM境界問題の定量化が可能に！

δρiatom δ v jatom =1
δρiatom δ v jatom = −1δρiatom δ v jatom = 0.0
RGB=(1-LRF, LRF, 0.5)/max LRF

Ueda et al. Int. J. Quantum Chem. 113, 336–341 (2013).

QM/MM法による構造最適化
計算コスト
QM > QM/MM >> MMmicroiterative scheme
QM部分とMM部分の構造最適化に異なるアルゴリズムや座標系を使うQM: macroiteration（内部座標＋擬ニュートンアルゴリズム）
MM: microiteration（カーテシアン座標＋共役勾配法）

構造最適化の手続き(1)断熱アプローチ
QMのmacroステップではMM環境は完全に
構造最適化される(2)交互スキーム
QMとMMの構造最適化は交互に行われる

溶媒分子のとりあつかい
シミュレーションの計算コストの要因①

タンパク質 1,870原子(112アミノ酸残基）
水 24,039原子（8,013分子）系のほとんどが水分子！
②

水分子の配置＋配向（高い自由度）構造の異方性
水素結合

多数の溶媒分子の配置・配向をどのように表現するか？

溶媒分子のとりあつかい
シミュレーションの計算コストの要因
多数の溶媒分子の配置・配向をどのように表現するか？
①
1. あらわにとりあつかう
・古典力学
MM
2. モデル化してとりあつかう②
構造の異方性
水素結合

・統計モデル
RISM-SCF, ER
・連続誘電体モデル
PCM, CPCM, SMx, COSMO

古典力学 (QM/MM, ONIOM)
H

H
H

O

H

O

H

H

O

O

O
H

O

H

CH3

C

H
H

O

O
H

H
O

CH3
O

H

H

H
O

O

H

H

O

O
O

H

H

O

H

H
H

H

CH3

QMCH3

C

H
H

H

利点：精度

H

O

O

H

O

H

H
O

H
H

O

H
O

MM

H

活性中心：量子力学、タンパク質：古典力学
溶媒環境：古典力学

欠点：計算量（時間（アンサンブル））

H

H

統計モデル (RISM-SCF, ER)
H

O

H

H

O

O
O

H

O

H

CH3

C

H
H

QM

H

H

O

O
H

H
O

CH3
O

H

MM
H
O

H

H
H

O

H

溶媒の分布関数H

利点：精度・計算量

活性中心：量子力学、タンパク質：古典力学
溶媒環境：統計理論

欠点：ダイナミクスのとりこみ連続誘電体モデル(PCM, CPCM, SMx, COSMO)
H
H

O

H

H

O

O
O

H

O

H

CH3

C

H
H

water
ε = 78.39

H

O

O
H

H
O

CH3
O

H

H
O

H

H
H

O

CH3

O Cavity
QM CH3
C
O
εMM
=1

H

H

利点：計算量活性中心：量子力学、タンパク質：古典力学
溶媒環境：連続誘電体

欠点：精度（水素結合のとりこみ）

PCMの形式
溶質分子が溶媒分子に囲まれたとき、ハミルトニアンH = H 0 +Vint

真空中の分子（溶質）のハミルトニアン

溶質−溶媒相互作用のポテンシャル

PCMでは、
1. 溶媒分子を誘電率εの連続媒質におきかえ、
ε = 78（水）PCMの形式
溶質分子が溶媒分子に囲まれたとき、ハミルトニアンH = H 0 +Vint

真空中の分子（溶質）のハミルトニアン

溶質−溶媒相互作用のポテンシャル

PCMでは、
2. 溶質の形にもとづいてcavity（空孔）をつくり、
ε = 78（水）ε0 = 1.0

PCMの形式
溶質分子が溶媒分子に囲まれたとき、ハミルトニアンH = H 0 +Vint
溶質−溶媒相互作用のポテンシャル

真空中の分子（溶質）のハミルトニアン

PCMでは、
3. cavity（空孔）の中に溶質分子をおく
ε = 78（水）OH

ε0 = 1.0

HPCMの形式
ε = 78（水）

溶質分子が溶媒分子に囲まれたとき、O

ハミルトニアンH = H 0 +Vint

真空中の分子（溶質）のハミルトニアンHH

ε0 = 1.0

溶質−溶媒相互作用のポテンシャル

溶媒効果を溶質の電荷密度と溶質によって分極した連続誘電体のつくる電場との相互作用として表す(Reaction ﬁeld法(RF法))

Vint = ∫ ρsolute ( r) Φ in,solvent ( r) dr
溶質の電荷密度

溶質によって分極した溶媒（連続誘電体）のつくるcavity内の電場

PCMの形式
誘電体のつくる電場は、cavity表面近傍に分布する分極電荷の与える電場に等しい
計算では、cavity表面を面積素片(tessera)にわけて、
面積素片kの nk
法線ベクトルΦ in ( rk )Φout ( rk )面積素片k（位置rk,電荷Qk, 面積Sk)

O
H

H

Qk 面積素片のもつ電荷Φ in,solvent ( r) = ∑面積素片の位置ベクトルk r − rk
面積素片のもつ電荷は古典電磁気学から求められる
比誘電率

面積素片の面積
面積素片の法線ベクトルε r −1Qk =
Sk ( ∇Φ in ( r) • n k )4πε rcavity内の電場

PCMの形式
cavity内の電場は溶質によるものと溶媒によるものにわけることができるΦ in ( r) = Φ in,solute ( r) + Φ in,solvent ( r)
溶質の電荷密度によりつくられるcavity内の電場

溶質によって分極した溶媒（連続誘電体）のつくるcavity内の電場

Φ in,solvent ( r) の計算には Qk が必要

Qk の計算には Φ in ( r) が必要Φ in ( r) の計算には Φ in,solvent ( r) が必要

くりかえし計算で解く
self-consistent reaction ﬁeld法（SCRF法）III. 金属タンパク質中の補欠分子

生物と金属バランスを保つことで生物の生命維持に必須のはたらきをする

多量
Na
K
Mg
Ca

微量
Fe
Zn
Cu
Mn
V

超微量
Ni
Cr
Co
Mo
W

呼吸・光合成・窒素固定
どのように金属と生体分子が相互作用（制御）しているか？

金属タンパク質「金属活性部位」とタンパク質が役割分担

金属をもつ生体分子の特徴（環境にやさしい触媒）
例：窒素（反応しにくい安定な気体）からアンモニアを生成

N2 + 3H2 → 2NH3・ハーバーボッシュ法

Fe3O4

N2 + 3H2

高温・高圧

2NH3

(500 ℃, 300気圧)・ニトロゲナーゼ（鉄：Fe, モリブデン：Mo)

N2 + 3H2

ニトロゲナーゼ

常温・常圧

2NH3

金属タンパク質の機能発現機構？

金属補欠分子タンパク質の生物活性において重要な、タンパク質に結合する非タンパク質

活性中心として機能を発現する

ヘム (Fe)スペシャルペア (Mg)

酸素発生複合体(Mn)[4Fe-4S]クラスター (Fe)Inside-out アプローチ
1. 金属補欠分子固有の電子構造を調べる
2. タンパク質中での電子構造を調べる
3. 1と2を比較する
static eﬀects (構造歪み、静電相互作用)
dynamic eﬀects (ゆらぎ)（例）シトクロムc酸化酵素のヘム a in CcO
タンパク質なし
タンパク質ありTakano and Nakamura,
J. Comput. Chem. (2010)

タンパク質は電荷移動を強める
PropAの電荷変化がプロトン移動のトリガー

遷移金属と生体分子の相互作用
酸化数の変化
電荷

強相関電子系
スピン

軌道

スピン状態の変化

stepwiseな電子移動

軌道・電荷・スピンが協同的に影響しあう
理論的取り扱いが非常に困難

強相関電子系
金属イオンの電子間にはたらく強いクーロン反発相互作用
さまざまな電子配置・電子相関を考慮にいれないといけない（例）酸素発生複合体 (Mn4CaO5クラスター)

可能な電子配置数
Mnの3d軌道
Oの2p軌道約100京通り！

ab initio DMRG法で可能になった
Kurashige, Y.; Chan, G. K.-L.; Yanai, T. Nature Chem. 2013, 447, 660–666.

Density Functional Theory (DFT)E = ET + ENE+ ENN + EEEApproximation

In DFTEDFT = ETs + ENE+ ENN + EJ + EXC
厳密な交換相関汎関数はわからない
ETs : kinetic energy of non-interacting system
EXC : exchange-correlation energy

Exchange-Correlation Energy
The exchange-correlation energy has the form:
Correction of electron-electron interaction

EXC = (ET – ETs) + (EEE – EJ)
Correction of kinetic energy

DFTの精度は交換相関汎関数項に依存する

DFTの利点
Excellent performance-to-cost ratio !
performance
DFT
HF
MP2
CC
cost
DFT
HF
MP2
CC

ex. Mean absolute error of heat of formation
for the G2 dataset (kcal/mol)
3.0 (B3LYP), 7.1(BLYP), 90.8 (SVWN)
149.2
23.8
11.5
ex. Benzene (C6H6) single point
O(N3)~O(N4)1 min 38 sec (SVWN), 2 min 20 sec (BLYP)O(N4)1 min 48 sec
O(N5)2 min 36 sec
O(N7)2 hours 29 min 20 sec
Xu, X. and Goddard III, W. A. , Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 2673-2677.
Jensen, F. Introduction to Computational Chemistry, 2nd ed. Wiley (2007)

DFTの欠点
Obscure systematic improvability
There is not clear path to the exact solution of the Schrödinger eq.,since DFT show very good accuracy for a wide range of applications
due to the error cancellation.

DFTの欠点
Obscure systematic improvability
There is not clear path to the exact solution of the Schrödinger eq.,since DFT show very good accuracy for a wide range of applications
due to the error cancellation.

DFTの交換相関汎関数のチェックが必要！
（特に遷移金属系）

III-1.酸素運搬タンパク質ヘムエリスリン

Takano, Y.; Yamaguchi, K.; Nakamura, H. Int. J. Quantum Chem. 2013, 113, 497–503

ヘムエリスリン
Hr is an oxygen transport protein found in several marine
invertebrates.
OxyHr

DeoxyHr[Fe(µ-OH)Fe]
His73
His101

His77 HOFeIIO2[Fe(µ-O)Fe]
His77

His54
FeII

O
O
O
O
Glu58
Asp106

Jab = –13 cm–1

His73

His25

His101hydroperoxide
O
H
O
His54

O
FeIII

FeIII

O
O
O
O
Glu58
Asp106

His25

Jab = –77 cm–1

In the reac+on with molecular oxygen, two electrons from the FeII–FeII center and a
proton from the bridging µ-OH group are transferred to O2 binding to the ﬁvecoordinated Fe center, yielding the µ-O bridged diferric site possessing a terminalhydroperoxide (oxyHr).

二核鉄活性中心
Diiron ac+ve sites with µ-O and/or µ-OH bridges are found in many
proteins.
Hemerythrin (Hr)
His
His
His

H

O

O

O
Fe

Fe
O
O

O

O

Glu Asp
Oxygen transport

Metane monooxygenase
Ribonucleo?de reductase
(MMO)
(RNR)
Glu
H2O
H2O
Glu
H2O
H
His
His
O
His
O
Asp
His
O Fe
Fe
Fe
Fe
O
Glu
Glu
O
His Glu
H
O
O
O
O
His
Glu
Oxygen ac?va?on(Conversion of methane
to methanol)

Glu
Oxygen ac?va?on(Conversion of ribonucleo+de
to deoxyribonucleo+de)

The structure of the ac?ve site are similar to each other,
but the func?ons are diﬀerent.

以前の研究
X-ray
Fe2

E
α

β

–0.8

–8.0

Opt
Fe2

O2
α

β

α

β

π∗, π⊥∗
dxy (δ⊥)dxy (δ⊥)

–9.0
–10.0
π∗, π⊥∗–11.0

dxy (δ⊥)–12.0

dx2–y 2 (δ)
dyz (π)dxy (δ⊥)dzx (π ⊥)dz2 (σ)dyz (π)dzx (π ⊥)dz2 (σ)dx2–y2 (δ)

–13.0
–14.0
–15.0
–16.0
タンパク質環境は活性中心の構造を歪めることで、金属のd軌道のエネルギー準位を制御し酸素結合を可能にしている

Takano, Y.; Isobe, H.; Yamaguchi, K. Bull. Chem. Soc. Jpn. 2008, 81, 91–102

目的
ヘムエリスリンのタンパク質環境による静電相互作用は酸素運搬機能にどのようにはたらくのか？

アイデア
タンパク質環境＝タンパク質環境を構成する原子電荷との直接の静電相互作用

MM＋

活性中心の電子状態によって誘起される原子電荷の揺らぎ（分極効果）

PCM

手法の選択How do we determine an appropriate broken-symmetry
theoretical treatment?"
1. Use a method which many people use.
→ UB3LYPUB3LYP provides very accurate results for organic systems.
But how about transi?on metal systems?

2. Use a method which can reproduce experimental data.
Magne?c coupling constants (Jab) were reported.
(OxyHr: –77 cm–1, DeoxyHr: –13 cm–1)

How to calculate magnetic coupling constants? "化学結合と磁気的相互作用◆Example (H2 molecule)Magnetic Coupling Constants, JabHeisenberg Hamiltonian"O

H = −2∑ Jab Sa ⋅ S b

Fe

Fe
O
O

O

O

◆Jab values!
LS

J ab = HSSaSb

E− HS E
2

S −

LS

S2

Jab > 0: ferromagnetic"
Jab < 0: antiferromagnetic"€

Each iron site is regarded as a spin site."
Magnetic coupling constants, Jab, can be calculated using the Heisenberg hamiltonian.

"

計算方法
◆Electronic Energy"
E = ET + EV + EJ + EXC!
ET ! Kinetic energy term"
EV! Nuclear-electron attraction"
EJ! Coulomb integral"
EXC! Exchange-correlation functional"

EXC = C1EHF + C2ESlater + C3EBecke + C4EVWN + C5ELYP !
UHF
UB2LYPUB3LYPUB3LYP*
UBLYP

C1
1.0
0.5
0.2
0.15
0

C2
0
0.5
0.8
0.85
1.0

C3
0
0.5
0.72
0.72
1.0

C4
0
1.0
1.0
1.0
1.0

C5
0
1.0
0.81
0.81
1.0

モデル
OxyHr (PDB ID: 1HMO)"
His77
His73
His101

H

O

O
His54

O
FeIII

FeIII

O
O
O
O
Glu58
Asp106

His25

DeoxyHr (PDB ID: 1HMD)"
His73
His101

His77 H
O
FeII

His54
FeII

O
O
O
O
Glu58
Asp106

His25

His"
Asp, Glu"

methyl imidazole"acetate"Magnetic Coupling Constants, Jab

DeoxyHr (Exptl.: –14 cm–1)
SolvationUHF

UBHandHLYPUB3LYP

UBLYP

真空

–3.49

–10.4

–145.1–396.1

PCM

–3.59

–10.6

–128.2–326.4

PCM + MM

–3.58

–10.6

–137.9–398.8

UBHandHLYPは妥当的な磁気的相互作用をあたえるMagnetic Coupling Constants, Jab

OxyHr (Exptl.: –77 cm–1)
SolvationUHF

UBHandHLYPUB3LYP

UBLYP

真空

–14.1

–49.0

–145.1–396.1

PCM

–17.2

–64.0

–128.2–326.4

PCM + MM

–16.7

–64.4

–137.9–398.8

PCM + MM、PCMは妥当的な磁気的相互作用をあたえるOrbital Energy Gap

タンパク質環境による静電相互作用は活性中心の、金属のd軌道のエネルギー準位を制御し酸素結合に寄与しているOrbital Energy Gap
E (eV)

O2 π*軌道Fe dδ軌道

‒0.8
‒6.71
‒8.27
‒8.84
Opt

X-ray

タンパク質による構造歪み

×○

X-ray+PCM+MM○

タンパク質による静電相互作用

×

×

○

O2 π*軌道との軌道相互作用

弱

強

まとめ
1. BHandHLYP法はOxyHrとDeoxyHrの双方で妥当的な磁気的相互作用を与えた
2. ヘムエリスリンの蛋白質環境による静電相互作用は活性中心の金属のd軌道のエネルギー準位を制御し酸素結合に寄与しているIII-2. 光合成反応中心スペシャルペア

Yamasaki, H.; Nakamura, H.; Takano, Y. Chem. Phys. Lett. 2007, 447, 324–329
Yamasaki, H.; Takano, Y.; Nakamura, H. J. Phys. Chem. B 2008, 112, 13923–13933

光合成
植物などの光合成色素をもつ生物がおこなう、
光エネルギーを化学エネルギーに変換する生化学反応

二酸化
炭素 + 水
(H2O)
(CO2)

hν

酸素
糖
(O2) + (C6H12O) +光合成電子伝達反応（明反応）+

炭素固定反応（暗反応）

ATP

光合成電子伝達反応（明反応）

数々のタンパク質複合体により構成されている
Essen+al細胞生物学

光化学系

反応中心: 光エネルギーを電子エネルギーに変換
Essen+al細胞生物学

光合成反応中心
光エネルギーを電子エネルギーに変換(量子収率 > 0.8)hνSpecial pair (SP)
e–

BchlM

BchlL

PheoM

PheoL
QM

Fe

chromophores

QL

光合成反応中心
光エネルギーを電子エネルギーに変換hνHis(L173)

Special pair (SP)
BchlM

SPL

BchlL

PheoM

PheoL
QM

M-branchFe

QL

L-branch

SPM

His(M202)

色素が疑似二回対称性を持つにもかかわらず
電子移動経路は一方向スペシャルペアカチオンラジカルのスピン密度
wild type
Leu(M160)His
Leu(M160)His + Phe(M197)His

スペシャルペア

SPM

SPL

Arg(L135)Glu
Arg(L135)Leu
Leu(M160)His + His(L168)PheLeu(M160)His + Leu(L131)His

SPL

Phe(M197)His
Arg(M164)Leu

SPM

Arg(M164)Glu
Leu(L131)His
His(L168)Phe
His(L168)Phe + Phe(M197)HisLeu(L131)His + Phe(M197)HisLeu(L131)His + His(L168)Phe
100% 90% 80% 70% 60% 50% 40% 30% 20% 10% 0%スペシャルペアカチオンラジカルのスピン密度は非対称
Lubitz et al., Acc. Chem. Res. 2002, 35, 313

目的スペシャルカチオンラジカルのスピン密度の非対称性の原因は何か？
光合成反応中心の非対称な電子移動タンパク場の効果
QM/MM

l スピン密度の偏りが再現された
l スペシャルペア単体でもスピン密度に偏りが生じている
l タンパク質の静電相互作用がスピン密度の非対称性が強める
Yamasaki, H.; Nakamura, H.; Takano, Y. Chem. Phys. Lett. 2007, 447, 324–329[QM]方法:UB3LYP 基底関数6-31G(d)(0.25) [MM] AMBER96

側鎖の影響Acetyl group

Nonpolar groups

13keto carbonyl
group
Methyl ester(Mes) group
Phytyl-I (Phyt) group
Phytyl-II (Phyt) group

Mes基とPhyt基の配向の違いによりスペシャルペアの電子構造の非対称性が生じた

側鎖の影響
Yellow: positive Δρ
Blue: negative Δρ

InvertedMes基とPhyt基の配向が異なる

Phyt

Mes
S244(L)

Cys247(L)Gly280(M)Gly283(M)配向の違いは保存されている

1aig

Rhodobacter
sphaeroides

2i5n

1eys

Rhodopseudomonas
viridisPurple bacteria

2axt

Themochromatium tepidum

Thermosynechococcus elongatus
Cyanobacteria

タンパク質環境が配向の違いを与える
Side view

Gly280Ile284

Met248

Ser244

Ser206

Rb. sphaeroides (red)Phe180 Rh. viridis (blue)
T. tepidum (green)
T. elongatus (yellow)

Ala176

Leu209Thr277
Leu156Top view

Leu131

Tyr210

Phe181M-sideL-side４つの結晶構造の重ね合わせでは、
Mes基とPhyt基のまわりはよく重なっている

配向の違いを与える
アミノ酸残基は保存されている
Leu/Val

Phe/Val Phe/Phe

Ala/Gly

Ile/IleSer/Thr配向の違いを与える
アミノ酸残基は保存されている
Leu or Phe/Val

Leu/Val Tyr/Leu?/Leu Ile/Ile

Ser/Gly
Gly or Ala/Leu
Mes基とPhyt基の配向を決めるアミノ酸残基はよく保存されている

NMRによる検証

Gupta, K. B. S. S. et al. J. Am. Chem. Soc. 2013, 135, 10382–10387

まとめ
1. Mes基とPhyt基の配向の違いがスペシャルペアの電子構造の非対称性を引き起こした
2. タンパク質の静電相互作用がスピン密度の非対称性を強めた
3. Mes基とPhyt基の配向を決める配列は保存されていた

IV. 生体高分子生化学的機能解析のための分子計算技術の開発

Inside-out アプローチ
1. 金属補欠分子固有の電子構造を調べる (QM)
2. タンパク質中での電子構造を調べる (QM/MM)
3. 1と2を比較する
static eﬀects (構造歪み、静電相互作用)
dynamic eﬀects (ゆらぎ)

高精度化した半経験的分子軌道法によるQM部分の高速化+
超並列ab initio QM/MM-MD連成計算プログラム
Platypusの開発

IV-1. 高精度化した
半経験的分子軌道法によるQM部分の高速化

Saito, T.; Kitagawa, Y.; Takano, Y. J. Phys. Chem. A 2016, 120, 8750–8760.
Saito, T.; Takano, Y. Chem. Lett. 2017, 46, 1567–1569.

QM/MM-MD
シミュレーションの計算コスト
time for QM/MM
single point calculation / s

600

UDFT/DZP levelsize of QM region

number of steps in MD rune.g.,) 1ns = 1×106 step

Provided that every QM/MM energy and force
computation takes 10 min.,it adds up to the 1,000,000 QM/MM computations→ one MD run takes 20 years !!

Semiempirical MolecularOrbital (SE-MO) modeltime for QM/MM

single point calculation / s

600

UDFT/DZP levelroughly 2 or 3 orders ofmagnitude faster0.1

Semiempirical MO (SE-MO)

size of QM region

number of steps in MD rune.g.,) 1ns = 1×106 step

1ns MD run requires 70 days.
M. J. S. Dewar and W. Thiel, J. Am. Chem. Soc. 99, 2338 (1977).
W. Thiel, WIREs. Comput. Mol. Sci. 4, 145 (2014).

SE-MO法の精度
GMTKN30-CHNOF database(430 reactions mostly involving closed-shell molecules)
L. Goerigk and S. Grimme, J Chem. Theory Comput. 6, 107 (2010).
L. Goerigk and S. Grimme, Phys. Chem. Chem. Phys. 13, 6670 (2011).

MAEs (kcal/mol)*

遷移金属への対応

MNDO

27.3

No

AM1

14.7

No

PM3

11.3

No

PM6

10.2

YES

OM3

6.3

No

＊mean absolute errors, excluding the MB08-165 subset

P. Dral et al., J. Chem. Theory Comput. 12, 1097 (2016).

Given the applicability to a wide range of molecules,
we chose PM6, which supports 70 elements in the periodic table.
J. J. P. Stewart, J. Mol. Model. 13, 1173 (2007).

Some problems using SE-MO?

■methylene (1,3 CH2) radicalC
H

H

bent : H-C-H = 102.3˚■ cyclobutadiene (1C4H4) radical
H

H

H

linear with the H-C-H angleof 180˚

C

C

C

C

H

non-planarsquare planar (D4h)■ transition state of C=C epoxidation with a Mn complex

×

fail to find a TS

The spin-unrestricted PM6 (UPM6) computation does not work.

目的
Our primary purpose is to improve the accuracy of PM6 in describing
electronic and geometric structures of open-shell species,
for which a slight modification was made.
reparameterized PM6(rPM6)
the original PM6

training set containing (di)radicals using
the UB3LYP/6-31G* model as reference data••••
part of GMTKN30 was also included to maintainthe performance for general-purpose applications
The number of optimizable electronic and core−core parameters adds up to
43 (P) for the basic elements (H, C, N, O).
T. Saito, Y. Kitagawa, Y. Takano, J. Phys. Chem. A, 120, 8750 (2016).

Application to Methylene(included in the training set)
Method1CH

2

C-H H-C-H
UAM1
1.068 140.6
UPM3
1.067 131.3
UPM6
1.018 180.0
UrPM6
1.081 111.9
UB3LYP/6-31G* 1.090 113.5
Expt.
1.107 102.3
a With the AP correction

3CH

C-H
1.062
1.063
1.019
1.039
1.082
1.077

2

H-C-H
151.7
149.9
180.0
152.6
133.1
134.0

S-T gap[kcal/mol]a
34.1
42.7
31.5
11.2
6.0
9.0

It is not surprising that the present UrPM6 shows great improvement
over the original UPM6.
The standard SE-UMO calculations are likely to provide strong diradical
characters that entail the significantly large S-T gaps.

rPM6 shows great improvement
over the original UPM6

■methylene (1,3 CH2) radical
C
H

bent
H-C-H = 111.9˚H■ cyclobutadiene (1C4H4) radical
H

H

H

C

C

C

C
H

square planar (D4h)planar■ transition state of C=C epoxidation with a Mn complex
CPU time ratio (vs. DFT) "force calc.
1 : 4608"
Hessian calc. 1 : 3559"
Saito, Takano. Chem. Lett. 2017, 46, 1567–1569.

IV-2. 超並列ab initio QM/MM-MD
連成計算プログラム
Platypus

Takano, Y.; Nakata, K.; Yonezawa, Y.; Nakamura, H. J. Comput. Chem. 2016, 37, 1125–1132.

超並列ab initio QM/MM-MD連成計算プログラム
Platypusの開発(PLATform for dYnamics Protein Uniﬁed Simulation)

E = EQM + EMM + EQM/MM
Platypus-QM

Platypus-MM

HF, DFT, CASSCF, CIS, 並列性能に優れた新規
MP2プログラムの長距離ポテンシャルの開発
高並列化による高速化 (zero-multipole法)+
Chain-of-State法 (疎結合並列計算）

現在までの研究開発成果•

RDF, UHF, R-DFT, U-DFT, CIS, CIS(D), CIS-DFT, MP2, CASSCF の実装•

エネルギー・力計算における積分計算のハイブリッド並列化とSIMD化•

2電子積分計算の原子軌道基底から分子軌道基底への変換部分の高度化•

CASSCF計算におけるdirect CI 計算の高速化（CAS(16,16)まで可能）

QM部分の高速化・超並列化
・

粗結合 (chain-of-state法)を利用した超並列化

超並列化されたQM/MD計算による反応自由エネルギー

「京」での並列性能(QM)
並列化率 99.9888%
総数32,768コア

!"#

RHF/cc-pVDZ

実行効率





並列化率 99.9888%
実行効率 7.27〜2.24％









 



 

$
&

$%

(8192コア：4.26%)
CASCI(16,16)/6-31G**
総数16,384コア
並列化率 99.9728%

 

光合成反応中心スペシャルペア(280原子） !"

実行効率 15.43〜6.80％





(8192コア：7.85%) 

 # 



#
%

$

「京」での並列性能(QM/MM-MD)
実行効率

98,304コアのハイブリッド並列の実行 !"

RHF/cc-pVDZ

Strong scaling
総数2,048コア 





並列化率： 99.7893%
実行効率 10.60〜1.94％







 








#





総数98,304コア
Ace-ALA-NMe+
1,080 水分子

128プロセス x 768レプリカ

 

Weak scaling (Chain of State)

実行効率 3.42〜3.20％
Weak scaling効率：93.36%

 !" 

 #$ 

#$#
%

$

励起状態QM/MDシミュレーション
sirius発色団-643.2Ground state
1st excited state-643.2

Total energy / a.u.

-643.3-643.3-643.4-643.4-643.4-643.5-643.5

0

5

10

15

20

25

Time / pssirius発色団(1,374原子) QM: 発色団
SA-CASSCF(8,8)/MINI-4 MM: 水∆t = 0.25 fs
120,000ステップ

30

V. まとめタンパク質のQM/MM計算
1) タンパク質は巨大で、ヘテロで、ダイナミックな系
2) タンパク質の機能の発揮は局所部位の電子状態に依存する「どのようにしてモデル化するか?」
功：計算量の減少 vs. 罪：情報の欠損
・QM/MM法の適用高精度化した半経験的分子軌道法によるQM部分の高速化
超並列ab initio QM/MM-MD連成計算プログラムPlatypusの開発

・DFT汎関数の検証（遷移金属タンパク質）
・溶媒分子の取扱→PCMなど

Protein as ampliﬁer
・光合成反応中心のスペシャルペア
タンパク質なし
タンパク質あり
実験値

タンパク質が電子非対称性を強める
・ヘムエリスリンタンパク質がFeとO2の相互作用を強める

タンパク質環境の電気的・立体的摂動によって、補欠分子の固有に持つ性質が強められ、制御されている

Protein as ampliﬁer
タンパク質
溶媒環境

料理人

活性中心
h0p://www.michiba-shunsara.jp/menu/sz044.html

包丁共同研究者
ヘムエリスリン光合成反応中心・中村 春木 （大阪大）
・山崎 秀樹 （大阪大）
rPM6開発・齋藤徹 （広市大）・北河康隆 （大阪大）

Platypus開発
・中村 春木 （大阪大）
・中田一人 （NEC）
・米澤康滋 （近畿大）
・山中秀介 （大阪大）

外部資金
新学術領域「コンピューティクス」
新学術領域「３Ｄ活性サイト科学」
JST-CREST「ライフサイエンスの革新を目指した
構造生命科学と先端的基盤技術」「次世代生命体統合シミュレーション
ソフトウェアの研究開発」