暗黙の溶媒和

暗黙的溶媒和（連続溶媒和と呼ばれることもある）は、溶媒を個々の「明示的な」溶媒分子ではなく連続媒体として表す方法であり、分子動力学シミュレーションやその他の分子力学の応用で最もよく使用されます。この方法は、タンパク質、DNA、RNA、多糖類の折り畳みや立体配座転移、生体高分子とリガンドの会合、生体膜を介した薬物の輸送など、構造的および化学的プロセスにおける溶質 - 溶媒相互作用の自由エネルギーを推定するためによく適用されます。

暗黙の溶媒和モデルは液体では正当化され、平均力のポテンシャルを適用して、多くの高度に動的な溶媒分子の平均的挙動を近似することができます。しかし、生体膜やタンパク質の界面や内部も、特定の溶媒和特性や誘電特性を持つ媒体と見なすことができます。これらの媒体は、脂質二重層の「極性プロファイル」などのさまざまな解析関数によって特性を記述できるため、必ずしも均一ではありません。^[1]

暗黙的溶媒法には、歴史的に最初のモデルであるアクセス可能な表面積(ASA) に基づくモデルと、より最近の連続体静電気モデルという 2 つの基本的なタイプがありますが、さまざまな方法のさまざまな修正と組み合わせが可能です。アクセス可能な表面積 (ASA) 法は、移動のギブス自由エネルギーと溶質分子の表面積との間の実験的な線形関係に基づいています。^[2]この方法は、自由エネルギーのエンタルピー成分のみを含む分子力学法や静電気法とは異なり、溶媒和の自由エネルギーを直接操作します。溶媒の連続体表現により、計算速度も大幅に向上し、溶媒の立体配座の不完全なサンプリングから生じる統計的平均化の誤差も減少します。 ^[3]そのため、暗黙的溶媒と明示的溶媒で得られるエネルギーランドスケープは異なります。^{[4]暗黙的溶媒モデルは生体分子のシミュレーションには便利ですが、これは近似法であり、パラメーター化と}イオン化効果の処理に関連する特定の制限と問題があります。

最も単純な ASA ベースの方法における 溶質分子の溶媒和自由エネルギーは次のように表されます。

${\displaystyle \Delta G_{\mathrm {solv} }=\sum _{i}\sigma _{i}\ AS​​A_{i}}$

ここで、 は原子i のアクセス可能な表面積、 は 原子iの溶媒和パラメータ、すなわち単位表面積あたりの特定の原子 i の溶媒和の自由エネルギーへの寄与である。異なるタイプの原子（炭素（C）、窒素（N）、酸素（O）、硫黄（S）など）に必要な溶媒和パラメータは、通常、一連の有機化合物の移動自由エネルギーの計算値と実験値の最小二乗近似によって決定される。実験値のエネルギーは、溶質の標準モル濃度を使用して、異なる溶液または媒体間のこれらの化合物の分配係数から決定される。 ^{[5] [6]}${\displaystyle ASA_{i}}$${\displaystyle \sigma_{i}}$

特に、溶媒和エネルギーは、溶質分子を溶媒から真空（気相）に移動させるために必要な自由エネルギーです。このエネルギーは、分子力学で計算された真空中の分子内エネルギーを補うことができます。したがって、必要な原子溶媒和パラメータは、最初は水-ガス分配データから導出されました。^{[7]しかし、タンパク質と}脂質二重層の誘電特性は、真空よりも非極性溶媒の誘電特性に非常に似ています。そのため、より新しいパラメータは、オクタノール-水分配係数^[8]または他の同様のデータから導出されています。このようなパラメータは、実際には2つの凝縮媒体間の移動エネルギー、または2つの溶媒和エネルギーの差を表します。

ポアソン-ボルツマン方程式(PB) は、イオンを含む溶媒中の溶質の静電環境を表します。 cgs単位 では次のように表すことができます。

${\displaystyle {\vec {\nabla }}\cdot \left[\epsilon ({\vec {r}}){\vec {\nabla }}\Psi ({\vec {r}})\right]= -4\pi \rho ^{f}({\vec {r}})-4\pi \sum _{i}c_{i}^{\infty }z_{i}q\lambda ({\vec {r}})e^{\frac {-z_{i}q\Psi ({\vec {r}})}{kT}}}$

または（mksの場合）：

${\displaystyle {\vec {\nabla }}\cdot \left[\epsilon ({\vec {r}}){\vec {\nabla }}\Psi ({\vec {r}})\right]= -\rho ^{f}({\vec {r}})-\sum _{i}c_{i}^{\infty }z_{i}q\lambda ({\vec {r}})e^{\frac {-z_{i}q\Psi ({\vec {r}})}{kT}}}$

ここで、は位置依存の誘電率、は静電ポテンシャル、は溶質の電荷密度、は溶質から無限遠にあるイオンiの濃度、はイオンの価数、qは陽子の電荷、kはボルツマン定数、Tは温度、は溶液中のイオンに対する位置rの位置依存のアクセス可能性の係数（多くの場合、一様に1に設定されます）。ポテンシャルが大きくない場合は、方程式を線形化してより効率的に解くことができます。 ^[9]${\displaystyle \epsilon ({\vec {r}})}$${\displaystyle \Psi ({\vec {r}})}$${\displaystyle \rho^{f}({\vec{r}})}$${\displaystyle c_{i}^{\infty }}$${\displaystyle z_{i}}$${\displaystyle \lambda ({\vec {r}})}$

この方程式は理論的に確固とした根拠があるが、近似なしで計算するには計算コストがかかる。様々な汎用性と効率性を持つポアソン・ボルツマン方程式の数値ソルバーが開発されており、^{[10] [11] [12]}特殊なコンピュータハードウェアプラットフォームを備えたアプリケーションも含まれている。^[13]しかし、PBソルバーの性能は、より一般的に使用されている一般化ボルン近似の性能にはまだ及ばない。^[14]

一般化ボルン(GB) モデルは、正確な (線形化された) ポアソン-ボルツマン方程式の近似です。これは、内部の誘電率が外部の溶媒と異なる球の集合として溶質をモデル化することに基づいています。このモデルの関数形式は次のとおりです。

${\displaystyle G_{s}=-{\frac {1}{8\pi \epsilon _{0}}}\left(1-{\frac {1}{\epsilon }}\right)\sum _{i,j}^{N}{\frac {q_{i}q_{j}}{f_{GB}}}}$

どこ

${\displaystyle f_{GB}={\sqrt {r_{ij}^{2}+a_{ij}^{2}e^{-D}}}}$

そして ${\displaystyle D=\left({\frac {r_{ij}}{2a_{ij}}}\right)^{2},a_{ij}={\sqrt {a_{i}a_{j}}}}$

ここで、は自由空間の誘電率、はモデル化される溶媒の誘電率、は粒子iの静電荷、は粒子iとjの間の距離、 は有効ボルン半径と呼ばれる量（長さの次元を持つ）である。^[15]原子の有効ボルン半径は、溶質内での原子の埋没度を特徴付けるものであり、定性的には原子から分子表面までの距離と考えることができる。有効ボルン半径の正確な推定は、GBモデルにとって重要である。^[16]${\displaystyle \epsilon_{0}}$${\displaystyle \epsilon }$${\displaystyle q_{i}}$${\displaystyle r_{ij}}$${\displaystyle a_{i}}$

一般化ボルン (GB) モデルに疎水性溶媒アクセス可能表面積 (SA) 項を加えたものが GBSA です。これは最も一般的に使用される暗黙の溶媒モデルの組み合わせの 1 つです。このモデルを分子力学の文脈で使用することを MM/GBSA と呼びます。この定式化は、明確に定義された三次構造を持つ短いペプチドのネイティブ状態をうまく特定できることが示されていますが、^[17]他の研究で GBSA モデルによって生成されたコンフォメーション アンサンブルは、明示的な溶媒によって生成されたものと大きく異なり、タンパク質のネイティブ状態を特定しません。^[4]特に、塩橋は、おそらく静電遮蔽が不十分なために過剰に安定化され、ネイティブよりも高いアルファ ヘリックス集団が観察されました。膜の静電環境を近似するために GB モデルのバリアントも開発されており、膜貫通タンパク質の膜貫通ヘリックスの折り畳みにある程度成功しています。^[18]

もう一つの可能​​性は、アドホックな迅速な戦略を使用して溶媒和自由エネルギーを推定することです。高速暗黙溶媒の第一世代は、原子あたりの溶媒アクセス可能な表面積の計算に基づいています。原子タイプのグループごとに、異なるパラメーターが溶媒和への寄与をスケールします（上記の「ASAベースのモデル」）。^[19]

CHARMM 19力場には別の戦略が実装されており、EEF1と呼ばれています。 ^[20] EEF1はガウス形の溶媒排除に基づいています。溶媒和自由エネルギーは

${\displaystyle \Delta G_{i}^{solv}=\Delta G_{i}^{ref}-\sum _{j}\int _{Vj}f_{i}(r)dr}$

iの基準溶媒和自由エネルギーは、グループ i が本質的に完全に溶媒にさらされている適切に選択された小分子に対応する。積分はグループjの体積V _jにわたって行われ、合計はi の周りのすべてのグループjにわたって行われる。EEF1 はさらに距離依存（非定数）誘電体を使用し、タンパク質のイオン性側鎖は単純に中和される。これは真空シミュレーションより 50% 遅いだけである。このモデルは後に疎水効果で拡張され、Charmm19/SASA と呼ばれた。^[21]

溶質の周りに水分子の層または球を含め、暗黙の溶媒でバルクをモデル化することが可能。このようなアプローチは、MJ Frischと同僚^[22]および他の著者 [23] によって提案されている。^{[ 24} ]たとえば、文献^[23]では、バルク溶媒は、一般化ボルンアプローチと、クーロン粒子のペアワイズ相互作用に使用されるマルチグリッド法を使用してモデル化されている。これは、静電計算の粒子メッシュエワルド総和(PME) 法による完全な明示的溶媒シミュレーションよりも高速であると報告されている。溶媒和に関する情報にアクセスして取得できるさまざまなハイブリッド方法が利用可能である。^[25]

PB や GB のようなモデルでは平均静電自由エネルギーを推定できますが、水分子や溶媒分子の組織化に対する溶質による制約から生じる(主に)エントロピー効果は考慮されません。これは疎水効果と呼ばれ、疎水性コアを持つ球状タンパク質の折り畳みプロセスにおける主要な要因です。暗黙的溶媒和モデルは、疎水効果を考慮する項で拡張される場合があります。これを行う最も一般的な方法は、溶媒アクセス可能表面積 (SASA) を疎水効果の程度の代理として採用することです。ほとんどの著者は、この効果の範囲を 5 ～ 45 cal/(Å ² mol) としています。^[26]この表面積は溶質に関係しますが、生理学的温度では疎水効果は本質的にエントロピー的であり、溶媒側で発生することに注意してください。

PB、GB、SASA などの暗黙的溶媒モデルには、水分子がランダムに衝突してファンデルワールス反発により溶質の動きを妨げることで生じる粘性がありません。多くの場合、これは構成と位相空間のサンプリングがはるかに高速になるため望ましいことです。この高速化は、明示的溶媒と比較して達成される CPU 高速化に加えて、シミュレートされた時間単位あたりに訪問される構成が増えることを意味します。ただし、速度論が関心事である場合は、誤解を招く結果につながる可能性があります。

ハミルトン力学の代わりにランジュバン力学を使用し、特定の溶媒に適切な減衰定数を選択することで、粘性を追加することができます。^[27]実際の二分子シミュレーションでは、衝突頻度を大幅に下げることで、コンフォメーション探索を大幅に（場合によっては100倍まで）高速化できます。^[28]また、最近では、溶媒を介した運動量移動と関連する熱変動を考慮するために、変動流体力学に基づくサーモスタットを開発する研究も行われています。 ^[29] ただし、タンパク質の折り畳み速度は、すべての領域で粘性に線形に依存するわけではないことに留意する必要があります。^[30]${\displaystyle \gamma}$

第一溶媒和殻の溶質-溶媒水素結合は有機分子、特にイオンの溶解度に重要である。その平均エネルギー寄与は暗黙の溶媒モデルで再現できる。^{[31] [32]}

すべての暗黙の溶媒和モデルは、溶質の非極性原子は一緒に集まったり非極性媒体を占有したりする傾向があるのに対し、溶質の極性および荷電グループは水中に留まる傾向があるという単純な考えに基づいています。ただし、異なる種類の原子からの反対のエネルギー寄与を適切にバランスさせることが重要です。長年にわたっていくつかの重要な点が議論され、調査されてきました。

湿潤した1-オクタノール溶液には約2Mの水が含まれているため、タンパク質や生体膜の近似としては不十分であり、シクロヘキサンの方がはるかに適した近似であることが指摘されている。^[33] 脂質二重層を横切るさまざまな化合物の受動透過性障壁の調査から、1,9-デカジエンは二重層内部の良い近似値となり得るという結論に至ったが、[ ^34] 1-オクタノールは近似値が非常に低かった。^[35]タンパク質工学データ からタンパク質内部に導出された溶媒和パラメータのセットもオクタノールスケールとは異なっており、非極性原子の場合はシクロヘキサンスケールに近かったが、極性原子の場合はシクロヘキサンスケールとオクタノールスケールの中間であった。^[36] したがって、タンパク質の折り畳みとタンパク質膜結合のモデリングには、異なる原子溶媒和パラメータを適用する必要がある。この問題は依然として議論の余地がある。この方法の元々のアイデアは、有機分子の実験的分配係数からすべての溶媒和パラメータを直接導き、溶媒和自由エネルギーを計算できるようにするというものでした。しかし、最近開発された静電モデルの中には、すべての種類の原子に対して20 または 40 cal/(Å ² mol)というアドホックな値を使用するものもあります。このようなモデルでは、極性原子の存在しない「疎水性」相互作用は、大きな静電エネルギーペナルティによって無視されます。

厳密に言えば、ASA ベースのモデルは溶媒和、つまり液体または均一な媒体間の移動のエネルギーを記述するためにのみ適用されるべきです。固体状態でのファンデルワールス相互作用エネルギーを表面エネルギー単位で表すことが可能です。これは、タンパク質工学やリガンド結合エネルギーを解釈するために時々行われ、^[37]脂肪族炭素の「溶媒和」パラメータは約 40 cal/(Å ² mol)となり、^[38]これは水から液体炭化水素への移動で得られる約 20 cal/(Å ² mol) の 2 倍の大きさです。なぜなら、このようなフィッティングによって導出されるパラメータは、疎水性エネルギー (つまり 20 cal/Å ^{2 mol) と固体状態の脂肪族基のファンデルワールス力のエネルギーの合計を表し、}アルカンの融解エンタルピーに対応するためです。^[36] 残念ながら、簡略化されたASAベースのモデルでは、タンパク質構造や分子結晶において類似の極性を持つ原子のクラスター化の原因となる固体状態の異なる種類の原子間の「特定の」距離依存相互作用を捉えることができません。このような原子間相互作用のパラメータは、タンパク質内部の原子溶媒和パラメータとともに、タンパク質工学データから近似的に導出されています。^[36] 溶媒分子がタンパク質の結合空洞と強く結びつき、タンパク質と溶媒分子が連続した固体を形成すると、暗黙の溶媒和モデルは破綻します。^[39] 一方、このモデルは、水から流動性脂質二重層への移動を説明するためにうまく適用できます。^[40]

さまざまな暗黙的溶媒和モデルとパラメータセットの性能を評価するには、さらなるテストが必要です。これらは、疎水性および両親媒性アルファヘリックス（α）などの非常に単純な構造を持つ少数の分子に対してのみテストされることがよくあります。この方法は、数百のタンパク質構造に対してテストされることはほとんどありませんでした。^[40]

荷電基のイオン化は、暗黙的溶媒和の連続静電モデルや標準的な分子力学および分子動力学では無視されてきた。水から誘電率が約 3 (脂質二重層) または 4 ～ 10 (タンパク質内部) の非極性媒体へのイオンの移動には、ボルン方程式および実験からわかるように、かなりのエネルギーがかかる。しかし、荷電タンパク質残基はイオン化可能であるため、非極性環境では電荷を失うだけで、中性pHではエネルギーコストは比較的小さく、ヘンダーソン・ハッセルバルヒ方程式ΔG = 2.3RT (pH - pK)によれば、Asp、Glu、Lys、および Argアミノ酸残基の場合、約 4 ～ 7 kcal/mol である。このようなイオン化効果のエネルギーコストが低いことは、イオン化残基が埋め込まれたタンパク質変異体で実際に観察されている。^[41]および膜中の疎水性αヘリックスペプチドには、中央にイオン化可能な残基が1つあります。^[42]しかし、PB、GB、GBSAなどのすべての静電的方法は、イオン化可能な基が非極性環境で帯電したままであると想定しているため、静電エネルギーが大幅に過大評価されます。最も簡単に利用できる表面積ベースのモデルでは、この問題は、荷電原子のさまざまな溶媒和パラメータを使用するか、いくつかの修正を加えたHenderson-Hasselbalch方程式を使用して処理されました。^[40] しかし、後者のアプローチでも問題は解決されません。荷電残基は、分子内イオン対と水素結合に関与している場合、非極性環境でも帯電したままになる可能性があります。したがって、Henderson-Hasselbalch方程式を使用しても、エネルギーペナルティが過大評価される可能性があります。このようなイオン化効果を説明するより厳密な理論的方法が開発されており、^[43]このような方法を暗黙の溶媒和モデルに組み込むための取り組みが進行中です。^[44]