エネルギー表示法による溶解度の計算、および光毒性予測/VCDによる構造決定古賀 良太
株式会社クロスアビリティ 代表取締役
rkoga@x-ability.jp
共同研究者：
エネルギー表示法： 松林伸幸、櫻庭俊（京都大学）
光毒性予測/VCD（計算、ソフト開発）： 山口徹（TSテクノロジー）、古川祐貴（クロスアビリティ）、安田耕二（名古屋大学）
2012年1月17日 第323回CBI学会研究講演会

• エネルギー表示法(ER)– 京都大学・松林准教授が2000年に考案した高速な自由エネルギー計算手法。応用論文あり– 2011年10月から弊社と北浦教授によるJST(ASTEP)『フラグメント分子軌道法を中核とした量子創薬計算システムの構築』で共同研究中。
• GAMESS FMO-Amber-ERのAmber-ERのI/Fを開発• 光毒性予測／VCD– 某社のブラインドテストを契機として、弊社開発のGPGPU QMが活用できる分野であると考えている。応用論文なし• 本発表は背景と現状説明

2

『創薬段階における物性評価の重要性』
日薬理誌 (Folia Pharmacol. Jpn.) 127, 213-216(2006)• 薬理活性の増強だけに注力して探索研究を進めてしまうと、化合物特性に物性情報を反映させるような展開が困難。
• 迅速な溶解性評価がほしい。計算時間かかると意味がないのでlogP（cheminfoがベター）– 迅速な溶解性評価と分配係数logDによる経口吸収性の脂溶性が必
要とされる– 精査段階で結晶からの溶解度が極めて低くなる場合がある– 結晶性や結晶形を収集することが重要で、結晶多形の事前評価が望まれる※ 面による速度論の違いもあり平衡論だけでは難しい
3

エネルギー表示法によるアミノ酸側鎖アナログ計算
amino acidanalog soluteamino acidanalog soluteamino acidanalog soluteAla

methaneVal

propaneLeu

iso-butaneIlen-butaneSermethanolThrethanolPhetolueneTyrp-cresolCysmethanethiolMet

methyl ethyl sulfide

AsnacetamideGlnpropionamideTrp

3-methylindoleHid

4-methylimidazoleHie

4-methylimidazole

平均誤差

free energy ∆µ of hydration (kcal/mol)

4

実験値

2

現実問題への適用可能性

00.7 kcal/mol-2

0.7 kcal/mol-4エネルギー表示法
力場（ポテンシャ
ル関数）の性能

-6
-8

-10

Experimental
FEP
ER（近似計算、短時間）

0.5 kcal/mol
近似（汎関数）の性能

自由エネルギー摂動法-12
Ala Val Leu Ile Phe Ser Thr Tyr Cys Met Asn Gln Hid Hie Trp（厳密計算、長時間）厳密法に遜色無い精度を数十分の一の時間で達成する自由エネルギー計算手法

Karino et al., CPL
2010

自由エネルギー計算～計算法の発展と理論の発展
溶液・界面・生体分子系
統計力学に基づく計算
構造を表す分布関数やダイナミクス
高速で正確な自由エネルギー計算

量子計算におけるエネルギー計算が難しい（時間がかかる）

MDシミュレーションと溶液理論の結合手法
溶媒和自由エネルギー∆µ = ○○ kcal/mol

現実の使用
に耐える
正確さ

[

近似汎関数

エネルギー分布関数~

ρ̂ e (ε ) = ∑ δ (ε − v(ψ , x i ))i

熱力学積分法（厳密手法）

0-5 ε (kcal/mol)
ε (kcal/mol)

容易に
計算可能

・・・

純溶媒

溶液系

純溶媒系と溶液系を結ぶ（仮想的）中間状態を（多数）用意]∆µ = ∫ dεερe (ε ) − F ρ e (ε )MD
（典型的に）
分子内・分子間相互作用

動径分布関数とエネルギー分布関数
1000

2

親／疎

~
~
100
80

1

溶液のとらえ方005r (Å)r (Å)
原子間距離

10

原子毎

分子丸ごと

e

n

遠／近

500

3

gO-O
OO
gO-H
OH

ρ (ε)/ρ (Å mol/kcal)
エネルギー分布関数

動径分布関数
g(r)
g(r)

3

60
40

親溶媒

20
0
-5

0

ε (kcal/mol)

分子間2体エネルギー ε (kcal/mol)

1 g/cm3・25 ºCにおける
水の動径分布関数

1 g/cm3・25 ºCにおける
水のエネルギー分布関数

自由エネルギー
動径分布関数 g(r)exp(-E/kBT)

距離1次元情報の抽出
Kirkwood, 1930年代
本格的発展, 1970年代

エネルギー表示エネルギー1次元情報の抽出

構造・配置の多次元情報（アンサンブル）

cytochrome c 系の水和解析
solute: horse heart cytochrome c (1HRC),
104 residues and heme, 1748 atoms,
total charge of +7
solvent: 20000 water molecules (TIP3P)
temperature: 300K
タンパク質の分子構造
ゆらぎ (RMSD)

タンパク質固定構造（配座）における
水和自由エネルギー∆µの収束挙動
block average over 100 ps
cumulative average-1490∆µ (kcal/mol)

RMSD (Å)

2.0
1.5
1.0
0.5
0
0

5

10
15
MD time (ns)

20

-1500
-1510
-1520
0

0.5

1.0
1.5
Calculation time (ns)2.0

 固定構造につき、1 nsで ~3 kcal/molの誤差で、水和自由エネルギーを計算
 タンパク質構造（配座）と水和効果の関連
Karino et al.,
 タンパク質（溶質）－水（溶媒）相互作用の引力・斥力成分の役割
JCP 2011

溶媒を扱うタンパク質の自由エネルギー計算

水和自由エネルギーµ (kcal/mol)
solvation∆free
energy ∆µ (kcal/mol)

平衡ゆらぎにおけるタンパク質構造と水和の関係

タンパク質分子内エネルギーと
水和自由エネルギーの補償

-1350
-1400

Eintra + ∆µ = const

-1450
-1500

 スクリーニング／モデリングでは・・・
serine protease (1PPE)

-1550
-1600
-1650
-1700

-300
-200
-100
soluteintramolecular energy Eintra (kcal/mol)
タンパク質分子内エネルギー

Eintra (kcal/mol)

x-ray

rmsd 5.5 Å

rmsd 9.5 Å 平衡ゆらぎにおける引力効果（主に水素結合）の支配
 大振幅運動における斥力効果
少数のスナップショット配置でスクリーニング（主に排除体積効果）の効果

タンパク質の水和自由エネルギー解析が全原子レベルで可能に

「溶質」「溶媒」「溶媒和」概念の拡張～統一的な概念構成に向けて「溶媒」＝溶液系に最初からあるもの「溶質」＝溶液系に後から入ったもの

e通常の意味での溶媒和

脂質膜・ミセルへの分子結合タンパク質への
基質結合

電子の付加
（還元）

溶液理論構成の要件
 内部自由度のある分子
 不均一系
 混合溶媒e通常の意味での溶媒和

脂質膜・ミセル タンパク質への基質結合への分子結合 気体様低密度（1対1の「結合」）から
液体様高密度領域をカバー
 量子論との結合（QM/MM法）
 非物理的挙動を導かないような数学的に堅牢な定式化

電子の付加
（還元）

エネルギー表示溶液理論の構成
スナップショット配置における2体エネルギー分布（ヒストグラム）
ρ̂ (ε ) = ∑ δ (ε − v(ψ,x i ))ereference solvent

10ixi : i番目の溶媒分子の座標
v(y,x) : 溶質－溶媒間の2体相互作用
ε: 2体エネルギー値（分布関数の横軸）

 同じ溶質－溶媒相互作用を持つ配置（構造）をグループ化（射影）する

10
10

分布関数

ψ: 溶質の座標

溶媒和自由エネルギー∆µ = − k BT log exp (− β ∫ dε ερˆ e (ε ))
3
2

10

1

~
~

10

2.0

10

1.5
1.0

10

分布

0.5

状態密度

0

10

-5
0
5 105 1010
 どのような近似でも、同じ溶質－溶媒相互作用を2体エネルギー ε (kcal/mol)持つ配置（構造）が溶媒和自由エネルギーに同じ重みで効くという統計力学の原理を破らない ⇒ 精度の高い計算

 どのようなポテンシャル関数に対しても1次元の座標 ⇒ 計算スピード向上
 溶質および溶媒分子が、共に、全体として1つのものとして扱われる⇒ 内部自由度 ⇒ 外場 ⇒ 不均一系・QM/MM系⇒ 低密度／濃度極限 ⇒ 超臨界流体・混合溶媒
 DFT（密度汎関数理論）的な理論構成 ⇒ 系統的定式化

-2
-4
-6
-8
-10

Matubayasi et al.,JCP 2006, JCP 2008

ミセルとリン脂質二重膜
DMPC (1,2-dimyristoyl-sn-glycero-3SDS (sodium dodecyl sulfate)ミセル
phosphatidylcholine)二重膜
hydrophobic tailhead group
water

0.04

0.04

density (Å )

0.03

-3

-3

density (Å )

0.03

0.02

0.02

0.01

0.01

0

0
0

5

10

20
15
r (Å)

OSO3hydrophobic

headgroup

25

headgroup

glycerol

tail

30

water

0

5

10

15

20

25

30

(Å)
r z(Å)

H2O

regionaqueous

各領域における溶媒和自由エネルギーの算出

I

II

COIII

IV
V
NMe +

VI

3

O
O
CO

H2O
O

PO3-

hydrophobic glycerol headgroupaqueousDMPC二重膜20-2-4methanebenzeneethylbenzene-6-8∆µ (kcal/mol) (kcal/mol)
溶媒和自由エネルギー∆µ (kcal/mol) (kcal/mol)
溶媒和自由エネルギー

SDSミセル
2
0
-2

CO
CO2
benzeneethylbenzene

-4
-6
-8

I

II

hydrophobic

III

IV

OSO3

-V

headgroup

VI

H2O

bulk

aqueousI

IICO IIIIV 3+ V
NMe

O
O
CO

VI

bulk

H2O
O

PO3-

hydrophobic glycerol headgroupaqueous 疎水性領域における顕著な安定化
 疎水性・親水性領域の違いは、ミセルの場合により顕著 分子結合の強度（分配係数）とサイトを決定可能
 疎水性溶質の分布は、ミセル中でより局所的
結合寿命も評価可能
 ミセル・膜内部と水領域での、溶質の安定性の逆転

還元過程の解析－電子１個を溶質と見るー
溶質種：電子
混合溶媒種：電子が付加される分子（イオン）（普通の意味での）溶媒（水など）

溶液内

ΨN

ε QM = ΨN +1 H QM ΨN +1 − ΨN H QM ΨN

酸化状態

∆µ

e溶液内

ΨN+1
e-

還元状態
Takahashi et al.,CPL 2008, JCP 2008

ε MM = v(ΨN +1 , x ) − v(ΨN , x )exp(− β ∆µ ) = exp − β ε QM + ∑ ε MM,i  
i

 N


電子分布は、溶媒配置によって揺らぐエネルギー座標を構成する相互作用は多体的
ρ(εQM)とρ(εMM)（及び関連相関関数）から、∆µを構成
FAD (Flavin Adenine Dinucleotide) の活性部位（イソアロキサジン環） の還元への適用∆µQM = -39.5 kcal/mol∆µMM = -40.6 kcal/mol
大きな水の効果（波動関数部分に匹敵）ERのまとめ•

Gromacs・NAMDに対応、Amberは対応中••

GUIはフィアラックス社のMF Amberに移植予定

利点（積分方程式主体のRISMと比べて）•自由エネルギーが定量的に計算できる。

•

量子化学計算を用いた場合、電子雲を部分電荷に落とさず、電子雲のまま適用できる。

•

分子間サイトの関係がexplicitな構造座標として表す必要がないので、サイト間の相互作用の近似が不要で、rigidでも溶質flexibleでもどちらでもよい。また、均一系であるか不均一系であるかを問わないので生体膜や界面活性剤にも適用できる。

MM-ERの評価（１）： 概要
C1の溶媒和自由エネルギー（GAFF、溶質flexible、
Tip4p/amber99sb）: -22.2422±0.1073 kcal/mol※計算時間は１化合物４コアで４～５時間程度（５割がMD）

J.Chem.Inf.Model.
2011, 51, 69-82
MMPBSA -> ER

c1/1000water : -22.2422±0.1073
c2/1000water : -22.4811±0.1422
c3/1000water : -27.9471±0.2576
c4/1000water : -25.2569±0.1451
c5/1000water : -39.5630±0.3022
十分なサンプリングをしている
c6/1000water : -33.8582±0.2325①c1 + 1000water = 1ns, 10fs毎に c7/1000water : -33.5672±0.2735
c8/1000water : -31.5105±0.1425
snapshot ----> 105 frames
c9/1000water : -32.3265±0.2298
c10/1000water : -4.1927±0.1029②1000water = 10ps, 10fs毎にc11/1000water : -27.4158±0.1881
snapshot ----> 103 frames
c12/1000water : -32.0990±0.1704③c1単体 = 10ns, 10fs毎にc13/1000water : -31.8347±0.2275
c14/1000water : -31.0157±0.1804
snapshot ----> 106 frames
c15/1000water : -37.2129±0.2295
計算時間短縮は溶質rigidにする、
c16/1000water : -31.3534±0.1880
サンプリングを必要最低限に減らす、 c17/1000water : -32.7385±0.1674
c18/1000water : -34.5272±0.2568
MD & ERをGPUで動かす etc

MM-ERの評価（１）： MM-ERとMM-PBSAとの比較
検証：MM-PBSAのsolvationをERに変更することで、結果は改善するか？
0
-430

-420

-410

-400

y = 0.034x + 9.0
R² = 0.18

-1
-3
-5

J. Chem. Inf. Model. 2011, 51, 69–82のPB部分をERに変更してプロット

-25

-2
-4

DG( MM-ER)-390

DG( exp)

DG( exp)

-440

0
-20

-15

-10

-5

y = 0.067 - 4.1
R² = 0.09

-1
-2
-3
-4
-5

-6

-6

-7

-7

-8

DG(MM-PBSA)

J. Chem. Inf. Model. 2011, 51, 69–82

結果：MM-PBSAのsolvationをERに変更することで、結果は改善した
1.R2乗値は0.09から0.18に改善
2.実験的に結合自由エネルギーが上位3個の化合物の順位は、
MM-PBSAは5位、10位、12位、MM-ERは1位、9位、3位と改善している。

-80Winmostarでsdfファイル読込

#!/bin/bash

Gromacs-ERのバッチ処理（前半）mpd
acpype -i $1.pdb -n $2
cd ./$1.acpype
cp ../c1.acpype/water.xtc ../c1.acpype/water.log .
cp ../c1.acpype/c1.mdp .
mv c1.mdp $1.mdp
cp ../c1.acpype/c1w_*.mdp .
cp ../c1.acpype/c1w.top .
sed s/c1/$1/g c1w.top > $1w.top
rm -f c1w.top
mv c1w_nvt.mdp $1w_nvt.mdp
mv c1w_npt.mdp $1w_npt.mdp
mv c1w_run.mdp $1w_run.mdp
grompp -maxwarn 10 -f $1.mdp -c $1_GMX.gro -p $1_GMX.top -o $1.tpr
mpiexec -l -n 4 mdrun -s $1.tpr -c $1.gro -x $1.xtc -g $1.log
genbox -cp $1_GMX.gro -box 3.2 3.2 3.2 -cs tip4p.gro -maxsol 1000 -o $1w.gro -p $1w.top
grompp -f $1w_nvt.mdp -c $1w.gro -p $1w.top -o $1w_nvt.tpr
mpiexec -l -n 4 mdrun -s $1w_nvt.tpr -c $1w_nvt.gro
grompp -f $1w_npt.mdp -c $1w_nvt.gro -p $1w.top -o $1w_npt.tpr
mpiexec -l -n 4 mdrun -s $1w_npt.tpr -c $1w_npt.gro
grompp -f $1w_run.mdp -c $1w_npt.gro -p $1w.top -o $1w_run.tpr
mpiexec -l -n 4 mdrun -s $1w_run.tpr -c $1w_run.gro -x $1w_run.xtc -g $1w_run.log
20Gromacs-ERのバッチ処理（後半）（前ページの続き、、、）
mpiexec -l -n 4 mdrun -s $1w_run.tpr -c $1w_run.gro -x $1w_run.xtc -g $1w_run.log# ../../tools/gromacs/gen_structure -I/opt/gromacs-4.5.3/share/top -t $1w.top$HOME/enerepr-0.1.214/tools/gromacs/gen_structure_2 -I/opt/gromacs4.5.3/share/top -t $1w.top $1
cd soln$HOME/enerepr-0.1.214/tools/gromacs/gen_input -x ../$1w_run.xtc -log ../$1w_run.log$HOME/enerepr-0.1.214/ermod_gromacs
cd ../refs$HOME/enerepr-0.1.214/tools/gromacs/gen_input --traj ../water.xtc --log ../water.log--flexible ../$1.xtc$HOME/enerepr-0.1.214/ermod_gromacs
cd ../$HOME/enerepr-0.1.214/slvfe > slvfe_$1.log
cd ../※enereprは京大版ERソースコード（GPLで３月末公開予定）
21

MM-ERの評価（２）： SAMPLのchlorideの評価
7

20

6
15

5
4

10

3

ethane_calc

2ethane_expt1

ethane_calc-expt

5biphenyl_exptbiphenyl_calc-expt

0

0
-1

biphenyl_calc

0

2

4

6

0

8

2

4

6

8

10

12

-5

-2
-3

-10

ethaneのみOPLS-AA
12
10
8
6
dioxa_calc

4dioxa_expt2dioxa_calc-expt

0
-2
-4
-6

0

2

4

6

8

10

横軸： chlorideの数
縦軸： kcal/mol
Gromacs 4.5.3、力場Amber-GAFF※ethane単体（横軸0）だけOPLS-AA+ R.E.DとGAMESS-USで電荷を決めることで精度が改善した。
結果は、Clが無い場合に自由エネル
ギーがプラスに出ており、その後、Cｌの数にほぼ比例して更に増えていく。 これに対して、実測ではClの数への依存性は小さい。

MM-ERを使ってみて• MDをGromacsで行うことでトラジェクトリを得て、それを使い、複数のリガンドセットやSAMPL（オープンアイ社による溶媒和自由エネルギーコンテスト）の提供
化合物群を用いて、ERで溶媒和自由エネルギーを得た。
• ΔΔGの比較のためレンジも違うし全うな評価はできないが、そもそも化合物／蛋白質間の相互作用のところで誤差があり、精度向上には限界があるので、MM-PBSAよりMM-ERの方が改善されたといえる。
• SAMPLで活用したAmber-GAFFでは、Clが、単なる斥力源になっている。
原因について。Amber-GAFFのLJパラメータの元になったOPLS-AAのChlorobenzeneのパラメータを見ると、Clに非常に強い-0.180の電荷が付けられている(opls_264)。 つまり、OPLS-AAでは、Clに強い負電荷を載せることを前提に設定したためLJが弱くなっており、 そのままパラメータを持ち込んだので、電荷決定をRESPで行うAmberと噛み合っていないことが推測される。
力場の問題であり、ERの理論自体の問題によるものではない。

ERの欠点• FEPに比べて適用範囲が狭い– 現時点でリガンドバインディングはうまくいってない（サンプリングスキームの改善を検討中）• 現状、一番良い結果でR2=0.4程度– ρとρ0が重なっていないと使えない（粗視化モデルとの融合により解決可能？）• 溶媒和自由エネルギーしかやらない• MDの欠点を改善するものではなく、引き継ぐ– MDとの組み合わせによる計算手法である• 理論がちょっと難しい– ただ、理論を理解しなくても使える

24

固体の溶解度• 結晶多形の溶解をERをつかって計算できるかもしれないが、速度論の扱いが必要かもしれない。研究中。
• 低分子の溶解度よりも実用上では重要かもしれない。
• テラヘルツスペクトルと基準振動解析（量子化学計算）の比較によるペレット（混合結晶）の構造予測が可能。
25

（以下、参考）テラヘルツ領域の分光
物質の分子間相互作用に由来する未開拓領域
1010

1011

1012

マイクロ波

1013
赤外

テラヘルツ
分子間相互作用回転
振動（共有結合由来）
THz spectrum

帰属困難
Broad

テラヘルツ領域における基準振動解析
2-aminopyridine ペレット
50

17

56

平面構造の基準振動(cm-1)
Twist 19.8
Out-of-plane bend 59.5
Asymmetric stretch 66.0-> (×0.847)
16.8, 50.4, 55.9

量子化学計算による基準振動解析の結果（孤立系）と、実験値の比較

テラヘルツ領域における基準振動解析

ねじれ16.5 cm-1

面外ずり 67.7

面外変角 28.4

面内変角／
ずり 77.4

非対象伸縮
63.6

面内伸縮／
変角 109.4

テラヘルツ領域における基準振動解析

ねじれ16.5 cm-1

面外ずり 67.7

面外変角 28.4

非対象伸縮
63.6

In-plane
赤外領域の局所モード
shear/bend（アミノ基伸縮）
3905.577.4

面内伸縮／
変角 109.4

既存手法の補完
分子間相互作用および構造解析で使われているサービス• Biacore（表面プラズモン共鳴）
相互作用強度、結合速度（ノンラベリング）• プロテイン・チップ（マイクロアレイ）
ハイブリダイゼーションの程度を蛍光強度で測定•

ファーマコフォア・モデリング
化合物データベースを用いて、リガンド－タンパクの構造サンプリング•

リプレイス・補完例
既存のデータベース、もしくは計算で構築する化合物情報データベースに基づき、テラヘルツスペクトルを入力として、相互作用部位の挙動、曲率計算（二次微分）による相互作用強度計算、周辺環境情報を出力する。計算だけで安価。

光毒性予測/VCDによる構造決定• 光毒性予測は第一励起エネルギー、VCDによる構造決定はVCDスペクトルを計算する。
• 計算結果のコンフォメーション依存性が大きいために、多数のコンフォメーションを発生させ、最安定なコンフォメーションを見つける必要がある• そのためには配座探索と高速な量子化学計算が不可欠である• 計算化学の専門家ではなく、安全性・物性の専門家が担当するケースも多く、windowsのUIを含めたユーザビリティが重要になる

32光毒性もVCDもQMで同時計算可• HOMO-LUMO Gapでdescriptorを扱う報告あり– 光毒性には第一励起エネルギーが関係している– 計算コストの観点から第一励起エネルギーと相関値の出るHOMO-LUMO Gapを採用– ただ、圧倒的に第一励起状態はHOMO-LUMOから構成される化合物のみ有効
Journal of Molecular Structure: THEOCHEM 719 (2005) 191-200• 第一励起エネルギーが最小コストで計算でき
るTDDFT計算が必要と考えているが、芳香族系炭化水素ならAM1-ZINDOがベター。
33

半経験的による励起エネルギー計算• 最も高速でそれなりなのはoptをAM1、励起エネル
ギーをZINDOにする。Gaussianの入力例%chk=temp.chk#n AM1 Opt--Link1–%chk=temp.chk#n ZINDO Geom=Check td=(nstates=5)※おそらくMOPACの方が高速• 半経験的方法だとパラメータがない場合があるし、あってもNa以降の原子を含む系では構造最適化の精度もガクンと落ちる。
34

QMによる励起エネルギー計算
TDDFT計算（光毒性スクリーニング用ナミキライブラリの一部）# B3LYP/6-31G* td=(nstates=10)の第一励起エネルギー

4.5634(eV)

2.7711(eV)

3.7267(eV)Journal of Molecular Structure: THEOCHEM 719 (2005) 191200では、7.2±0.7 (eV)でHOMO-LUMO gapのdescriptorを作ったとあるので、励起エネルギーとの絶対値とはダイブ違う値。
35

VCDによる絶対構造推定の流れ
各安定配座のエネルギー計算（ボルツマン分布）

計算による
スペクトル

配座探索ソフト
MOPAC
G09, GAMESS

マッチング

機械による
スペクトル

36

VCDによる絶対構造推定• X線で絶対構造はわかるが、X線解析しづらい
サンプルに対してVCDは有効• Gaussian社のウェブサイトに記述あり
http://www.gaussian.com/g_whitepap/vcd.htm– 標準物質からの相対エナンチオマー純度の決定– 絶対構造決定– 溶液中の配座決定
Biotools社、VCD計測機
37

エナンチオマーとジアステレオマー• エナンチオマー間でエネルギー差はない、VCDで区別するしかない。
• ジアステレオマーと幾何異性体(cis/trans)はエネルギー差がある。
• IRでジアステレオマー、VCDでエナンチオマーを区別するため、QMは有効だとGaussian社
は言っている。
http://www.gaussian.com/g_whitepap/vcd.htm
38

某社のブラインドテスト（１）構造リジッドアロマターゼ阻害剤

光学異性体毎にMOEのlowmode MDで発生した8種類の配座を提供
上記のVCDスペクトルがどちらの光学異性体のものかを予想

N
N

N

N

N

Br
O

N

Br
N

(R)-1-[(4-cyanophenyl)(3-bromo-4-hydroxyphenyl)-methyl]-1H-[1,2,4]triazoleON(S)-1-[(4-cyanophenyl)(3-bromo-4-hydroxyphenyl)
39-methyl]-1H-[1,2,4]triazole下記計算でテスト（１）に合格
aromataseはrigidなので同じconfigurationの配座のvcd傾向は一致しているとみて、sとrのsdfの中の最安定な配座のみ
b3lyp/6-31g* freq=vcd一点計算したところ、1400cm-1付近のピークが明らかにr体とs体で逆。方向一致しているr体と判断。
J. Med. Chem., 2008, 51, 4226-4238

方向が一致

R体

S体

40

某社のブラインドテスト（２）構造フレキシブル＋tautomer

H

H
H N

N

H

GT-2331
N3-H type209 conformers

N

N

H
H

GT-2331
N1-H type670 conformers

2種類のTautomerについてMOEのlow mode MDで配座を発生（合計879配座）
それぞれのTautomerの存在確率は？
41

Values

テスト（２）の問題

備考

Ei [kcal/mol]

0.8 相対エネルギー

E0 [kcal/mol]

0.0 相対エネルギー

T [K]298.15 温度

g [-]

1

N(Ei)/N(E0)

状態数(singlet:g=1,
triplet:g=3)0.258986383 存在比

kB [J/K]

1.381E-23 ボルツマン定数

Avog

6.02E+23 アボガドロ定数

ToJou

4.185851821 ジュール変換
存在比

存在割合

N(E0)

1.000

79.4%

N(Ei)

0.259

20.6%

• 真空中の計算で構造決定して良いか？
• 量子化学計算の精度で十分か？

 ei − e j 
Ni gi= exp−Nj gj
 k BT – エネルギー差・温度・状態数が分かれば存在比が計算できる。0.8kcal/molの差で20%変わる。
– よって、存在比をエネルギーから正確に議論する場合は、1kcal/mol以下の精度が必要（hartreeで点6桁以上の精度）。doubleの計算が必要。
42

計算時間はかけられない• 溶媒効果の計算は時間かかるので省きたい• 構造最適化に時間がかかるから、できる限りエネルギー一点計算で済ませたい• 励起エネルギーもCCSD/6-31G**くらいでちゃんと計算したいが無理• ab initioコンフォメーションサーチするにしても、B3LYP/6-311G*くらいが限界と考えるべき43

報告によると• 下記論文では、実験値はN1:N3=65:35といっている。
Minick et al., Chirality, 19:731-740 (2007)• 真空中の計算では90-95:5-10という報告あり
de Vries AH, van Duijnen PT. Theoretica calculationof tautomer equilibria of 4-substituted imidazoles in
the gas phase and in solution. Biophys Chem
1992;43:139 147.
だが、ダイブ古い計算の文献なので信用できない 44※ 計算はTSテクノロジー社協力

１．課題
GT-2331 N1-H typeとN3-H typeの存在比を推算

H

H
H
N

N H
GT-2331
N1-H type

Optimized Structure at B2PLYP/6-31G*

H
H N

N
GT-2331
N3-H type

Optimized Structure at B2PLYP/6-31G*
45

２．準備

イミダゾイル基

tautomer(互変異性体)の存在比計算で考えられる可能性

NNH・ 物質の平衡状態として互変異性体が存在している場合
R
R⇒ 今回のGT-2331はこのタイプ
tautomer間のエネルギー差を計算し、ボルツマン分布を計算することで、存在比を推算できるかもしれない。

HN
N

・ 互変異性を起こさせるために何か化学的な操作（塩基を加える等）が必要な場合⇒ ケト-エノールなど
化学的な操作の過程（操作=反応過程の遷移状態）を計算、もしくは、原料からの合成過程でtautomerが発生する過程の遷移状態を計算し、その活性化エネルギー差を議論に加えなければならない。
（遷移状態の計算は非常に大変）
ボルツマン分布(分子iとjの比率)

 ei − e j 
Ni gi= exp−Nj gj
 k BT (e: 全エネルギー、g: 状態数、kB: ボルツマン定数、T: 温度)46

３．最もオーソドックスな手法(手法1)
手順
１．2D構造を参考に、3D構造を作成
２．MM2でフルコンフォメーションサーチ。最安定構造を抽出
３．QM(PM6 or B3LYP or B2PTLYP)で構造最適化（真空中）
４．ボルツマン分布より存在比を計算

Level
PM6
B3LYP/6-31G*
MP2/6-31G*//B2PLYP/6-31G*
B2PLYP/6-31G*

Type

Energy [a.u.]

N1-H1
N3-H1
N1-H1
N3-H1
N1-H1
N3-H1
N1-H1
N3-H1

0.10369
0.10367-654.95029-654.95470-654.18211-654.18687-653.47183-653.47635

ΔE
N(N1)/N(N3)
Ratio [-]
[kcal/mol][-]
0.01
0.9812
49.5%
0.00
50.5%
2.77
0.0094
0.9%
0.00
99.1%
2.99
0.0064
0.6%
0.00
99.4%
2.84
0.0083
0.8%
0.00
99.2%

結果
・計算レベルによって存在比が異なるが、全レベルでN3-H1 typeが安定であると計算されている。
・PM6はN3,N1でエネルギー差が殆どなく1:1で存在と推算。
・ab initio計算ではいずれもN1-H1は殆ど存在しないと推算。
・エネルギー差が3 kcal/mol 離れてしまうだけで、存在比に大きく影響してしまう。
・量子化学計算の精度の問題（3 kcal/mol以下の精度を議論できるか）、測定溶媒等の影響（計算結果は真空中でのもの。溶媒を考慮すればさらに精度が悪くなると考えられる。）
47

４．コンフォメーションを意識した手法（最適化なし）(手法2)
手順
１．MOEのlow mode MDで配座を発生(N1:670配座、N3:209配座)
２．MM(MMFF94s)とQM(PM6)で全ての配座に対してエネルギーを計算
３．最安定配座のエネルギーより存在比を計算
結果
・MM計算でのエネルギーはN1とN3間で23kcal/mol程度開きがあり、信頼できないと考えられる。
・PM6は80%がN1と推算。
・MMFF94sでは、エネルギーの開きが大きいため存在比の計算ができないが、総じてN3が安定であると推算。

48

4．コンフォメーションを意識した統計処理的な手法（最適化なし）(手法1)
Energyとして異性体中の最安定構造を採用した場合の結果

Level
MMFF94s 1SCF
PM6 1SCF (H.O.F.)

Type
N1-H1
N3-H1
N1-H1
N3-H1

EnergyΔE
N(N1)/N(N3)
Ratio [-]
[kcal/mol][kcal/mol][-]
40.58867
23.05
0.00000.0%
17.53611
0.00
100.0%
71.47418-0.81
3.947779.8%
72.2873
0.00
20.2%

N1、N3各タイプでのコンフォメーションにおけるエネルギー値の分散(縦軸は％で標準化)

GT2331 PM6 1SCF (H.O.F.)

GT2331 MMFF94s 1SCF

60.0%

80.0%
70.0%
60.0%

40.0%
30.0%

N1 type20.0%

N3 type

10.0%

頻度(%)

頻度(%)

50.0%

50.0%
40.0%

N1 type30.0%

N3 type

20.0%
10.0%

0.0%

0.0%
70.0 70.5 71.0 71.5 72.0 72.5 73.0 73.5 74.0 74.5

17.0 20.0 23.0 26.0 29.0 32.0 35.0 38.0 41.0

エ ネルギ ー (kcal mol-1 )エ ネルギ ー (kcal mol-1)

49

５．コンフォメーションを意識した手法（最適化あり）(手法3)
手順
１．MOEのlow mode MDで配座を発生(N1:670配座、N3:209配座)
２．MM(MMFF94s)とQM(PM6)で全ての配座に対して構造最適化計算を実施
３．最安定配座のエネルギーより存在比を計算
結果
・PM6での異性体の構造最適化計算の結果は、N1 : N3 = 76 : 24 となり、実測結果を支持していると考えられる。
（実験系（溶液中）のIRでの結果は、N1 : N3 = 65 : 35）
最適化なしの結果(N1 : N3 = 80 : 10)より実測値へ近づいた。
・MM計算でのエネルギーは依然としてN1とN3間で平均10kcal/mol程度の開きがありあまり信頼できないと考えられる。また、N3-H typeの方が安定と計算され、実測とは逆の結果。
・これらはあくまで真空中の結果であるので、正確に計算するためには、測定溶媒等の影響も考慮に入れる必要があると考えられるが、ラフなスクリーニングとして使うには充分だと考えている。
N1、N3の全配座の最安定構造(PM6)を初期構造として、
ab initio計算を行えば、結果は改善される！？
（云わば、PM6でのフルコンフォメーションサーチ）

50

５．コンフォメーションを意識した統計処理的な手法（最適化あり）(手法3)
Energyとして異性体中の最安定構造を採用した場合の結果

Level

Type
N1-H1
N3-H1
N1-H1
N3-H1

MMFF94s OPT
PM6 OPT (H.O.F.)

EnergyΔE
N(N1)/N(N3)
Ratio [-]
[kcal/mol][kcal/mol][-]
25.8438.97
0.00000.0%
16.8710.00
100.0%
63.97674-0.69
3.231376.4%
64.67128
0.00
23.6%

N1、N3各タイプでのコンフォメーションにおけるエネルギー値の分散(縦軸は％で標準化)

GT2331 MMFF94s OPT

GT2331 PM6 OPT (H.O.F.)
80.0%
70.0%
50.0%
40.0%

N1 type30.0%

N3 type

20.0%
10.0%
0.0%
63

63.5

64

64.5

エ ネルギ ー (kcal mol-1 )

65

65.5

頻度(%)

頻度(%)

60.0%

90.0%
80.0%
70.0%
60.0%
50.0%
40.0%
30.0%
20.0%
10.0%
0.0%

N1 typeN3 type

16

18

20

22

24

26

28

エ ネルギ ー (kcal mol-1)

51

６．PM6でフルコンフォメーションサーチ→ab initio (手法3アドバンス)
手順
１．手法3によりN1、N3の最安定構造(PM6)を計算(= PM6でフルコンフォメーションサーチ)（N1については211番目、N3については57番目の配座が最安定であった。）
２．QM(B3LYP or B2PTLYP)で構造最適化（真空中）
３．ボルツマン分布より存在比を計算

Level
PM6(手法3の結果)
B3LYP/6-31G*
MP2/6-31G*//B2PLYP/6-31G*
B2PLYP/6-31G*

Type

Energy [a.u.]

N1-H1
N3-H1
N1-H1
N3-H1
N1-H1
N3-H1
N1-H1
N3-H1

0.10195
0.10306-654.95215-654.95464-654.18428-654.18684-653.47418-653.47631

N(N1)/N(N3)
ΔE
Ratio [-][-]
[kcal/mol]-0.69
3.2313
76.4%
0.00
23.6%
1.56
0.0713
6.7%
0.00
93.3%
1.60
0.0669
6.3%
0.00
93.7%
1.34
0.1045
9.5%
0.00
90.5%

結果
・PM6での最安定配座を初期構造として、DFT計算を行った結果、N1-HとN3-H間でのエネルギー差
は手法1より小さくなったが、依然N3が安定と計算され、実測とは逆の結果となっている。
・MP2、B2PLYP(B3LYP + MP2)での計算結果も概ね、B3LYPの結果と同じであった。
・本手法では実験結果を支持する結果が得られなかった。
52

７．計算時間（PM6構造最適化計算）
PM6 構造最適化(EF PRECISE)
ソルバ: MOPAC2009 (アカデミックフリー)
http://openmopac.net/MOPAC2009.html※PM6 はGaussian09でも計算できます。
マシン：Intel(R) Core(TM)2 Quad CPU
計算時間：
合計(1コア使用時)
平均(1コア使用時)

Q6600 @ 2.40GHz

N1-H1
8386 秒 (2.3 時間)
12.52 秒N3-H1
3622 秒 (1.0 時間)
17.33 秒合計
12008 秒 (3.3 時間)

８．まとめ・オーソドックスな方法（N1, N3の構造最適化によるエネルギーを比較）では、ab initio の結果でも、実測値を再現できなかった。
・それぞれの配座(N1:670配座、N3:209配座)に対して、構造最適化を実施しながらエネルギーの平均値を計算した結果、PM6法で、実測値を再現する結果となった。（ただし、これは真空中での構造最適化計算であり、測定溶媒等の影響は考慮されていない。）
・PM6でラフスクリーニングができる可能性を示していると考えられる。
・ab initio計算では、1000～1万倍の計算時間が必要であると考えられるため、計算コスト的に不利であると考えられる。
・PM6での最安定配座→ab initioの計算では、実験結果が支持できなかった。
53

トータル計算時間
普通のデスクトップ１台で担当者が計算した概算時間※第四段階は最安定エネルギーから1kcal/molの配座を選択• 第一段階（構造構築, ISIS Draw）• 第二段階（配座探索, MOE）
– まず、水素付加、MMレベルの構造最適化 - 30sec– 次に配座探索(LowModeMD)を実施- <5min• 第三段階（semi emprical計算）– semi empricalの構造最適化を実施 - <5min• 第四段階（DFT計算, Gaussian09）– 全41配座のDFTレベルの構造最適化計算（低精度:B3LYP/6-31G*）• 1配座にかかる計算は1hour– 全41配座のDFTレベルのエネルギー計算（高精度:B3LYP/6-311+G**）• 1配座にかかる計算は1hour
54

テスト（２）の現状• 今回の計算– MM2フルコンフォメーションサーチ→最安定選択→abinitio（手法1） NG– PM6フルコンフォメーションサーチ→最安定で比較（手法3） OK– PM6フルコンフォメーションサーチ→最安定選択→ab initio(手法3AD) NG• 議論点– 測定溶媒の効果– 手法1と手法3ADで、ab initio時のエネルギー差が縮まっているので、
ab initioは初期構造にかなり依存するので、難しい。けれども、全配座ab
initio計算するのは時間的に大変なので、GPGPU QMが使えるかも– 全配座のPM6は意外と使えるかもしれない。

• 将来に向けて– 配座発生の方法を検討中– 全配座ab initioを実施中
55

某社担当者のコメント
ブラインドテストを実施した某社からのコメント• 配座探索アルゴリズムだけで最安定配座を絞りこんだり、半経験的分子軌道法では最安定配座を探索すると失敗する。むしろ、配座探索では可能性のある多数の配座を発生させるに留め、ab initioで最安定（＋α）の配座を探すことが大切である。
• 化合物構造の柔軟性が高かったり、Tautomerがある場
合には、可能性のある配座が多数になり、自社の計算リソースに与える負荷が高くなる。XA社のGPGPU QMソリューションは費用対効果に優れており、VCDによる構造決定には有効であると考えている。
56

GPGPU QMソリューション
• 弊社と名大の共同開発のプログラム
XA-CHEM-SUITE : XA-CUDA-QM• テラヘルツの基準振動解析、光毒性予測のためのTDDFT(HOMO-LUMO Gap)、VCDいずれも多量の量子化学計算を伴う。
• HF, DFTがマルチコアCPUに比して３～４倍高速• 現状、GAMESSまたはGaussianにアドオン可能• 価格対効果は３倍57balloonによる配座探索• 高速というよりもコストがかからない方法• 有償のものを含めていろいろあるが、今後はMOEがなくても使えるように、フリーで使える
balloonとの組み合わせを検討– MOPACのPM6やGAMESS/Fireflyなどを組み合わせ
れば完全にフリーになる（vcdはGaussian09のみ）– 100個の配座を生成するコマンド例# balloon -c 100 --randomSeed 51277 --output-file 3D.sdf -o sdf -input-file 2D.smi --noGA58

光毒性予測/VCD評価システム
H
H
N

N H

配座発生
（balloon）

GT-2331
N1-H type

再安定な配座から1kcal/molの構造を抽出し
HOMO-LUMO gapとVCDスペクトルを計算(e.g. TDDFT + g09 freq=vcd)

各配座について漏れなく
semi-empiricalエネルギー一点計算(MOPAC)
低精度ab initioエネルギー一点計算（GAMESS on GPGPU計算）

安定な配座について
高精度量子化学計算エネルギー＆最適化計算
同一ネットワーク内

Windows操作マシン（ユーザーはこちらのみの操作）
配座探索ソフトballoon
Winmostar + MF Amber

Linux演算サーバ（ユーザーは意識する必要がない）
ジョブスケジューラ
Gaussian09 B.01

MD/ER + QM(freq/vcd, tddft)• 京大版ERソースコードはGPLで３月末公開予定• プリビルド計算サーバ販売予定– GUIから溶質、溶媒、力場を決めれば計算可能– 弊社が一次サポート• 研究の議論はermod_user（アカデミックML）に代理投稿、もしくは松林研の許可を得てermod_userに参加– GPUでQM + MD(Amber / Gromacs)が加速
Windowsコンソールマシン
Winmostar+ MF Amber (ER GUI add-on)

Linux計算サーバ with GPU¥998,000（予価）MD/ER + QM
60