第2章 小分子
分子モデリング研究は,空間における原子の相対位置を直交座標で表し,コンピュータ内に分子のモデルを作り出すことから始まる。出発構造はその後の研究の成否を決定づけるので,合理的で信頼できるものを用いることが不可欠である。分子の3D座標は様々な方法で得ることができるが,基本的な方法は次の3種類である。
まず最初に,X線データを利用して分子の3D構造を組み立てる方法に焦点を合わせよう。小分子の結晶学的情報に関する最も重要なデータベースはCambridge結晶学データベースである[1]。このデータベースには,大きさが約500原子までの有機化合物と無機化合物を対象に,実験的に決定された原子座標が多数収録されており,その内容は絶えず更新されている。Cambridge結晶学データセンターは,データベースとそれを検索し結果を解析するためのソフトウェアをリースしている。データベースの検索は,問題の分子に関する3D構造情報を含む簡潔で読みやすいファイルを出力する。市販されている分子モデリング・パッケージのほとんどは,このデータファイルを読取ることができる[2,3]。
データベースに記録された原子座標は,ファイルをモデリング・プログラムへ読込む際に自動的に直交座標へ変換される。構造は次に分子グラフィックスを利用して3D表示され,研究に供される。
一般に,小分子のX線結晶構造はきわめて高い分解能で解析できるが,データの精度は保証の限りではない。水素原子はX線回折では観測が困難なため,その位置決めは常に問題を伴う。X線回折法は,X線が原子の回りの電子雲によって散乱される原理を利用している。水素原子は電子を1個しか持たないため,X線の散乱に及ぼすその影響は小さい。そのため,水素原子は構造決定の際に通常無視される。もちろん水素原子の位置は,標準的な結合長と結合角の値を利用すれば指定することができる。しかしこの手順に従った場合,水素原子が関係する結合の長さは,いずれもあまり精確なものにはならない。X線データベースからの情報を利用する場合には,原子座標,結合長および結合角に内部矛盾がないことを事前に確かめた方がよい。X線構造データを取り扱う際には,特に次の点に留意すべきである:
キラルな分子の場合には,さらに次の点も確認したい:
分子データファイルには,これらの項目がチェックされたX線構造データのみを保存すべきである。ファイルの編成,拡張名,書式および情報内容はプログラムに依存している。
結晶状態の分子構造は,結晶充填力の影響下にある。そのため,結合長と結合角は理論的な標準値と異なることがある。固体における配座は,柔軟な分子がとりうる多くの低エネルギー配座のうちの一つに対応している。それは,結晶単位胞内の隣接分子の影響を常に受けており,場合によっては,結晶内に入り込んだ溶媒分子の影響を受けることもある。生物学的に重要な立体配座を明らかにしたければ,エネルギー的に許容される他の配座も配座解析により調べる必要がある。大域的なエネルギー極小構造(global energy minimum structure)と呼ばれる最も安定な立体配座に関する知識は,エネルギー含量の高い配座体の存在確率を評価する際に有用である。この問題については, 2.3節で詳しく取り上げることにする。
分子を組み立てるきわめて一般的な第二の方法は,既に存在する部分構造ライブラリー(fragment library)を活用する方法である。この方法は,結晶学データベースにアクセスできないときや,目的の分子に対するX線結晶データが手に入らないときによく用いられる。最近の市販分子モデリング・プログラムは,そのほとんどが部分構造ライブラリーを利用した分子組立て機能を備えている。
部分構造ライブラリーは,分子の3D構造を簡単に組み立てるための電子的なキットとして利用することができる。部分構造プールの全エントリーは,最適化の済んだ標準的な構造で登録されているため,このライブラリーを利用して組み立てられた3D構造は,最初から十分許容できる妥当な構造になっている。そのため変更を必要とするのは,ほとんどの場合ねじれ角だけである。この変更は原子の重なり,すなわち非結合原子対の接近によるvan der Waals接触を避けるために必要である。縮合環系分子では,飽和環の連結様式に様々な可能性が考えられる。この問題を解決するためには,よく似た構造をもつ環系のX線データや実験データをできる限り参考にし,正しい環連結様式が選択されるよう心掛けなければならない。
分子を構成する各原子は,いわゆる原子タイプ(atom type)で定義される特有の属性を備えている。分子モデリングでは,原子はたとえば混成状態や体積などの違いにより区別される。原子タイプとは,これらのパラメータを考慮した原子の分類である。力場の原子論的部分はこれらの原子パラメータで作られている。ライブラリーから選択されたフラグメントの原子タイプはもちろん既に定義されており,一般にそれは正しい。しかし多くの場合,どの原子タイプに属するかについて判断を下すことはそれほど容易ではない。N‐アセチルピペリジンを例にとってこの問題を考えてみよう。
もし部分構造ライブラリーのピペリジン環とアセチル基を用いてN‐アセチルピペリジンを組み立てたならば,ピペリジンの窒素原子は四面体構造のsp3混成タイプになる。しかし,この窒素原子はアセチル基に結合しているので,平面型三角形のsp2混成をとるアミド窒素原子と見なすこともできる。このような場合,正しい判断を下すためには,組立てルーチンから得られる結果をX線データと比較するか,問題の構造要素に対して量子力学的な計算を行ない,信頼性のある構造を求めるかしなければならない。図1は,半経験的計算とab initio計算の結果を,力場計算による構造およびN‐アセチルピペリジン‐2‐カルボン酸の結晶構造[4]と比較したものである。
力場構造のsp3窒素原子は四面体構造を保持しているが,結晶構造や量子化学的に計算された構造では,窒素原子はほとんど平面型の構造をとっている。以後の計算における誤差を避けようとすれば,窒素原子には原子タイプとして平面型構造を割り当てなければならない。
置換飽和環系分子の組立てに際し,もう一つ問題になるのは,環が正しい配座になっているかどうかという点である。環の配座は置換基の影響を受けて変化するからである。シクロヘキサンは,有機化学で最も詳細に研究されている環状分子の一つである。可能な様々な配座とそれらを隔てるエネルギー障壁は,これまで多くの研究で取り上げられている[5,6]。この分子の最も安定な配座がいす形であることは疑う余地がない。このことは,モノ置換シクロヘキサンにおいても当てはまる。置換基の位置はエクアトリアル位が優位である。エクアトリアル位とアキシアル位の間のエネルギー差は,小さな置換基ではあまり明確ではないが,置換基が大きくなるにつれ増大していく[7]。部分構造ライブラリーを利用した分子の組立て作業は,常に実験データと照合しながら進めることが必要であり得策である。
分子の3D構造を生成する第三の方法は,いわゆるスケッチ・アプローチ(sketch approach)である。このアプローチでは,まずマウスポインタを鉛筆代わりに使用して,コンピュータ画面上に分子の2D構造式を描く。数はきわめて限られるが,小分子の標準的な部分構造データがライブラリーに入っている場合には,出発点としてそれを利用することもできる。作画が完了すると,画面上の2D構造式は3D構造へ変換される。うまく変換するためには,原子タイプを正しく指定する必要がある。スケッチ・アプローチはきわめて簡単な方法であり,得られた構造は一般にあまり満足できるものではない。そのため,スケッチ操作の各段階の終わりには,おおざっぱな構造最適化が自動的に行なわれ,分子構造の緩和が図られる。
前節で説明した手順に従って作製された分子構造は,構造の最適化を行なうことによって,エネルギー極小状態にすることができる。この最適化は通常分子力学法により行なわれる。「分子力学」という表現は,分子の構造とエネルギーを推定する手段として広く受け入れられている計算法を定義するのに用いられる。
分子力学では,量子力学的アプローチとは異なり,原子を構成する電子と原子核は計算の中に顕在的な形では含まれない。分子力学は,原子から構成される分子を,調和力によって相互に作用し合う質量の集合体として捕える。このような単純化により,分子力学は,小分子に止まらず,さらに大きな分子やオリゴ分子系に対しても適用できる比較的高速な計算法となった。
分子力学法の枠組では,分子内の原子は様々な長さのばね(結合)によって連結された様々な大きさをもつゴム球(原子タイプ)として取り扱われる。構成原子全体のポテンシャルエネルギーは,Hookeの法則を適用し次式から計算される。
Etot=Estr+Ebend+Etors+Evdw+Eelec+... (1)
ここで,Etotは分子の全エネルギー,Estrは結合伸縮エネルギー項,Ebendは変角エネルギー項,Etorsはねじれエネルギー項,Evdwはvan der Waalsエネルギー項,Eelecは静電エネルギー項をそれぞれ表している。全エネルギーは,計算の過程で原子座標に関して極小化される。
分子力学は,基準となるひずみのない結合長,結合角,ねじれ角および非結合相互作用からの偏差として,分子の全立体エネルギーを計算する。ひずみのないこれらの値と力の定数(実際には,経験的に当てはめられたパラメータ)は,ひとまとめにして力場(force field)と呼ばれる。式(1)の最初の項は,ひずみのない理想的な距離からの結合の伸縮により生じるエネルギーの変化を表している。原子の間に働く力は調和的で,結合伸縮エネルギー項は,式(2)で与えられる簡単な二次関数によって記述できると仮定される。
Estr=1/2kb(b−b0)2 (2)
ここで,kbは結合伸縮の力の定数,b0はひずみのないときの結合長,bは実際の結合長である。
より精密な力場では,三次項[1]や四次関数[2-4],Morse関数[5]が用いられる。
変角エネルギーに対しても同様に,ほとんどの場合,伸縮エネルギーと同じ調和型のばねのモデルが適用される。変角項を記述する関数は式(3)のようになる。
Ebend=1/2kθ(θ−θ0)2 (3)
ここで,kθは変角の力の定数,θ0は結合角の平衡値,θは実際の結合角である。
ねじれポテンシャルエネルギー項は,一般に式(4)のようなコサイン関数で表される。
Etors=1/2kφ(1+ cos(nφ−φ0)) (4)
ここで,kφはねじれ障壁,φは実際のねじれ角,nは周期性(完全な1周期内に現れるエネルギー極小点の数),φ0は基準のねじれ角(値はエネルギー極大点が0°にあるコサイン関数では通常0°,エネルギー極小点が0°にあるサイン関数では180°)をそれぞれ表している。
直接結合していない原子間のvan der Waals相互作用は,通常次式のようなLennard-Jonesポテンシャル[6]によって表される。
(5)
ここで,Aijは反発項の係数,Bijは求引項の係数,rijは原子iとjの間の距離である。
式(5)はLennard-Jonesポテンシャルの一つの形であるが,力場によっては,その変形が使われることも多い。よく追加されるのは,静電力を記述するための関数項である。この項に対しては,一般に次のCoulomb相互作用式が使用される。
(6)
ここで,εは誘電率,Q1,Q2は相互作用している原子の電荷,rは原子間距離である。
電荷は2.4.1.1節で説明する方法により計算できるが,力場によっては,経験的に導かれたパラメータの一つとして用意されている場合もある[2-4]。
交差項,面外項,水素結合項などの寄与も考慮し,さらに細かく分化したポテンシャルエネルギー関数を使用した力場も幾つか提案されている。関数の形は力場によって異なるので,ここでその詳細をすべて説明することはできない。関心のある読者は,優れた総説[7,8]が出ているので,それを参照していただきたい。
結合には自然な長さと角度があるというのが分子力学の基本的な考え方である。力場では,結合長と結合角の平衡値および対応する力の定数が定義され,力場パラメータ(force field parameter)と呼ばれる。これらはポテンシャルエネルギー関数の中で使用されるが,基準となる平衡値からの構造のずれは,分子の全エネルギーの増加を引き起こす。分子力学における全エネルギーは,理想的な幾何学的構造をもつ仮想分子に比べたときの分子内ひずみの度合を表しており,それ自体に物理的意味はない。
汎用性のある良い力場に必要な条件は,様々なクラスのできる限り多くの分子を合理的な精度で記述できることである。分子力学の信頼性は,使用されているポテンシャルエネルギー関数の形とその関数に組み込まれたパラメータの良否に依存している。したがって,もし重要な構造要素に対するパラメータが欠けているならば,その力場を利用して質の高い計算を行なうことは明らかに不可能である。このような状況を避けるためには,個々の研究に対し,その都度適切な力場を選択する努力が必要であろう。
広い範囲の有機化合物と小分子に適用できる力場が幾つか開発されており[1-4,9],また,特にタンパク質や他の生体分子の計算に適した力場を組み入れた分子力学プログラムも報告されている[10-12]。もし特定の原子タイプに対するパラメータが欠けているならば,その値を適当に推定し,力場へ追加するようにしなければならない[13-15]。
既に説明したように,前節の方法で作り出された分子の3Dモデルは,ほとんどの場合,そのままではエネルギー的に理想的な状態ではない。理想的な構造にするためには,次に構造の最適化(geometry optimization)を行う必要がある。最適化の過程において,分子構造はより安定な平衡構造へと緩和される。図1と表1は,アンギオテンシン変換酵素阻害薬ラミプリラート(ramiprilate)の結晶構造に一般的なTripos力場を適用し,構造最適化を行なったときの結果を示している。結晶構造における内部ひずみは,理想的な結合長からの偏差の影響を特に強く受ける。したがって最適化の過程では,それに対応したエネルギー項(結合伸縮項,変角項)がとりわけ大きく変化する。全エネルギーの顕著な変化にもかかわらず,ねじれ角はあまり影響を受けていない。このことは,結晶中では,分子はほとんどすべて低エネルギー配座で存在するという周知の観測事実とよく符合している。結晶構造中における分子の配座は,決してエネルギー的に不都合なものではない。図1によると,結晶構造と最適化された力場構造はほんのわずかしか違わない。この事実は,またパラメトリゼーションが適切になされた力場であればどれを使用しても,大差のない構造が得られることを示唆するものと考えられる。
構造最適化の際には,不適当なvan der Waals接触は前もって取り除いておくべきである。最適化により得られる極小エネルギー構造は,出発構造に依存するからである[7]。
力場計算が構造最適化の最も重要な標準的方法として地位を確立できたのは,その速度と精度,そして小分子だけではなく大分子系に対しても幅広く適用できるといった長所によるものである。量子力学的方法は,複雑で多大な計算費用を必要とするため,後節で取り上げる特別な問題を除いて適用されることはない。
我々は,次に分子力学で一般に使用されるエネルギー極小化法の幾つかに焦点を合わせることにしよう。極小化アルゴリズムは,単にポテンシャルエネルギー面上の局所的な極小点(local minimum)を見つけ出すだけであり,大域的な極小点(global minimum)を発見できるわけではないことに注意していただきたい。
エネルギー極小化のための方法は,2つのクラスに大別される。第一のクラスは最大傾斜法,共役勾配法,Powell法のような一次微分法であり,第二のクラスはNewton-Raphson法で代表される二次微分法である。
最大傾斜法は,エネルギー極小点へ近づくために,数値的に計算されたエネルギー関数の一次微分を利用する。この方法では,まず出発構造に対しエネルギーが計算される。次に,構成原子の1つが座標系の軸方向に沿って順番にわずかずつ動かされ,その都度エネルギーが再計算される。このプロセスは全原子に対して繰り返され,最後に,全原子がエネルギー表面の下り坂方向の新しい位置まで動かされる[7]。操作は,あらかじめ定めておいた極小化判定条件を満たした時点で打ち切られる。最大傾斜法は,極小点付近では収れんが遅くなるため,極小点から遠く離れた構造に対して適用されることが多い。この方法は,十分精密化されていない結晶データをもとに低エネルギー構造を生成させたり,2D構造式から組み立てられた3D分子構造を緩和させるのに特に適している。最大傾斜法は,ほとんどの場合,共役勾配法のようなより高等なアルゴリズムを用いる前に行なわれる,予備的なおおざっぱな極小化計算に利用される。
共役勾配法は,ある繰返しサイクルから次のサイクルへ至る関数の情報をすべて蓄積していく。そのためこの方法では,前のサイクルでなされた進捗は,次のサイクルで逆戻りすることはない。極小化の各ステップで勾配が計算され,新しい方向ベクトルを計算するための追加情報として利用される。ベクトルの方向は,ステップを進めるごとに極小点へ向かって絶えず改善されていく。最大傾斜法に比べ,計算に多大な時間と記憶容量を必要とする。しかし大きな系を扱う場合には,この共役勾配法が選択すべき方法である。共役勾配法により達成できる極小点への効率的な収れんは,繰返しサイクル当りの費用と計算時間の増加を補って余りあるものである。
Powell法は共役勾配法と非常によく似ている。この方法は共役勾配法に比べ収れんが早く,適用できる問題範囲も広い。しかしこのアルゴリズムは,ねじれ角の劇的な変化を伴うことがままあるので,その利用に際しては注意が必要である。Powell法は,配座解析後に行なうエネルギー極小化の手段として使用することはできない。せっかく捜し当てた低エネルギー配座が望まない形へ変化してしまう可能性があるからである。このような目的には共役勾配法を用いるのが得策である。
二次微分法としてのNewton-Raphson法は,探索方向を見極めるのに,勾配に加えて関数の曲率を利用する。二次微分はまた,関数が極小となる位置を予測する目的で用いられることもある。Newton-Raphson法の効率は,収れんが近づくにつれ向上する。この方法の欠点は,大きな系の計算に多大な時間と記憶容量を必要とし,ひずみの大きな構造では,極小化のプロセスが不安定になりやすいことである。そのため,このアルゴリズムは,おおざっぱに最適化された構造をさらに精確なエネルギー極小構造へ速やかに収れんさせたいような場合にのみ利用される。最適化の方法についてさらに詳しく知りたい読者は, 引用文献の16と17を参照されるとよい。
以上を要約すると,エネルギー極小化法の選択は,2つの因子,すなわち系の大きさと最適化の現在の状況に依存すると言えよう。極小点から遠く離れた構造を考えたとき,繰返しの最初の10〜100サイクルにおいて最良の結果を与える方法は一般に最大傾斜法である。この極小化のプロセスは,そのあと共役勾配法かNewton-Raphson法を用いて収れんさせることにより完結する。系が大きすぎて二次微分行列を保存したり計算したりできない場合には,共役勾配法が使用できる唯一の方法になる。極小化の操作は,結果が収れんするまで続けられる。
エネルギー極小化の際の収束判定基準は幾つかの方法で定義される。最大傾斜法のような非勾配極小化法では,分子系の現在の構造の適否は,エネルギーや座標の増分を用いて判定するしかないが,勾配極小化法では,この目的に原子勾配を利用することができる。この点に関して最良の結果を与えるのは,分子の各構成原子上における力の二乗平均勾配を計算する方法である。望ましくない構造領域は,最大微分を吟味することにより検出可能である。すべての微分が与えられた値よりも小さければ,エネルギー極小構造の質について疑う余地はない。収束判定値は,たとえば最大微分の場合,極小化の目的に応じて変化する。もし,ひずんだ分子を単に緩和させるだけで良いのであれば,最大微分に対し0.1 kcal mol-1 Å-1のような緩い収束条件でも十分である。しかし,最終的な極小点を見つけ出すためには,0.001 kcal mol-1 Å-1より小さな判定基準が最大微分に対して要求される。
収束条件の選択は,合理的な精度でエネルギー極小構造を決定することと,進展がそれ以上望めないとき,不必要な計算を極力避けることの間の釣合いを考えてなされるべきである[17]。
分子力学の計算は,通常真空条件下(ε=1)で行なわれる。非極性炭化水素類の場合には,計算に溶媒を含めても気相計算と同じ結果しか得られない。しかし,電荷や双極子をもつ分子の研究では,溶媒効果を考慮した扱いが必要である[7]。そうしなければ,静電相互作用の影響を強く受ける立体配座が,過大に見積もられることになるからである。力場は求引的な静電相互作用を最大化しようとするため,エネルギー的にはきわめて望ましいが,非現実的な低エネルギー配座が得られてしまうのである。この問題は,対応する溶媒の誘電率(dielectric constant,ε)を計算に含めることにより解決される[18]。たとえば,水のεは80である。小分子の静電場は,高分子のそれとは対照的に均一であると見なしてよい。このことは,一定の誘電率を使用することが原理的に許されることを意味している。文献には,実験的に決定された溶媒分子の誘電率が多数報告されている。それらを利用すれば,溶媒和分子のCoulomb項を正しく取り扱うことができるはずである。
電荷と溶媒和の問題は,最適化の際に電荷を無視して計算を行なうことにより,きわめて簡単ではあるが効果的に処理することができる。この方法は,たいていの場合満足な結果を与える。この方法が特に役立つのは,配座解析の結果を最適化するような場合である。電荷の使用は静電相互作用によって立体配座を大きく変えてしまう可能性があるからである。ただし,水素結合現象を記述するような場合には,電荷を考慮した取扱いが常に必要である。
静電相互作用の強度はr-1に比例して減少する。そのため,ある種の力場では,高分子表面へのリガンド分子の接近に伴う溶媒分子の置き換えの効果をシミュレートするため,距離依存型の誘電率が使用される。配座解析が薬理作用団探索の一部を成しているような場合,このような処置は特に重要である。
理論計算による結果は,可能な限り実験データによって検証されなければならない。NMRデータはこの目的に対してとりわけ有用である。報告されているNMRデータのほとんどは,クロロホルムや類似の有機溶媒中で測定されている。理論値と実験値の間に存在する隔たりは,対応する誘電率を力場のCoulomb項に明示的に含めることによって改善されるはずである。
誘電率の考慮は,溶媒効果をシミュレートする一つの方策である。これに代わる方法として考えられるのは,ばらばらの不連続な溶媒分子を含んだ箱で分子の回りを囲む方法である。しかし,この方法は重大な欠点がある。それは,計算に時間がかかりすぎること,また計算に含めることができる溶媒分子の数に制約が存在することである。
量子力学的方法についても,ここで簡単に議論しておかなければならない。この方法は計算化学におけるきわめて有用な付加手段であるからである。一般に,分子の幾何学的構造や相対配座エネルギーのような性質は,パラメトリゼーションが適切になされた一般力場を用いれば,広い範囲のさまざまな構造に対し,高い精度で計算することができる。しかし,目的の構造に対し力場パラメータが入手できない場合もある。量子力学的方法は,このような場合に構造最適化の手段として利用される。また,遷移状態や反応経路の計算も,分極や異常な電子分布の影響を受けた分子構造の解析と併せて量子力学の守備範囲である。量子力学的方法が他の方法に比べて不利な点は,計算にかかる費用が莫大で,しかも小さな分子しか扱えないという点である。したがって量子力学的方法は,この方法でなければ処理できない特別な問題に対してのみ適用されるべきである。本節の目的は,理論的な視点から量子力学的方法を議論することではなく,半経験的方法やab initio法の計算プログラムを利用する際に役立つ実際的な指針を提供することにある。これらの方法の理論的側面については,単行本や総説が数多く刊行されているので,関心のある読者はそれらをご覧いただきたい[19-22]。
ab initio法は,分子力学や半経験的分子軌道法とは異なり,経験的パラメータを使わずに実験データを再現する。したがってab initio法は,実験的な情報がほとんどないか,全く利用できないような状況下において特に適した方法である。
ab initio計算の質は,計算に使用される基底関数系(basis set)[23,24]と計算方法に依存している。基底関数系の選択が悪いと,幾ら長い時間をかけて計算を行なっても無意味な結果しか得られない。どの基底系を使用するかの判断は,計算の目的や研究される分子の性質により異なってくる。大きな基底系を用いても,必ずしも常に実験データと一致する結果が得られるという保証はない[25]。このことは,よく心に留めて置いていただきたい。
ここでは,最もよく知られている基底関数系についてのみ論ずることにしよう。STO-3G(3個のGauss関数により近似されたSlater型軌道)は,昔広く使用された最小の基底関数系である。この最小基底系は,基底状態の全電子を収容するに必要な最小数の原子軌道から成り立っており,原子の球対称性を仮定している。
最近,特に人気が出ている基底関数系は,原子価殻2倍基底関数系(split-valence basis set)である。この基底系では,原子価殻の軌道は内側と外側の2組の基底関数によって表される。このような改良を加えることにより,電子の滞留がより柔軟に記述できるようになった[26]。原子価殻2倍基底関数系はSTO-3G基底系を発展させたものであり,ab initio計算では,3-21G,4-31Gおよび6-31Gといった基底系がよく利用される。これらの基底系は,内殻軌道と一次短縮原子価殻軌道の展開に用いられる原始Gauss関数の数において異なっている[25]。たとえば,4-31Gは内殻軌道が4個のGauss関数,原子価殻の内側と外側の軌道がそれぞれ3個と1個のGauss関数から作られていることを示している。
基底関数系の次のレベルの改良は,分極関数の導入によってなされる。すべての非水素原子はd軌道を追加されるが,これによってp軌道の中心は核位置からずれ,軌道にひずみ(分極)が生じる。この補正は,小さな環をもつ化合物では殊に重要である[26]。分極基底系は6-31G*のように星印を付けて表される。この基底系(6-31G*)は内殻軌道に対して6個,sおよびp原子価殻の内側/外側軌道に対して3/1個の原始Gauss関数をそれぞれ用い,さらにd関数の6個1組を追加した形になっている。
基底関数系についてさらに詳しく知りたい読者は,そのための便覧や総説が刊行されているので,それらをご覧いただきたい[22,25]。
適切な基底関数系を選択するための一般的規則は残念ながら存在しない。行なうべき計算の水準は,必要とされる精度と計算の対象となる分子の性質に依存して変わる。普通の大きさの簡単な分子の場合,その構造最適化は3-21G基底系を使用すればまず十分である。しかし問題によっては,その程度の計算水準では,納得のいく結果が得られないこともある。たとえば,分極効果や電子の非局在化,超共役効果による影響を受けている分子構造を解析するような場合には,d軌道を含んだ上述の6-31G*やさらに高級な基底系を使用しなければならない(図2)。
コンピュータ技術の急速な発展にもかかわらず,高水準のab initio計算はまだ常に実行できるわけではない。計算上の法外な要求を克服する一般的方法は,3-21G基底系を用いて構造パラメータを最適化したのち,6-31G*レベルで波動関数を計算する方法である。この計算手続きは,しばしば6-31G*//3-21Gと呼ばれる。
問題の分子的性質に対する計算結果の精度は,より高級な基底系を使用することにより常に改善されるわけではない。計算の適切な水準を見つけ出すためには,実験値との偏差を調べ,各基底系から得られる計算結果の特徴を理解する努力が必要とされる。
分子力学とab initio計算の間に存在する深く大きな溝は,半経験的分子軌道法によって埋められる。この方法は,その本質において量子力学的方法である。ab initio法との主な違いは,計算時間を短縮するために,積分の評価に経験的パラメータを導入している点である。一中心反発積分と共鳴積分は,実験値をできる限りうまく説明できるパラメータによって置き換えられる。
半経験的アプローチは,また,分子的性質のほとんどが主に価電子によって決定されるという基本的な考え方に立っている。したがって,計算では価電子のみが考慮され,計算時間の更なる短縮が図られる。
半経験的方法はいろいろ提案されているが,それらが拠り所としている理論的仮定はいずれも同じであり,異なるのは近似のスタイルだけである[27]。AM1[28]とPM3[29-31]は,結果の精度と計算費用の問題をきわめてうまく妥協させた申し分のない方法である。AM1やPM3による計算は,小基底関数系を使用したab initio計算と同程度の精度で実験値を説明できることが知られている。半経験的方法がab initio法よりも優れている点は,計算速度が数桁速く,しかも200原子までの系であれば適用できる点である。しかしその結果については,慎重な点検を怠ってはならない。ab initio計算で不適切な基底関数系を選択した場合と同じように,半経験的研究においても,適切なパラメータの欠如は意味のない結果をもたらすからである。広い範囲の分子に対する半経験的方法の適用可能性と様々な分子的性質の計算をテーマとした総説が幾つか公表されている[28-31]。半経験的方法は,一般に第三周期の元素に対しては誤った結果を与えることに注意していただきたい。
分子は剛体ではない。室温における分子の運動エネルギーは,構成原子をすべて永久に動かせるほど十分大きい。このことは,分子とその構成原子の絶対的位置が固定されたものではなく,一重結合で結合した置換基の相対的位置も時間の経過と共に変化することを意味する。すなわち,一重結合を幾つか含んだ化合物は,各瞬間において異なる様々な回転異性体(配座体)の混合物で存在し,この混合物の定量的および定性的組成は絶えず変化している。もちろん,多く見出されるのは低エネルギー配座体である。
ある配座から別の配座への転換は,一重結合の回りのねじれ角の変化と主に関係している。結合長と結合角はわずかしか変化しない。分子配座の変化は,分子のポテンシャルエネルギーと幾何学的形の間の関係を記述する多次元面上の運動と見なすことができる。すなわち,ポテンシャル面上の各点は,個々の配座のポテンシャルエネルギーを表し,分子の安定配座はこのエネルギー面上の極小点に対応する。配座の相対的分布はその統計的重みに依存するが,この重みはポテンシャルエネルギーだけでなくエントロピーの影響も受ける。そのため,ポテンシャルエネルギー面上の大域的極小点―ポテンシャルエネルギーが最低となる配座―は,必ずしも統計的重みが最大の構造と一致するわけではない(詳細については,文献1を参照されたい)。
分子の多重配座の周知の例は,エタンのねじれ形と重なり形,n-ブタンのアンチ形とゴーシュ形,シクロヘキサンの舟形といす形である。エタン分子におけるCsp3‐Csp3結合の回りの回転は正弦形曲線のポテンシャル関数によって記述できる(図1)。
60°,180°および300°にあるエネルギー極小点はねじれ形に対応し,120°,240°および360°にある極大点は重なり形に対応している。ポテンシャルエネルギー関数(ポテンシャルエネルギー面)の極大点に対応する構造は常態では存在しないので,物理的性質や化学的性質を研究する場合,エタンではねじれ形のみを考慮すれば十分である。しかしエタンよりも大きく柔軟な分子では,そのような簡単な状況はもはや当てはまらない。室温においてエネルギー的に許容される回転異性体が幾つか存在するからである。たとえば,n-ブタンは室温において約70%がアンチ形で存在し,30%がゴーシュ形で存在する[2]。したがって,この柔軟な直鎖脂肪族化合物の物理的性質を議論する場合には,アンチ形とゴーシュ形の両方の配座を考慮しなければならない。同じことは,いす形と舟形をとるシクロヘキサンのような環式化合物に対しても当てはまる。
薬物分子の生物活性は,低エネルギー配座群の中に隠れている唯一つの活性配座に依存すると考えられる[3]。様々な化合物に対するこのいわゆる生物活性配座の探索は,医薬化学に課せられた主要な任務の一つである。受容体タンパク質の活性部位における特異的な高分子環境へ結合できるのは生物活性配座だけである。我々は活性配座に関する情報に基づいて,特定の受容体系に対する新しい薬剤を設計することができる。周知のように,活性配座は必ずしも最低エネルギー配座と同じではない。しかし,それはもちろん溶液配座の母集団から除外されるようなエネルギーの高い配座であることもない(この側面の議論については,文献4を参照されたい)。この意味で,低エネルギー配座の同定は,分子の構造と生物活性の間にある関係を理解するための重要なステップである。
NMRのような実験的方法は,分子がとりうる配座のうち1ないし数種の配座に関する情報を提供できるにすぎない。分子の配座ポテンシャルの完全な把握は理論的方法によってのみ可能である。配座解析のための理論的方法はいろいろ開発されており,応用例も数多く報告されている[5-12]。配座解析の最も一般的な方法は,ポテンシャルエネルギー面上にある極小点を系統的にすべて調べる方法である。しかしこの方法では,回転できる結合の数が増えると,極小点の数は劇的に増加してしまい,それらの極小点をすべて調べようとすれば,大変な時間と労力が必要になる。
配座解析に必要な時間は,エネルギーの計算に用いられる方法にも直接依存している。配座エネルギーは量子力学的方法か分子力学的方法のいずれかを使用して計算することができるが,量子力学的計算は非常に時間がかかるため,大きな分子や柔軟な分子に対しては適用することができない。そのため,配座探索用プログラムのほとんどは,通常エネルギーの計算に分子力学的方法を採用している。本節では以下,系統的な探索法と共に,配座解析へのモンテカルロ法と分子動力学法の応用について論じることにしよう。
系統的な探索法
[6,7,13]は,いろいろある配座解析法の中でも恐らく最も自然な方法である。この方法は,分子内のねじれ角の各々を系統的に変化させ,可能なすべての配座を発生させる。もし角の増分が妥当な値であるならば,この手続きは分子の完全な配座空間像を描き出してくれるはずである。系統的な探索で通常採用される角の刻み幅は30°である。これは360°完全に1回転させたとき,12個の配座体が発生することを意味する。刻み幅は最適値の近傍では,配座の極小位置を正確に決定するため,さらに5°まで小さくする必要がある。発生する配座の数は,角の刻み幅だけではなく回転できる結合の数にも依存している。すなわち,配座の数は回転可能な結合の数をnとしたとき,次式に示すようにn乗の割合で増加する。
配座の数=(360/刻み幅)n
たとえば,回転可能な結合を6個もつ分子に対して30°の刻み幅で系統的な配座探索を行なった場合,発生する配座体の数は2,985,984個(126)にも達する。これは到底処理することのできない膨大なデータ量であり,数を減らす以外その解決策はない。
データ削減の第一段階として利用できる方法はvan der Waalsスクリーニング―衝突試験―である。この試験は,配座のポテンシャルエネルギーが正確に計算される以前に行なわれる。互いに直接結合していない原子間にvan der Waals体積の重なりが検出される配座は,このスクリーニングによりすべて排除される。配座の妥当性を判定する数学的基準として用いられるのは,非結合原子対のvan der Waals半径の和である。van der Waals球の硬さは,いわゆるvan der Waals係数の値を調整することにより変えることができる。この乗法定数は原子の相互貫入性を司っており,その値の減少は非結合原子間の接触を和らげ,有効配座の数を増加させることになる。
衝突試験をパスした配座体は,分子力学法によりそのポテンシャルエネルギーが計算される。一般に配座エネルギーの計算では,静電相互作用は無視され電荷は考慮されない。また,配座解析は真空条件下で行なわれる。このことに対する根拠は次の通りである(2.2.4節参照)。もし静電相互作用を配座解析に含めるならば,計算の全プロセスははるかに複雑なものになる。質の高い原子電荷は,単に結合の様式だけではなく,原子の個々の空間環境に敏感である。それゆえ,初配座に対して計算された原子電荷は,ねじれ角が変化すればそれに合わせて常に更新されなければならない。また,電荷間に存在する強い静電相互作用は溶媒によって緩和されるが,計算にはこの効果も組み入れる必要があろう。しかしそれは,小分子に対してさえ莫大な計算時間を要求する。そしてさらに注目すべきことに,系の複雑さが増しても,溶液中の分子の配座的挙動に関し,分子内相互作用が消失するという事実以外,他に特に深い理解が得られるわけではない。同じ結果は,電荷を無視し,真空条件下で解析を行なっても得られるのである。リガンドの分子内接触も,受容体や酵素の活性部位においてはあまり問題にならない。
衝突試験をパスしたすべての配座体に対して配座エネルギーを計算する際,エネルギー窓(energy window)を適用することにより,配座数をさらに減らすことができる。エネルギー窓適用の拠り所は,極小点付近のそれよりもはるかに高いエネルギーをもつ配座が母集団の中に見出される確率はほとんどないという事実である。このことは,そのような高エネルギー配座が分子の生物活性に重要な役割を果たしてはいないことを意味する。このエネルギー窓の値は,問題分子の大きさと適用力場の両者に依存し,5〜15 kcal mol-1の範囲で変動する[11-15]。
すべてのフィルターを通過し残ったものは,その分子に対してエネルギー的に存在可能な配座群だけになっているはずである。しかし多くの場合,その数はまだ多すぎる。配座の多くは,たとえば回転箇所が1箇所異なるだけであり,互いにきわめて強い関連を示す。明らかに,これらは高い類似性を示す配座として1つのファミリーにまとめることができる。分子の配座的性質を記述するに際して,このようにして作られた配座ファミリーの各々から最低エネルギーをもつ配座体だけを取り出し,詳細な検討はそれらに対して加えるようにしても一般性を失うことはない。配座をこのようなファミリーへ分類するための方法は幾つか開発されている[15-17]。これまでに提案されている分類法は,いずれもねじれ角を分類のパラメータとして使用する。異なるのは配座体を各ファミリーへ割り付ける手順だけである。系統的な配座探索の過程で蓄積された大量のデータは,もう一つの可能性として,クラスター解析や因子分析のような統計的手法を利用し解析することもできる。これらの方法の詳細な議論については文献18をご覧いただきたい。
次に系統的な配座解析のプロセスを示す具体例として,我々のグループが行なった2種のH2-抗ヒスタミン薬,チオチジン(tiotidine)とICI 127032(図2)の研究を紹介しよう[19]。この研究は,分子モデリングパッケージSYBYLに組み込まれたSEARCHモジュールを利用して行なわれた[16]。
解析で使用した回転角の刻み幅は15°である。ただし,シアノグアニジンのメチル基に対しては,その対称性により0〜120°の範囲を30°の刻み幅で回転させるに止めた。配座の理論数,3.98×107個はvan der Waalsスクリーニングを適用することにより4.6×106個まで減少した。これは最初の数の約10%である。またエネルギー窓(15 kcal mol-1)の適用は,配座数をさらに90%減少させた。この時点でコンピュータに蓄積されている配座数は453,393個であった。これは,まだまともに処理するには多すぎる数である。そのため,残った配座は次のステップでファミリーへ分類された。この操作は我々のグループが開発したIXGROSプログラム[17]を利用して行なわれた(IXGROSの完全なソースコードは付録1に掲載されている)。この操作は最終的に227のファミリーを生成し,各ファミリーはエネルギーが最小の配座によって代表された。3.98×107個から227個への減少はきわめて印象的ではある。しかし,この少ないはずの227個でさえまだ多すぎるのである。この段階で,227の配座体のうちどれが活性配座体であるのかを判定することはできなかった。しかし,このことが,そしてこのことだけが我々の関心事なのである。解を見つけ出すためには,剛体か少なくとも半剛体でしかも生物活性な同族体を含める必要があった。しかもそれらの同族体は,類似の機構で同じ受容体部位へ結合しなければならない。一般に,柔軟な分子の活性配座を見つけ出すためには,同じ系列に属し強力でより硬い構造をもつ化合物との比較が必要である。H2-拮抗薬の場合,この硬く強力な薬物を代表するのはICI 127032であった。この剛体マトリクスの低エネルギー配座を考察し,IXGROSを繰返し適用した結果,残ったのは8つのファミリーであった。これらの配座は,チオチジンの生物活性配座を決定するために利用されたが,その結果は上首尾であった(図3)。
以上の例にも示されたように,柔軟な同族体系列の配座探索では剛体分子を含めることがポイントである。高い生物活性を示す剛体誘導体は,系列の他の成員に対するマトリクス(母型)として利用される。Marshallとその協力者[7]は,不活性な剛体化合物も解析に含めることにより,この方法の拡張を試みた。その結果,配座空間はさらに局限され,探索に要する時間も数桁短縮された。Marshallらの方法は「活性類似体アプローチ」として普及している。
モンテカルロ法―ランダム探索法―は,配座空間を全く異なるやり方で探索する。この手法は統計的性質を備えており[20],探索の各段階で配座は無作為に修正され新しい配座になる。
モンテカルロ法は最適化構造から出発する。繰返しの各ステップで,ねじれ角[11]または直交座標[8,9]に新しい値が無作為に割り付けられる。新しく生成した配座は分子力学を利用して最適化され,この無作為な処理は繰り返し実行される。エネルギーを極小化された配座はこれまでのステップで生成した構造と比較され,新しい構造である場合のみ保存される。ランダム探索法は配座空間の全領域を潜在的に包含している。しかし,このことが実際に成立するのは,処理が十分長時間続けられた場合だけである。それはきわめて長い時間かもしれない。なぜならば,新しいユニークな配座を検出する確率は,既に発見されている配座体の数が増えるに従って著しく低下するからである。また,たとえ非常に長い時間かけて計算を行なったとしても,配座空間が完全に包含されたかどうかは定かではない。それゆえ,解析の完全性を調べる手段の確立は非常に重要なテーマである。この完全性は,異なる初期配座から出発する処理を幾つか平行して行なうことにより,能率良く調べることができるはずである。もし結果が同じかほぼ同じであれば,探索は完全であると見なしてもよい。完全性はまた個々の低エネルギー配座に対する回復速度に基づいて判定することもできる。確率過程は,そのような配座を何度も再生するに違いないからである。
ランダム探索法の主な利点は,原理的にどんな大きさの分子でもうまく処理できる点である。しかし実際問題として,高度に柔軟な分子はそれぞれの配座空間で占める容積が大きすぎるため,その結果はしばしば収れんしない。ランダム探索法は,その他,系統的な探索法で取り扱うのがむずかしい環式化合物の研究にも応用され,成果を上げている。ランダム探索法の有効性を実際の例で示そう。Saundersら[12]は,ランダム探索法を含め様々な手法を用いてシクロヘプタデカンの研究を行なった。各種の方法を組み合わせた結果得られたエネルギー極小配座の数は全部で262個であった。使用した方法は,単独ではどれも262個すべての配座体を発見することはできなかったが,ランダム探索法の一つはそのうち260個を検出することができた。このことからも,ランダム探索法は様々なタイプの分子の配座解析に非常に適した方法であると言えよう。しかし,配座空間を完全に包含するためには莫大な計算時間が必要なことも残念ながら事実である。
系統的な配座探索法は,柔軟な分子のポテンシャルエネルギー面上にある多数の極小点を検出する有用な手段である。一般に,この方法は許容されるすべての配座を発生させることができるので,完全な配座探索を達成できる確率は高い。しかし,この方法の適用性に関しては明らかに限界が存在している。複数極小問題を解くことができるのは,回転可能な結合の数が少ないかなり小さな分子の場合だけである。
2.3.1節でも既に取り上げたが,回転可能な結合が6個ある分子を考えてみよう。系統的な配座探索では,発生する配座体の数が多すぎるため,データ処理に容易ならぬ問題が発生する。そのため,たとえば回転子が15個あるアラキドン酸(図4)のような柔軟な分子の研究は事実上不可能である。データを削減する方法を幾つか適用したとしても,この分子に対する系統的な配座探索は500,000個の異なる配座を生成する。配座空間は,この段階ではまだ完全にサンプリングされていないにもかかわらず,計算処理はデータのオーバーフローが原因でプログラムにより自動的に打ち切られる。
しかしランダム探索法による配座解析も,この分子に対しては,計算に要する時間の関係で同様に不可能である。たとえば,アラキドン酸よりも束縛の大きなシクロヘプタデカン分子の計算でも,Micro-VaxUコンピュータ上で約94日を必要とする[12]。
これと関連して,飽和環式化合物や部分飽和環式化合物に対する系統的な配座解析は,別のかなりめんどうな問題を提示する。探索の過程で,新しい環配座を作り出すために,結合を切断しなければならないのである。この手続きの有効性と信頼性は,幾つかの総説のテーマになっている[12,14]。
これらの問題を克服するためのきわめて一般的な戦略としては,配座空間の探索に分子動力学シミュレーションを利用する方法が考えられる。このアプローチの目的は,分子の時間依存的な運動特性を再現することである。分子動力学は,その基礎を分子力学に置いている。分子内の原子は,(既に2.2.1節で説明したように)使用した力場の法則に従って互いに相互作用していると仮定され,一定の時間間隔で,Newtonの第2法則と呼ばれる次の古典運動方程式が解かれる。
Fi(t)=miai(t) (1)
ここで,Fiは時間tにおいて原子iに作用する力,miは原子iの質量,aiは時間tにおける原子iの加速度である。原子に作用する力はポテンシャルエネルギー関数のグラジエントを用いて計算され,原子の初速度は動力学計算の初めに無作為にあてがわれる。系の初期原子座標に基づいて,時間tにおける原子の新しい位置と速度がまず計算される。次に,原子はこれらの新しい位置へ移され,新しい配座が作り出される。このサイクルはあらかじめ設定されたステップ数だけ繰り返され,その結果生成したエネルギー的に可能な配座の集まりはアンサンブル(ensemble)と呼ばれる。
Newtonの運動方程式は,あらゆる分子動力学アプローチにおいて同じように適用される。アプローチにより異なるのは,使用される積分のアルゴリズムである。運動方程式の積分で特によく使用されるのは,Verlet積分器[21]やその変形であるBeemanアルゴリズム[22]と馬跳び法(leap-frog scheme)[23]である。本書の枠組では,分子動力学理論をこれ以上詳しく議論することはできない。この問題については,さらに詳しく論じた総説が幾つか出ているので,関心のある読者はそれらをお読みになるとよいであろう[24-27]。
配座解析に分子動力学を適用しようとするならば,読者はその前にこの方法のもつ特徴をよく飲み込む必要がある。分子動力学は,伝統的な構造最適化法と異なり,配座間に横たわるエネルギー障壁に打ち勝つことができる。したがって,ポテンシャルエネルギー面において現在の位置から最も近い点以外の局所的極小点も検出することが可能なはずである。しかし,エネルギー障壁が高い場合や分子の自由度が非常に大きな場合には,系に含まれる配座体のうち幾つかは恐らく発見されることはないであろう。配座空間の大きさを考えると,決められたシミュレーション時間内に完全な配座探索を達成することは至難の業である。
配座のサンプリング効率を高めるため,分子動力学で広く用いられているのは,シミュレーションを高温条件下で行なう作戦である[27]。分子は,高温では,配座間に大きなエネルギー障壁があってもそれに打ち勝つことができ,完全な配座探索が達成される確率が高くなることを利用するのである。どのようなシミュレーション温度とシミュレーション時間が最も適当であるかは,当然のことながら問題にしている分子の属性に強く依存している。
次に,配座の柔軟性とシミュレーション温度との関係を調べた最近の研究を紹介しよう。使用したデータと資料はF.S. Jorgensen(Copenhagen,デンマーク)のご好意により提供されたものである。取り上げられた分子は実験的によく研究されているシクロヘキサンである。分子は,様々な出発配座と温度で分子動力学シミュレーションにかけられた*。図5はその結果を示している。
図によると,出発配座として用いられたシクロヘキサンのねじれ舟形配座(T1=0)は,400 Kではもう一つのねじれ舟形との間で振動しているにとどまるが,600 Kになると,いす形配座の一つ(T1=300)へ移行するに十分な運動エネルギーを獲得する。さらに温度を1000 Kまで上げると,いす形とねじれ舟形の両配座(T1=300→60)が併存し,いす形‐いす形間の相互転換も観測されるようになる。また,800ピコ秒(ps)後に観測されるのは,もっぱら一方のいす形配座(T1=60)だけである。
第二の研究で取り上げられた分子は,3種のメチル置換シクロヘキサン(1,1‐ジメチルシクロヘキサン,1,1,3,3‐テトラメチルシクロヘキサン,1,1,4,4‐テトラメチルシクロヘキサン)である。分子はいろいろな温度で分子動力学シミュレーションにかけられ,各温度で観測されるいす形‐いす形相互転換と実験的に測定された環反転のエネルギー障壁[28]とが比較された。表1によると,分子動力学の結果は,実験的に測定された環反転障壁の相対的な大きさを見事に反映している。高温分子動力学のこの実例は,系が配座空間の特定の領域に閉じ込められないようにするためには,十分高いシミュレーション温度を設定する必要があることを示している。
配座空間の探索へ分子動力学を適用するに当たって一般に用いられる戦略は,一定の時間間隔で配座を選び出し,それらを一番近い極小点まで移動させる方法である。この方法は,環式化合物系も含め小分子に対する配座解析研究の幾つかで利用されている[14,29]。この方面の研究で特に印象的な例は,セスキテルペンラクトン,タルプシガルギン(tharpsigargin)のポリヒドロキシ類似体(図6)に対する配座解析である。これもまたF.S. Jorgensenの研究室で行なわれた研究である。
________________________________________
*Tripos社(St. Louis,米国)のSYBYL(第6.0.3版)を使用。エネルギー極小化:Tripos力場,PM3局在電荷,誘電率ε=20,収束条件0.005 kcal mol-1Å-1。MDシミュレーション:各種の温度で1000 ps,全エネルギーを保存し,ps当り1配座をサンプリング。
この研究では,ポリヒドロキシ誘導体は1200 Kで分子動力学シミュレーションにかけられた
*。調べたかったのは環系の配座的性質である。シミュレーションの間,七員環は幾つかの配座をとり,それらの間で環の相互転換が幾度となく観測された。このことは,配座空間が広範囲にわたり探索されたことを示すものと考えてよい。抽出された配座の各々は次にエネルギーを極小化され,もっぱら七員環の配座に関して互いに比較された。二乗平均(rms)値が0.1Å以下の配座はすべて同一であると見なされた。以上の操作から,5種の異なる低エネルギー配座が得られた。理論的に見出されたこれらの三環式環系配座の一つは,幸いにもタルプシガルギンのNMRデータ[30]と一致した。図7はその配座を示している。
しかし場合によっては,抽出された配座を極小化するだけでは,最終的な極小配座に到達するのに十分ではないこともある。高温動力学シミュレーションの目的は,配座間のエネルギー障壁を越えるに十分な運動エネルギーを分子に供給することにある。しかしこのシミュレーションは,しばしば簡単な極小化の操作では緩和されない,極端にひずんだ幾何学的構造を与えることがある。
もしこのような事態が発生したら,高温徐冷分子動力学を行なうとよい
[31]。このアプローチでは,高温シミュレーションで抽出された配座は次に最適化されたのち,分子の内部ひずみを取り除くために,たとえば300 Kのような低温で再度揺り動かされる。最後にもう一度最適化を行なうと,高温シミュレーションのあとに単に構造最適化しか行なわない場合に比べて,より低いエネルギーの配座が得られる。この高温徐冷分子動力学シミュレーションをさらに変形した方法にいわゆる模擬徐冷法がある[32]。この方法では,系はシミュレーション温度を下げながら一定の時間間隔で冷却される。温度が0 Kに近づいたとき,分子はエネルギーの最も低い極小配座になっている。徐冷の最終サイクルで得られた構造は保存され,高温での次のシミュレーションの出発点として利用される。一連の低エネルギー配座を得るために,サイクルは数回繰り返される。得られた構造は既に極小配座に近いはずであるから,構造をさらに極小化する操作は必ずしも必要ではない。この方法の応用は,幾つかの研究のテーマになっている。引用文献の33と34は,さらに詳しく知りたい読者にとって参考になるであろう。
分子動力学シミュレーションは,他の配座探索法でうまくいかない場合に特に役立つ,配座空間をサンプリングするためのきわめて有用な補助手段である。ユーザーは適切な方法の選択とシミュレーション条件の設定に注意を払い,配座探索の完全性と結果の妥当性を確実なものにしなければならない。各アプローチにはそれぞれ長所と短所があることに留意し,実験的に誘導されたデータが入手できる場合には,極力それを確認に役立てるよう心掛けるべきであろう。
________________________________________
*Tripos社(St. Louis,米国)のSYBYL(第6.0.3版)を使用。エネルギー極小化:Tripos力場,PM3局在電荷,誘電率ε=20,収束条件0.005 kcal mol-1Å-1。MDシミュレーション:1200 Kで1000 ps,全エネルギーを保存し,ps当り1配座をサンプリング。
薬物‐受容体複合体の形成における最初の段階は認識事象である。受容体は,近付いてくる分子が特異的で堅固な結合を形成するに必要な性質を備えているかどうかを識別しなければならない。この認識過程は,最終的な相互作用複合体を形成する以前のかなり遠い距離から始まる。各分子を取り巻く3D静電場は,この認識に特に重要な役割を演じている。分極率や疎水性といったその他の分子的性質もまた,相互作用する分子の距離が近付くにつれ徐々に作用しはじめる。様々な化学的プローブを使用した分子間相互作用エネルギーの系統的計算とサンプリングによって決定される分子場が,複雑性を異にするあらゆるレベルの分子間相互作用を理解する上できわめて価値のあるデータセットであることは容易にお分かりいただけるであろう。本節では,これらの分子的性質の計算と解析に使用される方法を紹介し,その評価を試みることにしよう。
分子静電ポテンシャル(MEP)の知識は,分子相互作用と化学反応を研究する者にとってきわめて重要である。分子が互いに近付いたとき,最初に感応し合うのは遠距離の静電力である。一般に,相互作用力は静電力,誘起力および分散力の3成分から成り立っている。第一の型の相互作用は,電荷や永久双極子モーメントをもつ極性分子間に見られる。また,第二の型は,極性分子が無極性分子と相互作用する場合に現れる。すなわち,極性分子の双極子が電場を作り出し,無極性分子の電子分布を変化させて双極子モーメントを誘発するのである。さらに,2つの分子が共に無極性の疎水性物質であった場合でも,一方の分子の電子分布の時間的なゆらぎは,隣接するもう一方の分子に一時的な分子双極子モーメントを誘発する。この第三の型の相互作用は分散力(dispersion force)と呼ばれる。分散力は弱く,しかも相互作用する分子間の距離が増加すると急速に減衰する(2.2.1節参照)。しかしこの力は,中性の無極性分子間では引力の主要部分を形作っている。分散力はまたLondon力とも呼ばれる。
静電相互作用は引力か斥力のどちらかである。接近する分子の電気的に陽性な部分は,電気的に陰性な領域を捜してそこへ結合しようとするのに対し,同じ種類の電荷を帯びた部分は互いに反発し合う。明らかに,この非共有相互作用は分子の荷電領域の間で特に強く働く。分子内の電荷や永久双極子モーメントは,分子の回りに3D静電場を発生する。そのため分子から適当に離れた距離では,極性分子のみならず中性分子の場合でさえ,有意な分子静電ポテンシャルが存在する。このポテンシャルは,分子を取り巻く空間の3D格子点に置かれた正の点電荷と分子の電子分布間にはたらく相互作用エネルギーを表している。分子静電ポテンシャルを決定するためには,分子の電子的性質に対する正確な取扱いが要求される。したがってここでは,分子荷電密度を計算するための方法をまず優先し,最初に取り上げることにしよう。
分子の電子的性質は,正に荷電した原子核の回りの電子分布によって定まる。電子分布に関する詳細な情報は,X線回折のような実験や量子力学的方法を利用した計算から得ることができる。しかし,これらの方法がもたらすものは,三次元空間全体に広がった荷電密度の確率分布である。相互作用エネルギーを計算するためには,原子の中心に布置された点電荷が通常必要である。これは,分子電子分布のきわめて単純化された描像である。確率分布からこのいわゆる原子局在電荷(点電荷)への変換は,電子密度を原子中心に凝縮させることによって達成される。このような経験的な局在電荷は,その定義に幾らか任意性を残している。分子の電子分布が個々の原子の中心に割り付けられるため,分子の属性が原子の性質に置き換えられてしまうからである。局在電荷は観測可能な量ではない。点電荷による方法は,それが分子の物理的性質や化学的性質の相関や予測に活用できる場合のみ,科学的に正しい適切なモデルでありうる。前にも述べたように,全分子間相互作用エネルギーへの静電項の寄与は非常に大きい。そのため,一般に利用されている分子力学プログラムのほとんどは,原子局在電荷に依存するエネルギー項が計算できるようになっている。そのようなプログラムを使用すれば,数百個の原子からなる高分子でも,静電エネルギーの迅速な計算が可能である。
原子局在電荷を計算するために様々な方法[1]が開発されているが,方法論的に区別すると,それらは次の2つのアプローチに分類することができよう。
トポロジカルな電荷
トポロジカルな方法は,主に原子のタイプによる電気陰性度の違いに基づいた方法である。
この方法では,結合に直接関与している原子へ原子電荷を合理的に配分するため,結合に関する実験的な構造情報と原子の電気陰性度が適当な規則に従って結びつけられる。分子の幾何学的形や配座に関する情報は考慮されない。計算に必要なのは原子の結合度行列だけである。この方法を最初に提案したのはDel Re[3]である。それは飽和分子を対象としたものであった。共役系へ拡張したのはBerthod-Pullman[4]である。これらの方法は,現在もなおモデリング用プログラムの幾つかに組み込まれ利用されている。さらに新しいアプローチとしてはGasteiger-Hueckel法がある。これは,原子電荷のσ成分を計算するGasteiger-Marsili法[2]と古いHueckel理論を組み合わせた方法であり[12],前の方法と比較してより現実に即した結果を与える。原子電荷のπ成分はHueckel理論により迅速かつ効率良く計算される。全電荷はσ成分とπ成分の和で表され,π系に属する原子の形式電荷はπ系全体にわたって非局在化していると仮定される。計算はHueckel電荷,Gasteiger電荷の順に行なわれる。トポロジカルな方法の大きな長所は計算が速く,しかもその結果が実験データとよく一致する点である。しかし,大きな難点も同時に存在する。それは,特定の分子群を使った確認作業が必要であり,それがなければ結果を信頼できない点である。確認の手続きは省略されることが多いが,そのような研究はもちろん説得力をもたず意味をなさない。
トポロジカルな方法は,電荷を計算する標準的手段として市販のソフトウェア・パッケージにしばしば組み込まれている。
量子化学的方法
原子局在電荷の他の計算法は,すべて波動関数の量子力学的計算に基礎を置いている。波動関数は半経験的方法かab initio法のいずれかを利用して得られ,その選択は波動関数に要求される精度と使用できる計算装置の性能に依存している。波動関数から荷電密度を求める方法はいろいろある。最も古くかつ最も広く利用されている方法は,Mullikenのポピュレーション解析[7]である。この方法は,多くの量子力学プログラム[13-15]に標準的な方法として組み込まれている。ポピュレーション解析では,波動関数から得られる電子密度は各原子軌道の被占状態に基づき原子間に分配される。この方法は広く利用されてはいるが,以前から指摘されているように,結果が使用した基底関数系に強く依存する。そのため,たとえば表1に示されるように,しばしば非現実的な結果を与える[16,17]。自然ポピュレーション解析[18]は,Mullikenの方法のもつ問題点をほとんど取り除いた方法として知られる。しかし,この改良法はab initio波動関数の使用を前提としている。
量子力学的に計算された波動関数から原子電荷を生成する第二の方法は,分子静電ポテンシャルを利用して電荷を誘導する方法である[7-11]。これは,まだ最近開発された新しい方法であり,静電ポテンシャル(ESP)適合法とも呼ばれる。荷電密度は明確に定義される関数量である[19]。それは,分子に関する重要で詳細な情報を含んでいると考えられる。すべての電子は,空間における荷電の分布に幾ばくかの寄与をしているからである。荷電密度はX線回折により実験的に測定することもできる[20]。しかし,この方法は費用と時間がかかりすぎるため,標準的な方法として利用することは不可能である。分子の電子的性質に関するかなり正確なイメージは,一組の原子電荷を使用して3D電子密度を再現することによっても得ることができる。ESP適合法では,分子を取り囲む空間の各点で量子力学的に計算された荷電密度をできる限り正確に再現するため,原子電荷の最小二乗適合操作が試みられる。この方法は,Mullikenのポピュレーション解析よりもはるかに良い結果を与える[9,11]。
ある方法で求めた荷電分布が信頼性があり,実際の分子電子分布を正しく反映しているか否かは,実験データに基づいて判定されなければならない。比較的容易に測定できる物性の一つは分子双極子モーメントである。分子双極子モーメントは,原子点電荷に基づいて簡単かつ迅速に計算できる。しかも,多数の化合物に対する実測値が文献に報告されているので,その計算値を実測値と比較することができる[21]。ちなみに,比較は剛体分子を対象に行なうべきである。双極子モーメントは分子の配座に強く依存するからである。ある特定の電荷計算法の適用性を検討する場合,柔軟な分子全体ではなく,小さな硬いフラグメントだけを取り出し調べることも多い。表1は典型的な小剛体構造に対する双極子モーメントの計算値と実測値をまとめたものである。双極子モーメントは,これまでに紹介した幾つかの手順に従い,様々な方法と基底関数系を使用して計算されている。双極子モーメントは量子力学的に定義できる属性である。したがって,それは波動関数から直接計算することもできる。表1でSCFと記された欄に記載された値はそれである。6-31G**のように大きな基底系を使用して得られた結果は,実験値と特によく一致している。
特定の分子系を研究するに際しどの方法を使用すべきかの選択は,幾つかの因子に依存する問題である。分子の大きさは重要な役割を演じており,また使用できるコンピュータの性能も制限因子の一つである。
トポロジカルな方法は,電荷と関連のある物性の合理的な推定値をきわめて迅速に提供できる点で,量子化学的方法よりも優れている。トポロジカルな方法で求めた双極子モーメントは一般に実験値とよく一致するが,これは一部実験結果に対する校正の成果である。この方法の主な欠点は,分子の幾何学的形と配座を無視している点である。もちろん,パラメータ表にない原子タイプを含んだ分子を扱う場合には,正しい結果は期待できない。たとえば, 表1のメチルシランを見てみよう。Gasteiger-Hueckel法による結果が示されていないが,これはケイ素に対するパラメータが欠けているためである。
分子荷電密度から原子電荷を計算する方法は,もしその結果が相互作用エネルギーを計算するために,経験的エネルギー関数で使用されるのであれば,最良の選択肢であろう。表1からも推論できるように,大きなab initio基底系を使用することは必ずしも必要ではない。小さな基底系や半経験的なAM1法を用いても,実験値ときわめてよく一致する結果が得られるからである。しかし得られた双極子モーメントの質は,原子点電荷の生成に使用した方法に強く依存している。分子荷電分布から直接得られた結果(SCF)は,いずれもMullikenポピュレーション解析の結果よりも実験値とよく一致している(表1参照)。後者の方法は,どの基底系を用いても粗い間違った双極子モーメントしか与えない。
分子が100個以上の原子から構成されている場合,分子全体に対して波動関数の正確な計算を行なうことは不可能である。この障害は,分子を部分的に重なり合うフラグメントへ分割し,各フラグメントに対して計算を行なうことにより回避することができる。計算の終了したフラグメントは,その性質が元の分子の属性を正しく反映しているという期待の下に,復元され元の大きな構造に戻される。
しかし,一連の分子に対して質の高い点電荷が得られていても,分子の類似性が問題にされている場合には,これらの量はあまり議論の役に立たない。分子類似性は3D荷電分布に基づかなければ,十分な議論ができないからである。この重要で明確に定義された3D荷電分布は,MEPの形で表したとき初めて最も有効な利用が可能になる。
分子静電ポテンシャル(MEP)は,分子を構成する核と電子の集合が分子近傍の空間点に作り出した荷電密度と単位正電荷(プロトン)との間に生じる相互作用エネルギーとして定義される量である。実際の計算では,その対象となるMEP点の数を制限するため,一般に打切り値が設定される。MEPは分子モデリングの研究を支える非常に有用な手段である。それは分子の静電特性を記述する量として,分子相互作用を解析し予測するのに利用される。MEPを生成する方法は2つある。第一は,量子力学的に誘導された波動関数から直接MEPを計算する方法である。この方法は単刀直入かつ正確であり,理論的にも一番望ましいが,計算に時間がかかりすぎる。第二は,分子荷電分布を表す原子局在電荷に基づきMEPを計算する方法である。これは第一の方法に比べより簡単である。MEPは静電相互作用に関するCoulomb方程式を適用して計算される。もちろん,第一の方法がはるかに優れている。もし十分精確な波動関数が得られるのであれば,ぜひとも第一の方法を使用すべきである。
波動関数から直接誘導されたMEPの基底系依存性を調べた研究が数多く報告されている[22-25]。AM1法で求めた静電ポテンシャルは,ab initio法による結果と良く一致する[22] 。AM1法は,ab initioレベルで処理できない大きな分子に対して特に有効である。
MEPの可視化
分子静電ポテンシャル(MEP)の表示には様々な手法が使用される。分子の比較研究へのMEPの迅速かつ容易な適用を妨げている主な障害は,データ点の数が多すぎることである。分子静電ポテンシャルを特定の分子平面内で切断し,2D等高線図の形で表す方法は,MEPを可視化する手段として広く利用されている。図はグラフィックス画面上にカラー表示され,実時間で処理される。各等高線はそれぞれエネルギーが等しい点を結んでいる。原子核の影響が強い領域は,単位正電荷と反発的な相互作用を行なうため,正のポテンシャルを生成し,一方,電子密度が高い領域は,単位正電荷と求引的な相互作用を行なうため,負のポテンシャルを生成する。
表示モードを2Dから3Dへ切り替えることにより,さらに高度で複雑な表示を行なうことも可能である。一般に,分子は等ポテンシャルの殻によって完全に覆い包まれているため,このような切り替えを行なっても問題は生じない。殻の面上の各点における静電ポテンシャルは,すべて同じ符号と大きさをもっている。この表示法の助けを借りると,分子の回りの正と負の荷電領域の全体的な分布状態を明快に示すことができる。2D等高線図は分子静電ポテンシャルの姿を常に的確に表現できるとは限らないが,3D等ポテンシャル面を利用すれば,様々な分子の比較と定性的な解釈が常に可能である。
分子静電ポテンシャルの第三の表示法は,分子表面の計算と可視化の問題に密接に結び付いている。そのため,分子表面の定義についてここで少し詳しく触れなければならない。分子表面の正式な取扱いでは,原子核は点として扱われ,電子雲は原子核に中心を置く球で近似される。この電子雲の球をvan der Waals半径で表したとき,すべての球を加え合わせて作られる分子表面はvan der Waals表面と呼ばれる。van der Waals表面は分子の3D体積の境界面にほぼ等しい。分子モデリング研究でよく用いられるもう一つの表面は溶媒接触可能表面(solvent accessible surface)である。これはConnolly表面とも呼ばれる[26]。Connolly表面は,プローブ分子がvan der Waals表面の上を転がるとき,その中心が描く表面のことである。
第三の方法では,静電ポテンシャルはvan der Waals表面やConnolly表面に色分けして表示される。分子表面の各色は,静電ポテンシャルのエネルギー値を表している。この方法を用いると,分子の形とその静電的性質を同時に表示することができる。しかし,大きな分子の場合には,その像は非常に複雑なものになってしまう。この問題は,幾つかの表示モードを組み合わせることにより解決できる場合がある。ある表示モードでは隠れている領域が他のモードでは見えるようになるからである(図1参照)。
分子静電ポテンシャルは,原子点電荷よりもはるかに信頼できる静電的反応性の指標である。MEPとその3D表示は,リガンド‐高分子受容体間の相互作用を解析し予測するための有用な手段である。
同じ受容体部位へよく似た様式で結合する分子の静電ポテンシャルは,共通の特徴を備えているはずである。対応する原子を1つずつ突き合わせる操作が満足な結果をもたらさないような場合でも,MEPに基づいて重ね合わせを行なうと適切な解答が得られることはよくあることである(2.5.3節参照)。
一例として,ヒスタミンH2拮抗薬の静電ポテンシャルを調べた我々の研究結果[27]を図2に示した。シメチジン(cimetidine)のイミダゾール環とチオチジン(tiotidine)のグアニジノチアゾール環は,静電ポテンシャルの立場から見たとき,互いにうまく重なり合うことがお分かりいただけよう。
分子間の非共有相互作用は多くの生物学的過程に関与している。たとえば,受容体へのリガンドのドッキング,基質と酵素との相互作用,タンパク質の折りたたみなどはこの非共有相互作用によるものである。非共有力は,また結晶における分子配列の幾何学的形や対称性も支配している。一般に,結合は相互作用により生じるエネルギーがvan der Waals斥力に打ち勝った場合のみ生成する。分子相互作用場の概念は,互いに接近し合う分子間のエネルギー状態を調べる一つの手段として有用である。この相互作用場は,標的分子とその周囲に張り巡らされた3D格子の中を移動する化学的プローブの間に生じる相互作用エネルギーの変化を記述している。プローブは,結合の相手分子やその断片の化学的特徴が反映されたものでなければならない。分子相互作用場は,コンピュータ・グラフィックスを利用することにより3Dエネルギー等高線図として可視化することもできる。エネルギーが正の大きな値をとる領域はプローブと反発し合う領域であり,一方,負の大きな値をとる領域はエネルギー的に有利な結合性領域である。
分子相互作用場はGRID[28],CoMFA[29],HINT[30]のようなプログラムを用いて計算される。GRIDは,分子相互作用場の研究で最も広く利用されているプログラムの一つである。このプログラムは小分子だけではなく,酵素のように大きなタンパク質分子に対しても使用できる。必要な入力データは構成原子の直交座標値だけであり,選択できるプローブの種類は非常に多い。相互作用エネルギーは標的分子を取り囲む規則的な格子点上で計算されるが,関心の対象が標的分子の特定フラグメントにある場合には,計算格子をその範囲に局限することもできる。エネルギーの計算結果はデータファイルに保存される。このファイルは,グラフ表示や解析の機能をもつ一般的な分子モデリング用プログラム[31-33]のほとんどへ転送することができる。3D等高線図は任意のエネルギーレベルで作成することができ,グラフィックス画面上に標的分子と一緒に表示される。等高線は速やかに表示されるので,ユーザーはほとんど実時間で結果を検討することができる。
本章では,小分子に対する相互作用場の計算に焦点を合わせることにしよう。高分子に対する相互作用場は後章で議論する予定である(4.6.2節参照)。
GRID場の計算では,たとえば水分子,ヒドロキシル基,カルシウムイオンといった小分子,化学フラグメントおよび原子がプローブとして使用される。これらのプローブは,結合の相手―たとえば受容体タンパク質の潜在的な結合部位,結晶構造内部の隣接分子―の化学的特徴をシミュレートしたものである。プローブはGRID計算の間,標的分子を取り巻く規則的な3D格子点を系統的に移動する。そして各格子点において,プローブと標的分子間の相互作用エネルギーが,次の経験的エネルギー関数に基づいて計算される。
Etot=Evdw+Eel+Ehb (1)
ここで,Etotは全相互作用エネルギー,Evdwはvan der Waals相互作用エネルギー,Eelは静電エネルギー,Ehbは水素結合の形成による相互作用エネルギーをそれぞれ表している。
van der Waals相互作用エネルギーは,非結合原子間の分散引力と分散斥力を加え合わせたものである。プローブ原子は,原子間の反発と電子雲の重なりにより標的分子の内部へ入り込むことはできない。斥力は経験的なエネルギー関数によって評価できるが,二原子間の距離がvan der Waals半径の和よりも小さくなると,それは大きな正の値をとるようになる。分散相互作用の引力項は,核の回りの電子の相関的な運動により生じる誘起双極子相互作用によるものである。無極性分子における分散引力と近距離斥力の間の平衡は,Lennard-Jonesポテンシャル[34]を用いて記述することができる(2.2.1節の式5参照)。このポテンシャルはGRIDプログラムにも組み込まれている。
静電相互作用は,リガンドと高分子受容体の間のはたらく遠距離引力として重要である。
Coulomb方程式(2.2.1節の式6参照)は数学的形式が簡単であるため,分子力学プログラムにおいて静電項の計算に広く利用されている。この式の欠点は,誘電率の異なる複数の分子から成る不均一系を十分表現できない点である。溶質と溶媒の間の不連続性は,元の式を拡張したより包括的なCoulomb方程式[28]では考慮されている。GRIDに組み込まれているのは,この拡張されたCoulomb方程式の方である。
水素結合の指向的性質は,分子間相互作用の特異性を決定する上できわめて重要な役割を演じている。したがって,相互作用エネルギーを正確に評価するためには,この分子間引力項を正しく記述する努力が特に必要とされる。水素結合は,正に帯電した水素原子と電気的に陰性な受容原子の間にはたらく中距離相互作用である[35]。水素結合における受容原子と供与原子間の距離は,それらのvan der Waals半径の和よりも小さい。分散力や点電荷に基づいた他の非共有相互作用とは対照的に,水素結合相互作用は受容側ヘテロ原子の孤立電子対の性質と配向に依存した方向性をもっている。
GRID法は,このような角度からの要求にも応じられるよう,相互作用エネルギー式に水素結合に関するエネルギー項をはっきりとした形で組み込んでいる[36]。この項の関数形は実験データに適合するように展開され,相互作用の方向,様式,強度などのパラメータも,すべて実際に観測された結晶学的な実験データに基づいて決められている。
GRIDに組み込まれたプローブは,水素結合能,van der Waals半径,原子電荷といった様々なパラメータを用いて,その属性を規定できるようになっている。したがって,属性を詳細に指定すれば,高分子の活性部位で見出される重要な官能基の性質をきわめて的確に反映した,特異性の高いプローブを選び出すことが可能である。一例として,重要なプローブ3種の属性とパラメータを表2に示す。
GRIDは,また標的分子に含まれる原子のタイプを調べるためのパラメータ表を備えている。各原子のタイプを規定するパラメータは,van der Waals,静電および水素結合相互作用の強度である。GRIDプログラムは,綿密なパラメトリゼーションと組み込まれたプローブの多様性により,正確な相互作用場を知る手段として,小分子だけではなく高分子に対しても広く適用されている。
分子相互作用場は,分子モデリングの分野で広く利用されている[37-42]。使用される戦略は,リガンドと標的高分子に対して入手できる構造情報に依存する。高分子の3D構造が知られている場合には,我々は相互作用場を利用して,リガンドにとって有利な結合領域の位置を正確に突き止めることができる。このような領域は,次に同じ受容体へ結合する新しいリガンドを設計するための出発点として利用される。この戦略の詳細な手続きは, 第3章で説明する予定である。
しかし一般には,高分子受容体側の構造情報は全く得られず,リガンドの性質しか分からない場合がほとんどであろう。分子相互作用場は,このような状況下でも,受容体結合部位の大まかなモデルを作り出すのに役立てることができる。もちろん,モデルの作製に使われるリガンド分子は,すべて同じ受容体部位へ類似の機構で結合することが前提条件である。相互作用場のパターンを比較できるのは,そのような条件が満たされた場合のみである。リガンド分子は,またエネルギー等高線の相対位置と大きさもよく似ていなければならない。
比較に用いられる等高線のエネルギーレベルは,使用されたプローブの種類に強く依存する。一例として, 図3にCa2+チャンネル遮断薬ニフェジピン(nifedipine)に対する2種類の相互作用場を示した。
リガンドの相互作用場を比較することにより,我々は受容体の結合部位が備えていなければならない物理化学的属性を特定することができる。これらの属性は,結合部位のモデルを作製する際に化学構造に翻訳され利用される。もし標的高分子がタンパク質であるならば,そのモデルは,対応する相互作用領域に条件を満たすアミノ酸断片を適当に配置したものになろう。アミノ酸断片は,すべての活性リガンドに共通な結合部位となるよう選択され配置されなければならない。たとえば,疎水性の相互作用場はフェニルアラニン,トリプトファン,バリン,ロイシン,イソロイシンのような疎水性アミノ酸の位置を表していると考えてよい。もちろんアミノ酸の種類を正確に特定するためには,受容体マッピング(receptor mapping)と呼ばれるさらに進んだアプローチが必要である。このアプローチについては,次の2.5節で取り上げることになる。
検討する化合物の数が増えると,共通に存在する相互作用領域を通常の方法ですべて発見することは困難になろう。この問題は,同じプローブを使用して得られた様々なリガンド分子の相互作用場から,計算により共通する領域を抽出することで解決される。共通領域は,場の格子点を1つずつ比較する方法で数学的に検出することが可能である。ヒットした格子点はファイルに保存された後,共通相互作用場の作成に利用される[43]。
分子相互作用場をさらに詳しく比較し解析したい場合には,GOLPE[44]やPLS[45]のような化学計量学的アプローチが有用である。最近まで,分子場に基づいた構造活性相関研究のほとんどは,定性的な議論に終始する傾向があった。これは,化学計量学的アプローチで使用される統計的方法が数学的にも方法論的にも難解であり,化学計量学の専門家しか実際には使用できなかったことに原因があると考えられる。しかし,分子相互作用場が生物活性に寄与する因子の確認に役立つ価値ある手段であることは,定性的な解析による多くの研究例[36-43] からも明らかである。
CoMFA法(比較分子場解析法:Comparative Molecular Field Analysis)[29]は定量的な三次元構造活性相関(3D-QSAR)を研究するために開発された方法である。このアプローチでは,一群の化合物における生物活性(または化学活性)の変動と3D構造情報との関連を調べるために,統計的方法―化学計量学的方法(chemometrical method)―が利用される。分子間における生物活性の差異は,分子を取り巻く非共有相互作用場の形や強度の変化としばしば密接な関係にあるというのが,CoMFA法の根底をなす考え方である。言い換えると,化合物群の生物学的性質を理解するに必要な情報は,立体場と静電場の中にすべて入っていると考えるのである。CoMFA法では,GRIDアプローチのときと同じように,分子は立体格子の中央に置かれ,各格子点において,分子とプローブの間の相互作用エネルギーが計算される。分子の並置(alignment)はこの操作の前提となる重要な条件である。一群の分子は,薬理作用団の決定に使用される周知のアルゴリズムを利用して並置される(第3章参照)。開発者によれば,CoMFA法は薬理作用団を決定するために利用することもできる。典型的なCoMFA研究では,大まかに並置された化合物から出発する。あらかじめ定められたすべての格子点(3D格子の交点)で相互作用エネルギーが計算され,分子の立体場と静電場が作成される。場の共通領域の相対位置は,たとえばPLS[45] のような化学計量学的方法の助けを借りれば見つけ出すことができる。発見された共通場領域は,次にいわゆる「場適合」操作により被検化合物間の最適な重ね合わせを実現するために利用される。CoMFA法の出現によって,我々は,被検化合物が異なる構造系列に属するため,1原子ごとの突き合わせ操作では並置がうまくいかないような場合でも,分子相互作用場の共通領域を利用して,その最適化を図ることができるようになったのである。
CoMFA法は一般に広い範囲で使用できるが,その適用の過程で様々な問題が発生していることも事実である。解析の結果は,リガンドの配座,並置の状態,相互作用場の記述パラメータ,使用した統計的方法などに強く依存する[46]。この方法は,経験豊富なユーザーにとっては強力な手段となりうるが,初心者には恐らく理解しにくいむずかしい方法であることを指摘しなければならない。CoMFA法については,ここではこれ以上踏み込まないことにする。CoMFA法に関連したすべての特色と難点を詳細に記述することは,本章の範囲を越えているからである。このアプローチについて詳しく知りたい読者は文献47をお読みいただきたい。
既に議論したように,分子間にはたらく引力と斥力は様々な種類の相互作用によって制御されている。ここで取り上げる疎水相互作用も,そのような相互作用の一つである。分子間における疎水相互作用は,溶媒和殻を構成する溶媒分子の配向変化と結び付いたエントロピー効果や遊離型溶媒分子によってもたらされる複雑な過程である。有効な疎水相互作用が生じるためには,相互作用する疎水性表面は互いに密接に接触していなければならない[48,49]。疎水結合に対する理解を助けるために,一つの例を上げることにしよう。タンパク質の深い結合性空洞(キャビティー)にある無極性領域は直接溶媒和されることはない。近くの水分子はキャビティーをおおい隠し,互い同志の分子間水素結合によって安定な氷山状構造を形成している。キャビティーに入り込んだ基質はキャビティーの疎水表面と相互作用し,水分子の整然とした氷山状構造を破壊する。この破壊は,エントロピーと系全体の自由エネルギーを増加させる結果となる[49]。また基質分子の脱溶媒和は,遊離した溶媒分子を新たに生成する。
エントロピー効果は,簡単な計算法がないために,現状ではほとんどのモデリング研究において通常無視されている。しかし疎水結合すなわちエントロピー効果は,薬物‐受容体相互作用[50]やタンパク質の折りたたみ[51]において重要な役割を演じている。我々は,これらの過程のエネルギー収支に,是が非でも疎水相互作用を含めるようにしなければならない。
疎水性は,分子とその環境の間に生じる相互作用の熱力学的情報を数量化した分子の経験的属性の一つである。実験結果に基づいて疎水性効果を説明する試みは,これまでに幾つかなされている。実験的に測定できる疎水性の最も重要な尺度は,水と有機相の間における分子の溶媒分配係数(log P)である。これは,また幾つかのグループ[52-54]により開発された経験的方法を統括し改良する目的でも利用されている。log Pは,様々に置換された分子集合に対する溶媒分配の実測データから誘導された,疎水フラグメント定数を用いることにより予測することができる。これらのフラグメント定数は,元の分子集合中に見出される構造要素の相対的な疎水性を表しているので,化合物の全疎水性(log P)はそれを構成する要素のフラグメント定数をすべて加え合わせれば計算できるという訳である。フラグメント定数は,生物学的に重要な様々な有機分子の疎水性を高い精度で予測できることが確認されている。
logPは疎水性を簡単に一次元的に表現したものであり,単に分子全体の性質を反映しているにすぎないことに注意したい。リガンドと高分子の間の分子相互作用をさらに詳細に解析する必要がある場合には,この尺度だけでは不十分である。
上述の理由から,溶媒分配係数を利用して,疎水性の3D指標を作り出す試みが幾つかなされている。静電場との類比から編み出された疎水場アプローチは,そのような試みの一つである。この手法は,たとえばHINTプログラム[30]の中に組み込まれている。
疎水相互作用のHINTモデルは,溶質と溶媒の間の分子相互作用を写し出すもう一つの物理化学的性質として,溶解度が利用できるという事実に立脚している。HINT理論の枠組では,フラグメントレベルの溶媒分配データ(疎水フラグメント定数)は,さらに簡単な原子疎水定数[55]へ変形されたのち利用される。HINT法の基幹をなすパラメータはこの原子記述子である。原子疎水定数は,被検分子を構成する原子のすべてに割り付けることができる。しかも,分配係数と溶解度の実測データから誘導されているため,その値は疎水相互作用だけではなく,静電項やvan der Waals項のような他の分子相互作用成分も含んでいる。そのため,この方法で作り出された場は,疎水パラメータと親水パラメータの両者を織り込んでいる。この場は一般にヒドロパシー場(hydropathic field)と呼ばれる。ヒドロパシー場の計算は経験的関数を使用して行なわれる(関数形の詳細については文献55を参照されたい)。アルゴリズムに含まれるのは,原子疎水定数,溶媒接触可能表面の原子断面積および距離関数である。距離関数は,生体環境における疎水効果の距離依存性を適切に記述する上で必要である。HINTによる3D分子格子図の作成は,GRIDやCoMFAと同じ方式で行なわれる。
HINT研究から得られる結果は,分子の回りに疎水場と親水場が同時に表示された組合せ等高線図である。正符号をもつ格子点は疎水性領域を表し,負符号をもつ格子点は親水性(極性)領域を表している。データが経験的な性格をもつために,等高線を表示するエネルギーレベルを設定する作業はかなり面倒である。選択したエネルギーレベルの如何によって,場の可視領域の大きさが著しく変化するからである。一般的な目安として,親水効果の等高線を疎水効果のそれよりも2〜5倍高いエネルギーレベルにすれば,均衡の取れた等高線図が得られることが多い[56]。
一例として,Ca2+チャンネル遮断薬ニフェジピンにおける疎水場と親水場の外観を図4(a)に示した。
ヒドロパシー場の解析から得られた情報は,次に来る戦略に利用することができる。たとえば,特定分子系列の近傍における疎水性領域と親水性領域の分布に関する定性的な情報は,未知受容体高分子の3D地図を作成する際に活用することができよう。もし問題の分子系列が大きく複雑なものであれば,インターフェイスを利用して,生成したデータセットをCoMFAへ直接読み込ませ,さらに綿密な解析に委ねるようにすればよい[57]。
作成されたヒドロパシー場は,高分子受容体の構造が知られている場合には,リガンドの構造を最適化し,生物活性の増強を図る目的で利用することもできる。文献57は,ヒドロパシー場の応用について詳しく知りたい読者にとって参考になるであろう。
超分子空間に表示された疎水的および親水的性質の分布は,分子表面へ射影することもできる。たとえば,MOLCADプログラム[58]は,疎水性のような局所分子属性を色分けしてマッピングするのに,そのスクリーンとして分子のConnolly表面[26]を使用する。分子表面の各点における局所疎水性に原子やフラグメントの影響を正しく反映させるためには,距離依存型関数が必要であるが,この問題は,たとえば分子疎水ポテンシャル[59]を導入することにより解決される。このポテンシャルは分子静電ポテンシャル(MEP)に相当するものと考えてよい。MEPのところでも指摘したように,分子表面への局所属性の射影は,属性記述子の分布を視覚的に把握し解釈する上できわめて有効である。分子表面へ疎水性を表示することの主な利点は,ヒドロパシー場による方法と比較して,タンパク質のように大きな分子系の解析がはるかに容易に行なえる点である。2つの手法の理論的背景は同じであるから,得られた結果は定性的に同等である。いずれの方法でも,log Pの正確な実測値が知られている分子を利用すれば,信頼度の検定を行なうことができる。しかし分配係数は,荷電分布と同じように分子の配座による影響を強く受ける。また,水相から疎水相へ移行したとき,分子は配座を変える可能性があるから,状況はさらに複雑である。残念ではあるが,このような事情により,実際に検定に使用できるのは,分配係数が配座の影響を受けない剛体または半剛体の構造をもつかなり狭い範囲の分子に限られる。図4(b)はニフェジピンにおけるMOLCAD分子疎水ポテンシャル面の一例を示している。
我々は,これまで分子の物理化学的性質を計算し可視化する方法について説明してきた。本節では,化合物の薬理作用の理解と予測へこの知識をどのように活用するかについて論ずることにしよう。薬物分子の薬力学的効果は,ほとんどの場合,生理学的に重要な高分子タンパク質との間の相互作用によってもたらされる。その高分子は酵素か受容体である。いずれの場合も,薬物分子がその中に入り込んで結合を形成することのできる,高度に特異的な三次元のキャビティーが存在する。同じ酵素や受容体に結合し,定性的によく似た活性を発現する化合物は,類似の結合様式をもつに違いない。すなわち,これらの化合物は,結合相手の高分子から見て,立体的に対応した位置に同じ化学的性質を示す構造要素をもつはずである。同じ薬理作用を示す同族体群に共通するこのような薬理作用団(pharmacophore)を同定することは,分子モデリングを利用して解決すべき主要な課題の一つである。受容体のほとんどは,まだその3D構造が解明されていないので,この薬理作用団に関する情報は,薬物‐受容体相互作用を分子レベルで解明する上できわめて重要なよすがとなるものである。
では,薬理作用団を見つけ出すためには,どのようにしたらよいのであろうか。この疑問に答えるためには,我々はまず最初に薬理作用団を構成する要素を定義し,重ね合わせるべき官能基や原子団をはっきりさせなければならない。もちろん,これは自動化された客観的な手続きにより完全に答えることのできる性質の問題ではない。二分子間で突き合わせる原子対は前もって指定されなければならないからである。この操作は,もし構造活性相関に関する既知情報を活用しなければ,無用なデータを多数作り出すことになろう。構造活性相関情報の利用は,このような無駄を省き,検討を要する可能性の数を大きく減らすのに有効である。また,リガンド間の類似性は,分子全体を問題にしているのではないことに注意してほしい。リガンド分子のほとんどは,受容体へ結合したとき,結合部位の中に完全にすっぽり収まるわけではないからである。このことは,検討を要する可能性の数をさらに減らすことになる。
水素結合が薬理作用団の重要な要素であると思われる場合には,被検化合物の原子パターンに孤立電子対の方向と距離を付け加えるべきであろう。これは,たとえば,対応する位置に擬原子を置くことによって実現できよう。原子位置は,水素結合の受容部位か供与部位(ヘテロ原子へ結合した水素のみ)かに応じて異なるフラッグを付け標識しておくとよい。AUTOFIT[1]の重ね合わせモードでは,このような原子情報が最初の吟味に利用される。また,芳香環のような平面形の要素は特別の構造単位として取り扱い,環全体の代わりに,たとえば環の中心を突き合わせ点として指定してもよい。このような取扱いは,他の平面原子団に対しても同様に適用することができる。
被検分子群がきわめて柔軟な同族体のみで構成されている場合には,共通薬理作用団の探索は非常にめんどうで退屈なだけではなく困難を伴い,しかも意味のある結果が得られないことすらある。この探索作業は,もし剛体か少なくとも半剛体の化合物を中に含めるならばはるかに容易となり,その意義も著しく高まることになろう。もちろん,これらの化合物は高い活性を示す化合物でなければならない。そうでなければ,柔軟なリガンドに対する鋳型として利用できないからである。また,配座的に束縛されしかも高い活性を示すこのような分子を対象とすることで,我々は,被検分子が生物活性配座にあるか否かを確認する手間を省くことができる。
重ね合わせる分子の選択は,意味のある結果を得る上できわめて重要である。構造のよく似た化合物の重ね合わせは,簡単に行なえるけれどもあまり意味はない。得られる情報が少ないのである。効率良く情報を引き出すためには,骨格の異なる化合物を被検系列に含めるようにしなければならない。もちろんこのような系列では,重ね合わせの操作は単純な1原子ごとのものではなく,機能的に等価な構造要素や分子場を対象としたものになる。
さらにもう一つの問題点にも言及しておこう。不活性分子や低活性な分子をどのように取り扱うかという問題である。重ね合わせの操作は,最初は高い活性を示す分子だけで行なうのが得策であろう。しかし,そのようにして求めた薬理作用団は,次に不活性同族体や低活性な同族体を含めて検証され,改良の手が加えられなければならない。同じ方針は,同一受容体に作用する拮抗薬と作動薬の研究にも適用できよう。重ね合わせは,最初は作用の異なる2つの薬物群に対し別々に行なうべきである。しかし,このようにして得られた2つの薬理作用団モデルは,次の段階で恐らく1つに結合させることができよう。なぜならば,競合的拮抗薬は作動薬の受容体結合部位の少なくとも一部に結合することが多いからである。結合部位のこのような部分的共有は,作動薬と拮抗薬の間ではごく一般的である。しかし,必ずしも常にそうであるとは限らない。
重ね合わせの方法は幾つか知られている。よく用いられるのは,剛体回転やフレキシブルな適合操作により手動または自動的に突き合わせを行なう方法である。最小化は,適合原子対間の二乗平均(rms)変位と配座エネルギーの双方に対して行なわれる。その他重要な方法としては,分子表面や分子場の同等性に基づいて被検分子を並置する方法がある。
最も広く利用されているのは,対応する原子位置を最小二乗法を用いて重ね合わせる方法である。二つの分子は,対応する原子間の距離の二乗平均(rms)が最小になるように重ね合わされる。rms値は適合の良否を判定する尺度である。この方法は,一見よく似ている分子間の異同を検出する上できわめて有効であるが,突き合わせる原子対を前もって指定しなければならず,それが難点となっている。重ね合わせの結果は,指定した原子の如何によって明らかに異なってくるからである。原子の対応関係が前もって見出せない分子系列に対して適用できないのはもちろんである。さらに原子配置の厳密な類似性は,同じ受容体と相互作用するための前提条件ではない。薬理データと構造活性研究から,見掛け上類似性のない分子が共通の作用機構をもつと考えざるを得ない場合は非常に多い。最小二乗法による重ね合わせ操作は,このような場合に対しても明らかに無力である。
このような状況は,被検分子の数が多くない場合には,手動によるインタラクティブな重ね合わせを行なうことにより回避できるかもしれない。原理的には,分子が幾つあってもグラフィックス画面上で直接解析し,適合の状態を目で判断することは可能である。このやり方は確かに非常に創造的であり,実験的に見出された構造活性データの基礎をなす作用機構について,新しいアイデアを刺激する可能性を秘めている。しかし一方では,このような主観的な手続きは当然偏見を含まざるを得ない。しかも計算による最適化ができないため,その正確な再現はしばしば不可能である。
Marshallら[2,3]により開発された活性類似体アプローチ(Active Analog Approach)は,薬理作用団モデルの作製に利用できる迅速で効率的な探索手法である。この方法では,まず系統的なアルゴリズムを利用して,同族系列の分子に対し立体的かつエネルギー的に許容される代表的な配座が選び出される。そしてそれらの配座の各々について,受容体との相互作用に重要と考えられる薬理作用基間の相互距離が計算される。もし対応関係を見出すために,ある分子におけるこのような一組の距離を他のすべての分子のそれと比較しようとするならば,問題を解決することは小さな分子の場合を除いてまず不可能である。しかし幸いなことに,薬理作用団を特定するという目的には,分子の配座空間をすべて完全に調べる必要は全くない。活性なリガンドに共通な配座空間領域だけが我々の関心事なのである。前の2.5.1節でも指摘したように,柔軟な分子系列の配座解析は,剛体か半剛体の分子を含めて行なうのが常道である。活性類似体アプローチでも,配座探索は最も硬い分子から出発する。そしてこの分子に対する距離マップがまず作成され,次にこれらの距離を拘束条件として,より柔軟な分子の配座探索が実行される。すなわち,活性な剛体分子に対する探索結果は,被検系列の他のすべての成員がとりうる配座空間を探索する際に,その基盤として用いられる訳である。活性分子はすべて受容体モデルにフィットしなければならないから,探索は,そのモデルを満たす配座空間領域に限定して行なえば十分である。たとえば,モデルから導かれる拘束条件を満たすために,もし一対の原子間の距離が一定の範囲内になければならないとすれば,そのような条件を満足するねじれ角範囲だけを計算の対象とすればよい。活性類似体アプローチの威力は,28種のアンギオテンシン転換酵素(ACE)阻害薬を扱ったACE活性部位モデルの構築研究において実証されている[4]。この方法を適用した結果,探索時間は従来の系統的探索法を用いた場合と比較して3桁以上短縮されたのである。
一方SEAL[5]は,活性類似体アプローチとは対照的に,1原子ずつの重ね合わせ操作を表立って行なわないマッピング法である。この方法は,系列の異なる分子間で迅速な比較を行ないたいときに有用である。適合の良否は,二分子間のすべての原子対について,その変位を加え合わせた類似性スコア(similarity score)を用いて判定される。並置関数は,1原子ずつの重ね合わせアプローチにおけるような特定の原子対だけではなく,理論的に可能な原子対をすべて考慮したものになっている。そのため,重ね合わせの結果には,分子の全体的な形状がある程度反映される。プログラムは,また並置操作の際に物理化学的性質を含めるオプションも備えている。対求和項の中に生物活性に重要と思われる物理化学的性質を付け加えることができるのである。原報には,立体的性質としてvan der Waals半径,電子的性質として点電荷を追加し,並置の最適化を試みた事例が報告されている。
また,計算プロトコルの最初のステップで,重ね合わせの対象となる原子中心や部位点を偏りなく自動的に特定するためのマッピング手法も知られている。ここで言う部位点とは,水素結合の受容性や供与性といった性質を示す分子表面の点のことである。このような機能を備えている市販プログラム・パッケージとしては,APEX-3D[6],CATALYST[7]があり,その他にもDISCO[8],RECEPS[9],AUTOFIT[1]といったプログラムが報告されている。前にも述べたように,重ね合わせは,指定された原子と部位点のすべての可能な組合せに対して行なわれ,得られた並置結果はrms値に基づいてランク付けされる。
分子は,その原子骨格ではなく,van der Waals体積の表面や外側の特徴に基づいて互いを認識するので,分子類似性の評価は分子場を考慮したものでなければならない。したがって,この分子場のマッピングと比較は,重ね合わせアプローチにおける重要な研究課題の一つである。マッピングの操作では,分子は一定間隔で描かれた場の3D格子の内部に布置され,各格子点は荷電分布,疎水ポテンシャルあるいは体積情報といった性質尺度で特性付けられる。最適化操作によって有意で確実な結果が得られるよう,類似性の閾値を定義することもある。個々の格子点や関連格子点のクラスターは様々な重みをあてがわれ,構造活性相関へできる限りの配慮が払われる。適当な分子―配座的に束縛された分子が望ましい―が鋳型分子として選択され,その格子特性値が様々な性質に対する基準として使用される。すべての試行分子は回転と並進の操作により,格子点の特性値が鋳型分子のそれと最もよく適合するよう位置を調整される。配向を探索するこの計算は,きわめて多大な時間を必要とする。そのため,場の様々な性質を利用して計算の効率化を図る方法がいろいろ提案されている。Manautら[10]は分子静電場に基づき,分子表面の相互類似性を効率よく最大化する方法を考案した。またClarkら[11]やDeanら[12]のグループは,Lennard-Jonesポテンシャルから誘導される場の物理化学的性質を利用する試みや,通常の格子評価法をGauss型関数の積分で置き換え,静電ポテンシャルを近似する試みについて報告している。ちなみに適合度の指標は,たとえば共通に占有される格子点の数と格子点の総数の比として計算される。
分子表面を重ね合わせるための手段は既に存在している。それは分子間での原子の対応関係を要求しないため,系列を異にする分子を重ね合わせる目的にも有効に利用できる。しかし,この方法が日常的な研究手段として利用されるようになるためには,複雑な計算を処理し,大量の分子配座を迅速にこなすことのできる高性能なコンピュータの普及が不可欠である。
3D分子構造を生成する方法の一つに,分子の2Dまたは2.5D表現から出発し,その情報を3D形式へ変換する方法がある。2.1.3節で紹介したスケッチ・アプローチは,作画された平面構造式を3D情報へ変換することができたが,本節で紹介するCONCORD[1,2]とCORINA[3]は,たとえば製薬会社のデータベースに蓄積されている幾千もの化合物に対する構造情報を利用し,その2Dまたは2.5D結合表を3D分子構造へ変換する機能をもつプログラムである。この目的に適した大型商用データベースの一つは,Chemical Abstractsから入手することができる[4]。
CONCORDは,大型データベースにエントリーされた潜在的な生物活性分子の2D結合表から3D構造を生成するために特に開発されたプログラムである。3D構造を生成するために,CONCORDはきわめて詳細な結合長テーブルを使用する。そして結合長を割り当てるに先立って,原子番号,混成,結合形式のような情報の他,結合に含まれる原子の環境を細かく分析する。結合長のこのように厳密な選定手続きは,環構造を組み立てる際には特に重要である。正しい値からのわずかなずれが,環の配座に劇的な影響を及ぼすからである。
変換に当たっては,プログラムはいわゆる「小環の最小集合」をまず同定し,個々の環系に対して論理的解析を行なう。そして環に隣接する構造や束縛条件に基づき,環の組み立て方が決定され,さらに平面性や立体化学的束縛を考慮して,その大まかな配座が定められる。
多環式系の場合,縮合位置の原子が指定されなければ,CONCORDは最低のエネルギー含量のもつ異性体を生成する。環を組み立て接続した後,プログラムは内部ひずみを取り除くために全体の配座を修正し,環内のすべての原子に均等にひずみを分配する。この操作は十分に緩和された幾何学的形をもつ環状構造をもたらす。
構造生成の次のステップでは,非環式部分構造が付け加えられる。結合長と結合角は,ここでも前もって用意されているテーブルの値が使用される。組み立てた構造の内部に密接なvan der Waals接触が生じるのを避けるため,ねじれ角はエネルギー的に許される配座になるまで変更される。CONCORDの主な長所は,計算速度の他,分子のトポロジーが組立ての各段階で考慮されるという点である。その結果として,CONCORDは少ない計算時間で質の良い3D構造を生成することができる。2D情報の大型データベースを3D空間へ変換しようとする場合,これは採否を決める重要な判断基準である。
CORINAの動作仕様はCONCORDのそれと非常によく似ており,環系を生成するときの出発点はCONCORDと同じである。しかし,CORINAは環系を接続する際に異なるアプローチを使用する。まず環が融合され,おおざっぱな力場を用いて環配座のエネルギーが計算される。もし実際に選択した環の接続様式がエネルギー的に不利であることが判明したならば,エネルギー的に可能な他の新しい環配座が試みられる。安定な環構造が生成したら,続いて構造の最適化が行なわれる。
エネルギーに基づいたCORINAアプローチは,論理的規則を使用するCONCORDアプローチに比し効率が悪いため,環接続の問題を解くのに時間が余計にかかる。
環系が完成した後,非環式部分構造を構築する方式はCONCORDの場合と同様である。鎖は通常完全な伸張配座で環に付加される。もちろん,これはさらに改善を必要とする構造である。ねじれ角は,密接なvan der Waals接触のない配座が得られるまで回転を施される。このおおざっぱな配座探索の結果として,プログラムはまずまずの構造を生成する。
得られた配座は,結晶環境での配座や低エネルギー配座の一つに対応しているかもしれない。しかし,それは単に偶然そうなったにすぎないのである。可能なすべての低エネルギー配座を検出するためには,最終的に生成した構造に対し,さらに配座解析を行なわなければならない。
ここで紹介した2つのプログラムは,どちらも3D構造を生成する手段として有用であり,その最も重要な用途は,大きな2Dデータベースを3Dデータベースへ変換することである。生成した3Dデータベースは,続いて薬理作用団の3D探索に利用されることになる。
3D探索が価値をもつ理由は幾つかある。それは,たとえば配座解析や構造最適化のようなさらに進んだ研究のための出発点として役立つ,近似的な3D分子構造群を発生させることができる。もし得られた近似構造が十分質の高いものであるならば,それらは対応する分子構造の様々な性質を計算するためにそのまま使用することもできよう。
共通の受容体へ結合するためには,生物活性リガンド類は一定の化学的および幾何学的条件を満たさなければならない。受容体を構成する原子の3D座標がすべて分かっている場合には,受容体と相補的な関係にあり,結合部位と効率良く相互作用するリガンドを探し出すことは容易である。しかし,受容体の完全な構造は分からない場合がほとんどである。結合部位についての正確なデータが欠如している場合には,我々は次善の策として問題の薬物系列に共通する薬理作用団を定義してかかる必要がある。この薬理作用団モデルは,受容体についての限られた情報とこの受容体へ結合するリガンド類(活性,不活性)に関する追加情報を利用すれば作り出すことができる。得られるモデルはおおざっぱなものではあるが,大型データベースの3D探索における探索基準として十分有用である。
3D探索は,指定した化学的および幾何学的条件を満たす既知化合物だけではなく,有利なリガンド‐受容体相互作用に必要な特徴を備えた未知の化合物をも明らかにしてくれる。
3D構造の大型データベースと適当な探索用ソフトウェアは,有効な3D探索を行なう上で必要な前提条件である。一般に,3D探索に必要な3D分子構造は2つの方法によって得ることができる。第一の一番簡単な方法は,ケンブリッジ結晶学データベース[5]やCASレジストリー[6,7]のような既存の3Dデータベースを探索する方法である。これらの商用データベースは,GSTAT(ケンブリッジ結晶学データセンターから入手可能)のような3D探索プログラムやChemical AbstractsのCAS ONLINEサービスを介して利用することができる。3Dサーチ[8],ALADDIN[9]といった他の3D探索プログラムを使用してもよい。3D構造データベースを得る第二の方法は,前節で紹介したCONCORDやCORINAを利用する方法である。ユーザーは,これらのプログラムにより2Dデータベースを3D構造情報へ変換し,独自の3Dデータベースを作成することが可能である。
3D探索を行なうに当たって特に慎重を要する問題は,どのような3D探索基準を使用するかという問題である。3D探索の方向と成否は,採用した基準の如何によって大きく影響されるからである。既に述べたように,受容体タンパク質の正確な構造データはまだほとんど知られていない。そのため探索基準は通常,生物活性化合物から誘導されたSARに基づいて設定される。
生物活性構造の「訓練集合」は,活性化合物だけではなく不活性化合物も含むのが望ましい。不活性化合物を含めることにより,たとえば活性化合物に対して立体的に許されない領域を明らかにすることができるからである。活性化合物に不可欠な構造条件の定義は,2D分子構造に基づいて簡単に行なうことができる。たとえば,訓練集合に含まれる活性成員がすべて酸性部分構造や水素結合受容基をもつならば,その情報は探索基準に反映されなければならない。しかし,これらの部分構造相互の間の3D配置を定めることはそれほど容易ではない。もちろん訓練集合にかなり硬い活性化合物が含まれている場合には,幾何学的な必要条件を定義することは比較的簡単である。部分構造の3D配向を決める際にそれを鋳型として直接役立てることができるからである。しかし,訓練集合が柔軟な分子しか含まない場合には,配座解析を行なわなければ,共通の薬理作用団を特定することは不可能である(2.5節参照)。
通常,探索クエリーは3点を用いて記述される。薬理作用団の構成原子や部分構造はオブジェクトとして扱われ,対応するオブジェクト間の距離はそれぞれの3D配向を定める。3D探索の処理手順は,どのプログラムでもほとんど同じである。すなわち最初の段階では,おおざっぱな探索が試みられる。そして,必要な化学構造がいかなる相対位置でも欠如し,3D探索基準を満たさない化合物がまず取り除かれる。残った化合物は,次の段階で化学的要求と幾何学的要求の双方を満たしているか否かを吟味される。プログラムは,3D探索の結果として,探索の過程で選び出された化合物をすべて列挙した「ヒットリスト」を最終的に生成する。
3D探索は,効果的に行なわれた場合には,適正な数の活性化合物からなるヒットリストを生成するはずである。しかし,あまり特異性のない探索基準を使用した場合には,化合物の数が多すぎる特大のリストが得られるかもしれない。また,ヒットした化合物のほとんどは,オブジェクトの3D配置が不適当であるために,実際には受容体へ結合することはできないであろう。探索基準の設定に慎重を期すこと,それが妥当な大きさのヒットリストを得るための条件である。
しかし,探索基準が幾何学的にあまりに厳密すぎると,探索に望ましくない制約を課すことになる場合がある。この点を説明するために,次のような状況を想定してみよう。訓練集合のすべての成員に共通する薬理作用団のパターンが,化合物の低エネルギー配座に必ずしも常に対応していないような場合である。このような場合,もし薬理作用団の3D配置を探索クエリーとして使用すれば,データベースに収録されている活性化合物の低エネルギー配座の幾つかは,それらが厳密な幾何学的条件を満たさないために,探索の網に掛かることはないであろう。結果として,活性化合物が見落とされてしまうのである。このような状況に対しては,もっと融通のきく探索手順を適用することが望ましい。
この問題を処理するために幾つかの試みがなされている。化合物の様々な低エネルギー配座をすべて保存することは現実的ではない。2.3.1節で既に述べたように,配座的な柔軟性―異なる可能な配座の数―は回転可能な結合の数と共に増加する。小さなデータベースでも,通常幾十万もの化合物が既に収録されている。個々の化合物に対し,さらに数百ないし数千の低エネルギー配座を保存しようとすれば,膨大なディスク容量が必要となるばかりでない。探索に要する時間を考えると,それは完全に実行不可能である。
この問題に対する非常に有効な一つの解決策は,フレキシブルなオブジェクト間距離を導入することである[10]。すなわち,薬理作用団の構成原子対に対し,距離の最小値と最大値を設定し,両者の範囲内にあればいかなる距離も許容されるとするのである。もし被検分子が十分に柔軟で,指定されたフレキシブルなすべてのオブジェクト間距離を同時に満たすことができるならば,その分子は取り上げられ,ヒットした化合物として保存される。
今日,3D探索プログラムは,潜在的に活性な化合物の包括的なヒットリストを生成するきわめて有効な手段として広く利用されている。プログラムは,配座的な柔軟性の問題にも十分配慮し,しかもできる限り効果的に探索を行なうことのできる3D探索基準を詳細に設定することで更に改善されるであろう。ヒットリストを合理的に評価する手段の開発も現在進められている。
本書の枠組では,いよいよ重要性を増しつつある3D探索の問題について,ごく簡単にしか紹介することができなかった。これは,重要な展開が現在も進行しつつある比較的新しい手法である。この問題については総説が幾つか出ているので,さらに詳しく知りたい読者はそれらをお読みいただきたい[2,11]。