OpenMMの使い方を順番にたどる記事のPart 9です。
前回の記事（Part 8 ）では「低分子 - タンパク質複合体 」についてシミュレーション系の構築 と計算の実行 を行いました。この記事では結果の解析 を行います。
一般的な解析手法についてはPart 6 、Part 7 で確認しました。この記事では「低分子リガンドとタンパク質の相互作用解析 」に焦点を当てたいと思います。
コードは基本的にTeachOpenCADD T020：Analyzing molecular dynamics simulations とMaking it rain : Protein-Ligand smulations からお借りしています。
相互作用フィンガープリント / 水素結合解析をします
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
md_universe = mda.Universe("complex_simulation_data/complex_prep.pdb" , "complex_simulation_data/trajectory.dcd" )
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
追加の計算対象は「引数 groupselectionsselectで選んだメイン のAtomGroupについて重ね合わせた状態で、追加 のAtomGroupについて計算が行われます。
from MDAnalysis.analysis import rms
md_universe.trajectory[0 ]
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
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
pairwase_protein = diffusionmap.DistanceMatrix(md_universe, select = "protein" )
pairwase_protein.run()
dist_matrix_protein = pairwase_protein.results.dist_matrix
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に指定します。
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" )
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" )
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 ）というツールを使います。タンパク質とリガンドとの相互作用をフィンガープリント形式 にして解析してくれる素敵なツールです。
ProLIFとは別のツールで恐縮ですが、動作のイメージをつかむためにPLIP *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
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のLigNetworkLigProt ２Dダイアグラム が作成できます。
引数にはProLIFの「① 相互作用フィンガープリント解析結果（ifp ）のデータフレーム」と「②リガンドをMoleculeオブジェクトにしたもの」を渡します。ただし、データフレーム（①）は「引数return_atoms=True
from prolif.plotting.network import LigNetwork
df_lpt = fp.to_dataframe(return_atoms = True )
lmol = plf.Molecule.from_mda(lig)
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 」は以下のように指定できます。
atomgroup_M793N = md_universe.select_atoms("resid 93 and name N" )
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)
戻り値は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()
全体では「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角度は一つのプロットにまとめて描画します。
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()
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"
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()
hb_anal_OD1_df =pd.DataFrame(np.round(hb_analyzer_OD1.results.hbonds, 2 ), columns = column_names)
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()
hb_anal_OD2_df =pd.DataFrame(np.round(hb_analyzer_OD2.results.hbonds, 2 ), columns = column_names)
まずはD-A距離 についてAsp の2つの酸素原子（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 角度についてもプロットしてみましょう。
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線 共結晶では捉えられなかった相互作用
3. おわりに
以上、今回は低分子ータンパク複合体のMDシミュレーションについてトラジェクトリ解析を試してみました。
特に（イメージしやすい古典的な）水素結合 形成の解析に焦点を当てて、ProLIF やMDAnalysis を使った解析を追いかけてみました。・・・といってもほとんどTeachOpenCADD T20とMaking it rainのコードのCopy&Pasteですが。。。
個人的にはProLIF の解析で当初想定していたリガンド TAK-285とタンパク Hinge領域 との相互作用が検出できなかった ものの、もともとの共結晶構造 （PDB ID：3POZ）を眺めるだけでは想像し難いリガンド OH基 とタンパク D855 側鎖カルボキシ基 との間に水素結合を形成 する可能性が示唆された点が非常に興味深かったです。
・・・といっても、MD計算で初期に与えるランダムな速度分布の偶然の結果 かもしれません。
より突っ込んだ解析、考察は門外漢の私には手に余るので専門家の皆様にお任せします（丸投げ）。
「MD計算の解釈で元の構造の分解能 をどう考慮すれば良いか？」など、色々とわからないことが多かったです。今回も間違いが多いと思うのでご指摘をいただければ幸いです。
ではでは！
