GROMACS で膜シミュレーションを実行する

膜シミュレーションの実行

ユーザーは、特にタンパク質が関与する場合、脂質二重膜のシミュレーションを実行する際に、頻繁に問題に遭遇します。膜タンパク質のシミュレーションを試みるユーザーは、この「チュートリアル <https://tutorials.gromacs.org/membrane-protein.html>」が役立つ可能性があります。

膜タンパク質のシミュレーションのための1つのプロトコルは、以下の手順で構成されます:

  1. タンパク質と脂質に対して、パラメータが利用可能な力場を選択してください。

  2. タンパク質を膜に挿入します。(たとえば、あらかじめ形成された二重膜を使用するか、粗い粒子の自己組み立てシミュレーションを実行し、その後、原子レベルの表現に戻す。)

  3. システムを中和し、過剰な電荷を中和し、最終的なイオン濃度を調整するために、イオンを添加します。

  4. エネルギー消費を最小限に抑える。

  5. タンパク質に膜が適応するようにします。通常、タンパク質のすべてのヘビー原子に対して、約5〜10nsのMDシミュレーションを実行し、制約（1000 kJ/(mol nm2)）を適用します。

  6. 制約なしで平衡状態に達する。

  7. 本番環境でMDを実行する。

genbox を使用して水を加える

水分子が、あらかじめ形成された脂質膜の周囲に存在する場合、solvate を使用して水を取り込むと、水分子が膜の隙間に導入されることがあります。 これらを除去するためのいくつかの方法があります

  • 短いMDシミュレーションを実行し、これらの水を除去することで疎水効果を抑制します。一般的には、分子が迅速かつスムーズに排除され、全体的な構造を損なうことなく、水を含まない疎水相に到達するのに十分です。もし、あなたの設定が完全に水を含まない疎水相を初期状態として必要とする場合は、以下のアドバイスを試してみてください：

  • gmx solvate での -radius オプションを設定して、水の排除半径を変更します。

  • ``vdwradii.dat``を ``$GMXLIB``フォルダからワークディレクトリにコピーし、編集して、リン脂質の原子の半径を（特に炭素の場合は0.35～0.5nmの間）調整してください。これにより、 solvate が水挿入に適した隙間を見つけにくくなります

  • 手動で構造を編集して削除する（:ref:`gro`ファイルの原子数を調整し、および:ref:`topology <top>`の変更を考慮することを忘れない）

  • 別の人が書いたスクリプトを使用して、それらを削除してください。

外部資料

  • 膜シミュレーションのスライド、膜シミュレーションのビデオ - (エリック・リンダール)。

  • 膜タンパク質のシミュレーションに関するチュートリアル – これは、脂質二重層の中に埋め込まれたタンパク質のシミュレーションにおいて、どのような質問や問題が発生するかを示すように設計されています。

  • OPLS-AA力場とBergerの脂質を組み合わせる：詳細な説明、その理由、方法、およびテストについて。

  • 複数の膜タンパク質のトポロジーと異なる力場（GAFF、CHARMM、バーガー）：Shirley W. I. Siu、Robert Vacha、Pavel Jungwirth、Rainer A. Böckmannによる「膜の分子シミュレーション：異なる力場からの物理的特性 <https://doi.org/10.1063/1.2897760>」。

  • Lipidbook は、膜と膜タンパク質のシミュレーションで使用される脂質、界面活性剤、その他の分子の力場パラメータの公開リポジトリです。詳細は、J. Domański, P. Stansfeld, M.S.P. Sansom, and O. Beckstein. J. Membrane Biol. 236 (2010), 255—258. に記載されています。doi:10.1007/s00232-010-9296-8。

新規分子のパラメータ化

ほとんどのパラメータ設定に関する質問や問題は、以下の2つのルールを覚えておくことで、非常に簡単に解決できます

  • 異なる力場を混在させるべきではありません。:ref:`力場 <gmx-force-field>`は（最良の場合）、互いに整合性を保つように設計されており、他の力場と併用することは通常うまくいきません。もし、最初の力場で一部のシステムをシミュレーションし、別の力場で残りの部分をシミュレーションする場合、その結果は疑わしくなる可能性があり、レビュー担当者が懸念する可能性があります。適切な力場を選び、その力場を使用してください。

  • もし、新しいパラメータを開発する必要がある場合は、残りの力場が最初にどのように導出されたかに従って、それらを導出する必要があります。つまり、元の文献をレビューする必要があります。力場パラメータを導出する方法は唯一の方法はありません。必要なのは、残りの力場と一貫性のあるパラメータを導出することです。どのように行うかは、使用したい力場によって異なります。たとえば、AMBER力場の場合、非標準のアミノ酸のパラメータを導出するには、さまざまな量子計算を行う必要がある場合があります。一方、GROMOSまたはOPLSパラメータを導出するには、(a) さまざまな流体および液体状態の特性を調整し、(b) 経験、化学的直感、または類似に基づいてパラメータを調整することが含まれる場合があります。自動化されたアプローチに関するいくつかの提案は、こちらで確認できます：here。

GROMACS のシミュレーション経験をある程度持つことが、新しい力場や既存の力場のための新しい分子をパラメータ化する前に賢明な判断です。これらは専門的なトピックであり、特に研究プロジェクトとして学部生に提供するには適していません（ただし、高価な擬似乱数ジェネレーターを使用する場合は別）。GROMACS のリファレンスマニュアルの「第5章：相互作用関数と力場 <ff>」について、非常に詳細な知識が必要です。警告が十分に与えられていない場合は、以下に記載されている、異質な種に関するパラメータ化についてご確認ください。

別のアドバイス：パラメータの取得は、高級な宝石を購入するのと同じように、慎重に行うべきです。通りで「ダイヤモンド」のネックレスを10ドルで売っている人が現れたからといって、そこから購入する必要はありません。同様に、特に、パラメータの提供元が誰であるか、また、パラメータの取得方法が説明されていない場合、興味のある分子のパラメータを、聞いたことのないウェブサイトからダウンロードするだけでは、必ずしも最適な戦略とは限りません。

`PRODRG <http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg>`_のトポロジーを使用する際には、内容を確認せずに使用すると、その結果が`公開 <http://pubs.acs.org/doi/abs/10.1021/ci100335w>`_され、GROMOSファミリーの力場パラメータを適切に導出するためのヒントも含まれます。

珍しい種

したがって、タンパク質/核酸系のシミュレーションを行いたいが、様々な異種金属イオン（ルテニウムなど）に結合したり、機能に必要な鉄硫黄クラスタが存在したり、同様の状況である。しかし、（残念ながら）使用したい力場には、これらのパラメータが提供されていない。どうすればよいか？ |Gromacs|の`ユーザーフォーラム`にメールを送信し、FAQを参照する。

もし本当にこれらの分子動力学シミュレーションを試みるのであれば、文献から、または独自のパラメータ化によって、それらのパラメータを入手する必要があります。しかし、それを行う前に、一度立ち止まって検討することが重要です。なぜなら、そのような原子/クラスに対してパラメータがすでに存在しない理由は、時として存在するからです。特に、これらの原子/クラスに対して標準的なパラメータを開発/入手し、分子動力学で使用することが妥当かどうかを確認するために、以下の基本的な質問を自問することができます。

  • 量子効果（つまり電荷移動）が重要である可能性はありますか？（つまり、酵素の活性部位に二配位子を持つ金属イオンがあり、酵素の機能を研究することに関心がある場合、これは非常に重要な問題となる可能性があります）。

  • 私の選択した力場に対して、標準的な力場パラメータ化手法が、この種の原子/分子に対してうまく機能する可能性はありますか？（つまり、例えば Hartree-Fock 6-31G* が遷移金属を適切に記述できない場合など）

もし、これらの質問のいずれかに「はい」と答えた場合は、古典的な分子動力学以外の方法でシミュレーションを行うことを検討してください。

たとえこれらの質問に対する答えが「いいえ」であっても、興味のある化合物の専門家と相談してから、ご自身でのパラメータ化を試みることをお勧めします。さらに、これらの複雑なパラメータ化を試みる前に、より簡単なものをパラメータ化してみることをお勧めします。

平均力ポテンシャルの可能性

平均力場ポテンシャル（PMF）とは、特定の系のすべての構成における平均的な力を与えるポテンシャルを指します。 GROMACS でPMFを計算する方法はいくつかありますが、最も一般的な方法の1つは、プルコードを使用することです。アンブレラサンプリングを使用してPMFを取得する手順は次のとおりです。アンブレラサンプリングは、統計的にありえない状態のサンプリングを可能にします：

  • 反応経路（誘導MDシミュレーション、通常のMDシミュレーション、または任意に作成された設定から）に沿って、設定の系列を生成する。

  • 傘サンプリングを使用して、これらの設定をサンプリングウィンドウ内に制限します。

  • gmx wham を使用して、WHAM アルゴリズムを利用して PMF (Probability Matrix Function) 曲線を再構築します。

より詳細なチュートリアルは、こちらから参照できます: https://tutorials.gromacs.org/umbrella-sampling.html (傘状サンプリングに関する)。

集中エネルギー

特定の構成のエネルギーを計算することは、場合によっては便利な操作です。GROMACS でこれを実行する最良の方法は、mdrun の -rerun メカニズムを使用することです。このメカニズムは、tpr に含まれるモデル物理を、mdrun に提供された軌跡または座標ファイルにある構成に適用します。

mdrun -s input.tpr -rerun configuration.pdb

注意：提供された設定は、:ref:``tpr`ファイルを:ref:``grompp <gmx grompp>`を使用して生成した際に使用したトポロジーと一致している必要があります。 grompp <gmx grompp>`に提供した設定は、原子名を除いて、ほとんどの場合無関係です。この機能は、エネルギーグループ（参照マニュアルを参照）や、複数の設定の軌跡（この場合、デフォルトでは:ref:``mdrun <gmx mdrun>`が各設定に対して近傍検索を実行します。これは、入力が類似しているという前提を立てられないためです。）にも使用できます。

ゼロステップのエネルギー最小化は、エネルギーを報告する前に1ステップを実行しますが、ゼロステップのMD実行は、制約の存在下での可能な再始動に対応するために（回避可能な）問題を抱えています。そのため、これらの手順は推奨されません。

カーボンナノチューブ

ロバート・ジョンソンのヒント

引用自 Robert Johnson の gmx-users メーリングリストアーカイブでの投稿。

  • 必ず、「ターミナル」炭素原子が、ファイル内の構造において結合を共有していることを確認してください。

  • periodic_molecules = yes を、mdp ファイルに設定して、gmx grompp での入力に使用してください。

  • たとえトポロジーが正しいとしても、ナノチューブを不適切なサイズの箱に入れると、折り曲がりが生じる可能性があります。そのため、`VMD`_を使用して、ナノチューブと周期的な画像を視覚化し、画像間の間隔が正しいことを確認してください。間隔が小さすぎたり大きすぎたりすると、チューブに大きな応力が生じ、折り曲がりや伸びを引き起こす可能性があります。

  • ナノチューブの軸方向に圧力を加えないでください。 実際には、デバッグ目的で、問題が起こっているかどうかを特定するまで、圧力を完全にオフにしておく方が良いかもしれません。 そして、もし問題がある場合は、それが何かを特定してください。

  • x2top を特定の力場とともに使用する場合、分子の結合に関する仮定がなされます。ナノチューブの末端炭素原子は、周期的な場合、最大で2つの炭素原子に、非周期的な場合、または水素でキャップされた場合は1つの炭素原子にのみ結合されます。

  • 「-pbc」オプションを使用することで、無限長のナノチューブを生成できます。これにより、x2top は、ターミナルの C 原子が実際に化学結合を共有していることを認識します。したがって、grompp を使用しても、単結合の C に関するエラーは発生しません。

アンドレア・ミノアのチュートリアル

GROMACS を使用してカーボンナノチューブをモデリングする (また、http://chembytes.wikidot.com/grocnt にもアーカイブされています) には、OPLS-AA パラメータを使用して CNT の簡単なシミュレーションを設定するためのすべてが含まれています。 buildCstruct (Python スクリプトでターミナル水素も追加) や TubeGen Online (PDB 出力をコピーしてファイルに貼り付け、cnt.pdb という名前を付ける) などのツールを使用して、単純な CNT の構造を簡単に生成できます。

|Gromacs|の最新バージョンで動作させるには、通常、次の手順を実行する必要があります：

  • ディレクトリ cnt_oplsaa.ff を作成する

  • このディレクトリに、チュートリアルページのデータを使用して、以下のファイルを作成してください：

    • forcefield.itp を、itp セクションにあるファイル

    • atomnames2types.n2t ファイル（セクション n2t にあります）

    • aminoacids.rtp は、rtp セクションにあるファイル

  • カスタムフォースフィールドを使用してトポロジーを生成します（:ref:`gmx x2top`コマンドを実行するディレクトリに`cnt_oplsaa.ff`ディレクトリが存在するか、またはGMXLIBパスにある必要があります）。-noparam`オプションを使用すると、`:ref:`gmx x2top`がコマンドライン（-kb、-ka、-kd）から結合/角度/二面角の力定数を使用せず、代わりにフォースフィールドファイルを使用するように指示されます。ただし、この場合、次のステップ（二面角関数の修正）が必要になります。

gmx x2top -f cnt.gro -o cnt.top -ff cnt_oplsaa -name CNT -noparam

「ジヘドラルの関数タイプは、gmx x2top によって '1' に設定されていますが、フォースフィールドファイルでは '3' が指定されています。したがって、トポロジーファイルの [ジヘドラル] セクションにある関数タイプ '1' を '3' に置き換えてください。簡単な方法は sed を使用することですが、オペレーティングシステムに合わせて調整する必要がある場合があります。また、トップファイルを手動で確認し、ジヘドラルの関数タイプのみが変更されていることを確認してください。」

sed -i~ '/\[ dihedrals \]/,/\[ system \]/s/1 *$/3/' cnt.top

一度、構成が完了したら、システムをセットアップできます。たとえば、真空環境でのシミュレーション（お好みのパラメータを em.の mdp と md.の mdp に設定して実行する）：

少し大きめの箱に入れる:

gmx editconf -f cnt.gro -o boxed.gro -bt dodecahedron -d 1

真空中のエネルギーを最小化する:

gmx grompp -f em.mdp -c boxed.gro -p cnt.top -o em.tpr
gmx mdrun -v -deffnm em

真空中のMD:

gmx grompp -f md.mdp -c em.gro -p cnt.top -o md.tpr
gmx mdrun -v -deffnm md

軌跡を確認する:

gmx trjconv -f md.xtc -s md.tpr -o md_centered.xtc -pbc mol -center
gmx trjconv -s md.tpr -f md_centered.xtc -o md_fit.xtc -fit rot+trans
vmd em.gro md_fit.xtc