Zenn
💊

GROMACSのシミュレーション解析であるRMSD, RMSF, RoG, hydrogen bond, RDF解析を論文に沿ってやってみた

2024/12/02に公開
gromacs
insilico創薬
mdsimulation
タンパク質
tech

本記事はリガンド-タンパク質 MDsimulation（分子動力学シミュレーション）を論文に沿って行なった結果に対し、
論文で使われている解析方法であるRMSD, RMSF, RoG, hydrogen bond, RDF解析を行なっています。

動作検証済み環境

Windows 11 Home, 13th Gen Intel(R) Core(TM) i7-13700, 64 ビット オペレーティング システム、x64 ベース プロセッサ, メモリ：32GB

本記事は以下の論文のフォローになります。
※本論文のMD simulationはAMBERを用いて行なっておりますが、本記事ではGROMACSでフォローしています。

https://www.nature.com/articles/s41598-022-10364-z

本論文は標的の選定方法からスクリーニング、スクリーニングした薬物候補化合物の毒性、物性予測、MDシミュレーションまでしており、ほぼ無料のツールですることができます。in silico創薬の基本が詰まっているので、一緒に勉強していきましょう。

宣伝

本記事を見てくださり、ありがとうございます。

インシリコ創薬についてより学びたい方は

拙著
https://zenn.dev/labcode/books/xvsw0lszpqir77
で学び、さらに色々な方法で新薬探索を楽しんでいただければと思います！

また化合物の評価を行いたい場合は
https://zenn.dev/labcode/books/ar7sn7jlhhafgm
を見ていただければと大変嬉しいです。

In silico創薬とその流れ

In silico 創薬、in silicoスクリーニングとは

インシリコ創薬（in silico drug discovery）は、コンピュータを駆使し、シミュレーションやデータ解析を用いて新薬を設計・発見する手法です。

これにより、従来の実験的な方法に比べてコストと時間を大幅に削減できるとされています。

その中でもin silicoスクリーニングは、膨大なデータベースから有望な候補物質を迅速に特定し、その効果や副作用を予測することで、創薬の初期段階での効率化が図られます。

一見難しそうなin silico創薬ですが、現在では様々なアプリケーションやwebサイトがあり、それらを駆使すれば、誰でも簡単に創薬をすることもできます。本記事では、それらのアプリケーションを駆使し、in silico創薬を行った論文をもとに、手法をわかりやすく説明していきます。

In silico創薬の流れ

  1. ターゲット選定、準備：疾患の原因となる分子（ターゲット）を特定し、その3D構造を準備します。
  2. 化合物ライブラリの構築：in silicoスクリーニングに使用する化合物集団（＝ライブラリ）を作成します。
  3. in silicoスクリーニング：コンピュータ上で数百万の化合物を対象に、ターゲット分子との結合親和性を評価します。
  4. 分子動力学シミュレーション（MD simulation)：候補化合物とターゲットの相互作用を動的に解析し、安定性や効果を予測します。
    4.1.環境構築:Windowsでの環境構築を示しています。
    4.2.MD simulation:MD simulationを行います。結果の可視化を行います。
    4.3.RMSD, RMSF, RoG, hydrogen bond, RDF 解析(本記事):構造変化や分子間相互作用を評価するための解析手法を行います。
    4.4.MM-GBSA解析とMM-PBSA解析、decomposition解析、アラニンスキャン:分子間相互作用や結合エネルギーへの各残基の寄与を解析します。
  5. 物性、毒性評価：過去のデータを用いて、新たな候補物質の特性や副作用を予測するモデルを構築します。
  6. （実験的検証）：in silico解析で得られた候補物質を実験室で合成し、実際の効果や安全性を検証します。

すでに公開したものについては、リンクを貼っています。今回はMaterial and Methodの以下の部分をフォローしています。

この記事では上記の論文のFigure5, 7の作り方を説明します。分子ドッキングの結果が違うので、正確に同じグラフにはなりませんが、作成の仕方を参考にしてみてください。

Indexファイルの作成

それでは実際に解析していきます。

RMSDではタンパク質のリガンドの重原子の位置がどれくらい変化するかを計測します。

リガンドの重原子については、デフォルト設定では選択できないので、それを選択できるようにします。以下を実行してください。

gmx make_ndx -f md_0_10.gro -n index.ndx

gmx make_ndx コマンドは、GROMACSでインデックスグループを作成または編集するために使います。以下はそのコマンドの使い方です。

  • f em.gro: 構造ファイル（ここではem.gro）を指定します。このファイルからインデックスグループを生成します。
  • n index.ndx: 既存のインデックスファイルを読み込んで編集する場合に使います。ここでは index.ndx を読み込んで編集できます。

以下のように出力されます。

                      :-) GROMACS - gmx make_ndx, 2024 (-:

Executable:   /usr/local/gromacs/bin/gmx
Data prefix:  /usr/local/gromacs
Working dir:  /home/shizuku/protein_ligand_complex_ver3
Command line:
  gmx make_ndx -f em.gro -n index.ndx

Reading structure file
Going to read 1 old index file(s)

  0 System              : 169799 atoms
  1 Protein             : 13165 atoms
  2 Protein-H           :  6560 atoms
  3 C-alpha             :   890 atoms
  4 Backbone            :  2670 atoms
  5 MainChain           :  3562 atoms
  6 MainChain+Cb        :  4342 atoms
  7 MainChain+H         :  4414 atoms
  8 SideChain           :  8751 atoms
  9 SideChain-H         :  2998 atoms
 10 Prot-Masses         : 13165 atoms
 11 non-Protein         : 156634 atoms
 12 Other               :    29 atoms
 13 UNL                 :    29 atoms
 14 NA                  :    23 atoms
 15 Water               : 156582 atoms
 16 SOL                 : 156582 atoms
 17 non-Water           : 13217 atoms
 18 Ion                 :    23 atoms
 19 Water_and_ions      : 156605 atoms
 20 Protein_UNL         : 13194 atoms

 nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list groups
 'a': atom       '&': and  'del' nr         'splitres' nr   'l': list residues
 't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help
 'r': residue              'res' nr         'chain' char
 "name": group             'case': case sensitive           'q': save and quit
 'ri': residue index

>

今回新しくリガンドの重原子を作成したいため、

13 & a H*　と打ち、Enterしてください。

次に以下を打ち、21番目にリガンドの重原子を選択できるようにします。

name 21 UNL_Heaby

最後にqを打ち、インタラクティブモードから出てください。

RMSD解析

RMSD（Root Mean Square Deviation）は、分子動力学シミュレーションにおいて、分子の構造変化を時間に渡って追跡するための解析手法です。RMSDは、シミュレーション中の各フレームの原子の位置を、参照構造と比較して、どれだけ変化しているかを測定します。

今回は先ほど作成した UNL_Heaby　がシミュレーション中でどのくらい変化しているかを見ます。

まずはcdでMD simulationで得られたデータが入っているディレクトリに移動してください。

md_0_10.tpr、　md_0_10_center.xtc 、index.ndxの三つが入っているか解説してください。

RMSD解析

gmx rms -s md_0_10.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd.xvg

このコマンドは、GROMACSのgmx rmsツールを使って、分子動力学（MD）シミュレーションのRMSD（Root Mean Square Deviation、二乗平均平方根偏差）を計算するためのものです。特定の構造（例: エネルギー最小化後の構造）とシミュレーション中の構造との違いを時間経過に伴って測定します。RMSDは、シミュレーションの安定性や系の構造変化を評価するための重要な指標です。

gmx rmsは、GROMACSのRMSDを計算するためのコマンドです。このツールは、リファレンス構造（通常はエネルギー最小化後の構造など）に対する各時点のシミュレーション構造の座標を比較し、RMSDを計算します。

  • gmx rmsは、GROMACSのRMSDを計算するためのコマンドです。このツールは、リファレンス構造（通常はエネルギー最小化後の構造など）に対する各時点のシミュレーション構造の座標を比較し、RMSDを計算します。
  • sはリファレンス構造を指定するオプションです。ここでは、エネルギー最小化後の構造情報が入ったem.tprファイルを指定しています。このファイルは、RMSD計算時の基準となる構造です。通常はシミュレーション開始時の最適化された構造を使用します。
  • fは、RMSDを計算するために使用するシミュレーションの軌道（トラジェクトリ）ファイルを指定します。ここでは、軌道ファイルmd_0_10_center.xtcを指定しています。このファイルにはシミュレーション中の構造の座標情報が含まれており、特定の基準構造との比較に使われます。
  • nはインデックスファイルを指定するオプションです。index.ndxファイルには、特定の原子や分子グループが定義されています。このファイルを使うことで、RMSDを計算する際にどの部分（タンパク質全体、特定のサブユニット、リガンドなど）を比較対象とするかを指定できます。たとえば、タンパク質の主鎖のみや、リガンドのRMSDを計算したい場合に便利です。
  • tuは時間単位を指定するオプションで、ここではns（ナノ秒）を指定しています。これにより、出力されるRMSDデータの時間単位がナノ秒になります。デフォルトではピコ秒（ps）ですが、ナノ秒の方が長時間シミュレーションでは扱いやすいです。
  • oは、RMSD計算結果の出力ファイル名を指定します。ここでは、rmsd.xvgという名前のファイルにRMSDのデータが出力されます。このファイルは、Xmgraceなどのプロットツールで可視化できる形式（.xvg）です。時間とRMSD値の変化をグラフ化する際に使用します。

以下のように出力されます。

Command line:
  gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd.xvg

Reading file em.tpr, VERSION 2024 (single precision)
Reading file em.tpr, VERSION 2024 (single precision)
Select group for least squares fit
Group     0 (         System) has 169799 elements
Group     1 (        Protein) has 13165 elements
Group     2 (      Protein-H) has  6560 elements
Group     3 (        C-alpha) has   890 elements
Group     4 (       Backbone) has  2670 elements
Group     5 (      MainChain) has  3562 elements
Group     6 (   MainChain+Cb) has  4342 elements
Group     7 (    MainChain+H) has  4414 elements
Group     8 (      SideChain) has  8751 elements
Group     9 (    SideChain-H) has  2998 elements
Group    10 (    Prot-Masses) has 13165 elements
Group    11 (    non-Protein) has 156634 elements
Group    12 (          Other) has    29 elements
Group    13 (            UNL) has    29 elements
Group    14 (             NA) has    23 elements
Group    15 (          Water) has 156582 elements
Group    16 (            SOL) has 156582 elements
Group    17 (      non-Water) has 13217 elements
Group    18 (            Ion) has    23 elements
Group    19 ( Water_and_ions) has 156605 elements
Group    20 (    Protein_UNL) has 13194 elements
Group    21 (      UNL_Heaby) has    14 elements
Select a group:

4を選択してください。

Select group for RMSD calculation
Group     0 (         System) has 169799 elements
Group     1 (        Protein) has 13165 elements
Group     2 (      Protein-H) has  6560 elements
Group     3 (        C-alpha) has   890 elements
Group     4 (       Backbone) has  2670 elements
Group     5 (      MainChain) has  3562 elements
Group     6 (   MainChain+Cb) has  4342 elements
Group     7 (    MainChain+H) has  4414 elements
Group     8 (      SideChain) has  8751 elements
Group     9 (    SideChain-H) has  2998 elements
Group    10 (    Prot-Masses) has 13165 elements
Group    11 (    non-Protein) has 156634 elements
Group    12 (          Other) has    29 elements
Group    13 (            UNL) has    29 elements
Group    14 (             NA) has    23 elements
Group    15 (          Water) has 156582 elements
Group    16 (            SOL) has 156582 elements
Group    17 (      non-Water) has 13217 elements
Group    18 (            Ion) has    23 elements
Group    19 ( Water_and_ions) has 156605 elements
Group    20 (    Protein_UNL) has 13194 elements
Group    21 (      UNL_Heaby) has    14 elements
Select a group:

21を選択してください。

上記の設定は、シミュレーション中にリガンド（UNL_Heavy）がタンパク質の主鎖構造に対してどのように動いたか、その相対的な変化を計算するためのものです。RMSD値が時間とともに安定している場合、リガンドは基準構造（タンパク質の主鎖）に対して安定していることを意味し、大きく変動する場合はリガンドがシミュレーション中に動き回ったことを示します。

lsコマンドを押して、以下のrmsd.xvgあるか確認してください。

rmsd.xvgをテキストで開いた後、エクセルなどで張り付けてグラフを作成してください。

以下のようにRMSDができます。

このRMSDの変動パターンから、シミュレーションの最初の25 nsでは大きな構造変化が起こっていたが、以降は徐々に安定し、最終的には構造が収束していると考えられます。これにより、リガンドが結合部位に安定した可能性が示唆されますが、いくつかの時間点で再び変動が見られるため、完全に安定した結合ではない可能性もあります。

RMSF解析

RMSF（Root Mean Square Fluctuation）は、分子動力学シミュレーションで各原子や残基の時間に伴う変動量を測定する方法です。具体的には、各原子がシミュレーション中に平均位置からどれくらい揺れ動いたかを示します。これにより、分子のどの部分が柔軟で、どの部分が安定しているかがわかります。

gmx rmsf -s md_0_10.tpr -f md_0_10_center.xtc -o rmsf.xvg -res
  • gmx rmsf: GROMACSのRMSF解析を行うコマンド。
  • s topol.tpr: シミュレーションのトポロジーファイル（.tpr）を指定します。
  • f center.xtc: トラジェクトリーファイル（.xtc）を指定します。このファイルにはシミュレーションの座標データが入っています。
  • o rmsf.xvg: RMSFの結果を出力するファイル（.xvg形式）を指定します。
  • res: 残基ごとにRMSFを計算するオプションです。

以下のように出力されます。

(base) shizuku@DESKTOP-5I5GHRA:~/protein_ligand_complex_ver2$ gmx rmsf -s md_0_10.tpr -f md_0_10_center.xtc -o rmsf.xvg -res
                        :-) GROMACS - gmx rmsf, 2024 (-:

Executable:   /usr/local/gromacs/bin/gmx
Data prefix:  /usr/local/gromacs
Working dir:  /home/shizuku/protein_ligand_complex_ver2
Command line:
  gmx rmsf -s md_0_10.tpr -f md_0_10_center.xtc -o rmsf.xvg -res

Reading file md_0_10.tpr, VERSION 2024 (single precision)
Reading file md_0_10.tpr, VERSION 2024 (single precision)
Select group(s) for root mean square calculation
Group     0 (         System) has 169799 elements
Group     1 (        Protein) has 13165 elements
Group     2 (      Protein-H) has  6560 elements
Group     3 (        C-alpha) has   890 elements
Group     4 (       Backbone) has  2670 elements
Group     5 (      MainChain) has  3562 elements
Group     6 (   MainChain+Cb) has  4342 elements
Group     7 (    MainChain+H) has  4414 elements
Group     8 (      SideChain) has  8751 elements
Group     9 (    SideChain-H) has  2998 elements
Group    10 (    Prot-Masses) has 13165 elements
Group    11 (    non-Protein) has 156634 elements
Group    12 (          Other) has    29 elements
Group    13 (            UNL) has    29 elements
Group    14 (             NA) has    23 elements
Group    15 (          Water) has 156582 elements
Group    16 (            SOL) has 156582 elements
Group    17 (      non-Water) has 13217 elements
Group    18 (            Ion) has    23 elements
Group    19 ( Water_and_ions) has 156605 elements
Select a group:

3を選択してください。

rmsf.xvgが出てくるので、これをエクセルで書いてください。

以下のようになります。

このグラフは、タンパク質の残基ごとのRMSF (Root Mean Square Fluctuation) を示しています。RMSFは、各残基の平均的な位置からの振動（フレキシビリティ）を評価する指標です。以下のように解釈できます：

  • X軸 (Residue): 残基番号を表しており、312番目の残基から開始しています。グラフに表示されている範囲は、312番目から約700番目の残基に対応します。
  • Y軸 (RMSF): 各残基の振動の大きさ（Å単位）を示しています。数値が大きいほど、その残基がシミュレーション中に大きく動いている、すなわちフレキシブルであることを示します。
  • グラフの左側（残基番号が小さい部分）では、RMSFが比較的低く（1-2Å程度）、この領域の残基はあまり動いていない（安定している）ことを示しています。
  • 逆に、グラフの右側（残基番号が大きい部分）では、RMSFが高く（最大4Å以上）、この領域の残基は比較的動きやすいことを意味します。
  • 全体的に、特定の残基が局所的に大きなフレキシビリティを示している部分（例えば600番目付近）が見られ、この領域がタンパク質の可動性が高い部分であることを示唆しています。

これにより、動きが少ない安定した領域と、動きが大きく柔軟性が高い領域を特定することができます。リガンドの結合したタンパク質とリガンドが結合していないタンパク質をRMSF解析で比べて、リガンドの結合による影響をみることが多いです。

RoG解析

RoG解析（Radius of Gyration、回転半径）は、分子がシミュレーション中にどれだけコンパクトにまとまっているかを測定する指標です。これは、分子全体の重心から各原子の距離を考慮して計算され、タンパク質や分子の折りたたみ状態や構造の安定性を評価するのに使われます。

gmx gyrate -s md_0_10.tpr -f md_0_10_center.xtc -o gyrate.xvg

コマンド gmx gyrate は、分子の重心からの距離を基準に、シミュレーション中のタンパク質や分子の回転半径（gyration radius）を計算するために使用されます。

以下のように出力されます。

Command line:
  gmx gyrate -s md_0_10.tpr -f md_0_10_center.xtc -o gyrate.xvg

Reading file md_0_10.tpr, VERSION 2024 (single precision)
Reading file md_0_10.tpr, VERSION 2024 (single precision)
Available static index groups:
 Group  0 "System" (169799 atoms)
 Group  1 "Protein" (13165 atoms)
 Group  2 "Protein-H" (6560 atoms)
 Group  3 "C-alpha" (890 atoms)
 Group  4 "Backbone" (2670 atoms)
 Group  5 "MainChain" (3562 atoms)
 Group  6 "MainChain+Cb" (4342 atoms)
 Group  7 "MainChain+H" (4414 atoms)
 Group  8 "SideChain" (8751 atoms)
 Group  9 "SideChain-H" (2998 atoms)
 Group 10 "Prot-Masses" (13165 atoms)
 Group 11 "non-Protein" (156634 atoms)
 Group 12 "Other" (29 atoms)
 Group 13 "UNL" (29 atoms)
 Group 14 "NA" (23 atoms)
 Group 15 "Water" (156582 atoms)
 Group 16 "SOL" (156582 atoms)
 Group 17 "non-Water" (13217 atoms)
 Group 18 "Ion" (23 atoms)
 Group 19 "Water_and_ions" (156605 atoms)
Specify a selection for option 'sel'
(Select group to compute gyrate radius):
(one per line, <enter> for status/groups, 'help' for help)
>

1を選択してください。

gyrate.xvg　が出力されます。一行目と二行目のデータを使って、散布図を書いてください。

以下のように書けます。

  • シミュレーションの最初の部分では、RoGが35Å付近で変動していますが、徐々に減少し、約30Åで安定しています。
  • 250 ns以降は、RoGが一定の範囲で安定しており、これはシミュレーションが収束し、タンパク質の全体的な構造が安定したことを示しています。

水素結合数測定

MDシミュレーションにおける水素結合数の測定は、分子間および分子内での結合の強さや安定性を評価する重要な解析手法です。水素結合（H-bond）は、分子が互いに相互作用する際に形成される強力な非共有結合の一種で、タンパク質-リガンド相互作用や分子の安定性において非常に重要です。

以下を実行して、タンパク質とリガンドの水素結合の数を計算してください。

※gmx hbond コマンドではうまく解析できないことがあるそうです。

gmx hbond-legacy -s md_0_10.tpr -f md_0_10_center.xtc -num hydrogen-bonds.xvg -tu ns

gmx hbond-legacy コマンドを使用して、水素結合を解析し、-num オプションで結合の数を hydrogen-bonds.xvg に出力するコマンドの説明に移ります。

以下のように出力されます。

Command line:
  gmx hbond -s md_0_10.tpr -f md_0_10_center.xtc -num hydrogen-bonds.xvg -tu ns

Reading file md_0_10.tpr, VERSION 2024 (single precision)
Reading file md_0_10.tpr, VERSION 2024 (single precision)
Available static index groups:
 Group  0 "System" (169799 atoms)
 Group  1 "Protein" (13165 atoms)
 Group  2 "Protein-H" (6560 atoms)
 Group  3 "C-alpha" (890 atoms)
 Group  4 "Backbone" (2670 atoms)
 Group  5 "MainChain" (3562 atoms)
 Group  6 "MainChain+Cb" (4342 atoms)
 Group  7 "MainChain+H" (4414 atoms)
 Group  8 "SideChain" (8751 atoms)
 Group  9 "SideChain-H" (2998 atoms)
 Group 10 "Prot-Masses" (13165 atoms)
 Group 11 "non-Protein" (156634 atoms)
 Group 12 "Other" (29 atoms)
 Group 13 "UNL" (29 atoms)
 Group 14 "NA" (23 atoms)
 Group 15 "Water" (156582 atoms)
 Group 16 "SOL" (156582 atoms)
 Group 17 "non-Water" (13217 atoms)
 Group 18 "Ion" (23 atoms)
 Group 19 "Water_and_ions" (156605 atoms)
Specify a selection for option 'r'
(Reference selection, relative to which the search for hydrogen bonds in target selection will develop.):
(one per line, <enter> for status/groups, 'help' for help)
>

1を選択してください。

Available static index groups:
 Group  0 "System" (169799 atoms)
 Group  1 "Protein" (13165 atoms)
 Group  2 "Protein-H" (6560 atoms)
 Group  3 "C-alpha" (890 atoms)
 Group  4 "Backbone" (2670 atoms)
 Group  5 "MainChain" (3562 atoms)
 Group  6 "MainChain+Cb" (4342 atoms)
 Group  7 "MainChain+H" (4414 atoms)
 Group  8 "SideChain" (8751 atoms)
 Group  9 "SideChain-H" (2998 atoms)
 Group 10 "Prot-Masses" (13165 atoms)
 Group 11 "non-Protein" (156634 atoms)
 Group 12 "Other" (29 atoms)
 Group 13 "UNL" (29 atoms)
 Group 14 "NA" (23 atoms)
 Group 15 "Water" (156582 atoms)
 Group 16 "SOL" (156582 atoms)
 Group 17 "non-Water" (13217 atoms)
 Group 18 "Ion" (23 atoms)
 Group 19 "Water_and_ions" (156605 atoms)
Specify a selection for option 't'
(Target selection, relative to which the search for hydrogen bonds in reference selection will develop.):
(one per line, <enter> for status/groups, 'help' for help)
>

13を選択してください。

hydrogen-bonds.xvg が出力されているので、こちらの一行目と二行目のデータを使ってグラフを書いてください。

このグラフから、最初の25 nsではタンパク質とリガンドの間に比較的安定した水素結合が存在していますが、シミュレーションが進行するにつれて、水素結合の数が減少し、不安定な状態になっていることがわかります。これは、リガンドが結合ポケット内で再配置されているか、結合から外れ始めていることを示している可能性があります。

RDF 解析

Figure7のRDF plotを作成します。

タンパク質－リガンドの相互作用におけるRDF（Radial Distribution Function、放射分布関数）は、リガンドの原子とタンパク質の特定の原子や残基との間の相互作用を距離に基づいて解析する重要な指標です。RDFプロットを使うことで、リガンドとタンパク質の結合部位の相互作用の強さや距離依存性を視覚的に理解することができます。

上記のSer239とArg399に当たるアミノ酸とリガンドとのRDF解析をしていきます。

まずindexファイルを作成します。

gmx make_ndx -f md_0_10.gro -o index.ndx

以下のように出力されます。

Command line:
  gmx make_ndx -f md_0_10.gro -o index.ndx

Reading structure file
Going to read 0 old index file(s)
Analysing residue names:
There are:   890    Protein residues
There are:     1      Other residues
There are: 51780      Water residues
There are:    23        Ion residues
Analysing Protein...
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...

  0 System              : 168557 atoms
  1 Protein             : 13165 atoms
  2 Protein-H           :  6560 atoms
  3 C-alpha             :   890 atoms
  4 Backbone            :  2670 atoms
  5 MainChain           :  3562 atoms
  6 MainChain+Cb        :  4342 atoms
  7 MainChain+H         :  4414 atoms
  8 SideChain           :  8751 atoms
  9 SideChain-H         :  2998 atoms
 10 Prot-Masses         : 13165 atoms
 11 non-Protein         : 155392 atoms
 12 Other               :    29 atoms
 13 UNL                 :    29 atoms
 14 NA                  :    23 atoms
 15 Water               : 155340 atoms
 16 SOL                 : 155340 atoms
 17 non-Water           : 13217 atoms
 18 Ion                 :    23 atoms
 19 Water_and_ions      : 155363 atoms

 nr : group      '!': not  'name' nr name   'splitch' nr    Enter: list groups
 'a': atom       '&': and  'del' nr         'splitres' nr   'l': list residues
 't': atom type  '|': or   'keep' nr        'splitat' nr    'h': help
 'r': residue              'res' nr         'chain' char
 "name": group             'case': case sensitive           'q': save and quit
 'ri': residue index

>

Ser239は550番目のアミノ酸なので、

ri 550と記入し、Enter,

name 20 Ser239と記入し、Enterを押してください。

続いて、

Arg399は710番目のアミノ酸に当たるので、

ri 710 と記入し、Enter

name 21 Arg399と記入し、Enterを押してください。

q で戻ってください。

続いて、以下のコードを実行してください。

gmx rdf -s md_0_10.tpr -f md_0_10_center.xtc -n index.ndx -o rdf_Ser239_Arg399.xvg -b 0 -e 10000 -bin 0.01 -tu ns  -norm rdf
  • gmx rdf: GROMACSでRDFを計算するコマンド。
  • b 0: シミュレーションの最初のフレームを使って解析を開始します（0 nsから開始）。
  • e 10000: 解析を終了するフレーム（ここでは10,000 ns）。
  • bin 0.01: RDFを計算する際のビン幅。0.01 nmごとに距離を区切ってRDFを計算します。
  • tu ns: 時間単位をナノ秒（ns）に設定。
  • norm rdf: 正規化方法を標準的なRDF（全体の平均密度で正規化）に指定。

以下のように出力されます。

                        :-) GROMACS - gmx rdf, 2024 (-:

Executable:   /usr/local/gromacs/bin/gmx
Data prefix:  /usr/local/gromacs
Working dir:  /home/shizuku/protein_ligand_complex_ver6
Command line:
  gmx rdf -s md_0_10.tpr -f md_0_10_center.xtc -n index.ndx -o rdf_Ser239_Arg399.xvg -b 0 -e 10000 -bin 0.01 -tu ns -norm rdf

Available static index groups:
 Group  0 "System" (168557 atoms)
 Group  1 "Protein" (13165 atoms)
 Group  2 "Protein-H" (6560 atoms)
 Group  3 "C-alpha" (890 atoms)
 Group  4 "Backbone" (2670 atoms)
 Group  5 "MainChain" (3562 atoms)
 Group  6 "MainChain+Cb" (4342 atoms)
 Group  7 "MainChain+H" (4414 atoms)
 Group  8 "SideChain" (8751 atoms)
 Group  9 "SideChain-H" (2998 atoms)
 Group 10 "Prot-Masses" (13165 atoms)
 Group 11 "non-Protein" (155392 atoms)
 Group 12 "Other" (29 atoms)
 Group 13 "UNL" (29 atoms)
 Group 14 "NA" (23 atoms)
 Group 15 "Water" (155340 atoms)
 Group 16 "SOL" (155340 atoms)
 Group 17 "non-Water" (13217 atoms)
 Group 18 "Ion" (23 atoms)
 Group 19 "Water_and_ions" (155363 atoms)
 Group 20 "Ser239" (17 atoms)
 Group 21 "Arg399" (15 atoms)
Specify a selection for option 'ref'
(Reference selection for RDF computation):
(one per line, <enter> for status/groups, 'help' for help)
>

13を選択してください。

続いて、以下が出力されます。

Available static index groups:
 Group  0 "System" (168557 atoms)
 Group  1 "Protein" (13165 atoms)
 Group  2 "Protein-H" (6560 atoms)
 Group  3 "C-alpha" (890 atoms)
 Group  4 "Backbone" (2670 atoms)
 Group  5 "MainChain" (3562 atoms)
 Group  6 "MainChain+Cb" (4342 atoms)
 Group  7 "MainChain+H" (4414 atoms)
 Group  8 "SideChain" (8751 atoms)
 Group  9 "SideChain-H" (2998 atoms)
 Group 10 "Prot-Masses" (13165 atoms)
 Group 11 "non-Protein" (155392 atoms)
 Group 12 "Other" (29 atoms)
 Group 13 "UNL" (29 atoms)
 Group 14 "NA" (23 atoms)
 Group 15 "Water" (155340 atoms)
 Group 16 "SOL" (155340 atoms)
 Group 17 "non-Water" (13217 atoms)
 Group 18 "Ion" (23 atoms)
 Group 19 "Water_and_ions" (155363 atoms)
 Group 20 "Ser239" (17 atoms)
 Group 21 "Arg399" (15 atoms)
Specify any number of selections for option 'sel'
(Selections to compute RDFs for from the reference):
(one per line, <enter> for status/groups, 'help' for help, Ctrl-D to en

20を選択し、Enter、

21を選択し、Enter

ctrl-Dで終了します。終了すると、計算が始まります。

rdf_Ser239_Arg399.xvg が生成しているので、グラフにしてみてください。

青線（Ser(239)）のピーク

  • ピーク位置: 青線のピークは 約4.0 nm の距離に位置しています。このピークはSer(239)に対する参照粒子が、約4.0 nmの距離で最も多く存在していることを示しています。
  • g(r)の値: g(r)の最大値は 約8.0 です。これは、Ser(239)の周りに4.0 nmの距離で、全体の平均密度に対して8倍の密度で粒子が集まっていることを示しています。これは、Ser(239)と他の粒子が4.0 nmの距離で強い相互作用を持っている可能性を示唆します。

赤線（Arg(399)）のピーク

  • ピーク位置: 赤線のピークは 約4.8 nm に位置しています。Arg(399)との相互作用がこの距離で最も強く、粒子が集まっていることを示しています。
  • g(r)の値: Arg(399)に対するg(r)の最大値は 約6.0 です。Ser(239)に比べると、少し密度が低いことを示していますが、それでも強い相互作用を持っていることがわかります。

2つのピークの比較

  • 距離の違い: Ser(239)は 約4.0 nm、Arg(399)は 約4.8 nm の距離で最も強い相互作用を示しています。このことから、これら2つの残基が異なる距離でリガンドや他の分子と相互作用していることが示されています。
    • これは、リガンドや他の分子が特定の距離でこれらの残基に結合する、もしくは相互作用することを示唆します。
  • 強度の違い: g(r)のピーク値から、Ser(239)の方がArg(399)よりもやや強い相互作用を示しています。これは、Ser(239)周囲の局所密度が高いことを意味しています。

以上でRMSD, RMSF, RoG, hydrogen bond, RDF解析は終了です！

参考文献

https://www.udemy.com/share/105vCm3@9eQi5TNCBbh2arjqlQq11eIj5xU4FyZWgdCUVdVrt0drqyT67M0pVeRk39zckzTBVw==/

宣伝

本記事を見てくださり、ありがとうございます。

インシリコ創薬についてより学びたい方は

拙著
https://zenn.dev/labcode/books/xvsw0lszpqir77
で学び、さらに色々な方法で新薬探索を楽しんでいただければと思います！

また化合物の評価を行いたい場合は
https://zenn.dev/labcode/books/ar7sn7jlhhafgm
を見ていただければと大変嬉しいです。

Discussion