結晶計算法によるアスピリン結晶の配座多形解析

小畑 繁昭, 後藤 仁志


Return

1 はじめに

医薬化合物の結晶形は,薬物の体内動態(ADME)や品質管理に大きな影響を与える.このため,複数の結晶形,すなわち結晶多形が存在する場合,それらを同定し,最適な結晶形に制御することが,医薬品の特許,および臨床製剤や製品化の段階において非常に重要な課題である[1].一般に,結晶多形は結晶中の分子の構造がほぼ同じで充填様式が異なるパッキング多形(packing polymorphism)と,配座が異なることによって充填様式が異なる配座多形(conformational polymorphism)に分類される[2].したがって,結晶多形の可能性を調べるためには,配座異性体,分子配向,空間群,格子定数など,すべての自由度の組み合わせを考慮しなければならず,その探索空間がいかに広大であるかは容易に想像できる.
最近,我々は,結晶計算法による結晶多形解析を目的の一つとして,配座空間探索が可能な分子計算プログラムCONFLEX[3, 4]に,これとは別に開発した結晶計算プログラムKESSHOU[5 - 7]のアルゴリズムを導入した[8, 9].この新しいCONFLEX / KESSHOU [10](以下,「CONFLEX/K」と表記する)を用いれば,広大な結晶多形の探索空間内の,少なくとも配座多形の可能性の一部を,網羅的に探索できる.そこで本研究では,新しい結晶構造が発見され話題になったアセチルサリチル酸 1,すなわちアスピリンの配座多形の可能性を調べるため,CONFLEX/Kを用いて単分子の立体配座解析と結晶多形解析を行う.また,この結晶計算に汎用高精度分子力場として知られるMerck Molecular Force Field(以下,「MMFF94」)[11 - 13]とMM3[14 - 16]を適用し,実験値やab initio計算結果との比較から,単分子の立体配座解析と,配座多形解析における結晶力場の精度について考察する.


Figure 1. Chemical structure of acetylsalicylic acid (Aspirin) 1. Three dihedral angles are determined as follows; t1: C1-C2-O-C8, t2: C2-O-C8=O, t3: C2-C1-C7-O(H)

2 研究背景

1は,1899年にドイツのLeverkusenを本拠地とするBayer社が「アスピリン」[17]という商品名で発売して以来,鎮痛剤,消炎剤,解熱剤として,また近年では抗血小板剤[18, 19]として,世界中で広く利用されているベストセラー薬である.1971年,Vaneの報告により,多くの生理作用を持つプロスタグランジンの生合成を阻害し,薬理作用を発現させることが明らかになった[20].このため,この最も古い合成薬の効能がさらに広がる可能性に,今も注目が集まっている [21, 22].
1のX線結晶構造は,長い間,唯一つの構造[23, 24]だけが知られていたが,その一方で,融点や溶解速度などの重要な物性の差異が実際に観測されている[25, 26].このため,複数の結晶形が存在するかどうかについて,以前から議論されていた[25 - 28].しかし,結晶多形の存在を示す明確な証拠が得られていなかったため,それらの違いは多形によるものではないと考えられるようになっていた[27, 28].ところが,2005年,Vishweshwar等が新しい結晶構造("form II"と命名された)を単離し,解析したことを報告したことから[29],1が少なくとも2種類の結晶多形,すなわち,従来から知られていた結晶構造form I(form IIと区別するため)とform IIの二つの結晶多形を形成できることが明らかとなった.form Iとform IIはどちらも空間群P21/c,z=4であり,分子の立体構造もほぼ一致し,結晶の空間充填様式もよく似ている(Figure 2).しかし,重なり合う二量体層がc軸方向に対して結晶座標で0.5だけずれているため,近接するアセチル基間の相互作用の組み合わせが異なり,それぞれ特徴的な分子間相互作用ネットワークを形成する.また,Bond等の詳細な追試により,form IIは観測できるが,それだけで構成された結晶を単離できないことから,form IIがform Iと相互変換可能な微量成分であることがわかってきた[30, 31].このことは,1がform II以外にも微量な結晶多形を含んだ結晶である可能性を示唆している,と考えることもできる.
結晶計算法による 1の多形予測は,form IIが見つかる以前から行なわれており,いくつかの結晶構造が提案されてきた[32 - 34].その最初の報告として,GavezzottiとFilippiniがform Iとほぼ同じ配座で,空間群P-1の結晶構造を提案している[32].また,Payne等は,ベンゼン環とアセチル基が同一平面上にある平面配座が空間群P21/cの結晶構造を形成する配座多形の存在を示唆した[33].そして,2004年,OuvrardとPriceは,DFT計算(B3LYP/6-31G**)から得られたいくつかの安定配座と複数の空間群から様々な結晶多形を生成し,彼等が開発した独自の多極子結晶力場[35, 36]を用いて結晶エネルギーを計算した[34].その結果,form Iとほぼ同じ配座,同じ空間群で,充填様式が異なる新しい結晶多形を,最も安定な結晶構造として予測した.そして,とても興味深いことに,この結晶構造は,その一年後,前述したVishweshwar等が観測したform IIと酷似していたのである[29].このように,計算化学に基づく予測が最も有名な 1の結晶多形を予言するという大きな成功を収めたことによって,結晶計算を含む分子シミュレーションによる結晶構造解析が,実践的な結晶多形スクリーニングの研究手法の一つとして注目を集めるに至っている.


Figure 2. Molecular packing arrangements and positions in the crystal polymorphs of 1. The neighboring dimers (blue and green) were related by two-fold screw axes along b axes.

3 方法

3. 1 結晶計算法

CONFLEX/Kによる結晶計算は,非対称単位に含まれるすべての分子(以下,「オリジナル分子」)と,空間群で定められている対称操作と結晶の並進対称性にしたがって有効結晶半径Dmaxまで展開した非対称単位(以下,「レプリカ単位」)内の分子(以下,「レプリカ分子」)で構成される球状の結晶を構築する.そして,この球状結晶の結晶エネルギーを次のように評価する:結晶エネルギーEcrystalは,非対称単位内の全相互作用エネルギーEintraと,非対称単位内のすべてのオリジナル分子と,Dmax以内に展開されたレプリカ分子との間に生じる相互作用エネルギーElattice(これを格子エネルギーと呼ぶ)の総和として次式で定義される.

ここで,Nは非対称単位内の分子を構成する総原子数, Mは結晶空間に展開されたレプリカ単位の数を表す.また,Min(RiJ)はオリジナル分子とレプリカ分子の間の最近接原子間距離であり,本研究では,RiJDmax以内にあるレプリカ分子に関わるすべての非結合相互作用をElatticeに含めた(その他に,原子対単位や非対称単位で相互作用を考慮することも可能である).なお,1のX線結晶構造解析において非対称単位内の分子数は1(z' =1)であり,本研究の結晶計算においても同様にz' =1として結晶構造を構築した.また,Dmaxは結晶計算の精度を決定する主要因の一つであり,扱う分子に依存するが,1については既に詳細な検討を行なった[9].ここでは,十分に計算誤差が小さくなることが分かっているDmax=80Aを採用した.この場合,計算対象となる結晶構造におよそ12,000分子が含まれることになる.
CONFLEX/Kでは,3種類の結晶構造最適化法を実行することができる.一つは,分子の構造を固定して空間配置(単位格子と分子配向)のみを緩和する「結晶格子最適化法(LOPT: Lattice Optimization)」である.また,それとは逆に,格子を固定して,レプリカ分子に囲まれたオリジナル分子の構造を緩和する「分子内構造最適化法(MOPT: intra-Molecular Optimization)」がある.さらに,これら二つの最適化法を自己無撞着に至るまで繰り返すことによってすべての結晶パラメータを緩和する「完全結晶構造最適化法 (FOPT: Full Optimization)」は,最も高価な最適化法として実行することができる.ただし,いずれの最適化法を使っても,空間群は維持され,結果としてより高い対称性の空間群に到達する可能性はある.また,これらの最適化法は分子配向の全方位探索をしないため,最適化によって得られた結晶構造は,初期状態に依存した局所解である可能性がある.
本研究で行う配座多形解析は,非対称単位内のオリジナル分子の配座を変えることにより生じる結晶構造の変化を調べることであり,必然的に適用する最適化法は,LOPTかFOPTになる.FOPTは空間群を維持するという制約の下で,結晶エネルギーにおける極小点まで最適化することが可能であるが,計算時間も要することから,多くの結晶多形を評価するには適していない.したがって,本研究では,より高速に最適化結晶構造を求めることが可能なLOPTのみを適用した.この場合,結晶中におけるオリジナル分子の構造変位は起こらないため,Eintraは孤立系の配座エネルギーと同じ値になる.

3. 2 結晶力場

比較的高精度な結晶計算法では,結晶構造を再現するようにパラメータを最適化し,多くの原子が近接することによって誘起されるある種の多体効果を評価することができる,結晶計算専用の力場を適用することが多い.しかし,結晶計算の精度は,計算対象となる結晶の大きさ,積算する相互作用の範囲,あるいはEwald法などの近似計算や階層化・疎視化技術の導入など,結晶をどのようにモデル化するかということにも強く依存している.また,孤立真空中や気相中のように高精度な参照データを入手し難い.したがって,汎用性の高い,高精度な結晶力場の開発は,一般に困難であるとされている.
CONFLEXに導入したKESSHOUは,汎用結晶力場を構築することも目的の一つとして開発してきたが,現時点では専用の力場はなく,現在開発中である.そこで本研究では,主に孤立真空中,あるいは気相中の分子構造やエネルギーを再現するために開発された汎用分子力場を結晶計算に適用した.ここでは,分子内,および分子間ポテンシャルとして,CONFLEX/Kに装備されているMMFF94 とMM3 を適用した.
MMFF94は,非結合相互作用をvdW相互作用項,水素結合相互作用項,および電荷-電荷相互作用項に分けて計算する[12, 13].vdW相互作用と水素結合相互作用のポテンシャル関数に,Halgrenが希ガスの原子間ポテンシャルの精密な解析に基づいて定式化したBuffered-14-7型ポテンシャル関数(Eq. 3)[11]を,また,電荷-電荷相互作用に,クーロン力型ポテンシャル関数(Eq. 4)を採用している[13].

ここでRiJは原子iと原子J間の距離,qは結合電荷増分(Bond-charge increment)法で求められた原子電荷である.また,Dは有効誘電率であり,孤立真空中も結晶中も同じ値1.0を用いた.その他の記号の意味は,原則として,Halgrenの論文の定義に従う[13].MMFF94の水素結合相互作用項は,vdW相互作用項と同じ関数(Eq. 3)を採用しているが,平衡距離パラメータR*に対して0.8,ポテンシャル井戸の深さeに対して0.5のスケーリング係数をかけることにより,平衡距離を短く,ポテンシャル井戸を浅くなるように修正している[13].このように,MMFF94は,水素結合相互作用を原子間距離のみに依存する等方的な相互作用として扱う.
一方,MM3は,非結合相互作用を,vdW相互作用項(Eq. 5)と水素結合相互作用項(Eq. 6),そして,静電相互作用項として,結合双極子法に基づく双極子-双極子相互作用(Eq. 7)の和として定義している.

ここで,Eq. 5のRiJは原子iと原子J間の距離である.Eq. 6のRX-HRYH,およびbの各内部座標の定義をFigure 3 (a) に,Eq. 7のRab,およびcの定義はFigure 3 (b) に示した.Eq. 6とEq. 7におけるDは有効誘電率であり,一般には気相中を想定した値1.5を用いる.本研究では,孤立系(単分子)のMM3計算では有効誘電率に1.5を用いるが,結晶計算では1.0を適用した.その他の記号の意味は,Allinger等の論文の定義に従う[14, 16].MM3計算において,Eq. 6は水素結合対HY間の結合数が4つ以上離れた場合に適用され,3つの場合には,通常のvdW相互作用関数Eq. 5が適用される.このように,MM3は,MMFF94とは異なり,水素結合ポテンシャル関数の引力項にcosine関数を適用し,また,静電相互作用を結合双極子間の相互作用として評価することにより,非結合相互作用の異方性をある程度考慮した多極子力場と考えてもよいだろう.
なお,本研究では,共役p電子をPPP法で評価するMM3計算は行わない.また,1の計算においてMM3(92)の標準パラメータに含まれないポテンシャルパラメータは,Allinger等が提唱するParameter Estimation法を使って生成した[37].


Figure 3. Definition of internal coordinates: (a) hydrogen bonding interaction and (b) dipole-dipole interaction.

3. 3 配座空間探索と配座多形解析

配座多形解析を行なうため,まず,CONFLEXを用いて 1の配座空間探索を行った.この際,網羅する配座空間を決める探索上限値(SEarch Limit, SEL)と,探索過程で保持する配座エネルギー上限値(ESAV)を最安定配座から21 kcal/mol以内とし,また,回転異性体が考えられる全ての側鎖のねじれ角,つまり,t1t3の回転と,カルボキシル基の水酸基(O=C−O-H)の回転に対してStepwise Rotation法による構造変位(Local Perturbation)を適用した.配座空間探索における構造最適化にはMMFF94を採用し,さらに,得られた全ての配座異性体を初期構造としてMM3による構造最適化を行なった.また,同様に,Gaussian03[38]を用いてMP2/6-31+G**による構造最適化も行なった.
配座空間探索で得られた配座異性体から配座多形を構築するため,CONFLEXのGUIであるBARISTA[39]を用いた.非対称単位内における各配座異性体の位置と配向は,X線結晶構造解析で実測されたform Iとform IIの非対称単位を基準として,そのベンゼン環炭素と,これに結合するカルボキシル基炭素とエーテル酸素の原子位置を重ね合せることによって決定した.そして,各配座異性体を導入した非対称単位を用い,実測の空間群と格子定数に基づいて配座多形を生成した.次に,生成した配座多形を初期構造として,前述したLOPTによって結晶構造最適化した.得られた最適化結晶構造に対してPLATON[40]によりreduced cellを計算し,重複する結晶構造の有無を調べた.

4 結果

4. 1 孤立系の立体配座解析

CONFLEXによる配座空間探索の結果,およそ20 kcal/mol以内に10種類の配座異性体を得た.これらを初期構造として,MM3,MP2による構造最適化,基準振動解析を行い,MP2については,すべてエネルギー極小構造であることを確認した.MM3においては,後述するように,一つの配座が構造最適化によって別の配座と同じ構造に収束したため,最終的に9種類のエネルギー極小構造を確認した.各計算法による配座エネルギー(最安定配座からのエネルギー差)をTable 1に,最適化構造の三つのねじれ角t1~t3Figure 1)をTable 2に示した.また,計算法の違いによる最適化構造の変位を比較するため,同じ初期構造に由来する最適化構造に対して,ベンゼン環炭素,およびそのC1に結合するカルボキル基炭素と,C2に結合するエーテル酸素(Figure 1)の位置で重ね合わせた結果をFigure 4に示した.これらの図表において,各配座異性体を容易に区別するため,同じ初期構造に由来する配座異性体はMMFF94の配座エネルギーの低い順に1から10,次式に示すt2t3のねじれ角範囲に従って"S","+","-","A",そして,カルボキシル基がs-cisのときを"c",s-transのとき"t"として,4つの記号で表記した.

配座探索で得られた配座異性体の中で最安定配座と二番目に安定な配座は,いずれの計算法においても同じであり,それぞれ1SAcと2SScであった.2SScはX線構造解析において結晶中で観測される構造に最も近い.1は結晶中でカルボキシル基が向かい合った水素結合二量体を形成することから,結晶中の2SScの二量体が水素交換をすることによって1SAcが生じると考えてもよいだろう.両者のエネルギー差は,MMFF94において1.9 kcal/molと比較的大きいが,MP2で0.7 kcal/mol,MM3では0.3 kcal/molであり,それほど大きな差はないが,孤立真空中(気相中)における最安定配座が,最も安定な結晶構造を形成するとは限らない,という例の一つと言える.
ここで適用した分子力場の特徴は,安定な二つの配座異性体以外に現れている.Table 1より,MMFF94の配座異性体のエネルギー順位は,6SStを除けば,MP2のそれと一致するが,配座エネルギー値は1.2-3倍ほど不安定になっていることがわかる.6SStの配座は,MP2では3番目に安定な構造であり,MMFF94とMP2の最適化構造を比較すると(Table 2 and Figure 4),t3に大きな違いが見られる(MMFF94: 34.30°, MP2: 18.86°).同様な構造の違いは8AStでも観測される.この違いは,Figure 4から明らかなように,カルボキシル基の水酸基がs-transになることで,その水酸基水素と隣接するエーテル酸素との間で分子内水素結合を形成することに起因する.このH...O間距離をMMFF94とMP2の最適化構造で比較すると,6SStではそれぞれ2.089 Aと1.844 A,8AStで2.022 A と1.885 Aとなり,MMFF94はMP2よりもこの間の水素結合を弱く評価していることがわかる.これら以外の分子構造に関しては, MMFF94とMP2の最適化構造はほぼ一致している.
一方,MM3の配座異性体のエネルギー順位は,MP2のそれと部分的に一致せず,配座エネルギー差も小さくなる傾向にある(Table 1).また,多くのMM3の最適化構造において,t1は130°付近,t3がおよそ0°か180°の値になる(Table 2).特に,t3については,MM3のポテンシャルにおいて,カルボキシル基がベンゼン環と同一平面を形成するとより安定になるため,9A+tと10A-tは同じ最適化構造に収束した.これら以外にも,Figure 4より明らかなように,MM3の最適化構造は,他の計算法による最適化構造とは異なる傾向にある.その原因の一つは,本研究で適用した追加パラメータにあり,Allinger等が提唱するParameter Estimation法の問題であるが,今後,それらのパラメータを最適化する必要があると考えている.
CONFLEXが創出した10種類の配座異性体を,Glaserが報告した9種類のエネルギー極小構造(RHF/6-31G*とB3LYP/6-31G*)[41]と比較したところ,8種類の配座がほぼ一致する構造であった.CONFLEXが見逃した一つの異性体は,MP2/6-31+G**においてもエネルギー極小構造であるが,MMFF94で構造最適化を行うと3S-tに,MM3では7SAtに収束した.一方,CONFLEXが創出した配座異性体のうち2つの配座異性体(7SAtと9A+t)はGlaserの報告にはなく,少なくとも我々の知る限りにおいて,本研究で初めて見つかった配座異性体である.

Table 1. Comparison of conformational energies of conformers optimized by MMFF94, MM3 and ab initio methods.
DE(kcal/mol)
LabelaMMFF94MM3MP2c
1SAc0.0000.0000.000
2SSc1.9230.2860.735
3S-t5.5733.5534.092
(S-t)--5.263 d
4AAc7.6212.0875.408
5ASc8.1702.0845.563
6SSt9.8963.6663.668
7SAt10.2183.9086.157
8ASt14.8924.0987.503
9A+t16.7236.23011.158
10A-t20.124- b14.427
a Four characters mean that the first is the conformation number ordered by conformational energies in MMFF94 calculations, the second and third indices (S/A or +/-) are t2 and t3 dihedral angles expressions, and the last small letter (c/t) represents either s-cis or s-trans conformation of carboxylic group. b This conformer is not the energy minimum on the MM3 energy potential. c MP2/6-31+G**. d This is corresponding to one of the conformers reported by Glaser [41]. This is not the energy minimum on the MMFF94 energy potential.


Figure 4. Superimpositions of conformers optimized by MMFF94 (green), MM3 (blue), and MP2/6-31+G** (red) methods.

Table 2. Three dihedral angles of the optimized conformers
Dihedral angles (°) c
Label aMethodt1t2t3
1SAcMMFF9483.095.24-176.16
MM3131.58-7.43-179.71
MP264.618.33-163.95
2SScMMFF9483.812.657.79
MM3132.87-10.360.18
MP263.903.4723.95
3S-tMMFF9479.15-6.90-71.01
MM388.23-22.91-4.54
MP287.93-7.77-58.13
(S-t)MP2 b57.216.09-134.45
4AAcMMFF9499.08168.93158.55
MM3133.38130.79-179.90
MP288.62-175.99148.82
5AScMMFF94100.87174.15-30.09
MM3132.35135.300.24
MP291.86-174.97-35.43
6SStMMFF9499.142.0334.30
MM3125.51-12.48-0.81
MP2105.443.7118.86
7SAtMMFF9497.929.66134.22
MM3128.84-5.55-179.67
MP288.9611.61140.82
8AStMMFF94100.69170.9730.32
MM376.21-144.21-0.18
MP298.15176.4821.95
9A+tMMFF94101.26169.77126.40
MM3133.57128.65179.95
MP293.23-179.19130.99
10A-tMMFF9497.68167.09-131.59
MM3133.57128.65179.95
MP271.76-177.85-139.54
a See the text and also the footnote a in Table 1. b See the footnote d in Table 1. c Definitions are referenced at the footnote in Figure 1.

4. 2 CONFLEX/K結晶計算の適用

CONFLEX/Kによる結晶計算の精度を評価するため,既知のX線結晶構造form I とform IIを構築し,MMFF94とMM3の結晶力場を用いてLOPTによる結晶構造の最適化を行なった.この際に,水素位置をそれぞれの力場パラメータのC-H,O-H平衡結合長に補正した構造を,非対称単位内のオリジナル分子として適用した.LOPTによる最適化後の結晶構造の格子エネルギーと格子定数,単位格子の体積,および既知のX線結晶構造との差をTable 3に示した.また,X線結晶構造と最適化構造を重ね合わせた結果をFigure 5に示した.これらの図表より明らかなように,MM3で最適化した結晶構造はMMFF94のそれよりも, form Iとform IIをよく再現している.また,他の結晶計算法と比較しても,MM3が最も良い結果を示した.このように,MM3の非結合相互作用項は,少なくとも 1の既知の結晶構造を再現するために必要な精度で分子間ポテンシャルを与える.一方,MMFF94は,MM3ほど実験値を再現することはできなかったが,実用上,問題が生じるほど,大きな破綻はなかった.ただし,結晶構造を大きく評価する傾向にある.

Table 3. Comparisons of crystal structures (form I and II) among the experimental and various computational methods.
ObservedLattice ConstantsRMSD
PolymorphElatticeVolume ea / Ab / Ac / Ab / °Average
Method/ kcal mol-1/ A3(D / %)(D / %)(D / %)(D / %)/ %
form I
Expl. (r.t.) a854.211.4306.59111.39595.68
Gavezzotti b-29.40801.611.05 (-3.3)6.61 (0.3)10.98 (-3.6)90.3 (-5.6)3.7
Ouvrard c-25.36869.011.389 (-0.4)6.758 (2.5)11.350 (-0.4)95.86 (0.2)1.3
This work
MMFF94-50.32898.011.751 (2.8)6.460 (-2.0)11.862 (4.1)94.22 (-1.5)2.8
MM3-53.73864.911.465 (0.3)6.703 (1.7)11.346 (-0.4)97.24 (1.6)1.2
form II
Expl. (180 K) d835.812.1526.50611.368111.57
Ouvrard c-25.41866.112.124 (-0.2)6.696 (2.9)11.475 (0.9)111.60 (0.0)1.5
This work
MMFF94-50.58891.112.730 (4.8)6.395 (-1.7)11.952 (5.1)113.67 (1.9)3.7
MM3-54.38855.012.137 (-0.1)6.621 (1.8)11.351 (-0.1)110.39 (-1.1)1.0
a Ref. 24. b Ref. 32. c Ab initio optimized structures were used. See also ref. 34. d Ref. 31. e Volumes were calculated from lattice constants.


Figure 5. Superimpositions of the X-ray crystal structures (red) and the corresponding optimized crystal structures calculated by the lattice optimizations with MMFF94 (green) and MM3 (blue). RMS differences between the observed and the optimized ones were calculated by reference to the positions of the heavy atoms in unit cell. (a) RMSD (form I vs. MMFF94) = 0.185 A, (vs. MM3) = 0.093 A. (b) RMSD (form II vs. MMFF94) = 0.214 A, (vs. MM3) = 0.090 A. These crystal illustrations are generated by Mercury (copyright CCDC 2001-2007, http://www.ccdc.cam.ac.uk/mercury/).

4. 3 配座多形解析

MMFF94とMM3を結晶力場として配座多形を評価するため,それぞれの力場による最適化構造を非対称単位として,実測のform Iとform IIの結晶構造情報から,それぞれ20種類と18種類の配座多形を構築した.そして,これらを初期構造としてLOPTによる結晶構造の最適化を行った.最適化構造における分子回転と単位格子に対するエネルギー勾配の平均二乗差は,それぞれ1×10-3 kcal/mol/rad,1×10-3 kcal/mol/A以下まで到達した.各配座多形のEcrystalを,EintraElatticeのそれぞれの最小値からのエネルギー差(DEintraDElattice)の積算としてFigure 6に示した.また,DEcrystalが15 kcal/mol以内の配座多形の結晶構造情報をTable 4Table 5に示した.これらの図表において,最適化した配座多形を区別するために,非対称単位に導入した配座異性体のラベルに,form Iとform IIをそれぞれ意味する"/I"と"/II"を付加して表記することにした.
Figure 6から明らかなように,どの力場を適用しても,X線で観測された結晶構造に近い構造を示した2SSc/Iと2SSc/IIが最も安定な結晶構造であった.ただし,X線結晶構造解析では,form IIよりもform Iの方がわずかに安定である[31]が,これを再現できたのは,配座の立体構造に関して問題があると思われるMM3だけで,MMFF94ではform IIの方が0.3 kcal/mol優位であった.Ouvrard等の結晶計算[34]においても,同様にform IIの方が0.1 kcal/mol程度安定である(Table 3)ことから,現在,実行可能な結晶計算法ではこのわずかなエネルギー差を再現することは難しいといえる.


Figure 6. Conformational polymorphic energies calculated by lattice optimizations with (a) MMFF94 and (b) MM3. The green rectangles show DEintra, and the blue ones show DElattice.

Table 4. Conformational polymorphic structures calculated by lattice optimization with MMFF94
Energy / kcal mol-1Densitya / Mg m-3Lattice ConstantsEuler Angle Differenceb
PolymorphDEintraDElatticeDEcrystala / Ab / Ac / Ab / °Dq1 / °Dq2 / °Dq3 / °
2SSc / II1.927.050.001.33312.8366.39511.973114.074.14-1.461.26
2SSc / I1.927.350.301.31711.7276.56011.82593.88-2.47-0.11-1.18
1SAc / II0.0012.553.581.27512.5256.49113.108118.331.95-10.680.86
1SAc / I0.0013.094.121.27111.0676.60312.88092.290.056.200.24
6SSt / II9.908.149.071.14912.8069.5649.409115.40-34.72-4.06-30.33
6SSt / I9.908.949.871.13312.6778.9219.42097.7922.41-5.25-14.95
10A-t / I20.120.0011.150.94612.67211.3898.84998.21-20.25-18.36-12.12
10A-t / II20.121.9513.110.96115.61010.5429.276125.4112.010.87-2.10
a Density = z.W/(NA.VUC), where z is the number of molecules in the unit cell, W is the molecular weight of a molecule, NA is the Avogadro's number, and VUC is the volume of the unit cell. b These values were described by z-y-x angle system [42].

Table 5. Conformational polymorphic structures calculated by lattice optimization with MM3
Energy / kcal mol-1Densitya / Mg m-3Lattice ConstantsEuler Angle Differenceb
PolymorphDEintraDElatticeDEcrystala / Ab / Ac / Ab / °Dq1 / °Dq2 / °Dq3 / °
2SSc / I0.290.940.001.43413.3035.41111.59187.431.93-4.8110.82
2SSc / II0.291.200.261.39714.4245.54211.920116.03-3.406.876.61
1SAc / I0.004.303.081.41813.0835.42211.90087.381.77-2.0212.35
1SAc / II0.006.014.791.32514.3655.44312.859116.20-5.391.928.91
9A+t / I6.230.005.011.35810.5708.43110.23474.9773.5728.29-15.04
4AAc / II2.096.547.401.2387.3876.84219.12292.630.46-28.25-4.92
5ASc / II2.089.6610.521.30310.4346.50813.51590.76-2.721.29-3.21
7SAt / II3.919.6812.371.36611.8007.9549.42598.38-7.768.30-11.71
6SSt / II3.6711.7414.191.33611.0898.3349.993104.24-20.90-0.87-17.78
3S-t / II3.5512.2414.561.27811.6277.50612.545121.30-40.34-18.66-52.41
a See the footnote a in Table 4. b See the footnote b in Table 4.

前述したように,結晶中において,水素結合による二量体を形成しているカルボキシル基間で水素交換が起こるとすると(Figure 5),2SSc/Iと2SSc/IIは,それぞれ1SAc/Iと1SAc/IIに変換すると考えられる.これらは,どちらの結晶力場においても比較的安定な結晶構造として評価されたが,1SAc/Iと1SAc/IIのDElatticeは2SSc/Iや2SSc/IIと比べ,MMFF94でそれぞれ3.8と3.6 kcal/mol,MM3で3.1と4.5 kcal/molも高くなる.つまり,孤立真空中で優位であった1SAcは,結晶中では2SScよりも不利な分子間相互作用エネルギーを示す.その最も大きな要因は,1SAc二量体の水素結合距離が,2SScのそれよりもMMFF94で0.03-0.07 A,MM3で0.01-0.03 A離れ,この間の水素結合が弱くなっていることにある.ただし,分子が密に充填した結晶中においては,他にも多くの相互作用が複雑に関係しているため,必ずしも一つの要因でエネルギー差を説明できるわけではない.
その他の配座多形において,カルボキシル基がs-transである配座多形が,比較的安定な結晶構造として評価されたことは興味深い.これらはs-cisの配座多形のように安定な水素結合二量体を形成できないが,隣接する分子との間にChain型の水素結合ネットワークを形成する可能性がある.このため,s-transの配座異性体を用いた配座多形は,LOPTによって結晶中の分子の向きが大きく回転し,その影響は格子定数(格子長と格子角)の大きな変化にも現われている(Tables 4, 5).ここで適用したLOPTは非対称単位内の分子の内部座標をすべて固定する比較的制約の多い結晶構造最適化法ではあるが,Chain型の水素結合ネットワークを形成できるほど,分子配向や格子定数を驚くほど大きく変化させたことになる.その結果として,例えば, MMFF94による10A-t/Iや,MM3の9A+t/I(Figure 7)などは,配座多形の中で最小のElatticeを示した(Tables 4, 5).特に,後者は,Ecrystalにおいて1SAc / IIよりもわずか0.2 kcal/molの差しかない.では,このような配座多形が新しい結晶多形として存在する可能性はあるのだろうか?


Figure 7. Chain motif of hydrogen bonding network appearing in conformational polymorphic structures "9A+t/I" optimized by MM3.

5 考察

新しい結晶多形の存在の可能性を検討する前に,本研究で,結晶力場として適用したMMFF94とMM3の特徴(問題点)について,整理してみる.
孤立系の立体配座解析(前節)において述べたように,MMFF94はMP2に匹敵するほど配座の立体構造を再現することができたが,分子内水素結合の評価に問題があった.おそらく,このことが原因で,既知の結晶構造の再現性は,それほど高くなかった.一方,MM3はX線結晶構造を高い精度で再現することができたが,配座異性体の立体構造はうまく再現することはできなかった.このように,二つの力場は,結晶力場として用いるには一長一短があり,特に,既知の結晶構造とは大きく異なる配座多形を評価する場合,注意が必要である.しかしながら,MM3を分子間ポテンシャルとして適用した場合,少なくとも1の結晶構造の再現性は驚くほど高くなることから,配座異性体の立体構造をうまく再現することさえできれば,本研究で検討してきた他の配座多形についても,より正確に評価できることが期待される.
そこで,MMFF94とMP2の配座異性体から構築した配座多形について,MM3を結晶力場としてLOPTによる最適化を試みた.ただし,配座多形の構築にMMFF94構造を用いた場合には,MMFF94の配座エネルギーをEintra, MM3で評価したElatticeとしてEcrystalを算出した.これをMM3//MMFF94と表記する.MP2構造から結晶多形を構築した場合にも,同様に,Eintra にMP2の配座エネルギーを適用してEcrystalを算出し,これをMM3//MP2と表記する.Figure 6と同様に,DEintraDElatticeを積算した,EcrystalFigure 8に示した.
Figure 8から明らかなように,4つの安定な配座多形以外は,すべて10 kcal/mol以上も不安定な結晶構造であった.気相中や溶液中の配座変換や化学反応の熱力学を考えれば,10 kcal/molのエネルギー差は存在比としてほぼ0 %であることから,配座多形として存在する可能性は,極めて低いと思われる.ただし,その存在を考える上で,気相中や溶液中では,反応の遷移状態を比較的容易に超えられることが分布則を満たす条件であり,これが結晶中にも当てはまるとは限らないことに留意する必要があるだろう.仮に,結晶多形間の相転移に必要な活性化エネルギーが,気相中や溶液中のそれよりも大きいとすると,その相転移は一般的な条件では起こり難いが,一度起こってしまうと戻り難いはずである.残念ながら,現時点で我々は,結晶多形の存在と結晶エネルギーの関係を決定する実験,および計算結果を持っていないため,どの程度の結晶エネルギー差のある配座多形までが存在できるかどうかを判断することはできないが,今後,広範囲に調べる必要があると考えている.
以上の結果から,本研究で調べた配座多形が存在する可能性についてまとめる.まず,どの結晶力場を適用しても,圧倒的に優位な結晶構造である2SSc/Iと2SSc/IIは,それぞれ,既知の結晶構造form Iとform IIとほぼ一致していた.つまり,これらは実在する.また,これらの二量体間で水素交換が起こるとすれば,1SAc/Iと1SAc/IIもまた,微量成分として存在すると考えてよいだろう.ただし,実験でform Iとform IIを単離できないように,2SScと1SAcで構成された結晶多形を単離することはできないと思われる.それ以外の,本研究で検討した配座多形は,その存在の可能性は極めて低いと思われる.
最後に,MM3//MMFF94とMM3//MP2の結晶エネルギーの評価がよく似ていることは(Figure 8),他の化合物の配座多形解析への応用や,その評価の高速化を考えた場合,とても重要である.しかし,MMFF94の配座エネルギーにおける過大評価や分子の立体構造のわずかな違いによる結晶エネルギーの変化は,明らかに存在している.こうした分子の立体構造を変化させないLOPTによる配座多形解析では,当然のことながら,配座エネルギーと立体構造をいかに再現できる分子ポテンシャルを用意できるかが,より正確な評価を期待するためには重要であるといえる.


Figure 8. Conformational polymorphic energies calculated by lattice optimizations with (a) MM3//MMFF94 and (b) MM3//MP2.

6 まとめ

本研究では,1の配座多形の可能性を調べるため, CONFLEX/Kを用いて,立体配座解析と配座多形解析を行った.徹底的な配座探索の結果,10種類の配座異性体を見つけ,この中には,これまで報告のなかった新しい2つの配座異性体の存在を確認した.得られた配座異性体をMM3,およびab initio法(MP2/6-31+G**)で構造最適化を行ったところ,MMFF94構造はMP2構造とよく一致するが,MM3構造は側鎖のねじれ角が異なっていた.配座エネルギーについて,MMFF94はMP2よりも高く,MM3は低く評価する傾向にあった.一方,既知のX線結晶構造をLOPTで最適化したところ,MM3が極めて高い精度で結晶構造を再現できることがわかった.配座多形解析において,既知のX線結晶構造に近い2SSc/Iと2SSc/II,および,その水素結合二量体の水素交換で生じる1SAc/Iと1SAc/IIが安定な結晶構造として評価され,それ以外の配座は10 kcal/mol以上も不安定であると評価された.
本研究では,LOPTを使って既知のX線結晶構造の情報に基づいた極めて限定的な配座多形解析を行ったが,より広範囲な結晶多形解析を行うためには,より高速で網羅的なスクリーニング技術が必要になる.そこで注目される指標の一つとして,結晶密度があげられる.一般に,分子性結晶の構造は,より高密度に充填することによってより安定な結晶を形成する傾向にあると考えてもよいだろう.Table 4Table 5に示したように,各配座多形の密度と格子エネルギーを比較すると,高密度な配座多形の格子エネルギーが低くなる傾向にあり,両者は相関関係にあることに気づく.このことは,今後,配座多形だけでなく,空間群のスクリーニングを含む結晶多形解析の高速評価法を開発する上で重要である.
LOPT以外の最適化法が興味深い情報を与えることについても,ここで触れておく.例えば,空間群の制約以外にすべての自由度を最適化するFOPTを適用すると,基準振動解析を行って熱力学諸関数(自由エネルギー)を求めることができる.また,振動数と振動モードを算出し,例えば,相転移をうながす結晶振動を明らかにすることができる.一方,MOPTは,結晶場が分子構造に与える影響を定量的に評価し,分子構造を緩和する.このため,X線結晶構造解析において,解析精度を向上させるために利用することも可能である.いずれの最適化法も,目的に合わせて正しく選択することによって,新しい知見が得られるだろう.これらの研究成果についても,近いうちに報告する予定である.

7 Supporting Information

Table 6. Additional parameter of angle bending for MM3
Atom Types
IJKkqq0
375500.695110.000

Table 7. Additional parameters of torsion for MM3
Atom Types
IJKLV1V2V3
1375500.0002.470-0.600
50375240.0002.470-0.600
78375500.0002.470-0.600
35050750.00011.6000.000

本研究を進めるにあたり,アスピリンに関する有意なご助言をいただいた徳島大学大学院ヘルスバイオサイエンス研究部の山内あい子准教授に感謝いたします.本研究の一部は,日本学術振興会科学研究費補助金の援助を受けて行われました.また,本研究で行った結晶計算の一部は,東京工業大学学術国際情報センターの西川武志准教授の御協力のもと,同センターのスーパーコンピュータTSUBAMEを用いて行いました.すべての関係者に対してここに感謝の意を表します.

参考文献

[ 1] I. Amato, Chem. Eng. News, 85, 27-28 (2007).
[ 2] 佐藤清隆 著, 松岡正邦 監修, 結晶と多形-序論, 結晶多形の最新技術と応用展開, シーエムシー出版, 東京 (2005), pp. 3-17.
[ 3] H. Goto, E. Osawa, J. Am. Chem. Soc., 111, 8950-8951 (1989).
[ 4] H. Goto, E. Osawa, J. Chem. Soc., Perkin Trans., 2, 187-198 (1993).
[ 5] E. Osawa, P. M. Ivanov, H. Goto, M. Yamamoto, J. Ruzinski, A. Aoki, BS/KESSHOU, JCPE, P111.
[ 6] E. Osawa, H. Goto, T. Sugiki, K, Imai, Study of Intermolecular Interactions using Crystal Structure Database as Reference: A Preliminary Report on the Adjustment of van der Waals Constants, Intermolecular Interactions, ed. by Gans and Boeyens, Plenum Press, New York (1998), pp. 121-134.
[ 7] 大澤映二, 後藤仁志 著, 松岡正邦 監修, 結晶多形の予測アルゴリズム, 結晶多形の最新技術と応用展開, シーエムシー出版, 東京 (2005), pp. 118-135.
[ 8] K. Ohta, H. Goto, Mol. Sci., 1, NP001 (2007).
[ 9] S. Obata, H. Goto, J. Comput. Aided Chem., 9, 8-16 (2008).
[10] H. Goto, S. Obata, T. Kamakura, N. Nakayama, K. Ohta, CONFLEX version 6.2, Conflex Corp., Tokyo (2007).
[11] T. A. Halgren, J. Am. Chem. Soc., 114, 7827-7843 (1992).
[12] T. A. Halgren, J. Comput. Chem., 17, 490-519 (1996).
[13] T. A. Halgren, J. Comput. Chem., 17, 520-552 (1996).
[14] N. L. Allinger, Y. H. Yuh, J.-H. Lii, J. Am. Chem. Soc., 111, 8551-8566 (1989).
[15] N. L. Allinger, Z. -q. S. Zhu, K. Chen, J. Am. Chem. Soc., 114, 6120-6133 (1992).
[16] J.-H. Lii, N. L. Allinger, J. Phys. Org. Chem., 7, 591-609 (1994).
[17] ASPIRIN, Bayer, http://www.aspirin.de/.
[18] D. L. Bhatt, E. J. Topol, Nat. Rev. Drug Discov., 2, 15-28 (2003).
[19] S. P. Jackson, S. M. Schoenwaelder, Nat. Rev. Drug Discov., 2, 1-15 (2003).
[20] J. R. Vane, Nat. New. Biol., 231, 232-235 (1971).
[21] A. T. Chan, E. L. Giovannucci, J. A. Meyerhardt, E. S. Schernhammer, G. C. Curhan, C. S. Fuchs, J. Am. Med. Assoc., 294, 914-923 (2005).
[22] A. T. Chan, S. Ogino, C. S. Fuchs, N. Engl. J. Med., 356, 2131-2142 (2007).
[23] P. J. Wheatley, J. Chem. Soc., 6036-6048 (1964).
[24] Y. Kim, K. Machida, T. Taga, K. Osaki, Chem. Pharm. Bull., 33, 2641-2647 (1985).
[25] A. G. Mitchell, D. J. Saville, J. Pharm. Pharmac., 21, 28-34 (1969).
[26] M. de Bisschop, J. Pharm. Belg., 25, 330 (1970).
[27] R. R. Pfeiffer, J. Pharm. Pharmac., 23, 75-76 (1971).
[28] A. G. Mitchell, B. L. Milaire, D. J. Saville, R. V. Griffiths, J. Pharm. Pharmac., 23, 534-535 (1971).
[29] P. Vishweshwar, J. A. McMahon, M. Oliveira, M. L. Peterson, M. J. Zaworotko, J. Am. Chem. Soc., 127, 16802-16803 (2005).
[30] A. D. Bond, R. Boese, G. R. Desiraju, Angew. Chem. Int. Ed, 46, 615-617 (2007).
[31] A. D. Bond, R. Boese, G. R. Desiraju, Angew. Chem. Int. Ed, 46, 618-622 (2007).
[32] A. Gavezzotti, G. Filippini, J. Am. Chem. Soc., 117, 12299-12305 (1995).
[33] R. S. Payne, R. C. Rowe, R. J. Roberts, M. H. Charlton, R. Docherty, J. Comput. Chem., 20, 262-273 (1999).
[34] C. Ouvrard, S. L. Price, Cryst. Growth Des., 4, 1119-1127 (2004).
[35] D. J. Willock, S. L. Price, J. Comput. Chem., 16, 628-647 (1995).
[36] D. S. Coombes, S. L. Price, D. J. Willock, M. Leslie, J. Phys. Chem., 100, 7352-7360 (1996).
[37] N. L. Allinger, X. Zhou, J. Bergsma, J. Mol. Struct. (THEOCHEM), 312, 69-83 (1994).
Additional parameters used in this work are available from the supporting information.
[38] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery, Jr., T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez, and J. A. Pople, Gaussian 03 Revision D.01, Gaussian, Inc., Wallingford CT (2004).
[39] S. Sugatsuki, Y. Takata, T. Shinohara, T. Fujita, H. Goto, K. Ohta, BARISTA version 1.2, Conflex Corp., Tokyo (2007).
[40] A. L. Spek, PLATON: A Multipurpose Crystallographic Tool, Utrecht University, Utrecht, Netherlands (2008).
[41] R. Glaser, J. Org. Chem., 66, 771-779 (2001).
[42] M. G. Rossmann, D. M. Blow, Acta Crystal., 15, 24-31 (1962).


Return