OpenMMの使い方を順番にたどる記事のPart 9です。

前回の記事（Part 8）では「低分子 - タンパク質複合体」についてシミュレーション系の構築と計算の実行を行いました。この記事では結果の解析を行います。

一般的な解析手法についてはPart 6、Part 7で確認しました。この記事では「低分子リガンドとタンパク質の相互作用解析」に焦点を当てたいと思います。

コードは基本的にTeachOpenCADD T020：Analyzing molecular dynamics simulationsとMaking it rain : Protein-Ligand smulationsからお借りしています。

1. キナーゼと阻害剤
2. 低分子 - タンパク質複合体のトラジェクトリ解析
3. おわりに

1. キナーゼと阻害剤

具体的な解析を眺める前に、計算対象の構造の特徴について前提知識のおさらいをしておきます。

今回の構造はPDB ID: 3POZでEGFRタンパクのキナーゼドメインとリガンド TAK-285のX線共結晶構造です。

キナーゼは「ATPなどの高エネルギーリン酸結合を有する分子からリン酸基を基質あるいはターゲット分子に転移する（リン酸化する）酵素の総称」（Wikipedia - キナーゼ）です。なのでATP結合部位があることが多く、EGFRタンパクの場合以下のようなイメージになるでしょうか？*1

ATPのアデニンと水素結合を作っているのがヒンジ領域（Hinge region）と呼ばれる部分で、EGFRでは「Gln(Q) 791, Leu(L) 792, Met(M) 793」に相当します（あとで解析の際に残基番号を使います）。

キナーゼ阻害剤は結合様式によって分類され「ATP結合部位とその近縁のみ」に結合するものはType-Iと呼ばれるそうです。

木下誉富キナーゼを標的とした構造生物学および創薬の現状日本結晶学会誌 59，174-181(2017)

計算の対象としている複合体の構造では結合部位はこのようなイメージです（見やすさ重視で加工してます）。

アデニンに類似したピロロピリミジン環（？）がヒンジ領域のMet 793としっかり水素結合を作っていそうです。ただし、水素結合に関与しているピリミジン部位の窒素原子の位置は異なるようです。ATPとTAK-285の構造式を２次元でならべるとそっくりに見える部分構造でも、タンパク質との複合体では完全には重ならないというあたりに構造生物学の面白さを感じます。

なお、ピロロピリミジン環の先に伸びるフェニル環部分はEGFRポケットのゲートキーパー、疎水性ポケットへと伸びています。初期のキナーゼ阻害剤に対する獲得耐性変異としてT790Mなどが知られているそうです。確かに結合部位内で変異が起きてしまったらリガンドが結合しづらくなりそうですね。*2

ざっとですがキナーゼ阻害剤のおさらいは以上です。

2. 低分子 - タンパク質複合体のトラジェクトリ解析

2-1. トラジェクトリの動画

最初にPyMolで確認したトラジェクトリの動画を貼っておきます。（前回の記事の最後に貼ったもの）

リガンドが結構動いてますね。

2-2. MDAnalysisでRMSD解析

2-2-1. RMSDの推移

それでは具体的な解析を進めていきましょう！

まずはシミュレーションが崩壊していないかの確認も含めてトラジェクトリ全体のRMSDの推移を見てみます。

MDAnalysisを使うので、データを読み込んでUniverseオブジェクトを作成しアラインメントをとりましょう。

import MDAnalysis as mda
from MDAnalysis.analysis import align

# topology（.pdb）とtrajectory（.dcd）からUniverseオブジェクトを構築
md_universe = mda.Universe("complex_simulation_data/complex_prep.pdb", "complex_simulation_data/trajectory.dcd")

# Cα原子でアラインメントを取る
alignment = align.AlignTraj( mobile = md_universe, reference = md_universe, select="name CA", in_memory=True)
alignment.run()

トラジェクトリ全体のRMSDはMDAnalysis.analysis.rms.RMSD()で確認できます。

今回は主鎖（backbone）のRMSDを取ると同時に、タンパク質全体（protein）とリガンド（resname LIG）も追加で計算します。*3

追加の計算対象は「引数 groupselections」で指定します。selectで選んだメインのAtomGroupについて重ね合わせた状態で、追加のAtomGroupについて計算が行われます。

from MDAnalysis.analysis import rms

# フレームのリセット
md_universe.trajectory[0]

# RMSD解析（メイン "backbone", 追加 ["protein", "resname LIG"]）
RMSD_analysis = rms.RMSD(md_universe, reference = md_universe, select = "backbone", groupselections =  ["protein", "resname LIG"], ref_frame=0)
RMSD_analysis.run()

結果をプロットしましょう。

results.rmsdに「[[frame, time(ps), RMSD(A), RMSD(A), RMSD(A)], [...], [...]]」形式で入っているので、PandasのDataFrameにしてからmatplotlibでプロットします。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 結果のDataFrameへの変換
column_names = ['frame', 'time(ps)', 'backbone', 'protein', 'ligand']
rmsd_df = pd.DataFrame(np.round(RMSD_analysis.results.rmsd, 2), columns = column_names)

# 結果のプロット
plt.plot(rmsd_df['time(ps)'], rmsd_df['backbone'], label="backbone")
plt.plot(rmsd_df['time(ps)'], rmsd_df['protein'], label="protein")
plt.plot(rmsd_df['time(ps)'], rmsd_df['ligand'], label="ligand")
plt.legend(loc="upper left")
plt.xlabel("Time (ps)")
plt.ylabel("RMSD [$\AA$]")
plt.title("RMSD of protein and ligand")

plt.show()

主鎖が青色、タンパク全体がオレンジ色、リガンドが緑色です。タンパク質に比べるとリガンドの振れ幅が大きいですね。ただし、これはbackbone基準で合わせた場合なので、リガンドだけを考慮した場合とは異なることにご注意ください。

2-2-2. フレーム間のRMSD

先に見たRMSDは最初のフレーム（トラジェクトリ 0 ps）に対するものだったので、フレーム間のペアワイズのRMSDを見てみましょう。diffusionmap.Distancematrixを使います。

from MDAnalysis.analysis import diffusionmap

# proteinについてペアワイズRMSDを計算
pairwase_protein = diffusionmap.DistanceMatrix(md_universe, select = "protein")
pairwase_protein.run()
dist_matrix_protein = pairwase_protein.results.dist_matrix

# ligandについてペアワイズRMSDを計算
pairwase_ligand = diffusionmap.DistanceMatrix(md_universe, select = "resname LIG")
pairwase_ligand.run()
dist_matrix_ligand = pairwase_ligand.results.dist_matrix

結果は2Dのマトリックスです（ndarray）。ヒートマップで可視化しましょう！

proteinとligandでヒートマップの色のレンジを合わせるため、最初にフレーム間の距離の最大値をとっておいてvmaxに指定します。

# コードはTeachOpenCADD T20より
# フレーム間の距離の最大値を取得（heatmapのvmaxに使用）
max_dist = max(np.amax(dist_matrix_ligand), np.amax(dist_matrix_protein))

# プロットを作成
fig, ax = plt.subplots(1, 2)
fig.suptitle("RMSD between the frames")

# タンパクについてsubplot
img1 = ax[0].imshow(dist_matrix_protein, cmap="viridis", vmin=0, vmax=max_dist)
ax[0].title.set_text("protein")
ax[0].set_xlabel("frames")
ax[0].set_ylabel("frames")

# リガンドについてsubplot
img2 = ax[1].imshow(dist_matrix_ligand, cmap="viridis", vmin=0, vmax=max_dist)
ax[1].title.set_text("Ligand")
ax[1].set_xlabel("frames")

fig.colorbar(img1, ax=ax, orientation="horizontal", fraction=0.1, label="RMSD (Å)")

タンパク質のRMSDは時間推移で近いフレーム間では小さく、離れたフレーム間では大きくなっているようです。時間発展に伴って少しづつ全体の構造が変わっているようですね。

対して、リガンドではフレームが離れた箇所でもペアワイズのRMSDが小さくなる部分がありました。時間経過のプロットでも大きくRMSDが動いている箇所があったので、一定の範囲内で構造の揺り戻しが起きているのかもしれません。

2-3. Protein-Ligand Interaction Fingerprint（ProLIF）で相互作用解析

2-3-1. ProLIF？

つづいて具体的に「どの残基がリガンドと相互作用をしているか？」解析したいと思います。

ここではProLIF（Protein-Ligand Interaction Fingerprints）というツールを使います。タンパク質とリガンドとの相互作用をフィンガープリント形式にして解析してくれる素敵なツールです。

論文：Bouysset, C., Fiorucci, S. ProLIF: a library to encode molecular interactions as fingerprints. J Cheminform 13, 72 (2021).
ProLIF GitHub
ProLIF Documentation

ProLIFとは別のツールで恐縮ですが、動作のイメージをつかむためにPLIP（Protein-Ligand Interaction Profiler）の解析図を貼っておきます。*4

上図のようにPLIPではタンパク質とリガンド間にある（非共有結合の）相互作用を種類別に分けて解析してくれます。「ある残基との間にある種類の相互作用が有るか？無いか？」をフィンガープリントとして捉えるというイメージです。

PLIPの入力はPDBファイルで１点の静止画ですが、ProLIFを使うとMDトラジェクトリ（動画）を対象とした解析ができます。

ProLIFのインストールにはいくつか必要なライブラリがあります。ドキュメントのInstallationを参考にインストールしましょう。

推奨はcondaで専用の仮想環境を作った上で、githubからpipでインストールするみたいです。

# 仮想環境を作る
conda create -n prolif

# 仮想環境のアクティベート
conda activate prolif

# 主な依存ライブラリのインストール（これら以外は自動でインストールされる）
conda config --add channels conda-forge
conda install rdkit cython

# ProLIFのインストール
pip install git+https://github.com/chemosim-lab/ProLIF.git

※ 注意「pip install prolif」は使えないそうです。

私はMDお試し用の仮想環境に上記でインストールしました。

2-3-2. ProLIFでフィンガープリント解析

それではProLIFで相互作用フィンガープリントを解析してみましょう。

ProLIFはMDAnalysisで読み込んだトラジェクトリを対象として解析できます。リガンドおよびタンパク質の指定が必要なのでMDAnalysisでslect_atomsしてからProLIFのFingerprintを使いましょう。

import prolif as plf

# MDAnalysisでリガンドとタンパクを選択
lig = md_universe.atoms.select_atoms("resname LIG")
prot = md_universe.atoms.select_atoms("protein")

# デフォルトの相互作用解析
fp = plf.Fingerprint()

# トラジェクトリを渡して解析実行
fp.run(md_universe.trajectory[::10], lig, prot)

上記はトラジェクトリを10フレームごとに解析する場合です（スライスで指定）。今回は50フレームのトラジェクトリなので５セットの結果が出ます。

デフォルトではリガンドの周囲6.0Åの範囲にある残基を対象として解析を行います。

結果はPandasのDataFrameに出力することができます。

df_fp = fp.to_dataframe()
df_fp

ProLIFのデフォルトの解析で見つかった相互作用（True）は「0フレーム目でASP155.Aが水素結合ドナー（HBDonor）」として働いているというもののみでした。残りのフレーム、また他の残基との間に相互作用は検出できなかったようです。

今回のシミュレーション構造は「残基番号 700 〜 1017」なのでASP155.Aは「鎖AのD855」に相当します。ちょうどDFGループのDですね。

期待していたのはヒンジ領域（QLM）との相互作用でしたが出なかったものは仕方ない...残念。

2-3-3. ProLIFでLigand Interaction Network（LigPlot）

DataFrameのみだと味気ないのでもう少し凝った描画をしてみましょう。

ProLIFのLigNetworkを使うと、LigProtのようなお洒落な２Dダイアグラムが作成できます。

引数にはProLIFの「① 相互作用フィンガープリント解析結果（ifp）のデータフレーム」と「②リガンドをMoleculeオブジェクトにしたもの」を渡します。ただし、データフレーム（①）は「引数return_atoms=True」を設定している必要があるのでもう一度作り直しています。

from prolif.plotting.network import LigNetwork

# ① LigPlot用のifpのDataFrame
df_lpt = fp.to_dataframe(return_atoms = True)

# ② リガンドをMDAnalysisのオブジェクトからProLIFのMoleculeオブジェクトにしておく
lmol = plf.Molecule.from_mda(lig)

# ifpのDataFrameから2Dネットワークを作成
# フレーム 0から2D描画を作る場合
net = LigNetwork.from_ifp(df_lpt, lmol, kind='frame', frame=0)
net.display()

できました！アスパラギン酸（ASP155）の酸性側鎖とリガンド（TAK-285）の水酸基との間に水素結合ができているようです。シンプルでみやすいですね。

今回は最初のfram 0でのみ相互作用が検出（True）されていたので LigNetwork.from_ifpの「引数 kind='frame', frame=0」としました。

全てのフレームを合算した結果を描画したい時は「kind='aggregate'」とすればOKです。また、「threshold=0.3」とすれば合算した３割以上（閾値 0.3以上）の割合で見えている相互作用に限定して描画することもできます。

ちなみにPyMolでトラジェクトリの最初の方を確認するとこんな感じです。

ヒンジ領域も結構近いのに検出してくれないのかー。厳しい。。。*5

以上でProLIFを使った解析はおしまいです。

2-4. MDAnalysisで原子間距離の解析

つづいてMDAnalysisを使ってもう少し細かい解析を行いましょう。ProLIFで全体をざっとみただけでは拾いきれなかったものが見えるかもしれません。

まずは原子間距離を測ってみたいと思います。

シミュレーション時にリガンドの各原子につけた名前は以下のようになっています。

また、タンパク質は以下のように前処理をしているため、シミュレーション上の残基の番号は元の構造と700ずれます。

MDAnalysisでは各原子を「残基ID（resid） or 残基名（resname）」と「原子の名前（name）」の組み合わせで指定することができます（select_atoms）。

先の注意点を踏まえると「Met 793（→resid 93）の主鎖N原子」と「リガンド（resname LIG）のN4」は以下のように指定できます。

# EGFRのMet793のN
atomgroup_M793N = md_universe.select_atoms("resid 93 and name N")

# LigandのN4
atomgroup_LIG_N4 = md_universe.select_atoms("resname LIG and name N4")

この２原子の距離はdistances.distを使って求められます。

from MDAnalysis.analysis import  distances

dist_M793N_LigN4 = distances.dist(atomgroup_M793N, atomgroup_LIG_N4)
print(dist_M793N_LigN4)

#  [[93.        ]
#  [ 1.        ]
#  [ 3.30475404]]

戻り値はnumpy.ndarrayで「[[一つ目のAtomGroupの残基ID], [２つ目のAtomGroupの残基ID], [原子間距離(Å)]] 」となっています。

従ってMet残基とリガンドN4原子の距離は3.3Åとなります。トラジェクトリに沿った変化を眺めてみましょう。

# トラジェクトリの各フレームごとに原子間距離を求めてリストに格納
distlist_M793N_LigN4 = []

for _ in md_universe.trajectory:
    dist_t = distances.dist(atomgroup_M793N, atomgroup_LIG_N4)[2][0]
    distlist_M793N_LigN4.append(dist_t)

# 最小値、最大値
print("shortest length(A) : ", np.round(min(distlist_M793N_LigN4), 2))
print("longest length(A) : ", np.round(max(distlist_M793N_LigN4), 2))

# プロット
plt.plot(distlist_M793N_LigN4)
plt.ylabel("distance (Å)")
plt.xlabel("frame")
plt.ylim(0.0, 4.0)
plt.title("Distance between Met793's N and LIG's N4")
plt.show()

# shortest length(A) :  2.92
#  longest length(A) :  3.92

全体では「2.9Å〜 3.9Å」の範囲を行ったり来たりしているようです。

ちょっと長めですが最小値周りでは水素結合を形成する距離に入っていると言って良いかもしれません。

ところでこういう場合、結晶構造の分解能(PDB ID : 3POZ, Resolution: 1.50 Å）はどう考慮すれば良いのでしょうか？MD計算したから解消されていると考えるのか？計算後も誤差として分解能を考慮すべきか？・・・よくわかりません。ご教示いただけると嬉しいです。

2-5. MDAnalysisで水素結合の解析-1

さて、原子間距離の解析はできましたが水素結合を考慮するには角度も合わせて考える必要があります。

MDAnalysisには水素結合解析の機能も用意されています（hydrogenbonds.hbond_analysis））。

設定する引数が多いですが１つ１つはシンプルです。

① ドナー（D）原子（donors_sel(str)）
② 水素原子（H）本体（hydrogens_sel(str)）
③ アクセプター（A）原子（acceptors_sel(str)）
④ ドナーと水素原子（D-H）のカットオフ距離（d_h_cutoff(float)、デフォルト 1.2）
⑤ ドナーとアクセプター（D-A）のカットオフ距離（d_a_cutoff(float)、デフォルト 3.0）
⑥ ドナー - 水素 - アクセプター（D-H-A）の角度（d_h_a_angle_cutoff、デフォルト 150）

水素結合の定義そのものという感じですね。

それではドナーを「Met 793（→resid 93）の主鎖N原子」、アクセプターを「リガンド（resname LIG）のN4」として解析してみましょう。

from MDAnalysis.analysis.hydrogenbonds import hbond_analysis

# ① ドナー原子
donor_MetN = "resid 93 and name N"
# ② 水素原子
donor_MetH = "resid 93 and name H"
# ③ アクセプター原子
acceptor_LigN = "resname LIG and name N4"

hb_analyzer = hbond_analysis.HydrogenBondAnalysis(md_universe,
                                              donors_sel = donor_MetN,
                                              hydrogens_sel = donor_MetH, 
                                              acceptors_sel = acceptor_LigN,
                                              d_h_cutoff = 1.2,
                                              d_a_cutoff = 3.0,
                                              d_h_a_angle_cutoff = 150)
hb_analyzer.run()

結果はresults.hbondsにnumpy.ndarrayとして格納されていて「[<frame>, <donor atom id>, <hydrogen atom id>, <acceptor atom id>, <distance>, <angle>] 」となっています。

DataFrameに変換しましょう。

column_names = ["frame", "donor_index", "hydrogen_index", "acceptor_index", "distance", "angle"]
hb_anal_df =pd.DataFrame(np.round(hb_analyzer.results.hbonds, 2), columns = column_names)
hb_anal_df.set_index("frame")

全50フレームのうち4フレームでのみ（デフォルトの）水素結合形成条件を満たしました。1割以下なのでProLIFで検出されなかったのも仕方ないかもしれませんね。

このままプロットしてもつまらないのでカットオフを緩めて再解析してみましょう（d_a_cutoff = 3.5）。解析結果のD-A距離とD-H-A角度は一つのプロットにまとめて描画します。

# ドナー - アクセプター カットオフ距離を3.5にして再解析
hb_analyzer_35 = hbond_analysis.HydrogenBondAnalysis(md_universe,
                                              donors_sel = donor_MetN,
                                              hydrogens_sel = donor_MetH, 
                                              acceptors_sel = acceptor_LigN,
                                              d_h_cutoff = 1.2,
                                              d_a_cutoff = 3.5,
                                              d_h_a_angle_cutoff = 150)
hb_analyzer_35.run()

# DataFrame化
hb_anal_df_35 =pd.DataFrame(np.round(hb_analyzer_35.results.hbonds, 2), columns = column_names)
hb_anal_df_35.set_index("frame")

# プロットを作成
fig = plt.figure(figsize=(10, 6))

dist_plot = hb_anal_df_35["distance"].plot(style="o", x="frame")
dist_plot.set_ylabel("distance (Å)")
dist_plot.set_xlabel("frame")
dist_plot.set_ylim(0.0, 4.0)

angle_plot = hb_anal_df_35["angle"].plot(secondary_y=True, style="*", x="frame")
angle_plot .set_ylabel("angle (°)")
angle_plot .set_xlabel("frame")
angle_plot .set_ylim(0.0, 180)

fig.legend(loc="lower right", bbox_to_anchor=(1, 0), bbox_transform=dist_plot.transAxes)
plt.show()

D-A距離（distance）とD-H-A角度（angle）は一見連動して動いているように見えますが所々逆の動き（short length & wide angle）になっているのが面白いですね。

TeachOpenCADD T20では同様の解析を経て、距離のばらつきに対して角度は180度近くで安定化していることからリガンドの結合に重要な相互作用となっているのではないかと考察されていました。・・・なるほど、そう解釈するのか！

2-6. MDAnalysisで水素結合の解析-2

MDAnalysisによる水素結合解析の方法がわかったので、ProLIFで見つかった「ASP 855（→resid 155）の側鎖OD1/OD2原子」と「リガンド（resname LIG）のO3」との相互作用も解析してみましょう。

コードは同じことの繰り返しですが、アスパラギン酸側鎖のカルボキシ基には酸素原子が２つあるので両者で違いがでるか？興味があります。

水素結合のカットオフを緩めて解析してみました（d_a_cutoff = 4.0, d_h_a_angle_cutoff = 120）。

# ① ドナー原子
donor_LigO = "resname LIG and name O3"
# ② 水素原子
donor_LigH = "resname LIG and name H22"
# ③ アクセプター原子−１
acceptor_AspOD1 = "resid 155 and name OD1"
# ③ アクセプター原子−２
acceptor_AspOD2 = "resid 155 and name OD2"

# Asp OD1原子について解析
hb_analyzer_OD1 = hbond_analysis.HydrogenBondAnalysis(md_universe,
                                              donors_sel = donor_LigO, hydrogens_sel = donor_LigH, 
                                              acceptors_sel = acceptor_AspOD1,
                                              d_h_cutoff = 1.2, d_a_cutoff = 4.0, d_h_a_angle_cutoff = 120)
hb_analyzer_OD1.run()

# Asp OD1原子の解析結果データフレーム
hb_anal_OD1_df =pd.DataFrame(np.round(hb_analyzer_OD1.results.hbonds, 2), columns = column_names)

# Asp OD2原子について解析
hb_analyzer_OD2 = hbond_analysis.HydrogenBondAnalysis(md_universe,
                                              donors_sel = donor_LigO, hydrogens_sel = donor_LigH, 
                                              acceptors_sel = acceptor_AspOD2,
                                              d_h_cutoff = 1.2, d_a_cutoff = 4.0, d_h_a_angle_cutoff = 120)
hb_analyzer_OD2.run()

# Asp OD2原子の解析結果データフレーム
hb_anal_OD2_df =pd.DataFrame(np.round(hb_analyzer_OD2.results.hbonds, 2), columns = column_names)

まずはD-A距離についてAspの2つの酸素原子（OD1/OD2）の結果を一緒にプロットしてみましょう。

# D-A距離についてOD1/OD2を同時にプロット
plt.scatter(hb_anal_OD1_df["frame"], hb_anal_OD1_df["distance"], color = 'red', marker="o", label="with OD1")
plt.scatter(hb_anal_OD2_df["frame"], hb_anal_OD2_df["distance"], color = 'blue', marker="*", label="with OD2")

plt.xlabel("frame")
plt.ylabel("distance (Å)")
plt.ylim(0.0, 4.0)

plt.legend()
plt.show()

Asp855側鎖のカルボキシ基のうち、OD2原子（red circle）はリガンド水酸基と4.0Åの距離に入ったのは3/50フレームのみでいずれもギリギリ4.0以下でした。

一方で、OD1原子（blue star）は26/50フレームで4.0Å以下の距離に入り、3.0Å以下の点も多くみられます。

同様にD-H-A角度についてもプロットしてみましょう。

# D-H-A角度についてOD1/OD2を同時にプロット
plt.scatter(hb_anal_OD1_df["frame"], hb_anal_OD1_df["angle"], color = 'red', marker="o", label="with OD1")
plt.scatter(hb_anal_OD2_df["frame"], hb_anal_OD2_df["angle"], color = 'blue', marker="*", label="with OD2")

plt.xlabel("frame")
plt.ylabel("angle (°)")
plt.ylim(0, 180)

plt.legend()
plt.show()

角度についても、OD2原子（red circle）は3フレームとも150°よりも小さく、水素結合の基準を満たしそうにありませんでした。

また、OD1原子（blue star）も半数近くは150°以下となっており、距離、角度ともに基準を満たすものは10フレーム前後（2割程度）に落ち着きそう、という結果になりました。

以上の比較から、今回の複合体トラジェクトリではAsp855側鎖カルボキシ基の一方の酸素原子(OD2)とは水素結合を形成するものの、もう一方の酸素原子（OD1）に結合が切り替わることは観察されませんでした。

もう少しつっこむと、30フレーム（600 ps）以降では距離、角度ともに（デフォルトの）水素結合要件を満たすものはありませんでした。つまり、1nsの短いシミュレーションの初期でたまたま見られた偶然の結合の可能性もあります。

ご参考までに計算対象の元のPDB（3POZ）の構造を貼っておきます。

PyMolのデフォルトではリガンド（TAK-285）の水酸基はArg841の主鎖カルボニル酸素およびAsn842の側鎖アミドカルボニル酸素と相互作用しています。一方、タンパク（EGFR）のAsp855はK745と水分子を介したネットワークを形成しているようです。

この共結晶構造解析だけからリガンド OH基とタンパクのAsp855が水素結合を形成する可能性を考慮するのは難しそうです。今回のMDシミュレーションで見られた相互作用をランダムに与えられた初速の産んだ偶然と捉えるか？、X線共結晶では捉えられなかった相互作用と捉えるか？、繰り返しのMDシミュレーションやより長時間のシミュレーションを行えばわかるのかもしれませんが、私の手には負えません。。。。

3. おわりに

以上、今回は低分子ータンパク複合体のMDシミュレーションについてトラジェクトリ解析を試してみました。

特に（イメージしやすい古典的な）水素結合形成の解析に焦点を当てて、ProLIFやMDAnalysisを使った解析を追いかけてみました。・・・といってもほとんどTeachOpenCADD T20とMaking it rainのコードのCopy&Pasteですが。。。

個人的にはProLIFの解析で当初想定していたリガンド TAK-285とタンパク Hinge領域との相互作用が検出できなかったものの、もともとの共結晶構造（PDB ID：3POZ）を眺めるだけでは想像し難いリガンド OH基とタンパク D855 側鎖カルボキシ基との間に水素結合を形成する可能性が示唆された点が非常に興味深かったです。

・・・といっても、MD計算で初期に与えるランダムな速度分布の偶然の結果かもしれません。

より突っ込んだ解析、考察は門外漢の私には手に余るので専門家の皆様にお任せします（丸投げ）。

「MD計算の解釈で元の構造の分解能をどう考慮すれば良いか？」など、色々とわからないことが多かったです。今回も間違いが多いと思うのでご指摘をいただければ幸いです。

ではでは！

*1:PDB ID: 2TIXのリガンドは正確にはATPそのものではなくアナログの「Phosphoaminophosphonic acid-adenylate ester（ANP）」と呼ばれるものです。あくまでイメージとしてご容赦ください。

*2:Structure Based Drug DesignでTAK-285をベースにした構造から展開し、T790M/L858R変異に対しても活性を有する化合物を得る文献もありました。ACS Med Chem Lett. 2012(4)201。PubMedで読めます。

*3:リガンドの名前（LIG）は前回シミュレーション系構築の際につけました

*4:PLIPの文献はこちら Melissa F Adasme, Katja L Linnemann, Sarah Naomi Bolz, Florian Kaiser, Sebastian Salentin, V Joachim Haupt, Michael Schroeder, PLIP 2021: expanding the scope of the protein–ligand interaction profiler to DNA and RNA, Nucleic Acids Research, Volume 49, Issue W1, 2 July 2021, Pages W530–W534,

*5:なお、ProLIFはトラジェクトリ以外のデータも解析できます。Fingerprintオブジェクトを作成後、run()ではなくgenerate()を使えば１対１、run_from_iterable()を使えば複数を対象にした解析ができます。ドッキングの結果の複合体に対して適用したりすることができるようです。参考：Working with docking poses instead of MD simulations