OpenMMをステップバイステップで 〜Part 10 全体のまとめと感想、あと失敗例〜

OpenMMの使い方を順番にたどるシリーズのまとめです。

この記事はシリーズ全体の総括と各記事へのリンクまとめが主目的です。

また、備忘録としてうまくできなかったことものせておきます。ProLIFの踏み込んだ設定やParmEdのファイル変換機能にふれています。

1. 全体の総括

「OpenMMをステップバイステップで」というタイトルでMDシミュレーションの計算の設定や実行、結果の解析方法について複数の記事にまたがって眺めてきました。

前回で「タンパク質ー低分子複合体のシミュレーション解析」まで確認し、（Part 1で設定した）当初の目標 「TeachOpenCADDのT019 · Molecular dynamics simulationとT020 · Analyzing molecular dynamics simulations の内容」を全てたどり終わりました。区切りが良いのでこのシリーズは一度終わりとして、全体の総括をしたいと思います。

各記事でとりあげた内容とツールを以下の表にまとめています。最初からタンパク質-低分子複合体のシミュレーションを行うのはハードルが高かったので、タンパク質単独の前処理方法から始めて、順番に計算の設定・実行、解析を眺めました。

記事リンク	おもな内容	ツール
Part 1	GUIツールを使ったタンパク質構造の前処理	OpenMM Setup
Part 2	Part 1の前処理操作をPythonで実行	PDBFixer
Part 3	GUIツールを使ったシミュレーション条件の設定	OpenMM Setup
Part 4	Google Colab上でMDシミュレーションを実行	OpenMM
Part 5	シミュレーション実行スクリプトの仕組みと書き方(Python)	OpenMM
Part 6	トラジェクトリ解析① (RMSD,RMSF,回転半径)	MDAnalysis, NGLView
Part 7	トラジェクトリ解析② (PCA, cross correlation)	MDAnalysis, pytraj
Part 8	低分子-タンパク質複合体の計算系の設定と計算	OpenMM, OpenFF, MDTraj
Part 9	低分子-タンパク質複合体のトラジェクトリ解析	ProLIF, MDAnalysis

メインに利用したツールはOpenMMでしたが、他にも様々なオープンソースのツールが必要でした。出てくるツール群の概要はだいたい以下の記事にまとめてあります。

magattaca.hatenablog.com

また、追加の情報として「OpenMMのよくある質問（FAQ @ Github）」を以下に翻訳しました。MD計算の設定や結果について書かれていて参考になりました。

magattaca.hatenablog.com

2. きっかけと感想

さて、今回MDシミュレーションを試してみようと思ったのは言うまでもなくTeachOpenCADDという素晴らしい教材が用意されていたからですが、それに加えてTwitterでMaking it rainというプロジェクトを教えていただいたことも大きかったです。

Making it rainについてはこちらの記事で取り上げました。

magattaca.hatenablog.com

MD計算にあこがれはあったものの、難しそうだし高性能なコンピュータがないとどうしようもないと思っていましたが、Google Colabを利用して計算できること、使いやすいオープンソースのプログラムが色々と開発されていると知って「これは遊べるぞ！！」となりました。

OpenMMにGUIツール（OpenMM Setup）が提供されていることを知ったのも大きかったです。「① GUIで使い方の流れと出力を把握 → ② Pythonのコードで同じことができるか確認」という手順を踏めたことで、私でも凡その流れを把握することができました。

実際に計算・解析してみて、「おおっ本当にタンパク動くんだな！」「X線だと綺麗にはまっているリガンドもMDだと割とパタパタ動くんだな！」となって楽しかったです。

特に、複合体のMDトラジェクトリで元の複合体構造にはなかった水素結合がリガンドとタンパクの間に形成されていたのは驚きでした。「相互作用も動的なもので、共結晶構造だけが唯一の正解の結合ポーズではないんだなー。」ってことが実感できたのが何よりでした。

さて、優秀な先人の皆様のおかげで学習する環境は整っているものの、私にとっては難しく、予想以上に長期のシリーズ記事になってしまいました。ここまでお付き合いいただきありがとうございました。

取り上げている計算内容はTeachOpenCADDとほとんど一緒なのに長さは約5倍！！・・・この無駄な長さの中に少しでも役立つ情報があればとても嬉しいです。

3. このシリーズでできなかったこと

MD計算は非常に範囲が広く、深い分野なのでこのシリーズで扱えなかった内容はとてもたくさんあります。

例えば、タンパク質の構造前処理一つとっても@Ag_smith先生のこちらの記事ほどは踏み込めていません。

qiita.com

ジスルフィド結合や金属イオンなど、難しいポイントはたくさんありそうです。

また、低分子の力場についてもOpen Forcefieldで開発されている新しい力場をうまく適用することができませんでした。

低分子-リガンド複合体のシミュレーションということで、「相互作用のエネルギーはどうなるのか？」「結合親和性の見積もりにはどのように活かせば良いのか？」といった点も非常に興味深いところではありますが、私の能力では手に負えず扱うことができませんでした。

AmberではMMPB-SA法（？）といったツールがあるようです。（以下の記事ではシミュレーション時間 500nsくらいは欲しい、といった記載がありました。Google ColabのGPUでも1nsに1時間くらいかかったのに・・・）

blog.livedoor.jp

OpenMMを使った結合自由エネルギー計算ツールとしてはChodera Labの開発したYANKといったものがあるようです。

まだまだ他にもトピックがたくさんありますね！もっといい方法があるよー、などなどご教示いただけると嬉しいです。

4. （おまけ）試したけどうまくいかなかったもの

総括は以上として、以下にはうまくいかなかった操作を備忘録として置いておきます。

4-2. ProLIFの相互作用パラメータを変えた再解析

4-2-1. 課題

記事（Part 9）では「低分子 - タンパク質複合体」のMDシミュレーションのトラジェクトリ解析で水素結合の解析を行いました。その中で、ProLIFとMDAnalysisに水素結合検出の差がありました。この点についてProLIFの相互作用の定義を変更することで解析結果が変わるか検討しました。

4-2-2. おさらい

まずおさらいですが計算対象の構造はPDB ID: 3POZでEGFRタンパクのキナーゼドメインとリガンド TAK-285のX線共結晶構造です。

図の通り共結晶構造からはEGFRのヒンジ領域「Met793の主鎖NH」と「リガンドのN原子」との間に水素結合ができても良さそうです。

実際、MDAnalysisによるトラジェクトリ解析では、デフォルトの定義「ドナー(D) - アクセプター(A) 距離 < 3.0Å」、「D - H - A角度 > 150°」を満たすフレームがいくつか見つかりました。

一方で、ProLIF を用いたInteraction Fingerprint解析（IFP）では、この領域に水素結合は検出されませんでした。

以上の原因を探るためProLIFの相互作用検出の定義を確認し、期待した相互作用を検出できるように定義を変更できるか検証します。

4-2-3. ProLIFの相互作用の定義

まず、ProLIFの相互作用の定義について、元論文の記述を確認しましょう。

• 論文：Bouysset, C., Fiorucci, S. ProLIF: a library to encode molecular interactions as fingerprints. J Cheminform 13, 72 (2021). （オープンアクセス、CC BY 4.0）

判定基準の距離と角度はTable 2に書かれています。

図中に枠で示した通り、水素結合（HBAcceptor、HBDonor）は「距離 ≦ 3.5 Å、角度 130° 〜 180°」です。MDAnalysisの定義よりも緩いですね。こちらは問題なさそうです。

それでは判定基準の原子パターンはどうでしょうか？こちらはSMARTSパターンで定義されているようです。

HBAcceptorは「（N or O or F or アニオン）and (カチオンではない)」です（多分）。「中性か負電荷を帯びているN原子、O原子、F原子」と考えると水素結合のアクセプターとして納得です。

一方でHBDonorは原子番号を使って定義をされていて「窒素(#7)-H or 酸素(#8)-H or 硫黄(#16)-H」です。硫黄が入っているところに生体分子をターゲットとしたツールであることを感じますね。

ここで注目したいのはHBAcceptorの原子の種類の指定方法です。「N」のように大文字のアルファベットで書いた場合、SMARTSでは「アリファティックな窒素原子」を意味するはずです*1。ですので、芳香属性の窒素原子はこのアクセプターの定義では検出できない可能性があります。

4-2-4. RDKitでSMARTSパターンの認識をチェック

ProLIFの定義をふまえて「芳香族性の原子が水素結合アクセプターからはじかれている」可能性が出てきました。

「SMARTSパターンでどのような構造が認識されるか？」、RDKitで検証してみましょう。ProLIFも構造式の取り扱いはRDKitをベースとしているようなので、RDKitで検証しておけば挙動が予測できるはず・・・多分。

今回計算対象のリガンド（TAK-285）の構造を作って、HBAcceptorのSMARTSパターンと部分構造検索でマッチするか検証します。

from rdkit import Chem
from rdkit.Chem import Draw

# SMILESからTAK-285を構築
smi_tak285 = "CC(C)(CC(=O)NCCn1ccc2c1c(ncn2)Nc3ccc(c(c3)Cl)Oc4cccc(c4)C(F)(F)F)O"
rd_tak285 = Chem.MolFromSmiles(smi_tak285)
rd_tak285H = Chem.AddHs(rd_tak285)

Draw.MolToImage(rd_tak285H)

SMARTSからMolオブジェクトを作ってGetSubstructMatchesでマッチする部分構造を見つけ、ハイライトして描画します。

# SMARTSのMolオブジェクト化
HBA_query = Chem.MolFromSmarts('[N,O,F,-{1-};!+{1-}]')

# 部分構造検索と該当箇所の取り出し
HBA_tak285 = rd_tak285H.GetSubstructMatches(HBA_query)
matched_HBA = [a[0] for a in HBA_tak285]

# 描画
Draw.MolToImage(rd_tak285H, highlightAtoms = matched_HBA)

少し小さいですが赤い丸で示されているのがSMARTSに適合した部分です。構造式の左から時計回りにOH基、アミド、アニリンN、エーテルのO、CF3基のFと言ったものが認識されています。

ここで芳香族ヘテロ環（ピロロピリミジン環？）の窒素原子はクエリにかかっていないことに注目です。

MDAnalysisで解析した対象のリガンド窒素原子は、ProLIFではそもそも水素結合アクセプターを取りうる原子として認識されていなかった可能性が示唆されました。

ちなみに、HBDonorの定義のように原子番号でSMARTSを定義するとどうなるでしょうか？

# N → #7に変更して再検証
# SMARTSのMolオブジェクト化
HBA_query_2 = Chem.MolFromSmarts('[#7,O,F,-{1-};!+{1-}]')

# 部分構造検索と該当箇所の取り出し
HBA_tak285_2 = rd_tak285H.GetSubstructMatches(HBA_query_2)
matched_HBA_2 = [a[0] for a in HBA_tak285_2]

# 描画
Draw.MolToImage(rd_tak285H, highlightAtoms = matched_HBA_2)

今度はリガンドの全ての窒素原子がハイライトされました。水素結合ドナーになりそうにない窒素原子（ピロール様のN？）まで含んでしまっているところが問題ではありますが、、、適切なクエリを作るのは難しいですね。

4-2-5. ProLIFの相互作用パラメーターを変えよう

RDKitによる検証で水素結合ドナーの定義を変えてやれば、前回の宿題の水素結合も認識される可能性が出てきました。

ProLIFでは検出する相互作用のパラメーターを変えたり、新しい相互作用を定義することができます。

• ドキュメント①：2.1 Changing the parameters for an interaction
• ドキュメント②：2.2. Writing your own interaction

ProLIFのデフォルトの各相互作用のクラスはドキュメントのprolif.interactionsに記載があり、水素結合アクセプターは「prolif.interactions.HBAcceptor(donor='[#7,#8,#16][H]', acceptor='[N,O,F,-{1-};!+{1-}]', distance=3.5, angles=(130, 180))」となっています。

アクセプターのパラメーターを変えた新しいクラス（HBAcceptor_allN）とを作ってみましょう！

import prolif as plf

class HBAcceptor_allN(plf.interactions.HBAcceptor):
    def __init__(self):
        super().__init__(acceptor="[#7,O,F,-{1-};!+{1-}]")  

元々のクラスを継承した新しいクラスを作成し、super()で親クラスの__init__メソッド実行時のパラメータを書き換えています*2

ここではアクセプターのSMARTSのNを#7に変更しました。

同様にして水素結合ドナー（prolif.interactions.HBDonor(donor='[#7,#8,#16][H]', acceptor='[N,O,F,-{1-};!+{1-}]', distance=3.5, angles=(130, 180))）についても新しいクラスを作っておきましょう。

class HBDonor_allN(plf.interactions.HBDonor):
    def __init__(self):
        super().__init__(acceptor="[#7,N,O,F,-{1-};!+{1-}]")

これらの新しい水素結合クラスを使って、ProLIFの水素結合認識範囲はひろがるでしょうか？？？

4-2-6. MDトラジェクトリの再解析

それではMDトラジェクトリをもう一度ProLIFで解析してみましょう。

まずはMDAnalysisで読み込んで解析の準備をします。

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()

リガンドとタンパクを選択します。

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

新しくつくったProLIFの相互作用（['HBAcceptor_allN', 'HBDonor_allN']）を指定したFingerprintを作ります。

# 解析したい相互作用のFingerprint
fp = plf.Fingerprint(['HBAcceptor_allN', 'HBDonor_allN'])

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

df_fp = fp.to_dataframe()
df_fp

検出されたのは前回と同じくASP155との相互作用だけでした。解析対象をトラジェクトリの全フレームにしてもASP155以外は見つかりませんでした。

定義を変えてacceptorの窒素原子を全て（アリファティック以外も含む）に変更しても、目的のヒンジ領域との相互作用は検出できないようです。残念。

4-3. OpenMMからAmberのファイルに変換してMMPBSA

先にAmberToolsにはエネルギー計算（MMPB-SA.py）するためのツールがあると書きました。

Google検索しているとOpenMMで計算した結果をAmberの対応ファイルに変換してエネルギー計算しようとしている例などが出てきました。

そこでまずはファイル変換方法を試しました。

4-3-1. ParmEdのインストール

ファイルの変換ができるツールとしてParmEdがあります。（Making it rain でも使われていました）

• ParmEd ドキュメント
• ParmEd GitHub

インストールはcondaでOKです。

conda install -c conda-forge parmed

4-3-1. ParmEdのファイル変換でエラー

ParmEdでファイル変換を行うには、目的の構造をParmEdのStructureオブジェクトに変換したのち.save()で望みのファイル形式を選択します。

OpenMMからParmEdのStructureを作るには、load_fileで「①トポロジーと②系（SystemのXML）」を読み込みます。

import parmed as pmd

# トポロジーのためにPDB読み込み
test_pdb = pmd.load_file("complex_simulation_data/complex_prep.pdb")

# OpenMMで作ったsystem.xml情報の読み込み
test_system_xml = pmd.load_file("complex_simulation_data/system.xml")

# ParmEdのStructure作成
test_struct = pmd.openmm.load_topology(test_pdb.topology, test_system_xml)

ところがこれをAmberのファイル形式（prmtop or parm7）に変換しようとするとエラーが出ました。

test_struct.save("test_amb.prmtop")

"""エラ-メッセージ
〜省略〜
       1579 for bond in self.bonds:
    -> 1580     bond.type.used = True
       1581 self.bond_types.prune_unused()
       1582 data['BOND_FORCE_CONSTANT'] = [type.k for type in self.bond_types]

    AttributeError: 'NoneType' object has no attribute 'used'
"""

以上のように、「AttributeError: 'NoneType' object has no attribute 'used'」というエラーが出ます。bond_typeがNoneTypeになっているのがいけないようです。

4-3-2. エラーの原因と回避策

同様のエラーの報告例があり、回避策も２通りありました。

• GitHub Issue: Bond-related error saving AMBER prmtop from ParmEd structure #930
• Issues Antenna: Error writing out prmtop file from OpenMM system

上記の議論によると、OpenMMでforcefieldからcreateSystemで系を構築する際に、拘束条件（constraints）が設定されているとbondのパラメータを消してしまうそうです。そのためParmEdがbond typeをアサインできずNoneのエラーになるとのこと。

回避策①として、createSystemの引数にflexibleConstraints=Trueを与えて系を構築する方法が挙げられていました。

例えば、Part 8で使ったコードを修正すると以下のようにします。

# flexibleConstraints=TrueでSystemをつくる
system = FF.createSystem(complex_pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1.0*unit.nanometers, 
                         constraints = HBonds, rigidWater=True, ewaldErrorTolerance=0.0005, flexibleConstraints=True)

# SystemをXMLで出力する
with open("flexConst_system.xml", mode="w") as file:
    file.write(XmlSerializer.serialize(system))

以上のようにして作成したSystemのXMLをParmEdで読み込んでみます。

# 作り直したSystemのxmlを読み込み
test_system_xml2 = pmd.load_file("flexConst_system.xml")

# ParmEdのStructure作成
test_struct2 = pmd.openmm.load_topology(test_pdb.topology, test_system_xml2)

# Amber形式で出力
test_struct2.save("test_amb.prmtop")

作り直したSystemの情報を使うと無事「test_amb.prmtop」ファイルが作成されました！

また、OpenMMのSystemを作り直すのが面倒な場合、回避策②としてダミーのbond typeを割り当てる方法も挙げられていました。

# ダミーのbond typeを割り当てる
bond_type = pmd.BondType(1.0, 1.0, list=test_struct.bond_types)
test_struct.bond_types.append(bond_type)
for bond in test_struct.bonds:
    if bond.type is None:
        bond.type = bond_type

# Amber形式で出力
test_struct.save('bond_type_assinged.prmtop')

こちらの方法でも無事「bond_type_assinged.prmtop」ファイルが作成されました！

4-3-3. MMPB-SAの失敗

Amber形式にファイル変換する方法がわかったので、MMPB-SA計算ができるかどうかを試します。

計算に必要なファイルは複合体、タンパク質のみ、リガンドのみのトポロジーのようです。

そこで先のParmEdのStructureからそれぞれ必要な部分構造を抜き出して以下のようにファイルを作成しました。

# タンパク質を残基番号 1-317 で選択して出力
pmd_protein = test_struct[':1-317']
pmd_protein.save('protein.prmtop')

# リガンドを残基名 LIG で選択して出力
pmd_ligand = test_struct[':LIG']
pmd_ligand.save('ligand.prmtop')

# 複合体を出力
pmd_complex = test_struct[':1-317, LIG']
pmd_complex.save('complex.prmtop')

print(pmd_protein)
print(pmd_ligand)
print(pmd_complex)

#  <Structure 5120 atoms; 317 residues; 5176 bonds; PBC (orthogonal); parametrized>
#  <Structure 63 atoms; 1 residues; 66 bonds; PBC (orthogonal); parametrized>
#  <Structure 5183 atoms; 318 residues; 5242 bonds; PBC (orthogonal); parametrized>

原子数を見る限り目的の構造を選択できていそうです。

これらのファイルを使って、Making it rainを参考にante-MMPBSA.py``MMPBSA.pyを使おうとしました。

ところが構造を作れないといったエラーがでてうまくいきませんでした。エラーの中身を見てもよくわかりませんでした。Amberの勉強をしないとわからなさそうです。。。

以上、試してうまくいかなかったこと２点でした。

5. おわりに

以上、今回は関連記事のまとめとMD計算を試した感想でした。また、失敗例は使い道はありませんが、検討中にわかったこと（ParmEdの使い方）もあったので備忘録として残しておきました。

Part 1から3ヶ月以上もかかってしまいましたが、素人のチュートリアルお試し記事にお付き合いいただきありがとうございました！

色々とご指摘・ご教示いただければ幸いです。

ではでは！

*1:参考：Daylight SMARTSの説明

*2:参考：ProLIF ドキュメント 2.1.2. Reparameterizing an interaction with another name