Prediction of the Multiple Binding Poses in Classical System Using Metadynamics

Abstract

We have performed metadynamics to explore the ligand docking poses with an RNA aptamer in three collective variables (CVs) space, one distance and two dihedral angles. We showed that the free energy surface (FES) of the ligand binding has several local minima. Furthermore, we also demonstrate that each metastable structure can be deduced from CVs of each free energy minimum.

1 目的

分子シミュレーションによって標的分子(タンパク質，核酸など)に対する低分子の結合活性(結合標準自由エネルギー)を定量的に予測することは，基礎科学的側面のみならず創薬などへの応用面からも興味深い．特に，近年，計算機システムの進歩によって大規模な並列計算が可能になったため，計算例の報告が増加している．自由エネルギーの計算では標的分子と低分子の結合構造が必要である．しかし，自由度が多くなるに従い，その結合構造を予測することは非常に困難となる．したがって，実験によって結合構造が決定されている系を計算対象にすることが多かった．ここでは，高精度な結合構造予測の可能性を調べるために，flexibleな標的分子とrigidな低分子の系に対してメタダイナミクスを応用してその結合ポーズの探索を試みた．

2 方法

一般に，標的分子の結合サイトにおける自由エネルギー面は極めて複雑である．そのため，単純な熱平衡シミュレーションでは準安定な状態に落ち込んでしまい，そこから抜け出せなくなってしまう．ここでは，このような状況を克服する方法としてメタダイナミクス [1]を用いる．メタダイナミクスでは，一般に次式のようなガウス型のバイアスポテンシャルを堆積していくことで，準安定な状態からの脱出を図る．集団変数をs，その時間発展をσ(t)として，次のガウス型ポテンシャルを用いることが多い．   

$V(s,t)=\sum _{t<t\text{'}}W(t\text{'})\exp \left[\sum _i\frac{{(s_i-{\sigma }_i(t\text{'}\left)\right)}^2}{2{\delta }_i^2}\right]$(1)

ここで，Wはガウシアンの高さ，δは標準偏差である．well-tempered ensembleではWを初期値W₀から(2)式に従ってリスケールしていく．系の温度Tに対して仮想的な温度T+ΔTで集団変数をサンプリングすることから(T+ΔT)/Tはバイアスファクターと呼ばれ，バイアスのない場合が1で，∞の場合が通常のメタダイナミクスとなる．   

$W\left(t\right)=W_0{e}^{-\frac{V\left(\sigma \right(t),t)}{\Delta T}}$(2)

テオフィリン分子/RNAアプタマー複合体構造のモデル化(Figure 1)には1O15 (PDB) [2]を用いた．RNAの力場はff14SB [3]を使い，テオフィリン分子の記述には点電荷にRESP [4]，力場にGAFF [5]を用いた．ユニットセル中に3個のMg²⁺イオンをGoudaら [6]のモデルに従って配置し，セルを中性化するように26個のNa⁺イオンを導入した．メタダイナミクスの集団変数はFigure 2のように選んだ．C₁，C₂，C₃，C₄はそれぞれテオフィリン分子の重心，{A7, C8, G26}の重心，{C22, U24}の重心，およびテオフィリン分子の窒素原子TEP (N)で，1つの距離dと2つの二面角ϕ₁，ϕ₂で張る空間を探索した．NVTの条件下で2 fs刻みでサンプリングを行い，1 psごとにバイアスポテンシャル(W₀ = 0.29 kcal/mol，δ = 0.2 A，0.35 rad，0.35 rad)を堆積し，ほぼ収束する800 nsまで実行した．なお，バイアスファクター は5とした．実行中にRNAの構造が破壊されることを防ぐために，テオフィリン分子の周りでRNAの構造安定化に寄与しているbase tripleのRMSD上限がテオフィリン分子の存在しない場合の熱揺らぎ程度になるように調和ポテンシャルで拘束を行った．また，テオフィリン分子が平面的であることから180度回転したものも考えられるが，本報告では実験で決定された結合ポーズ以外に着目するために，PDBの向きのものだけに限定した．全てのシミュレーションは，分子動力学エンジンにGROMACS-4.6.7，メタダイナミクスのプラグインPLUMED-2.1.3 [7]を用いて行った．

Figure 1.

 Theophylline/RNA complex. A red circle depicts the binding site.

Figure 2.

 Collective variables {d,ϕ₁,ϕ₂}.

3 結果と考察

シミュレーションの結果，最安定構造から自由エネルギー差2.5 kcal/mol以内に4個のlocal minimumが発見された．特に，d = 8Å近傍には最安定状態を含む3個の極小状態が存在することが分かった(Figure 3)．ここで，図の等高線の間隔は2 kcal/molである．

Figure 3.

 Free energy surfaces for theophylline binding to an RNA aptamer at d= 8 A.

この自由エネルギー面でlocal minimumを与える座標から，それぞれの結合構造を抽出する．ここでは，各準安定構造を与える座標の構造を初期構造にして，エネルギー極小化を行い，最終的にNPT条件下での結合構造を得た(Figure 4)．各結合構造の室温での集団変数の時間平均はメタダイナミクスのものに近かった．また，構造P₁は最安定構造で，これは実験で決定されたものをよく再現した．Figure 3から分かるように，構造P₁から構造P₂へは，ϕ₁の線形変化がエネルギー障壁の小さな径路で，これはFigure 4でテオフィリン分子が回転して水素結合を切断し，その後に新たな水素結合を形成する過程に対応している．これらの結果は，メタダイナミクスで得られた自由エネルギー面から，一意的な結合構造を抽出できることを示唆している．今後，抽出した構造を標準自由エネルギー計算の初期構造に用いた場合の予測精度，自由度の大きな低分子の場合に同様な手法で結合構造を抽出できるかどうか，などを検討する．

Figure 4.

 Obtained (meta-) stable binding poses around the binding pocket. The most favorable pose is P₁.

参考文献

• [1]   A. Barducci, G. Bussi, M. Parrinello, Phys. Rev. Lett., 100, 020603 (2008).
• [2]   G. M. Clore, J. Kuszewski, J. Am. Chem. Soc., 125, 1518 (2003).
• [3]   J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser, C. Simmerling, J. Chem. Theory Comput., 11, 3696 (2015).
• [4]   C. I. Bayly, et al., J. Phys. Chem., 97, 10269 (1993).
• [5]   J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A. Case, J. Comput. Chem., 25, 1157 (2004).
• [6]   H. Gouda, I. D. Kuntz, D. A. Case, P. A. Kollman, Biopolymers, 68, 16 (2003).
• [7]   G. A. Tribello, et al., Comput. Phys. Commun., 185, 604 (2014).