Implementation of Pair Interaction Energy DecompositionAnalysis and Its Applications to Protein-Ligand Systems

Abstract

フラグメント分子軌道(Fragment Molecular Orbital; FMO)法はタンパク質などの巨大分子の電子状態計算を可能にする方法であり，FMO計算によって得られるフラグメント間相互作用エネルギー(Inter Fragment Interaction Energy; IFIE)がタンパク質-リガンド間相互作用などを理解する上で有用であるため創薬研究などに用いられている．本研究ではFMO計算プログラムABINIT-MP及びそのプリポストBioStation Viewerに，IFIEを静電相互作用エネルギー(ES)，交換反発エネルギー(EX)，電荷移動相互作用エネルギー(CT+mix)，分散エネルギー(DI)の4つのエネルギー成分に分割・解析できる機能PIEDA (Pair Interaction Energy Decomposition Analysis)を実装し，複数のタンパク質-リガンド系で実証計算を行った．その結果，ノイラミニダーゼ-オセルタミビル，EGFRチロシンキナーゼ-エルロチニブ，エストロゲン受容体-リガンド複合体の3つの系がそれぞれ異なる特徴的なタンパク質-リガンド間相互作用を持っていることが示され，IFIEを各エネルギー成分に分割して評価できるPIEDAの有用性が示された．

1 はじめに

近年の目覚ましいコンピューターとソフトウェアの発展に伴い，高精度で実用的なシミュレーションが創薬やものづくりの現場に普及しつつある．フラグメント分子軌道(FMO)法 [1,2,3,4]は，タンパク質などの巨大分子をアミノ酸などの単位でフラグメント(断片)に分割し，複数のフラグメントの組み合わせにおける電子状態を計算することで分子全体のエネルギー等を計算する方法である．並列コンピューティングとの相性が良く，現実的な時間内でのタンパク質の高精度量子化学計算を可能にする上，フラグメント間の相互作用エネルギー(Inter Fragment Interaction Energy; IFIE)が得られるため，医薬品化合物などの設計において有用である．FMO法はこの10年の間に主に生命化学分野において飛躍的に発展し，すでに創薬やワクチン開発の場で利用され始めている [5,6]．

ABINIT-MP [5,7,8,9]はFMO法の主要プログラムの1つであり，その専用プリポストBioStation ViewerとのパッケージであるBioStationとして，東京大学生産技術研究所を中核拠点とした文部科学省次世代IT基盤構築のための研究開発「イノベーション基盤シミュレーションソフトウェアの研究開発」をはじめとするプロジェクトにて開発，公開されてきた [10]．商用利用による普及も進んでいる [11]．本研究では，前述のIFIEを静電相互作用エネルギー(ES)，交換反発エネルギー(EX)，電荷移動相互作用エネルギー(CT+mix)，分散エネルギー(DI)の4つのエネルギー成分に分割し，その詳細を評価できる手法であるPIEDA (Pair Interaction Energy Decomposition Analysis) [12]をMIZUHO/BioStation [11]に実装し，これを3つのタンパク質-リガンド系に適用した．用いた系は，ノイラミニダーゼ-オセルタミビル(商品名タミフル®)，EGFRチロシンキナーゼ-エルロチニブ(商品名タルセバ®)，エストロゲン受容体-17β-エストラジオール複合体である．PIEDAを用いることでタンパク質と医薬品化合物などの相互作用をより詳細に解析することが可能になり，新規医薬品の設計に有用な知見が得られると期待される．

2 手法

FMO法ではフラグメント間の相互作用エネルギーであるIFIEを計算することができ，例えば医薬品化合物とタンパク質の各アミノ酸残基との間の引力的または反発的相互作用の大きさを評価することができる．北浦，Fedorovら [12]により提唱されたPIEDA (Pair Interaction Energy Decomposition Analysis)はIFIEを各種エネルギー成分に分解し，注目する相互作用がどのようなエネルギー成分から成り立っているかを評価できる方法である．PIEDAの理論的背景や計算手法については論文 [12]にて詳細に解説されているが，本稿でも主要な式に関して再度概説する．

2.1 エネルギー分解法

PIEDAは北浦・諸熊のエネルギー分解法(Energy Decomposition Analysis; EDA) [13]の考え方をFMO法のIFIEに適用した方法である．例えば水分子の二量体があった場合，これらの水分子間の結合エネルギーΔEは二量体の全エネルギーから各々の水分子の全エネルギーを引くことで求められる(Supermolecule方式)．その結果得られる水分子の結合エネルギーは数kcal/mol程度となるが，その内訳を静電相互作用，交換反発，電荷移動相互作用等の各種成分に分ける方法がEDA法である [13]．例えば，酸素原子間の距離が2.98 Åの水分子二量体ではHF/4-31G*レベルでΔEが-7.72 kcal/molとなるが，その内訳として静電相互作用が-8.98 kcal/mol，交換反発が4.19 kcal/mol，電荷移動が-2.11 kcal/molなどと求まり，静電相互作用が主なエネルギー成分であることがわかる [13]．

2.2 PIEDAの各種エネルギー成分

PIEDAによりフラグメントI，J間のIFIE (ΔE_IJ)を静電相互作用エネルギーΔE^ES，交換反発エネルギー(ΔE^EX)，電荷移動相互作用エネルギー(mixed termを含む, ΔE^CT+mix)，分散相互作用エネルギー(ΔE^DI)の4成分に分けることができる(式(1))．   

$\Delta E_{IJ}^{}=\Delta E_{IJ}^{ES}+\Delta E_{IJ}^{EX}+\Delta E_{IJ}^{CT+mix}+\Delta E_{IJ}^{DI}$(1)

ここで静電相互作用エネルギーΔE^ESはフラグメント間の電子-核，電子-電子，核-核の引力的または反発的な相互作用の和として計算される(式(2))．   

$\begin{array}{c}\Delta E_{IJ}^{ES}=Tr\left(D^I{u}^{I,J}\right)+Tr\left(D^J{u}^{J,I}\right)\\ \text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}+\sum _{pq\in I}\sum _{rs\in J}{D}_{pq}^I{D}_{rs}^J\left(pq|rs\right)+\Delta E_{nucler}^{IJ}\end{array}$(2)

ここで式(2)の第1，2項は電子-核間の引力に対応し，フラグメントI (またはJ)の電子とフラグメントJ (またはI)の核との間の核引力積分(u^I,Jまたはu^J,I)にフラグメントI (またはJ)の密度行列Dを掛けて求める．次に式(2)第3項は電子-電子の反発に対応し，フラグメント間の4中心2電子反発積分(pq|rs)に各フラグメントの密度行列Dを掛けて求める．第4項はフラグメントIとJのダイマーの核反発エネルギーから各モノマーの核反発エネルギーを引いたものである．

次に式(1)第2項の交換反発エネルギーΔE^EXについて述べる．交換反発はPauliの排他原理による占有軌道間の反発的相互作用であり，通常のHartree- Fock法では2電子交換積分の項によって計算される．PIEDAにおいてΔE^EXはHartree-FockレベルのIFIEで式(3)の関係が成り立つことを利用して求める．   

$\Delta E_{IJ}^{HF}=\Delta E_{IJ}^{ES}+\Delta E_{IJ}^{EX}+\Delta E_{IJ}^{CT+mix}$(3)

IとJのダイマーの全エネルギー計算において，通常はSCF計算が行われ電子密度が自己無撞着になるまで反復計算が行われるが，ΔE^EXは電子密度を各モノマーの電子密度に固定した状態で得られる全エネルギーE^EXを用いて式(4)に示すように計算する．   

$\begin{array}{c}\Delta E_{IJ}^{EX}=E_{IJ}^{EX}-Tr\left(D_{IJ}^{EX}{V}_{IJ}\right)-E{\text{'}}_I-E{\text{'}}_J\\ -\Delta E_{IJ}^{ES}\end{array}$(4)

ここでE^EXはフラグメントのダイマー形成による電荷移動の寄与を無視した全エネルギーなので式(3)における第3項が0となる．また，式(4)の第1から4項はIFIEを計算するための式に対応し，先に求めたΔE^ESとの差からΔE^EXを計算できる．

ΔE^ESとΔE^EXが計算できれば式(3)からΔE^CT+mixを求めることができる．ここでΔE^CT+mixは全体から他のエネルギー成分を差し引いた余りとして求めるため，電荷移動とその他の相互作用がカップリングした高次の相互作用エネルギー(mixed term)を含むが，その値の大きさは一般的な系では無視できる程度の大きさであることが示されている [13]．

最後にΔE^DIはMP2計算などによって求めるダイマーの電子相関エネルギーから各モノマーの電子相関エネルギーの和を引いて求めることができる．

2.3 フラグメント間の電荷移動量の評価

論文 [12]ではフラグメントIからJへの電荷移動量Δq_I→Jについても評価する方法が提唱されている．注目するフラグメントIの原子電荷について，モノマーの場合のMulliken電荷q_I^Iとダイマー形成時のMulliken電荷q_I^IJの差分を取る方法であり，以下の式で求められる．   

$\begin{array}{c}\Delta q_{I\to J}^{CT}\equiv q_I^I-q_I^{IJ}\\ q_I^X\equiv -\sum _{\alpha \in I}{Z}_{\alpha }^X+Z_I\left(X=I\text{or}IJ\right)\\ Z_{\alpha }^X=\sum _{\mu ,\nu \in \alpha }{D}_{\mu \nu }^X{S}_{\mu \nu }^X+\sum _{\mu \in \alpha }\sum _{\nu \notin \alpha }{D}_{\mu \nu }^X{S}_{\mu \nu }^X\end{array}$(5)

ここでαは原子のインデックス，μとνは基底関数のインデックスであり，DとSはそれぞれ密度行列と重なり行列である．

このフラグメント間の電荷移動量を評価する方法についてもBioStationに実装し，Δq_I→JをPIEDAの各種エネルギー成分と合わせて出力するようにした．電荷移動の相互作用の大きさをΔE^CT+mix項によって評価することができ，さらにΔq_I→Jを評価することによって電荷移動の方向を知ることができる．

3 実装

PIEDAの計算機能をABINIT-MPに，PIEDAの計算結果を可視化する機能をBioStation Viewerにそれぞれ実装した [11]．なお，PIEDAの検証については，先行論文 [12]でEDAとの数値比較がなされており，水クラスターn量体 (n = 2−16)において両者の各エネルギー成分が1％程度の誤差で一致することが示されている．本実装では，先行論文で用いられているGAMESSプログラムとの数値比較を行い，一致することを確認した．ABINIT-MPでのPIEDAの実行は非常に簡単であり，従来のFMO計算のインプットファイルにPIEDAの実行を指示するネームリストとキーワードEnergyDecomposition = 'YES'を追加するだけで良い．FMO計算が完了するとログファイルにはIFIEの出力に続いてFigure 1に示すようにPIEDAで計算されたエネルギー成分であるES，EX，CT+mix，DI及びΔq_I→Jがフラグメントの組み合わせ毎に出力される．これは水分子5量体の例であるが水素結合を形成する水分子間の静電相互作用エネルギーが最大−16.7 kcal/molと大きいことがわかる．これに交換反発エネルギーなどが加わり正味の相互作用エネルギーは−8から−9 kcal/mol程度となる．また，各エネルギー成分の絶対値はES，EX，CT+ mix，DIの順に大きいこと，水素原子から酸素原子に0.04 e程度の電荷が移動していることなどもわかる．

Figure 1.

 (a) Output of PIEDA analysis for the system of five water molecules (b) calculated at FMO-MP2/6-31G** level.

GUIを用いた可視化解析は計算完了時に出力されるチェックポイントファイル(cpfファイル)をBioStation Viewerで読み込むことで実施できる．ここではエストロゲン受容体の複合体を例にGUIの基本的な機能を紹介する．従来からの機能としてBioStation Viewerでは注目するフラグメント(ここではリガンド)と周囲のフラグメント(ここではタンパク質の各アミノ酸)の間のIFIE値に基づいて立体構造を色づけすることができたが，今回新たに実装した機能では， PIEDAの各エネルギー成分(ES，EX，CT+mix，DI)及びΔqに基づいて立体構造を色づけることできるようになった(Figure 2 (a-e))．それぞれのカラーグラデーションは，ESが赤/青(−/+)，EXが紫(+)，CT+mixがシアン(−)，DIが緑(−)，Δqが赤/青(−/+)となっており，色の強弱が相互作用の強さを表している．また，Figure 2 (f)に示すように，Figure 2 (b)~(e)のうち最も大きなエネルギー成分によって各フラグメントを色づけする"Main component”表示機能も実装した．

Figure 2.

 Visualization of PIEDA results using BioStation Viewer: (a-d) four energy components of interaction energy between the estrogen receptor (ribbon) and its ligand (yellow, ball & stick), (e) charge-transfer estimated by change of Mulliken atomic charge,(f) "main component" representation, in which each amino acid colored by the largest energy component in the four components.

4 解析事例

3つのタンパク質－リガンド系に対してPIEDAを適用した．具体的には，(1)ノイラミニダーゼ(neuraminidase; NA)とオセルタミビル(商品名タミフル®)，(2) EGFRチロシンキナーゼとエルロチニブ(商品名タルセバ®)，(3)エストロゲン受容体(ER)と17β-エストラジオールについて，リガンドと周辺アミノ酸残基間の相互作用成分解析を行った．尚，計算レベルはFMO2-MP2/6-31G*を用い，フラグメント分割に関しては，タンパク質はアミノ酸残基単位，リガンドは全体を1フラグメントとして扱った．

4.1 インフルエンザノイラミニダーゼ(NA)

抗インフルエンザ薬であるオセルタミビルは，増殖したウイルスの遊離過程において，ウイルス膜表面タンパク質であるNA (Figure 3 (a)) が宿主細胞表面のシアロ糖鎖を切断する機能を阻害する．阻害剤とNAの周辺アミノ酸残基との間の相互作用をFMO法によって定量的に評価することで，阻害メカニズムの解明や阻害剤の設計などに役立てることができると期待される．

Figure 3.

 (a) A complex of influenza surface protein NA (pink and green) and oseltamivir (yellow). (b) Active conformation of oseltamivir.(c) Inter fragment interaction energy (IFIE) between oseltamivir and its surrounding amino acids. The four energy components are shown by red (ES), pink (EX), cyan (CT+mix), and green (DI).

本解析ではPIEDAを用いてNA –リガンド間の相互作用を各エネルギー成分に分割し，相互作用の詳細を解析した．Figure 3 (c)にNA –リガンド間のIFIEを示す．この図において横軸はアミノ酸フラグメントであり，縦軸はそのアミノ酸残基とリガンドとの間のIFIEの各エネルギー成分の大きさ(kcal/mol)を表している．また，横軸のアミノ酸はリガンド－アミノ酸間のIFIEの総和が引力的に大きいアミノ酸を左から並べ，同じくIFIEの総和が反発的に大きいアミノ酸を右から並べ，中間のIFIEが0 kcal/mol付近のアミノ酸は省略している．この棒グラフから阻害剤とNAのアミノ酸残基間の相互作用はES成分が支配的であることがわかる．これはFigure 3 (b)に示すようにオセルタミビル活性体構造が両性イオンのリガンドであることから妥当な結果であるといえる．

Figure 4には従来のIFIE (Total IFIE)，各エネルギー成分(ES，EX，CT，DI)，およびmain component表示により各アミノ酸フラグメントを色づけした図を示す．ESによる相互作用が支配的であるため，Total IFIEとESの表示にほとんど違いがなく，Main component表示でもESによる引力的相互作用が支配的であることが視覚的にも確認できる．

Figure 4.

 (a) Conventional IFIE (Total IFIE), (b)-(e) the four energy components based on PIEDA (ES, EX, CT, DI), and (f) main component of the interaction energy between oseltamivir and its surrounding amino acids.

4.2 EGFRチロシンキナーゼ

EGFRチロシンキナーゼ(EGFRK, Figure 5 (a))は肺がんの分子標的薬のターゲットあり，既にゲフィチニブ(商品名イレッサ®)，エルロチニブ (Figure 5 (b)) などの阻害剤が上市されている．エルロチニブはπ電子を多く持つリガンドであり，Figure 5 (c)(d)に示すようにEGFRKとの結合状態では周囲のアミノ酸と多くのCH/π相互作用 [14,15,16]を示す．

Figure 5.

 (a) A complex of EGFRK (ribbon) and erlotinib (ball & stick). (b) Structural formula of erlotinib. (c),(d) CH/π interaction between erlotinib and surrounding amino acids by using CHPI program [14,15,16]. (e) IFIE between erlotinib and its surrounding amino acids.

Figure 5 (e)にEGFRK－リガンド間のIFIEの各エネルギー成分を示す．Figure 3 (c)のNA‐リガンドでは静電項(ES)が殆どであったが，ここでは分散項(DI)の寄与も大きいことがわかる．リガンドとのCH/π相互作用が示唆されているアミノ酸にFigure 5 (e)中に矢印を付けて示しているが，分散エネルギーが大きいアミノ酸と一致しており，CH/π相互作用の存在を裏付けている．

Figure 6に示した立体構造のエネルギー成分による色づけを見ると，Total IFIEで引力的な相互作用を示している残基でもES成分では反発的なものがあることがわかる．これらを安定化させているのはDI成分であり，main component表示と併せて見ることで各成分の寄与を明らかにすることができる．各残基の主要成分は主に，荷電・極性アミノ酸残基がES，疎水性残基がDIとなっている．

Figure 6.

 (a) Conventional IFIE (Total IFIE), (b)-(e) the four energy components based on PIEDA (ES, EX, CT, DI), and (f) main component of the interaction energy between erlotinib and its surrounding amino acids.

4.3 ER-17βエストラジオール

核内受容体の一種であるEstrogen Receptor (ER)は，女性ホルモンを内因性リガンドとしており，乳がんや骨粗しょう症のターゲットとなっている．これまでに，エストロゲン受容体への複数の化学物質の結合性をFMO計算を用いて評価できることが示されている [17,18]．

本解析では，PIEDAを用いてERと17β-エストラジオール(Figure 7 (a))との間の相互作用エネルギー成分を解析した．Figure 7 (b)にERとリガンド間のIFIEを示す．Total IFIEが引力的に大きかった上位2つのアミノ酸残基はES成分の寄与が大きく，それに続く残基ではDI成分の占める割合が大きいという結果になった．静電相互作用の寄与が特に大きかった2つの残基はGLU353とHIS524であり，Figure 7 (a)に示す通りリガンドとの間に水素結合を形成している．特にGLU353においてはCT成分も大きいことが判った．一方，その他の残基は大半が疎水性残基であった．Figure 8 (f)に示すmain component表示からはアミノ酸残基毎にESまたはDIが主要成分になっていることが示され，各残基の相互作用の特徴を視覚的に捉えることができる．

Figure 7.

 (a) Schematic picture of interactions between estrogen receptor (ER) and 17β-estradiol (EST). Direction and magnitude of charge-transfer (CT) are indicated by arrows. (b) IFIE between EST and its surrounding amino acids

Figure 8.

 (a) Conventional IFIE (Total IFIE), (b)-(e) the four energy components based on PIEDA (ES, EX, CT, DI), and (f) main component of the interaction energy between EST and its surrounding amino acids.

5 おわりに

本稿では，BioStationに新たに実装したPIEDAの概要と適用事例について報告した．PIEDAを用いることにより，これまで単一の値として得られていたIFIEを各種エネルギー成分に分けて計算することができるようになった．また，PIEDAで得られる各種エネルギー成分の可視化解析機能についても，main component表示などの新規の機能も含めて実装した．実証計算では，インフルエンザNA，EGFRチロシンキナーゼ，エストロゲン受容体の3つの系について，リガンドとの相互作用のエネルギー成分を解析し，それぞれの相互作用の特徴を明らかにした．PIEDA解析がリガンド結合を評価する上で有用であることが示され，創薬などの分子設計への活用が期待される．また，最近では複数のリガンドの結合性を統計的手法によって予測するための技術開発が行われており，FMO法との融合が模索されている [19,20,21,22,23]．PIEDAは各種エネルギー成分と電荷移動量という複数の説明変数を与えることができるという点で期待できる．

References

• [1]   K. Kitaura, T. Sawai, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett., 312, 319 (1999).
• [2]   K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett., 313, 701 (1999).
• [3]   D. G. Fedorov, K. Kitaura, CRC Press, Boca Raton, FL, 2009.
• [4]   D. G. Fedorov, T. Nagata, K. Kitaura, Phys. Chem. Chem. Phys., 14, 7562 (2012).
• [5]   S. Tanaka, Y. Mochizuki, Y. Komeiji, Y. Okiyama, K. Fukuzawa, Phys. Chem. Chem. Phys., 16, 10310 (2014).
• [6]   O. Ichihara, J. Barker, R. J. Law, M. Whittaker, Mol. Inf., 30, 298 (2011).
• [7]   T. Nakano, T. Kaminuma, T. Sato, Y. Akiyama, M. Uebayasi, K. Kitaura, Chem. Phys. Lett., 318, 614 (2000).
• [8]   Y. Mochizuki, K. Yamashita, T. Murase, T. Nakano, K. Fukuzawa, K. Takematsu, H. Watanabe, S. Tanaka, Chem. Phys. Lett., 457, 396 (2008).
• [9]   Y. Mochizuki, K. Yamashita, K. Fukuzawa, K. Takematsu, H. Watanabe, N. Taguchi, Y. Okiyama, M. Tsuboi, T. Nakano, S. Tanaka, Chem. Phys. Lett., 493, 346 (2010).
• [10]   ABINIT-MP and BioStation Viewer can be downloaded from the following web site. http://www.ciss.iis.u-tokyo.ac.jp/riss/dl/download/index.php#download_2．
• [11]  [ ABINIT-MP 5.0 and BioStation Viewer 13.01 (MIZUHO Rev.), T. Nakano, Y. Mochizuki, et al., Research and Development of Innovative Simulation Software (2011).
• [12]   D. G. Fedorov, K. Kitaura, J. Comput. Chem., 28, 222 (2007).
• [13]   K. Kitaura, K. Morokuma, Int. J. Quantum Chem., 10, 325 (1976).
• [14]   Y. Umezawa, M. Nishio, Bioorg. Med. Chem., 6, 493 (1998).
• [15]   M. Nishio, M., Hirota, Y. Umezawa, "The CH/π Interaction Evidence, Nature, and Consequences",1998, Wiley-VCH, New York
• [16]   M. Nishio, "New edition: Introduction to Intermolecular Forces in Organic Chemistry", Kodansha (2008)「」
• [17]   K. Fukuzawa, Y. Mochizuki, S. Tanaka, K. Kitaura, T. Nakano, J. Phys. Chem. B, 110, 16102 (2006).
• [18]   K. Fukuzawa, K. Kitaura, M. Uebayasi, K. Nakata, T. Kaminuma, T. Nakano, J. Comput. Chem., 26, 1 (2005).
• [19]   K. Mori, CICSJ Bull., 31, 108 (2013).
• [20]   K. Ohno, K. Mori, M. Orita, M. Takeuchi, Curr. Med. Chem., 18, 220 (2011).
• [21]   Y. Alexeev, M. P. Mazanetz, O. Ichihara, D. G. Fedorov, Curr. Top. Med. Chem., 12, 2013 (2012).
• [22]   S. Hitaoka, H. Matoba, M. Harada, T. Yoshida, D. Tsuji, T. Hirokawa, K. Itoh, H. Chuman, J. Chem. Inf. Model., 51, 2706 (2011).
• [23]   T. Yoshida, M. Shimizu, M. Harada, S. Hitaoka, H. Chuman, Bioorg. Med. Chem. Lett., 22, 124 (2012).