タンパク質からみる生命科学

計算生命科学のための量子化学基礎
東京大学 生産技術研究所
佐藤 文俊
2014/11/11
計算生命科学の基礎

2013年 ノーベル化学賞

http://www.nobelprize.org/nobel_prizes/chemistry/laureates/2013/press.html

2Agenda

 背景
 電子と原子核 分子軌道法 タンパク質量子化学計算における諸問題
 現状と将来展望

3

背景

4

量子化学と量子化学計算
 量子化学– 量子力学原理を化学の問題に適用する分野– 分子構造・物性・反応性の理論的解析

 量子化学計算– 量子力学方程式を解くことで量子化学の問題を取り扱う– 分子全体に広がる分子軌道を計算して分子の電子状態を求める→ 別名を分子軌道計算– 十分な信頼性があることが立証されている

7 November 2014

Diels-Alder 反応
 共役ジエン(ブタジエン)にアルケン(エチレン)が付加して6員環構造を生じる反応

RNA ポリメラーゼ

チューリングマシンと同じ働き（NHK, 1999）
平成26年11月7日

ATP 合成酵素

モーターと同じ働き

http://www.miraikan.jst.go.jp/focus/sciencetopics/06/index.html
平成26年11月7日

光合成反応中心タンパク質
BClp
スペシャルペア(P)

M鎖

L鎖

L鎖

M鎖
H鎖

光電素子と同じ働き
9酵素反応
 典型的な化学反応– 量子化学こそが直截的な解析法

平成26年11月7日

電子と原子核

11

量子力学の成立前後 I
 1900Planck の定数 1913

Bohr の原子モデル

 1926

Schrödinger 方程式

 1927

水素分子イオン

 1927原子価結合法 (Heitler, London)

7 November 2014

本講義における数式の規約
 表示– ボールドベクトル、行列、テンソル、演算子– ローマンスカラ 暗黙の型– 小文字

電子に関する物理量、インデックス– 大文字

原子核に関する物理量、インデックス

平成26年11月7日

定常状態のSchrödinger方程式
𝑯𝛹 = 𝐸𝛹
𝑁𝑎𝑡𝑜𝑚𝑯=−ℏ22𝑀𝐴𝑁𝑒𝑙𝑒𝑐∆ 𝐴 −𝐴=1𝑁𝑒𝑙𝑒𝑐 𝑁𝑒𝑙𝑒𝑐+𝑖=1
𝑗>𝑖𝑖=1𝑒2𝒓𝑖 − 𝒓𝑗

ℏ22𝑚𝑁𝑒𝑙𝑒𝑐 𝑁𝑎𝑡𝑜𝑚∆ 𝑖 −𝑖=1

𝑁𝑎𝑡𝑜𝑚 𝑁𝑎𝑡𝑜𝑚+𝐴=1𝐵>𝐴𝐴=1𝑍𝐴 𝑒 2𝒓𝑖 − 𝑹𝐴𝑍𝐴 𝑍𝐵 𝑒 2𝑹𝐴 − 𝑹𝐵

Born-Oppenheimer近似 (≃断熱近似)
𝛹 𝒓, 𝑹 = 𝛹 𝑒𝑙𝑒𝑐 𝒓, 𝑹 𝛸 𝑛𝑢𝑐 𝑹



原子単位

1 a.u. ==ℏ = 𝑒2 = 𝑚 = 1𝑁𝑒𝑙𝑒𝑐

𝑯𝑒𝑙𝑒𝑐 =𝑖=11− ∆ 𝑖 −2𝑁𝑎𝑡𝑜𝑚𝐴=1𝑍𝐴+
𝒓𝑖 − 𝑹𝐴27.21 eV627.51 kcal/mol
𝑁𝑒𝑙𝑒𝑐

𝑗>𝑖1𝒓𝑖 − 𝒓𝑗

Born-Oppenheimer 近似における原子核項
 原子核の動きは電子に比べてゆっくり– 原子核の座標変化ごとに電子の方程式は定常状態として扱う
𝑁𝑒𝑙𝑒𝑐

𝑯𝑒𝑙𝑒𝑐

𝛹𝑒𝑙𝑒𝑐

1− ∆ 𝑖 −2=𝑖=1

𝑁𝑎𝑡𝑜𝑚𝐴=1

𝑍𝐴+
𝒓𝑖 − 𝑹𝐴

𝑁𝑒𝑙𝑒𝑐

𝑗>𝑖1𝒓𝑖 − 𝒓𝑗

𝛹 𝑒𝑙𝑒𝑐= 𝐸 𝑒𝑙𝑒𝑐 𝑹𝐴 𝛹 𝑒𝑙𝑒𝑐
全電子エネルギー

𝐸 𝑇𝑜𝑡𝑎𝑙𝑹𝐴𝑁𝑎𝑡𝑜𝑚 𝑁𝑎𝑡𝑜𝑚= 𝐸 𝑒𝑙𝑒𝑐

𝑹𝐴+𝐴=1𝑁𝑎𝑡𝑜𝑚𝑯𝑛𝑢𝑐 𝛸 𝑛𝑢𝑐 =−𝐴=1

𝐵>𝐴

1∆ 𝐴 + 𝐸 𝑇𝑜𝑡𝑎𝑙 𝑹𝐴
2𝑀𝐴

𝑍𝐴 𝑍𝐵
𝑹𝐴 − 𝑹𝐵
𝛸 𝑛𝑢𝑐 = 𝐸𝛸 𝑛𝑢𝑐

Newton方程式と”ちから”
 Newton方程式
𝑁𝑎𝑡𝑜𝑚

𝐴=1

𝒑2𝐴+ 𝑈 𝑹𝐴
2𝑀𝐴𝑑2 𝑹𝐴𝑀𝐴= −𝛻𝑹𝐴 𝑈 𝑹𝐴2𝑑𝑡=𝐸

 原子核座標 (=分子構造) の変化に対して、電子はそれ以上速やかに定常状態へ
𝑁𝑎𝑡𝑜𝑚 𝑁𝑎𝑡𝑜𝑚𝑈 𝑹𝐴=

𝛹 𝑒𝑙𝑒𝑐∗ 𝑯𝑒𝑙𝑒𝑐 𝛹 𝑒𝑙𝑒𝑐 𝑑𝒓 +
𝐴=1

𝐵>𝐴

 Born-Oppenheimer 分子動力学法
𝑑2 𝑹𝐴𝑇𝑜𝑡𝑎𝑙 𝑹∴ 𝑀𝐴=−𝛻𝐸𝑹𝐴𝐴
𝑑𝑡 2
平成26年11月7日

= 𝑭𝐴

𝑍𝐴 𝑍𝐵
= 𝐸 𝑇𝑜𝑡𝑎𝑙 𝑹𝐴
𝑹𝐴 − 𝑹𝐵

Born-Oppenheimer MD・構造最適化の模式図
初期構造,制御パラメータ

位置、制御パラメータ
分子動力学
MD力、エネルギー…

𝑭𝐴 = −𝛻𝑹𝐴 𝐸 𝑇𝑜𝑡𝑎𝑙 𝑹𝐴

構造の時系列, 最適構造
エネルギー時系列…

分子軌道法
MO

M2タンパク質
H+H+流入→ 脱殻→ ウィルスのリボ核タンパク質が宿主核内に侵入

分子軌道法19

量子力学の成立前後 II
 1929

Dirac の提言

 1930Hartree-Fock 方程式 (Fock) 1931

Hückel 法

 1932“molecular orbital”の命名 (Mulliken)

7 November 2014

Dirac の提言
物理学の大部分と化学の全体の数学的理論に必要な基礎的物理法則は完全にわかっているということであり、困難は、ただ、これらの法則を厳密に適用すると複雑すぎて解ける望みのない方程式に行きついてしまうことにある。
したがって、量子力学を応用するための実用的な近似方法を発展させ、過度の計算を行うことなしに、複雑な原子集合体の主だった性質を説明できるようになることが望ましい。
P. A. M. Dirac, Proc. Roy. Soc. (London), A123, 714 (1929).

7 November 2014

分子軌道計算の曙













1938-55
1946
1950
1951
1951
1952
1953
1963
1965
1965
1965
1965

小谷テーブルコンピュータの誕生 (ENIAC)
Gauss型関数 (Boys)
Hartree-Fock-Roothaan 方程式 (Roothaan)Xα 法 (Slater)
ab initio HF 法
フロンティア軌道理論 (福井)Pariser-Parr-Pople 法
拡張 Hückel 法 (Hoffmann)
HITAC5020E 運転開始 (東大大型計算機センター)
エチレンの ab initio HF 計算 (Moskowitz, Harrison)
CGTO (藤永)Kohn-Sham 方程式
DFT (密度汎関数)法

7 November 2014

HITAC 8500

様々な分子軌道法
 経験法– Hückel, 拡張 Hückel 法

 半経験法– PPP, CNDO, INDO, MINDO 法– 多中心分子積分に近似を導入

 ab initio Hartree-Fock (HF) 法– もっとも標準的な分子軌道法、電子相関なし– 電子相関を入れるためには CI 法（莫大な計算量が必要）

 密度汎関数法, Density Functional Theory (DFT)法,
Kohn-Sham (KS) 法– 電子相関効果をポテンシャルに組み込んだ方法– ab initio HF 法と同程度の計算量

Roothaan 方程式• 微積分方程式 → 行列 (固有値) 方程式• 自己無撞着 (Self-Consistent Field; SCF) 計算
で固有値 𝝐 と固有ベクトル 𝑪 を求める方法

𝑭𝑪 = 𝑺𝑪𝝐𝐹𝑝𝑞 = ℎ𝑝𝑞 +𝐽𝑝𝑞 +𝛼 𝐾𝑝𝑞 +𝛽 𝜇𝑝𝑞
1 電子項
クーロン項

ℎ𝑝𝑞 =

𝑔𝑝 𝒓

𝐽𝑝𝑞 =

1− ∆−
2𝐴𝑍𝐴𝒓 − 𝑹𝐴𝑃𝑟𝑠 𝑝𝑞 𝑟𝑠𝑟𝑠

Fock 交換項

𝐾𝑝𝑞 =
𝜇𝑝𝑞 =

基底関数 (given) 𝑔𝑝 𝒓

𝑃𝑟𝑠 𝑝𝑟 𝑞𝑠

𝑆𝑝𝑞 =𝑔𝑝 𝒓 𝑔𝑞 𝒓 𝑑𝒓

𝑔𝑝 𝒓 𝜇𝑋𝐶 𝜌 𝒓 𝑔𝑞 𝒓 𝑑𝒓

𝑃𝑝𝑞 =𝐶𝑝𝑖 𝐶𝑞𝑖𝑟𝑠

交換相関項

𝑔𝑞 𝒓 𝑑𝒓𝑖

2 電子 (分子) 積分𝑝𝑞 𝑟𝑠 =1𝑔𝑝 𝒓 𝑔𝑞 𝒓𝑔𝑟 𝒓 𝑔𝑠 𝒓 𝑑𝒓𝑑𝒓′𝒓 − 𝒓′

𝜌 𝒓 =𝑃𝑝𝑞 𝑔𝑝 𝒓 𝑔𝑞 𝒓𝑝𝑞

インスリン最高占有軌道(HOMO)

等電子密度面上の静電ポテンシャル
最低空軌道(LUMO)

閑話休題：
実際の分子で分子軌道を可視化

NATURE, VOL. 432, 867

Molecular orbital wavefunction of N2a. tomographic image reconstruction
b. 2p sg orbital

平成26年11月7日

2010 年ノーベル化学賞
 クロスカップリング– 金属等を触媒にした化学基(base)の反応

現行の化学論文の例
Angewandte Chemie (IF=13.734 (2012)) より試作(合成)して→ 構造解析して→ 量子化学計算

タンパク質量子化学計算における諸問題

29

プログラム開発と大規模計算の始まり


1966

IBMOL の開発 (Clementi)



1970

Gaussian70 の開発 (Pople)



1972

TNF 錯体 (50原子) の ab initio HF 計算 (Clementi)



1974

MOLECULE 開発とポルフィンの計算 (Almlöf)

1976

水の分子動力学計算 (Lie, Clementi, 吉嶺)



1976

スーパーコンピュータ CRAY-1 の開発



1979

分子研電子計算機センター運転開始
HITAC M-180(IAP)x2



1981

福井ノーベル化学賞受賞



1992

Gaussian92 が DFT 法をサポート



1998

Pople, Kohn ノーベル化学賞受賞



2000

シトクロムc (1,738原子) の DFT 計算 (佐藤, 柏木)



2005

インスリン6量体 (4,728原子) の DFT 計算 (稲葉, 佐藤)



2013

Karplus, Levitt, Warshel ノーベル化学賞受賞

ナノスケールの分子機械
 特異な反応を省エネルギーで– (密に詰まった)複雑かつユニークな電子状態を利用

 大きく分けて3つのアプローチ– 重金属(レアメタル)• 大量生産技術• 資源確保・汚染– ナノ構造• 特徴ある構造• 生産・加工技術– タンパク質• タンパク質工学• 解析・設計技術

4

定量性ある高精度計算の必要性
体内温度における標準自由エネルギー変化数~十数 kcal/mol

𝐼𝑃 = 𝐸 𝑁 − 1 − 𝐸(𝑁)
電子を抜くのに・・・
鉄イオン ～1000 kcal/mol
ヘム～100ヘムタンパク質＜1032

生体分子量子計算の挑戦

QM/MM, ONIOM:• Gaussian

FMO:• GAMESS• ABINIT-MP

カノニカル分子軌道計算:• ProteinDF

33

Roothaan 方程式 (再出)
• 微積分方程式 → 行列 (固有値) 方程式• 自己無撞着 (Self-Consistent Field; SCF) 計算
で固有値 𝝐 と固有ベクトル 𝑪 を求める方法

𝑭𝑪 = 𝑺𝑪𝝐𝐹𝑝𝑞 = ℎ𝑝𝑞 +𝐽𝑝𝑞 +𝛼 𝐾𝑝𝑞 +𝛽 𝜇𝑝𝑞
1 電子項
クーロン項

ℎ𝑝𝑞 =

𝑔𝑝 𝒓

𝐽𝑝𝑞 =

1− ∆−2𝐴𝑍𝐴𝒓 − 𝑹𝐴𝑃𝑟𝑠 𝑝𝑞 𝑟𝑠𝑟𝑠

Fock 交換項

𝐾𝑝𝑞 =𝜇𝑝𝑞 =基底関数 (given) 𝑔𝑝 𝒓

𝑃𝑟𝑠 𝑝𝑟 𝑞𝑠

𝑆𝑝𝑞 =𝑔𝑝 𝒓 𝑔𝑞 𝒓 𝑑𝒓

𝑔𝑝 𝒓 𝜇𝑋𝐶 𝜌 𝒓 𝑔𝑞 𝒓 𝑑𝒓

𝑃𝑝𝑞 =𝐶𝑝𝑖 𝐶𝑞𝑖𝑟𝑠

交換相関項

𝑔𝑞 𝒓 𝑑𝒓𝑖

2 電子 (分子) 積分𝑝𝑞 𝑟𝑠 =1
𝑔𝑝 𝒓 𝑔𝑞 𝒓
𝑔𝑟 𝒓 𝑔𝑠 𝒓 𝑑𝒓𝑑𝒓′𝒓 − 𝒓′

𝜌 𝒓 =𝑃𝑝𝑞 𝑔𝑝 𝒓 𝑔𝑞 𝒓𝑝𝑞

タンパク質分子サイズと計算資源
 サイズに関するルール
残基数 : 原子数 : 電子数 : 軌道数 (行列次元)≃ 100: 2000 : 7000 : 10000
100残基タンパク質なら1億要素3GB/3行列
1000残基タンパク質なら 100億要素 300GB/3行列
 律速計算 分子積分計算： クーロン項、Fock交換項

 数値積分計算： 交換相関項 行列演算： 行列積、対角化

 並列計算 巨大行列を分散保持、律速計算毎にタスクを分割
 例えば 99.9x % → 数 1000~10000 並列規模

16

SCF の収束
 理想的な収束-2.0

Orbital energy / eV-3.0-4.0-5.0-6.0-7.0-8.00510
Iteration number1521

SCF 計算失敗例
 被占軌道と空軌道の入れ替わり
 ギャップが空いていない

22

計算シナリオと収束法
 タンパク質全電子計算のためのSCF初期値作成法– タンパク質計算を達成するために必須の技術– 計算シナリオ、QCLO法 SCF計算収束法 (本日は省略)– SCF計算を安全・安定に収束へと導く技術• 軌道重なり, プロジェクション, レベルシフト, ミキシング• 特にタンパク質のような巨大な系において重要– SCF計算収束を加速する技術• Anderson, DIIS, GEDIIS, …
• 正しく収束するならば、SCF計算は何回転してもかまわない→…が、収束加速はそのまま計算時間の短縮に貢献
23タンパク質全電子計算のシナリオ
GRGQFPREHTVIAWQNTTY ････ FPRETQH

アミノ酸シーケンス
ステップ1 (1残基)

ステップ2 (3残基)
フレーム分子ステップ3 (数残基)

平成26年11月7日

最終ステップ(全タンパク質)QCLO (Quasi-Canonical Localized Orbital) 任意の空間で局在化し、その空間内ではカノニカル軌道

HIS

THR

VAL

HIS

THR

VAL

HIS

THR

VAL

組み合わせるとより長いペプチド鎖の良い初期値にHIS

THR

VAL

ペプチド鎖 全エネルギー収束過程-2526.85

Total Energy [a.u.]-2526.90

本方法
従来の方法

-2526.95-2527.00-2527.05-2527.10-2527.15
1357

9 11 13 15 17 19 21 23 25 27
Iteration

約20倍解に近い初期値が得られた
41

現状と将来展望42

生化学と量子化学
 生化学– 生命現象を化学的に研究する分野– 「理論やるより試験管を振れ」
20～30年前までは、「化学」でも同じでした…

 量子化学 (再出)– 量子力学原理を化学の問題に適用する分野– 分子構造・物性・反応性の理論的解析– 一般的な低分子では十分な信頼性が立証

7 November 2014

量子化学を生化学に適用する際の壁 そもそも昔のハード・ソフトウェア資源ではできなかった– 構造のモデル化と計算方法の近似• 活性中心周辺のみ扱う (サイズの縮小)– アミノ酸残基(and/or側鎖)と基質– ヘムタンパク質→ ポルフィン＋上下配位子• 計算方法の近似 (信頼性の低い計算法)– 望みのプログラムが存在しない→ 1990年代の製薬会社の失敗

 構造を近似したり信頼性の低い方法を用いては生化学の諸問題に対応できない– 一般化学分子の標準法による生体分子の丸ごと計算
7 November 2014

生体分子量子計算の挑戦 (再出)QM/MM, ONIOM:• Gaussian

FMO:• GAMESS• ABINIT-MP

カノニカル分子軌道計算:• ProteinDF

45

タンパク質を含めた全電子計算の重要性
実験では当然の事実
側鎖や官能基の置換、アポタンパク質の付加により活性中心の電子構造が劇的に変化→全領域をQMで扱う必要性

MbCO差電子密度分布(上)側鎖を置換
Fe-ポルフィン→Fe-プロトポルフィリン(下)周辺ペプチド数の比較
3残基→14残基
2014/11/7

全電子計算例- ヘムタンパク質の分子軌道 -active site

HOMO
CO-binding to myoglobin(MbCO)

Mb：Fe(III)で失活
HOMOの成分94.7%： Feのd軌道99.9%： Heme+近位HisHOMO
cytochrome c

Cyt.c：Fe(II)⇔Fe(III)で電子移動
HOMOの成分85.2%： Feのd軌道主成分の同じHOMOでも分子軌道の広がりに違いが見られる– 機能の違いと相関(?)

47全電子計算例- インスリンの静電ポテンシャル タンパク質の相互作用に関わる重要な物理量

全電子計算から高精度な静電ポテンシャルが算出可能

ProteinDFによる全電子計算結果

古典計算結果(各原子の電荷はアミノ酸のRESP電荷を使用)

全電子計算と古典計算の差

主流の古典計算では、misleadingを招く恐れがある48

カノニカル分子軌道計算のサイズ
30000

密度汎関数法

基底関数の数

25000(306残基,4728原子)インスリン6量体

Hartree-Fock法

20000(104残基,1738原子)

15000
10000
5000

シトクロームcRNAフラグメント
クランビン
プロテインキナーゼC

01996

P53単量体
1998

2000

2002
西暦

2004
2006
2008
(注)論文の掲載年にて整理

機能性タンパク質 と ”モノづくり”
集光機能
広い波長範囲で光エネルギーを捕捉することで太陽エネルギー収集デバイスとして展開可能

光電変換機能高効率で電子を取り出すこと
で、実用的なバイオ発電システ
ムとして展開可能

アンテナタンパク質
CO2固定機能律速状態を打破する(1000倍でパーフェクト・ルビス
コ)ことで、CO2削
減や、食糧、バイオ燃料といった諸問題への展開可能

e-ルビスコ

光合成反応中心タンパク質

e-ヒドロゲナーゼ

H2発生機能酸素による失活を防止する
ことで、水素を大量製造するクリーンなデバイスとして展開可能

タンパク質研究のための道具の現状
 タンパク質工学– ほぼ完成 (力技も使える状況)– 自在なアミノ酸置換、大量生成

 遺伝子解析・システムバイオロジー– ターゲットタンパク質の絞り込みなどに有用

 立体構造解析– X 線や NMR 構造解析、古典分子 (動) 力学– なぜ機能するかという、後付的ビジョンを与える

 問題点– 加工技術はほぼ完成、計測技術も進んでいる– 制御・予測の技術が遅れているリアル系にフィードバックするシミュレーション 製造系の最適化プロセス

実機試作

フォーミュラカーの非定常空力解析(FrontFlow/Red)
床・車両相対速度の影響

 バイオ系の最適化プロセス？
??

酵素精製・遺伝子改変酵素試作

遺伝情報／ホモロジーサーチ 構造解析／MD

Diels-Alder 反応 (再出)共役ジエン(ブタジエン)にアルケン(エチレン)が付加して6員環構造を生じる反応
 この反応の機構は量子化学によって明確に説明される


共役ジエンの最高被占軌道(HOMO)とアルケンの最低空軌道(LUMO)の間の相互作用によって、2つのσ結合がそれぞれの分子のπ軌道面の同じ側で協奏的に生成

共役ジエンに HOMO のエネルギー準位を上昇させる電子供与性基、アルケンに LUMO のエネルギー準位を低下させる電子吸引性基を置換すると、軌道相互作用が大きくなり反応が加速

分子軌道、電子密度、状態密度、相互作用解析タンパク質の階層構造と電子状態との間で，詳細かつ直接的な解析

2014/11/754

文藝春秋 2010 年 9 月号より
 京大

山中伸弥 教授

ES細胞だけで作用する遺伝子を片っ端から全部見つけたいと思っていました。…(自分で一から作るしかないと思っていた)ところが、理研の仕事の中にES細胞由来のデータも入っていた。データを解析するプログ
ラムも、NCBIで無料で公開されていた。

 ソニー

大賀典雄 氏(1978年) 当時、ソニーではちょうどウォークマンが売れ出した時期でしたが、エンジニアとして世の中の誰も知らないこと(CD開発)をやるのが面白くないはずがない。

量子力学とタンパク質工学の融合が技術革新を生み出す

参考文献

ISBN-13: 978-4627241411
基礎理論・教科書

ISBN-13: 978-4627879119 ISBN-13: 978-4130624541
シミュレーションのシミュレーションソフトのノウハウ
作り方 (基礎編)