3D01
ELG-OPT 法によるタンパク質・ＤＮＡの効率的安定構造探索（九大院・総理工 1，JST-CREST2）○青木 百合子 1, 2，Liu Kai1，Xie Peng1，折本 裕一 1

An efficient searching of optimized geometries for proteins and DNA
by ELG-OPT method（Kyushu Univ.1，JST-CREST2）○Yuriko Aoki1, 2，Liu Kai1，Xie Peng1，Yuuichi Orimoto1【序】バイオ系やナノ系など大規模系の電子状態の効率的計算方法として、一次元系を想定した
Elongation(ELG)法 1,2)および三次元系のための 3D-Elongation (3D-ELG)法 2-4)を開発してきたが、DNA やタンパク質等生体系の諸現象解明のためには、効率的構造最適化が必須である。本研究では、まず一次元系を仮定した ELG 法にエネルギー勾配法を導入した ELG-OPT 法を開発してきた。
性能評価のため、非結合系、非局在化系、タンパク質や DNA モデルを対象に従来法と効率性の比較を行ったのでいくつかの応用例、および三次元系への 3D-ELG の展開について報告する。
【方法と結果】ELG 法では、ある程度の長さの高分子鎖に対して、攻撃分子を順次付加する過程
で、相互作用の必要な領域のみを対角化に含めることにより高速化している。各 ELG ステップにおいて、軌道を Frozen する部分と Active にする部分があるが、Active 部分のみに対して Gradient
計算の対象とすることができるため、局所的エネルギー勾配計算の繰り返しによって全系の安定構造を探索する。ELG 法の過程で、攻撃分子から遠く離れており、既に Frozen した部分に属する
構造パラメータは、たとえ最適化しても変化しないため固定し、攻撃分子との相互作用部分に関してのみ、領域局在化軌道(RLMO)基底のエネルギー勾配を計算している。以下いくつかのテスト
分子に対して、全系を直接扱った場合との比較で示す。
最も簡単な非結合系として、直線状(HF)n=48 分子鎖の構造最適化を行い、ELG-OPT 法と従来の構造最適化法による結果を比較した。最適化後の安定性から、6-31G (d,p)レベルでは ELG-OPT 法は従来法の OPT 法より安定となり、一部のパラメータでずれが認められた（図１左）。
Non-bonding system: (HF)48
Poly-(HF)48

Final step of optimization iteration

Both Conv and ELG QFMM off
Basis set: Sto-3G, 6-31G, 6-31G(d,p)(au)-4,800.0(HF)6 (HF)6 (HF)6 (HF)6 (HF)6 (HF)6 (HF)6 (HF)6
1 unit0

B unit=4FH

F

1.1Å 1.3Å 0.7Å

H
1.7Å

F

H

F

1.1Å 1.3Å 0.7Å

H

F

H

-4,800.2F

1.7Å 1.1Å 1.3Å 0.7Å

ELG vs CONV bond length difference

エネルギー差 Error= Eelg-Econv

CONVの構造の悪さからくる
21

level

Error/atom

STO-4.49E-07

6-31G

9.09E-09

6-31G(d,p)-3.26E-05

CONVの構造の悪さからくる21

図１-4,800.4

20

40

6-31G(d,p)

60

80

100

ELG-OPT
CONV-OPT-4,800.6-4,800.8

Eelg lower than Econv
Telg shorter than Tconv-4,801.0

ＥＬＧの方が低い

ELG-OPT は、従来法OPTより安定点を与えうる

Total energy

1 unit =(HF)6 =H

21-4,801.2

E: Energy T: Time

ELG-OPT 法と従来法の最適化結合長の差：分子内結合(∆RH-F)、分子間結合(∆RF…H)と収束状況

ELG-OPT 法により得られた構造に対して、従来法で１点計算を行うと全エネルギーが完全に一致することから、最適結合長の差は従来法の構造が本来の安定構造に到達していないことによるも

(a) Initial structure of extended polyalanineのである。
言い換えれば、同じ閾値下では、
図１右に示すように ELG-OPT 法が速く収束

1 unitAfter optimization at 6-31G basis set,

するのに対して、従来法では多くのステッ

ELG
Conv

プを要していることが原因である。その他、∆E= 3.85×10-8 au/atom

共役π電子系等の非局在化系への応用とし(b) Optimization of DNAて平面ポリアセチレンの構造最適化を行い、STO-3G basis set OPTTOL=1.0×10-2 (Default=1.0x10-4)∆E= -8.1 ×10-4 au/atom

最適化構造について従来法と良好な一致がELG-OPT locates a morestable point thanconventional calculation.

見られたため、タンパク質や DNA 等生体高分子系の構造最適化に応用した。
ELG-OPT 法をポリアラニンに適用したと

図2

ころ、ほぼ正確に従来法による最適化構造

ポリアラニンおよび DNA の ELG-OPT および
CONV-OPT による構造最適化

および全エネルギーを再現している（図 2 (a)）
。このように ELG-OPT 法によって得られた最適構
造が CONV-OPT 法による結果と一致する場合は、全エネルギー誤差は～10-8a.u./atom を達成しているが、図 2 (b)の DNA の例のように、ELG-OPT 法によって到達し得た最適構造の方が安定な場合がある(～－8×10-4 a.u./atomELG-OPT 法が安定)。その構造での CONV 法とは完全に一致しているため（～10-8a.u./atom)、マイナスの差は ELG 法から派生したものではなく、従来法が準安定構造
で最適化を停止したことからくる。これは、同じ初期構造から最適化を行なった場合に、ELG-OPT
法では従来法に比べてよりグローバルな安定構造探索が可能であることを示唆しており、従来の構造最適化では到達し得ない安定構造を見出すケースである。巨大生体系では、複雑な相互作用
が様々な局所構造を発生させ、全系をまともに最適化する従来法では、最適化の結果が初期構造に強く依存する。つまり、一旦 Local minimum に落ち込むと、エネルギー勾配法による電子状態のみの最適化ではそこから抜けにくくなるのに対して、ELG-OPT 法では、各部分空間で局所的な安定構造探索をしながら逐次的に最適化を行うため、全系の安定構造探索が容易になることが考えられる。複雑生体分子のグローバル最適化構造探索における本方法の有効性が期待できる。
本方法を一般的に、絡み合った複雑生体系にも適用可能とするため、三次元系への 3D-ELG 法を開発中である。

The accuracy of system with one
disulfide bond calculations at
6-31G basis set

本方法では、図３のインスリンモデル

SH

S-S

のように、高分子鎖が巻いて一旦
Frozen した部分に最接近した場合に、その部分に存在する RLMO を再活性化

PDB code: 1PAKConventional SCF
Step Number Iteration

る S－S 結合を精確に表現していることを確認済みであるが、構造最適化にEnergy deviation

Time to form

Time to solve

Number of atoms Iteration Time to form Time to solveof atomsof SCF

Fock (sec.)

Fock (sec.)

in active regionof SCF

Fock (sec.)

1

107

17

335.6

21.2

107

17

334.6

20.2

0.00E+00

2

160

17

719.7

154.6

160

15

540.5

301.6

-1.54E-10

3

214

17

1305.2

422.7

214

15

959.5

711.4

-1.22E-10

4

263

17

2137.6

701.5

156

15

1361.2687.7

8.66E-07

させて固有値問題に含める手順を自動化したものである。鎖の途中におけELG SCF

図3

Fock (sec.)

三次元系用 3D-ELG 法によるインスリンの計算

ついては現在開発中である。
1) A. Imamura, Y. Aoki, K. Maekawa, J. Chem. Phys., 95 (1991) 5419.
2) Y. Aoki, F. L. Gu, Phys. Chem. Chem. Phys., 14 (2012) 7640.
3) K. Liu, L. Peng, F. L. Gu, and Y. Aoki, Chem. Phys. Lett., 560(2013) 66.
4) K. Liu, Y. Yan, F. L. Gu, and Y. Aoki, Chem. Phys. Lett., 565 (2013) 143.

(au/atom)

3D02
転写因子 p53 の特異的 DNA 配列探索・認識機構：	
 マルチスケールシミュレーション研究	
 (京大院理 1,	 阪大蛋白研 2)	
 	
 
◯寺川	
 剛 1,	 肥後	
 順一 2,	
 高田	
 彰二 1	 Specific DNA sequence search and recognition mechanism oftranscription factor p53: Multi-scale simulation study(Kyoto Univ.1, Institute for protein research, Osaka Univ.2)◯Tsuyoshi Terakawa, Junich Higo, Shoji Takada【序】	  代表的な癌抑制性の転写因子である p53 は、4 つのドメイン (N 末端ドメイン、コ
アドメイン、4 量体形成ドメイン、C 末端ドメイン) で構成される。N 末端ドメイン、C 末端ドメイン、コアドメインと 4 量体形成ドメインをつなぐリンカー領域は天然変性領域であり、静的な構造情報を元にした機能解析を困難にしている。本研究の目的
は、粗視化 MD によって得られる動的な構造情報を元にして、p53 の機能に重要な特異的 DNA 配列の探索と認識に、ドメイン間の協調が果たす役割を明らかにすることである。
【モデルの構築】	  先行研究において、我々は、p53-非特異的 DNA 複合体の粗視化シミュレーションを行うことにより、C 末端ドメインが DNA 上をスライディングする一方で、コアドメインは探索中に DNA との結合、解離を繰り返すという p53 の探索メカニズムを明らかにした[1]。このメカニズムにおいてリンカー領域が認識プロセスにおいて重要な役割を果たすと考えられる。しかし、比較的単純な天然変性領域のためのモデル[2,3]をリンカー領域に適用した結果、過去に実験で得られた SAXS プロファイルを再現することができなかった。本研究では、まず、この領域の全原子マルチカノニカル MDを行い、その結果得られたコンタクトマップを再現するように、粗視化モデルのパラ
メータを決定した (Fig. 1)。このパラメータを用いて粗視化 MD を行うことによって得られた理論的な SAXS プロファイルは実験を再現した (Fig. 1)。特異的配列とコアドメインの相互作用エネルギーのパラメータは、過去の実験で測定された解離定数を再現するように決定した。
【シミュレーション結果と考察】	  次に、得られたパラメータを用いて、p53 と DNA の粗視化 MD[4]を行った (Fig. 2)。
結果として、p53 は C 末端ドメインによって DNA 上を走査することによって特異的配列を素早く探索することができることが明らかになった。また、C 末端ドメインは特異的配列に非特異的に結合することよってコアドメインの結合を妨げていることが示唆された。
[1] Terakawa, Kenzaki and Takada (2012) J Am Chem Soc 134:14555[2] Terakawa and Takada (2011) Biophys J 104:1450[3] Terakawa, Kameda, and Takada (2011) J Comput Chem 32:1228[4] Li, Terakawa, Wang and Takada (2012) Proc Natl Acad Sci U S A 109:17789

3D03
グアニン―シトシン塩基対アニオンにおけるプロトン移動ダイナミクスの理論的研究（埼玉大院理工）○杉岡 雄仁、高柳 敏幸

Theoretical study of proton transfer dynamics ofthe guanine−cytosine base pair anion(Saitama Univ.) ○Yuji Sugioka,

Toshiyuki Takayanagi

放射線が生体に照射されると DNA の損傷が起こることはよく知られている事実である。この損傷は、イオン化エネルギーより大きなエネルギーをもつ２次電子によって引き起こされると考えられてきた。しかし、2000 年に Boudaïffa らのグループが、イオン化エネルギーよりも低いエ
ネルギーを持つ電子が DNA 損傷を引き起こすことを明らかにした[1]。その結果、DNA の放射線損傷の機構が見直されるようになり、低エネルギー電子と生体関連分子の衝突実験が盛んに行われるようになった。これまでの研究結果から低エネルギー電子による DNA の損傷には、解離性
電子付着が重要な役割を果たしていることが分かっている[2]。つまり、電子の付着によって分子が一時的にラジカルアニオンになった後、分子の変形や結合の解離が起こる。我々の研究室では
核酸塩基対のラジカルアニオンに着目し、電子付着後の反応ダイナミクスを理論的に明らかにすることを目的として研究を行なっている。本研究では、グアニン―シトシン塩基対のラジカルア
ニオン(GC•を取り上げる。
Fig.1 はグアニン―シトシン塩基対に電子が付着したときの変化の概略を表している。Fig.1 の下の図は GC•の SOMO を表している。グアニン―シトシン塩基対の双極子モーメントは 6D 程度であるため、電子は最初に塩基対の双極子モーメントに束縛される。その後電子がシトシンの* 軌道に束縛されるに伴い、塩基対の平面構造が歪む。GC•の構造と、そのときのエネルギーと SOMO は BH&HLYP/6−311++G(2d,p)にさらに diffuse 関数を加えたレベルで行なった。この
過 程 の 遷 移 状 態 は ま だ 求 め ら れ て い な い が 、 電 子 の 垂 直 脱 離 エ ネ ル ギ ー (VDE : Vertical
Detachment Energy) はこれまでの研究結果と良い一致を示している[3,4]。
また、塩基対が* アニオンに変化すると塩基間でのプロトン移動反応が起こる。この反応のバ
リアは B3LYP/6−31+G(d)レベルで 3.54 kcal/mol であり、ゼロ点振動補正を行なうと 1.13
kcal/mol であった。我々は特にこのプロトン移動反応に着目し、分子動力学計算を行なった。ポ
テンシャル曲面は B3LYP/6−31+G(d)レベルの反応座標に沿ったデータをもとに内挿法で作成し、動力学計算には核を量子的に扱うリングポリマー分子動力学計算(RPMD)を選択した。さらに、
RPMD の結果を古典分子動力学計算(classical MD)の結果と比較した。RPMD と classical MD に
おいてそれぞれ 30 本のトラジェクトリーを計算した結果、classical MD に比べ、RPMD の方がプロトン移動反応が起こりやすいことが分かった。Fig.2 にその典型的な例を示す。Fig.2 はそれぞれの分子動力学計算における N−H の長さの時間変化を表している。これらの結果は振動の量

子化およびトンネル効果がプロトン移動反応に重要な役割を果たしていることを示唆している。
詳細については当日発表する。

Fig.1 電子付着後の GC•の変化

Fig.2 分子動力学計算における N−H の長さの時間変化[1] Boudaïffa, B. et. al. Science 2000, 287, 1658–1660.
[2] Gu, J. et. al. Chem. Rev. 2012, 112, 5603–5640.
[3] Smets, J. et. al. Chem. Phys. Lett. 2001, 342, 342–346.
[4] Gupta, A. et. al. J. Phys. Chem. B 2012, 116, 5579–5587.

3D04
ノイラミニダーゼに対するシアル酸誘導体の阻害メカニズムに関する理論的研究: LERE-QSAR(九大先導研) ○比多岡 清司、吉澤 一成

Theoretical Studies on Inhibition Mechanism of Sialic Acid Analogues
against Neuraminidases: LERE-QSAR Approach(IMCE, Kyushu Univ.) ○Seiji Hitaoka, Kazunari Yoshizawa【序論】
インフルエンザの治療には, ウイルスの増殖・遊離を担う酵素であるノイラミニダーゼ (NA)を選択的に阻害する目的で開発された zanamivir (Relenza)や oseltamivir (Tamiflu)が代表的な抗インフル
エンザ剤として広く使用されている. しかしながら, 特に Tamiflu 服用後の精神・神経症状の副作用が報告されており, これら NA 阻害剤のヒトに対する影響が懸念される. その原因の一つとして, ヒト・
ノイラミニダーゼ (human NA: hNEU)への影響が指摘され, 抗インフルエンザ剤の hNEU に対する阻害効果の検討が行われてきた. 現在, ヒトにおいては細胞内局在性や基質特異性の異なる 4 種類のヒト・ノイラミニダーゼ (hNEU1–4)が同定されているが, このうち 3 次元立体構造が解明されているのは hNEU2 のみである. その構造的特徴はインフルエンザウイルスの NA と比較的類似している
ため, NA 阻害剤の hNEU2 に対する作用メカニズムをインフルエンザ NA と対照させて理解することは, 高い選択性と副作用の少ない新規 NA 阻害剤の開発につながると考えられる.
本研究では, hNEU2 と抗インフルエンザ剤を含む一連のシアル酸誘導体の複合体について, 非
経験的フラグメント分子軌道 (fragment molecular orbital: FMO)法等による分子科学計算ならびにその結果に基づく自由エネルギー変化の線形表現 (linear expression by representative energy
terms-quantitative structure−activity relationship: LERE-QSAR)解析 [1]から, 複合体形成に伴う全自由エネルギー変化の変動を支配する相互作用様式ならびにその変動に対するシアル酸誘導体の各部分構造の寄与を原子・電子レベルで定量的に明らかにすることを目的とした.
【方法】
Chavas ら [2]により報告されている抗インフルエンザ剤を含む一連のシアル酸誘導体の hNEU2に対する阻害活性データに基づき, 本研究では Relenza や Tamiflu を含む合計 7 化合物を使用した. 欠落アミノ酸残基が最小である X 線結晶解析構造 (PDB code: 1VCU)を各化合物の初期構造として使用し, 分子動力学 (MD, AMBER)計算により各複合体構造を構築した. 得られた複合体構造に対して, FMO (MP2/6-31G)計算, Poisson−Boltzmann (PB)計算および量子化学的分極連続体
モデル (SCRF-CPCM/HF/6-31+G(d,p))計算を実行し, 結合自由エネルギー (Gbind = EbindHF Ecorr: EbindHF および Ecorr は, 静電的結合相互作用エネルギーおよび分散相互作用エネルギーをそれぞれ表す), 水和自由エネルギー変化 (Gsol)および解離自由エネルギー変化 (Gdiss)をそれぞれ得た. これら代表エネルギー項を用いて , LERE-QSAR 式 (Gobs  1 (EbindHF  Ecorr <Gsolpol>)  2 Gdisspol  const)を構築し, 実測値 (Gobs)の予測ならびにその変動を支配する相互作用エネルギーの同定を行った. Gsol およびGdiss については, 極性項 (Gsolpol, Gdisspol)に加えて,非極性項 (Gsolnonpol, Gdissnonpol)の寄与もあるが, これら非極性項の化合物間の変動は極性項のそれと比べてきわめて小さいため, Gsol  Gsolpol およびGdiss  Gdisspol としてそれぞれ評価した. また,Gsolpol は阻害剤の結合部位以外の微小な構造変化に対して敏感であるため, この効果を MD より得られるダイナミクストラジェクトリから統計平均値 (Gsolpol)として算出することで考慮した. なお,1 および2 は回帰分析より得られる線形定数である.
【結果および考察】
一連のシアル酸誘導体と hNEU2 の複合体形成に伴う実測のGobs の変動を, 物理化学的な意味づけとともに説明可能な統計的に有意である LERE-QSAR 式を得ることができた.
Gobs  0.0476 (EbindHF  Ecorr  Gsolpol)  0.0976 Gdisspol  2.25
n  7, r  0.985, s  0.237 kcal/mol, F  65.510図 1 は, Gobs の変動に対する各エネルギー項のEbindHF
Ecorr
Gdisspol

寄与を示す. 高い化学的精度において, 計算値(Gcalc)はGobs をきわめて良好に再現している.
ずかに大きいが, これは Tamiflu が有する amino基の結合と非結合状態における解離自由エネル
ギー差を反映している. Ecorr の寄与は小さく, この項を除外した相関式が良好に成立することを確認

5
Energy (kcal/mol)

Compound 7 (Tamiflu)のGdisspol は他と比べてわGsolpolGobsGcalc0−5

している. 一方で, EbindHF とGsolpolの各々の寄与は大きいものの, 両者の間には良好な逆相関関係が成立する結果として (r  0.969), 静電相
互作用エネルギー (EbindHF  Gsolpol)の寄与
は比較的小さくなる. しかしながら, この寄与はEcorr の寄与と比べて大きいため, 静電相互作用エ

−10

1

2

3

4

5

6

7

Compounds

図 1. シアル酸誘導体と hNEU2 の複合体形成に伴う各エネルギー項の変動 (compound
1 を基準として表す)

ネルギーがGobs の変動に対して支配的であることが示唆される. ここで, EbindHF の寄与がGsolpol
の寄与よりも相対的に大きいことから, 特にEbindHF の寄与が支配的となる. したがって, 一連のシアル酸誘導体と hNEU2 の複合体形成に伴う実測の全自由エネルギー変化の変動に対して, 両分子間における静電的結合相互作用が最も支配的な役割を果たしている.
以上の LERE-QSAR 解析の結果に加えて, Tamiflu のヒト (hNEU2)およびインフルエンザウイルス(N1-NA)のノイラミニダーゼに対する結合選択性の違いについても, シアル酸誘導体をフラグメント分割した詳細な結合相互作用エネルギー解析から当日議論を行う.
[1] Hitaoka, S; Chuman, H. J. Pestic. Sci. 2013, 38, 60–67.
[2] Chavas, L. M. G.; Kato, R.; Suzuki, N.; von Itzstein, M.; Mann, M. C.; Thomson, R. J.; Dyason, J.

C.; McKimm-Breschkin, J.; Fusi, P.; Tringali, C.; Venerando, B.; Tettamanti, G.; Monti, E.;
Wakatsuki, S. J. Med. Chem. 2010, 53, 2998–3002.

3D06

詳細釣り合いを課さないレプリカ置換法の提案（分子研，総研大） ○伊藤暁，奥村久士

Replica-permutation method without
the detailed balance condition(IMS, Sokendai) ○Satoru G. Itoh, Hisashi Okumura【序】 タンパク質の構造や機能を計算機シミュレーションを用いて調べるためには，タン
パク質に対する効率的な構造空間の探索が不可欠である．そのための手法として近年広く
レプリカ交換法 [1,2] が用いられるようになってきた．レプリカ交換法ではレプリカが温度空間をランダムウォークすることにより，シミュレーションが自由エネルギー極小状態に捕らわれることなく効率的な構造空間のサンプリングを実現する．最近我々はより効率的な構造空間のサンプリングを実現する手法としてレプリカ置換法を開発した [3]．
通常のレプリカ交換法ではレプリカ間の温度交換にメトロポリス法を用いるが，レプリカ置
換法ではレプリカ間の温度置換に諏訪・藤堂法 [4] を用いる．諏訪・藤堂法はメトロポリス
法とは異なり詳細釣り合いの条件を満足せずに状態遷移を行うモンテカルロ法であり，状態遷移のリジェクト率を最小化することができる．この方法をレプリカ置換に用いることで，従来のレプリカ交換法と比較してレプリカの温度空間の効率的サンプリングを実現することが可能となった．
【シミュレーション手法】 諏訪・藤堂法を用いて 2 個のレプリカ間で温度の置換を行うこととメトロポリス法を用いて 2 個のレプリカ間で温度の交換を行うことは等価なので，諏訪・藤堂法の利点を生かすためには 3 個以上のレプリカ間で温度の置換を行う必要がある．そこで，レプリカ置換法
では 2 個のレプリカ間での温度交 図１ レプリカ置換法の概略図．各レプリカが占有換だけではなく，全てのレプリカ間する温度の時間発展を示す．
での温度置換を考える．つまり，この方法では図１に示したように従来のレプリカ交換法では起こり得ないようなレプリカの温度遷移を実現することが可能となる．
一般に，M 個のレプリカに対するレプリカ置換法を実行することを考える．このとき，レプリ

カ置換法における状態は X   x1

i (1)

,xm

i ( m ),, xMi ( M )

 で表わされる．ここで，上付きの添

え字 i(m) は下付きの添え字 m に対する置換操作を行う関数であり，上付き・下付きの添え

字はそれぞれレプリカと温度を示す．また， xm   q  , p iii

 で， q, p は座標と運動量を表わmす．さらに，  は M 個のレプリカと M 個の温度の全ての組み合わせに対するラベルで，  1, , M ! である．例えば，レプリカが 3 個の場合は図２に示したように 3!  6 個のレプリカと 温 度 の 組 み 合 わ せ が 考 え ら れ る ． 状 態 X に 対 す る 重 み は

M





w  X     exp   m H xm
m 1i ( m ) で与えられる．ここで， m  1 kBTm で kB はボルツマン定数， Tm は温度， H はハ
ミルトニアンを表わす．
レプリカ置換法では，従来のレプリ
カ交換法と同様に，シミュレーション
中にレプリカ置換トライアルを行う．こ 図２ レプリカ 3 個に対するレプリカと温度の全組
み合わせのレプリカ置換トライアルによってレプリカと温度の組み合わせがからへ遷移する確率 P  X   X   は
P  X  X   v  X  X  
w  X 

で与えられる． v  X   X   は組み合わせからへの総確率流であり，諏訪・藤堂法ではv  X   X    max 0, min   , w  X    w  X     , w  X   , w  X     ,  S  S 1  w  X 1  , S   w  X  , S0  S M !
 1と計算される． ただし，ここで状態 X 1 に対する重み w  X 1  が全での重みの中で最大であるとしている．一般に， w  X   が最大の重みとなる場合については文献[3]を参照のこと．

【結果】 新しく開発したレプリカ置換法と従来のレプリカ交換法を二重井戸型ポテンシャ
ル，真空中の Met-enkephalin，水中の C-peptide の 3 つの系にそれぞれ応用し，その効率を比較した結果を紹介する．図はレプリカがある温度から別の温度への遷移する確率を示
す．図３から，レプリカ置換法では従来のレプリカ交換法と比較してレプリカが頻繁に温度遷移をしていることが分かる．

図３ レプリカの温度遷移の確率．RPM，M-RPM ではそれぞれレプリカ置換トライアルに対して諏訪・藤堂法，メトロポリス法を用いた．REM は従来のレプリカ交換を表わす．
【参考文献】[1] K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).
[2] Y. Sugita and Y. Okamoto, Chem. Phys. Lett. 314, 141 (1999).
[3] S. G. Itoh and H. Okumura., J. Chem. Theory Comput. 9, 570 (2013).
[4] H. Suwa and S. Todo, Phys. Rev. Lett. 105, 120603 (2010).

3D07

ヘリックス構造をもつペプチド分子の圧力依存性に関する理論的研究（1 分子研，2 総研大）○森義治 1，奥村久士 1,2Theoretical study on the pressure dependence of a helical peptide(1Inst. Mol. Sci., 2SOKENDAI) ○Yoshiharu Mori1 and Hisashi Okumura1,2【序】タンパク質の構造とその機能は密接に関連していることから，現在まで構造を探る研究が数多くなされている。最近ではタンパク質の構造を理解するために熱力学的なパラメータとして圧力も用いられるようになってきた。高圧力条件下におけるタンパク質の構造は活性状態の構造と関連があるという報告もある。このようなことから高圧力下における生体分子のふるまいを分子論的な観点から研究することは意味のあることである。
圧力による構造変化では，タンパク質の圧力変成とよばれる，タンパク質構造が高圧力により天然構造から変化するような現象があきらかになっている。一方ではペプ
チドに高圧力をかけることにより，その二次構造が増加するような場合もあることがわかっている（下図参照）。本研究では，そのような二次構造が増加するペプチドの中から AK16 ペプチドを選び，理論的な方法により考察した。

【方法】
本研究では主にアラニンから構成されているペプチドである AK16 ペプチドの構造
圧力依存性を分子動力学シミュレーションにより調べた。前述したように AK16 ペプ
チドは高圧力状況において，その二次構造形成率が上昇することが分かっている。本
研究では高圧力状態においても構造の緩和を促すため，拡張アンサンブル法のひとつである温度・圧力に関する焼き戻し法[1]を適用してシミュレーションを実行した。
系として AK16 ペプチド(Ace-YGAAKAAAAKAAAAKA-NH2)と水分子 4414 個を用意した。考察する圧力の範囲は 0.1 MPa（大気圧）から 1.4 GPa までとした。
この分子動力学シミュレーションの結果から，ヘリックス構造の形成率，部分モル
体積変化およびペプチドの構造の圧力依存性を計算し解析を行なった。
【結果】
シミュレーションの結果，ヘリックス構造の形成率は圧力の増加にともなって，はじめは減少するがその後増加した。また，ヘリックス構造からそうでない構造への部分
モル体積変化は，圧力の増加により負の値から正の値に単調に増加することが分かった（ΔV の図参照）。高圧力条件では，部分モル体積変化が正の値をもつという実験結果を再現することができた。さらにヘリックス構造をとった状態に着目し，この状態での慣性半径を計算すると，これは圧力とともに減少することがわかった（Rg の図参照）。この結果から，圧力によりヘリックス構造は縮んでいると考えられ，これは部分モル体積変化のふるまいを説明する。このヘリックス構造の縮みは，ペプチドの内の残基間距離の計算からも裏付けることができた。これらの解析結果から，圧力の増加にともなうヘリックス構造の正確なふるまいを分子動力学シミュレーションにより明らかにすることができた[2]。

ヘリックス構造をとった状態【参考文献】[1] Y. Mori and Y. Okamoto, J. Phys. Soc. Jpn. 79 (2010) 074003.
[2] Y. Mori and H. Okumura, J. Phys. Chem. Lett. 4 (2013) 2079.

3D08
QM/MM 法によるジオ ールデヒドラ ターゼのミューテーション解 析（1 九大先導研・2 岡山大工）○土井富一城 1・蒲池高志 1・虎谷哲夫 2・吉澤一成 1Computational Mutation analysis for Diol Dehydrataseusing QM/MM Method(1IMCS Kyushu Univ., 2Okayama Univ.) Kazuki Doitomi1, Takashi Kamachi1,
Tetsuo Toraya2, Kazunari Yoshizawa1【緒言】
ジオールデヒドラターゼ (DD)はアデノシルコバラミン由来の酵素であり、1,2-ジオールを相当するアルデヒドへと変換する反応を触媒する。この酵素は本来の基質ではないグリセロール(GOL)の脱水反応を触媒する際に、不活性化されやすいことが知られている。GOLが結合した構造では本来の基質である
1,2-プロパンジオールの場合と違い、GOLの3個目の水酸基がGln336とSer301の近傍に存在し、Ser301と水素結合するが、Gln336とは水素結合していない。虎谷らはこれらの残基が反応に与える影響を検討するために酵素に変異を導入しGln336Ala、Ser301Ala変異型を作製し、野生型よりも不活性化しにくいと報告した[1]。
本研究では野生型が不活性化する原因と変異型が野生型よりも不活性化しにくい理由についてQM/MM法
を用いて検討した。

図 1. 活性中心近傍の構造【計算方法】
酵素の全原子計算を行うために QM/MM 法を用いた。QM/MM 法は高精度で計算負荷が大きい量子化学(QM)計算と安価で高速な分子力場(MM)計算を組み合わせた方法である。QM 領域の原子は反応に重要な 10 個のアミノ酸残基，金属イオン(Ca2+)、GOL、アデノシルラジカルのリボース部位とした。変異型の計算においては Gln336、Ser301
をそれぞれアラニンに置換した。QM 領域では TURBOMOLE プログラムを用いて
B3LYP 法で計算を行い、基底関数は SV(P)を適用した。それ以外の領域を MM 領域
として力場に CHARMm を適用し、計算を DL_POLY プログラムで行った。計算のインターフェイスとして ChemShell プログラムを用いた。
【結果と考察】
図2に示すのはこれまで知られていた脱水反応と本研究により提案された不活性化に寄与する副反応の反応機構である[3]。始めにアデノシルラジカルがGOLの1位の水素を引き抜き、反応物ラジカルが生じる。この反応物ラジカルから水酸基転移反応

(脱水反応)と水素移動反応(副反応)が競争的に進行する。水酸基転移反応では2位の水酸基が1位に転移し、生成物ラジカルになる。一方、水素移動反応では3位の水酸基の水素が1位に移動する水素移動反応が進行し、アルコキシ
ラジカルが生じる。このアルコ
キシラジカルは不安定なのでホルムアルデヒドとグリコールラ
ジカルに分解される。紫外可視分光スペクトルはホルミル基とCoの相互作用の可能性を示して

図 2. 脱水反応(上)と副反応(下)の機構

おり、この副反応を支持する。
また、GOLがDDに結合した構造がふたつ存在すること(GS、GR)が知られており、GR構造では脱水反応のみが起き、GS構造
では不活性化が進行しやすい[2]。野生型のGSとGR構造それぞれにおける水酸基転移
反応と水素移動反応の活性化エネルギーを図3に示す。それらのエネルギー差はGS、GR構造でそれぞれ10.0、14.1 kcal/molである。GS構造はエネルギー差が小さく、水素移動反応が進行しやすく、不活性化しや

図 3. 不活性化と基質結合構造の関係

すいことが明らかにされた[3]。
変異型の反応についても野生型と同様の機構で検討した。Gln336Ala変異型では水素移動反応と水酸基転移反応の活性化エネルギーはあまり変化しなかった。一方、GR構造で結合するGOLの割合が増えることが予測された[4]。GR構造は水素移動反応が進行しにくいので、不活性化が抑制される。この結果は実験による提案[1]と一致している。しかし、Ser301Ala変異型ではGS構造で結合するGOLの割合が多くなった。ところが、GS構造の水素移動反応の活性化エネルギーが高くなることによって、副反応が進行しにくくなることが判明した。当日はこれまでに得られた知見を基に野生型と変異型の反応について議論する。
[1] Yamanishi, et.al. FEBS J. 2012, 279, 793.
[2] Bachovchin, et.al. Biochemistry 1977, 16, 1082.
[3] Doitomi, K.; Kamachi, T.; Toraya, T.; Yoshizawa, K. Biochemistry 2012, 51, 9202.
[4] Doitomi, K.; Kamachi, T.; Toraya, T.; Yoshizawa, K. to be submitted

3D09
Hsp90 の ADP 解離過程における自由エネルギープロフィールと相互作用解析
（金沢大・理工） ○川口一朋、齋藤大明、長尾秀実
Free energy profile and interaction analysis in the dissociation processof ADP from Hsp90(Inst. Sci. Eng., Kanazawa Univ.) ○Kazutomo Kawaguchi, Hiroaki
Saito and Hidemi Nagao【序】Hsp90 (Heat Shock Protein 90)は他のタンパク質（クライアントタンパク質）のフォール
ディングを補助する機能をする分子シャペロンの一つである。ガンやウイルス感染症など様々な病気に関わるタンパク質と関係があり、それらの治療ターゲットとしても注目されている。Hsp90
は二量体を形成し、二量体構造の変化の過程でクライアントタンパク質のフォールディングが起こる。この大規模な構造変化には ATP の結合、ATP から ADP への加水分解、ADP の解離が必須である。Hsp90 は三つのドメイン（N 末端ドメイン、M ドメイン、C 末端ドメイン）から構成
され、各ドメインの構造および二量体構造が X 線結晶構造解析により明らかになっている。
また、
Hsp90 と核酸の結合には二種類の水素結合（直接水素結合と、水分子を介した水素結合）が重要であることもわかっている。ATP の結合による Hsp90 の構造変化が明らかにされているが、ATPの結合および ADP の解離のメカニズムの詳細は分かっていない。そこで、我々は ADP の解離に着目し、Hsp90 と ADP の複合体構造を用いた分子動力学シミュレーション（MD）を実行した。
熱力学的積分法を用いて ADP の解離の自由エネルギープロフィールを求め、Hsp90-ADP 間の水素結合に関する解析を行った。
【方法】MD の初期構造には Hsp90 の基質結合部位である N 末端ドメイン（NTD）とそれに結合する ADP の複合体の X 線結晶構造（PDB ID: 1byq [1]）を用いた。基本セル内に水分子と Na+
イオンを配置し、計 40,673 原子の系となった（図 1）。力場には CHARMM と TIP3P を用いた。
カットオフは 10Å とし、長距離力の計算には PME を用いた。MD の実行には MODYLAS[2]を用いた。
熱力学的積分法を用いて、結合距離の関数として自由エ
ネルギープロフィールを求めた。ここで、Hsp90 と ADPの重心間距離を結合距離とした。熱力学的積分法では反応座標 r が r0 から r1 まで変化する時の自由エネルギー差 Gを以下の式で表すことができる。
r1

G(r0  r1 )   F (r ) r drr0

図 1. MD の初期構造

ここで、 F (r ) r は結合距離 r における平均力であり、MD のトラジェクトリから求める。結合距
離 r として 0.3～1.8, 2.0, 2.5 nm の 18 点を用意し、それぞれに対して結合距離を拘束した MD を実行した。また、各結合距離における水素結合距離の変化を求めた。
【結果と考察】熱力学的積分法により各点（各結合距離）での平均力 F (r ) r を求めた。その結果、結合距離 r
が 0.8 nm 以下の点で急激な斥力の増加が見られた。この領域は平衡距離より内側の領域であり、ADP がタンパク質内部に埋もれている状態である。また、r>1.5 で引力はゆるやかに減少していき 0 に近づく。遠方では引力が徐々に小さくなっていくことがわかる。図 2 に平均力の結果から求めた、結合距離に対する自由エネルギープロフィールを示す。自由エネルギーは結合距離が 0.8
nm で極小点をとることがわかった。また Hsp90-ADP 間の 10 個の水素結合に関して、その変化を調べた。詳細については当日報告する。

図 2. 平均力と自由エネルギープロフィール

図 3. 水素結合距離の変化【参考文献】[1] W. M. J. Obermann, H. Sondermann, A. A. Russo, N. P. Pavletich, F. U. Hartl, J. Cell Biol., 143(1998) 901.
[2] Y. Andoh, N. Yoshii, K. Fujimoto, K. Mizutani, H. Kojima, A. Yamada, S. Okazaki, K.
Kawaguchi, H. Nagao, K. Iwahashi, F. Mizutani, K. Minami, S.-i. Ichikawa, H. Komatsu, S.
Ishizuki, Y. Takeda, M. Fukushima, J. Chem. Theory Comput., 9 (2013) 3201.

3D10
脂質二重層膜におけるグラミシジン A の静電ポテンシャル特性と電場応答（金沢大院・自然）齋藤大明，川口一朋，長尾秀実

Electrostatic potential profile of lipid bilayer containing gramicidin A and its
response to the change of electrostatic field(Institute of Science and Engineering, Kanazawa University) Hiroaki Saito, Kazutomo
Kawaguchi, Hidemi Nagao【 序 】膜タンパク質は生体膜における物質の選択的透過，シグナル伝達，エネルギー変換等の生体機能に直接関わる重要な生体分子であり，これらの機能は膜タンパク質を介したイオン・分子
透過と密接な関係がある．例えば，抗菌性ペプチドとして知られるグラミシジン A は，膜内において二量体を形成することによりカチオンを選択的に透過させるイオンチャネルを形成する事が知られている[1]．これら膜タンパク質のイオン・分子透過機構の解明は，生体内における膜タン
パク質の機能理解のみならず，創薬や新規ナノデバイスの研究・開発における重要課題である．
	  生体膜は脂質分子の種類やその混合割合によって，極性基付近の静電ポテンシャル特性が大きく変わり，これによりイオン透過性も大きく変化することが知られている[1]．一方で，生体内における膜タンパク質を介したイオン伝導は，膜の上下のイオン濃度差による膜ポテンシャルにより起こる．実験では電極による電位差の膜ポテンシャルを人工的に発生させることによりイオン伝導が観測される．しかしながら，タンパク質−脂質二重層膜のような混合複雑系における実験観
測の難しさのために，これらイオン透過特性は未だ明らかではなく，分子シミュレーションによる詳細な解析が望まれている．そこで本研究では膜ポテンシャルを模した一様電場環境（膜厚方向）におけるグラミシジン A の分子動力学シミュレーションを行い、脂質膜内におけるグラミシ
ジンチャネルの静電ポテンシャル特性とその電場応答ついて解析した結果について報告する。
【 方 法 】本研究では，グラミシジン A の脂質二重層膜への添加効果の評価のために，膜タンパク
質−脂質二重層膜系の分子動力学シミュレーションを実行する．脂質分子は DMPC( diC14:0-PC)を用い、これら脂質分子で構成される脂質二重層膜へグラミシジン A を添加させ、MD シミュレー
ションを実行する。MD 計算は等温・等圧条件下で行い，分子力場は脂質/ペプチド系にはCHARMM36 を，水モデルには TIP3P を用いた．膜への電場効果は，原子の部分電荷 qi に膜厚方向への一様電場 E による力(F = qiE)を作用させることにより行った．いずれの系の計算も 50ns ま
でに構造が十分に平衡化している様子が示され，50ns 以降のデータを構造や圧力特性評価に用いた．解析には脂質二重層膜の膜面積(Alipid)や膜厚(dP-P)，脂質分子のオーダーパラメータ(-SCD)やア
シル鎖のゴーシュ構造比(Fgauche)を行った. 静電ポテンシャルの計算は系の電荷分布に対するポア
ッソン方程式を解くことにより評価した．分子動力学シミュレーションには NAMD2.7 を用いた．

(a)

(b)

図 1. 膜厚方向に対する静電ポテンシャルマップ	 (a) 膜ポテンシャル無し, (b) 膜ポテンシャル有り(500mV)

図 2. 膜厚方向に対する静電ポテ
ンシャル特性【 結 果 と 考 察 】 表１に各々の系における脂質膜の構造パラメータ(膜面積，膜厚，疎水鎖領域の厚さ，オーダーパラメータ，ゴーシュ構造比)を示す．グラミシジン添加効果の比較の為に，表の括弧内にグラミシジン無しの系の値も示した．解析の結果，グラミシジン添加により膜面積は減
少し，膜厚と疎水鎖領域の厚さは増加する結果が示され、実験値との良い一致も示された[2]．
	  図 1 にグラミシジンのイオンチャネルを通過する系の静電ポテンシャルマップ(YZ 面)を示す．
図 1 において，正の静電ポテンシャルを示す赤-黄-緑色領域は膜の疎水鎖領域に対応しており，
緑-青色領域は膜の極性基や水和領域に対応している．グラミシジンの重心は図の中心に対応しており，膜を通過する青色領域はイオンチャネル領域に対応している．これらの結果が示すように、カチオンはイオンチャネル内では静電的に安定であり，これらの静電特性がイオン透過の選択制に寄与していることが示された．膜ポテンシャル(500 mV)を加えた系では膜の上下の水和領域において電位差が生じ，チャネル内部を通じたイオン透過経路の様子が示された．図２にイオンチ
ャネル中心を通過する膜厚方向に対する静電ポテンシャルプロフィールの様子を示す．膜ポテン
シャルが無い系（青線）では静電ポテンシャルに大きな変化がないのに対し、膜ポテンシャルがある場合の静電ポテンシャル（赤線）では膜の水和領域やイオンチャネル内部の領域において平坦なポテンシャル値となる結果が示された．これは膜厚方向に対する電場に応答した、水分子や脂質の極性基の再配向が原因であると考えられる．詳細は当日報告する。
【参考文献】[1] Kelkar, D. A.; Chattopadhyay, A. BBA-Biomembranes 2007, 1768, 2011-2025.
[2] de Planque, M. R. R.; Greathouse, D. V.; Koeppe, R. E.; Schafer, H.; Marsh, D.; Killian, J. A.
Biochemistry 1998, 37. 9333-9345.

表 1．膜構造パラメータ. 括弧内は実験値