薬物開発の過程において、設計されたリガンドが標的となる生体高分子(以下受容体と呼ぶ)に特異的に結合することは活性発現の必須条件である。受容体の3次元構造情報が手に入る場合には、それに基づいて計算機上でリガンドを設計するという手法が非常に有効である。こうした情報の利用により新規リガンドの発見やその最適化過程において構造的示唆が得られるだけでなく、受容体に対する結合定数K_aを計算予測することでリガンド候補選択の労力が大幅に節約できる。K_aは結合の自由エネルギーG_bindと式(1)の関係にあるため、K_aの予測はG_bindの予測と等価である。

　現在までのところG_bindの予測を目的とする手法は数多く発表されているが、そのいずれも一般性・簡便性・精度などの点で問題があり、実際のリガンド設計に応用されるに至っていない。よって私はこれらの問題点を解消したG_bind予測法の開発を目的として研究を行った。

Fig.1 リガンド・受容体結合過程

【G_bindの定式化】

　リガンド・受容体の結合の過程をFig.1のように分解して考える。まず分子が接触面で脱溶媒和(G_desolvES)を受ける。続く会合の過程でリガンド重心の並進と回転の自由度及び分子内自由度が消失(G_tr-rot,G_conf)し、その一方で分子間立体反発(G_rep)・疎水性相互作用(G_Hphob)・分子間静電相互作用(G_interES)が形成される。そこでG_bindを(2)式のように、各自由エネルギー項にc_dE,c_tr,c_cf,etc.の係数をかけた和で表現することにした。

　式中の各エネルギー項の中身についても独自の定式化を行った。

(1)脱溶媒和エネルギー静電項G_desolvES

　分子表面の一部での脱溶媒和エネルギー静電項G_desolvESを直接測定する方法はなく、また計算に用いる理論式も一般的で簡便な形のものは確立していない。そこで本研究では部分表面についてのG_desolvESに対し独自の定式化を行い、実測可能な分子全体の溶媒和自由エネルギーG_solvのデータを利用して式中のパラメターを決定することにした。

　G_desolvESは、溶質分子表面の誘起表面電荷が分子同士の接触によって取り除かれる時のエネルギーと考える。これを計算するためまず分子表面を微小表面要素Sに分割し、要素i上での複合体形成により取り除かれる誘起表面電荷密度を今回独自に考案した(3)式で計算する。

Fig.2 実験値G_solvと計算値G_solvESの比較

　

　ここでex(i)は要素iが属する原子の複合体形成時の露出率で、疎水コアに近いパッキングをうけている時は0、溶媒に露出している時は1に近い値をとる。₀,_wは各々真空中・溶媒中での誘電率、は溶質分子の原子電荷由来の電場、は表面要素に垂直な単位ベクトルである。

　G_desolvESはこの誘起表面電荷と溶質分子の原子電荷の相互作用から(4)式により計算する。

　

　(3)式で分子全体についてを計算し(4)式の符号を逆転させれば、溶質分子全体の溶媒和自由エネルギー静電項G_solvESを計算できる。溶媒和自由エネルギーG_solvのうち最も大きな寄与はこの静電項G_solvESであると考えられるので、G_solv≒G_solvESとしてよいであろう。この計算法は3個のパラメターを必要とするが、これらは計算値G_solvESが実験値G_solvを再現するように決定した。Fig.2にパラメター決定に用いた120個の低分子化合物について実験値G_solvと計算値G_solvESを比較した結果を示す。

　受容体側のG_desolvESは最終的な結果に影響しないことが確かめられたので、この項はリガンド側についてのみ計算する。

(2)リガンドの並進・回転自由度消失G_tr-rot

　リガンド重心周りの自由度消失の効果はリガンドによらないので定数とした。

　

(3)リガンド分子内の自由度消失G_conf

　リガンド分子内のsp3型bond数に比例すると考える。メチル基などの末端のbondや露出率が大きいbond、環に含まれるbondは自由度消失の効果が小さいので数えない。

　ただし7員環以上の環についてはflexibilityを考慮に入れた。

　

(4)リガンド・受容体間の立体反発G_rep

　vdW potential(AMBER力場を使用)のうち12乗反発項を用いるが、induced fitを考慮に入れるため受容体側の原子の露出率でスケーリングする。

　

(5)疎水性相互作用G_Hphob

　これまでの方法では分子間接触表面積に単純に比例する形の定式化が一般的であった。しかし本研究では分子間接触表面積のうち相手分子のvdW表面と密着している部分をhydrophobic contact surfaceと定義して、これに比例する形での定式化を行った。

　

(6)分子間静電相互作用G_interES

　形式上は古典的なクーロン相互作用を用いる。誘電率については実験的に求められた実効誘電率関数を参考にしつつ、露出率に応じて変化する関数形を独自に考案した。

　

【G_bind計算式中のパラメター決定】

　このG_bind計算法には(2)式中の6個の係数、(9)式中のc_EPS及び露出率計算の半径R_envの合計8個のパラメターが必要である。これらは全てG_bind計算値が実験値を再現するように重回帰計算を用いて決定した。決定はArabinose binding protein9複合体、Avidin28複合体、Sialidase 17複合体、Trypsin40複合体、Dihydrofolate reductase37複合体の合計131複合体について実行した。パラメター決定の結果をFig.3に示す。また比較のため、分子力場のエネルギー項(静電項・vdW項)のみで重回帰を行った結果をFig.4に示す。図中の破線は誤差±2.0kcal/molを示す。今回開発したG_bind計算法は、従来の分子力場のエネルギー項のみによる計算法に比べ実験値からのr.m.s.誤差を大幅に改善した。

Fig.3本研究での結果

Fig.4従来の分子力場

　上でパラメターを決定した(2)式を用いて、パラメター決定に用いていないThymidylate synthase24複合体、HIV-1protease42複合体、FKBP35複合体についてG_bindを計算し、本方法を検証した(Fig.5)。その結果、本方法はこれらの系についてもパラメター決定に用いた系と同じ程度のr.m.s.誤差(〜1.5kcal/mol)でG_bindを予測できた。結合定数K_a自体が測定法の違いなどで1桁程度(〜1.4kcal/mol)の誤差を生じうることを考えれば、充分満足すべき結果といえる。同様に分子力場計算についても上で決定したパラメターを用いてこれら3つの系で計算を行ったところ、FKBPについてはr.m.s.誤差1.2kcal/molという良い結果を与えたものの、Thymidylate synthase,HIV-1 proteaseでは各々3.0kcal/mol,3.3kcal/molという大きなr.m.s.誤差を生じた。このことからも本方法の優位性は明らかである。

Fig.5 検証例左からThymidylate synthase,HIV-1 protease,FKBP

　多様な性質の受容体および多様な骨格を持つリガンドに対して、広い範囲で精度良くG_bind、従って結合定数K_aを予測できる計算法を開発した。本方法は計算コストの点でも負担が少なく簡便である。この方法を予測的に用いることにより、創薬におけるリード創製・改良などの過程を大幅に効率化できると期待される。